为什么在指数随机图模型 (ergm) 中引入“blockdiag”约束参数会大大减慢计算速度？

问题描述

我正在尝试在加权网络 (network_ex) 上运行一个有价值的指数随机图模型 (ergm)。该网络显示了四个不同群体中个体之间的相互作用。无法发生组之间的交互，因此模型中需要包含块对角约束。但是，当我包含这个块对角约束时，计算会卡在 Monte Carlo 最大似然估计的特定迭代上至少 24 小时（之后我停止了该过程）。当我删除约束时，计算只需要大约 30 秒。我很困惑，因为我希望块对角线约束的计算速度更快，因为可能变化的边数较少。这是一个类似于我的模型的玩具示例：

library(ergm)
library(ergm.count)
library(tergm)

# define response matrix
mat_ex=matrix(
c(0,1,5,3,12,4,7,10,6,11,2,13,9,8,0),25,25)

nb=1:25
ID=rep("ID_",25)
names=paste0(ID,nb)
rownames(mat_ex)=names
colnames(mat_ex)=names

# define useful nodes attributes
SeAg=c("F_Adult","F_Adult","M_Adult","M_Sub-Adult","M_Adult")
Group=c(1,4)
jjj=c(0.9195853,0.5267635,1.0521188,1.4936954,83.8427083,11.3156841,0.6956844,0.2678993,51.5156250,0.7075267,26.5343467,0.7918371,0.3461749,1.5377506,283.7933438,3.2720311,0.9384706,192.6997677,19.0026389,61.6710180,69.0940213,126.9722367,0.5726383,393.5911824,0.7805982)

# define useful edges attributes
mat_iii=matrix(
  c(0.0000000000,0.6413537173,1.7611808824,1.1045365579,0.0610428348,0.0000000000,0.6214644129,1.0141235877,4.9439358178,0.2849059121,0.0013686734,0.0374005963,0.4951618629,0.5858325592,3.8541140580,3.7315606330,0.0009005127,0.3271637800,0.2980769040,0.9923034427,1.4379618200,3.4747962774,0.2162076270,0.6839987193,0.1628661860,1.2295260710,1.1504634987,1.8073008815,0.5011267262,0.0506384403,1.6109005375,4.0166662575,0.4233161826,0.1060254672,0.3408496603,0.4853374582,0.8949786399,0.0698993674,0.6698120875,1.3822143417,0.3009462716,2.3082969479,0.1018110048,0.0326341498,0.3828057647,5.5739028034,0.2237564859,0.0323276664,3.8922582745,0.5733125772,0.6159767242,1.1723116376,2.3784384336,0.7850055818,1.5181554486,0.4208603649,1.7763122536,0.1936917389,1.7158463101,0.0289622499,0.8168095460,0.9342629045,0.7772395917,0.0115138397,1.4325614287,0.0000000000),25)

mat_count_ex=matrix(c(0,16,22,14,17,25)

# define response network
network_ex=as.network(x = mat_ex,directed = TRUE,loops = FALSE,matrix.type = "adjacency",ignore.eval=FALSE,names.eval='weight')
set.vertex.attribute(network_ex,"Sex_Age",as.character(SeAg))
set.vertex.attribute(network_ex,"Group",as.character(Group))
set.vertex.attribute(network_ex,"JJJ",as.numeric(jjj))
set.edge.attribute(network_ex,"III",as.numeric(mat_iii))
set.edge.attribute(network_ex,"Nbobs",as.numeric(mat_count_ex))

# run ergm 
ergm_ex <- ergm(network_ex ~ 
                       sum
                      +nonzero
                      +mutual(form="min",threshold=0)
                      +nodemix("Sex_Age",levels=c(1,2),form="sum")
                      +nodeicov("JJJ",form="sum")
                      +diff("JJJ",pow=1,dir="t-h",sign.action="identity",form ="sum")
                      +edgecov(mat_iii,form="sum") 
                      +edgecov(mat_count_ex,form="sum"),response="weight",reference=~Poisson,constraints = ~blockdiag("Group"),control = control.ergm(MCMC.interval = 1000,MCMLE.maxit = 200,init.method = 'CD',MCMC.samplesize = 1000,MCMC.prop.weights="random",MCMC.burnin=100,seed=123456,parallel=4))

这是我得到的：

最佳有效提案“distRLE”不能考虑提示“稀疏”。通过 CD-MCMLE 开始对比散度估计：最多 60 次的迭代 1：收敛测试 P 值：0e+00 对数似然提高了 0.9066。最多 60 次的迭代 2：收敛测试 P 值：1.5e-248 对数似然提高了 0.4615。最多 60 次的迭代 3：收敛测试 P 值：1.5e-176 对数似然提高了 0.6357。最多 60 次的迭代 4：收敛测试 P 值：6.1e-147 对数似然提高了 0.4542。最多 60 次迭代 5：收敛测试 P 值：1.1e-92 对数似然提高了 0.2675。最多 60 次迭代 6：收敛测试 P 值：2.4e-18 对数似然提高了 0.04343。迭代 7，最多 60：收敛测试 P 值：1.9e-05 对数似然提高了 0.01517。最多 60 次的迭代 8：收敛测试 P 值：4e-04 对数似然提高了 0.01191。最多 60 次迭代 9：收敛测试 P 值：6e-03 对数似然提高了 0.009803。最多 60 次迭代 10：收敛测试 P 值：1.7e-03 对数似然提高了 0.01049。最多 60 次迭代 11：收敛测试 P 值：6.2e-02 对数似然提高了 0.006888。最多 60 次迭代 12：收敛测试 P 值：6e-02 对数似然提高了 0.007122。最多 60 次迭代 13：收敛测试 P 值：2.9e-01 对数似然提高了 0.004809。最多 60 次迭代 14：收敛测试 P 值：1.5e-02 对数似然提高了 0.00833。最多 60 次迭代 15：收敛测试 P 值：1.6e-01 对数似然提高了 0.005631。最多 60 次迭代 16：收敛测试 P 值：4.3e-01 对数似然提高了 0.004001。最多 60 次迭代 17：收敛测试 P 值：5.8e-01 检测到收敛。停止。对数似然提高了 0.003348。完成的CD。启动蒙特卡洛最大似然估计 (MCMLE)：最多 200 次的迭代 1：优化步骤长度 0.0863421108600662。对数似然提高了 2.006。估计方程不在容差范围内。迭代 2 最多 200：优化步长 0.0875588561417017。这对数似然提高了 2.751。估计方程不在容差区。最多 200 次的迭代 3：

这运行得很快，但不再进展。

解决方法

一些建议：

开发版本默认使用自适应 MCMC 并自动执行许多其他操作。我建议在（但不包括）seed= 之前删除所有控制参数。（也许 seed 也是如此，因为并行处理下的再现性不是很好。）
要更好地了解函数的作用，请传递 verbose=TRUE 或 verbose=2 或更高版本。（更高的数字意味着更多的输出，但它们也会减慢速度。）
要查看采样器在做什么，请传递 MCMC.runtime.traceplot=TRUE 控制参数。
不幸的是，此采样器效率不高，因此实际上您可能需要运行很长时间。（有一张票：https://github.com/statnet/ergm.count/issues/5。）

r r statnet