Why introducing a 'blockdiag' constraints argument in exponential random graph models (ergm) drastically slow down the computation?

Question

I am trying to run a valued exponential random graph model (ergm) on a weighted network (network_ex). This network shows the interactions between individuals in four different groups. Interactions between groups cannot occur so a blockdiagonal constraint needs to be included in the model. However when I include this blockdiagonal constraint the computation gets stuck on a specific iteration of the Monte Carlo maximum likelihood estimation for at least 24h (I have stopped the process afterward). When I remove the constraints the computation takes only around 30 seconds. I am puzzled as I would expect the computation to be faster with the blockdiagonal constraints as the number of edges that may vary is lower. Here is a toy example similar to my model:

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

# define response matrix
mat_ex=matrix(
c(0,1,5,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,3,5,12,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,4,0,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,12,1,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10,3,6,11,6,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,1,2,6,13,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,2,0,11,9,10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,6,3,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,3,4,8,0,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,8,1,5,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,6,9,11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,3,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,3,0,9,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10,2,3,0,11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,5,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,3,0,2,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,7,0,3,3,1,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,6,0,7,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,3,0,1,0,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,5,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,2,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,1,4,0,0,2,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_Adult","M_Sub-Adult","F_Adult","M_Adult","F_Adult","M_Sub-Adult","M_Adult","F_Adult","M_Adult","M_Adult","M_Adult","F_Adult","F_Adult","F_Adult","F_Adult","F_Adult","M_Sub-Adult","M_Adult","M_Adult","M_Adult","F_Adult","M_Adult")
Group=c(1,1,1 ,1,1 ,2, 2, 2 ,2, 2, 2, 2, 3, 3, 3, 3, 3, 4 ,4 ,4, 4, 4, 4,4, 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.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.6413537173,0.0000000000,0.6214644129,1.0141235877,4.9439358178,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,1.7611808824,0.6214644129,0.0000000000,0.2849059121,0.0013686734,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,1.1045365579,1.0141235877,0.2849059121,0.0000000000,0.0374005963
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0610428348,4.9439358178,0.0013686734,0.0374005963
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.4951618629,0.5858325592,3.8541140580,3.7315606330,0.0009005127
,0.3271637800,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.4951618629,0.0000000000,0.2980769040,0.9923034427,1.4379618200
,3.4747962774,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.5858325592,0.2980769040,0.0000000000,0.2162076270
,0.6839987193,0.0000000000,0.1628661860,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,3.8541140580,0.9923034427,0.2162076270
,0.0000000000,1.2295260710,0.0000000000,1.1504634987,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,3.7315606330,1.4379618200
,0.6839987193,1.2295260710,0.0000000000,1.8073008815,0.5011267262,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0009005127
,3.4747962774,0.0000000000,0.0000000000,1.8073008815,0.0000000000,0.0506384403,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.3271637800,0.0000000000,0.1628661860,1.1504634987,0.5011267262,0.0506384403,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,1.6109005375,4.0166662575,0.4233161826,0.1060254672,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,1.6109005375,0.0000000000,0.3408496603,0.4853374582,0.8949786399,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,4.0166662575,0.3408496603,0.0000000000,0.0698993674,0.6698120875,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.4233161826,0.4853374582,0.0698993674,0.0000000000,1.3822143417
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.1060254672,0.8949786399,0.6698120875,1.3822143417
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.3009462716,2.3082969479,0.1018110048,0.0326341498,0.3828057647
,5.5739028034,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.3009462716,0.0000000000,0.2237564859,0.0323276664,3.8922582745
,0.5733125772,0.6159767242,1.1723116376,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,2.3082969479,0.2237564859,0.0000000000,2.3784384336
,0.7850055818,1.5181554486,0.4208603649,1.7763122536,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.1018110048,0.0323276664,2.3784384336
,0.0000000000,0.1936917389,1.7158463101,0.0000000000,0.0289622499,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0326341498,3.8922582745
,0.7850055818,0.1936917389,0.0000000000,0.8168095460,0.9342629045,0.7772395917,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.3828057647
,0.5733125772,1.5181554486,1.7158463101,0.8168095460,0.0000000000,0.0115138397,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,5.5739028034,0.6159767242,0.4208603649,0.0000000000,0.9342629045,0.0115138397,0.0000000000,1.4325614287
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000,0.0000000000
,0.0000000000,0.0000000000,1.1723116376,1.7763122536,0.0289622499,0.7772395917,0.0000000000,1.4325614287
,0.0000000000)
,25,25)

mat_count_ex=matrix(c(0,1,8,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,4,6,13,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,4,0,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,4,0,0,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10,13,2,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,3,13,13,12,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,6,9,22,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,2,0,14,11,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,10,3,0,4,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,3,5,10,0,14,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,10,1,6,10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,0,4,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,7,12,16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,4,10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,3,0,13,13,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,14,2,4,0,14,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,5,3,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,3,1,3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,7,1,4,4,3,11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,8,0,8,2,0,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,2,7,0,2,2,3,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,4,5,5,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,7,0,1,2,0,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,2,0,1,2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,5,2,6,1,0,5,0),25,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,"III",form="sum") 
                      +edgecov(mat_count_ex,"Nbobs",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))

And here is what I get:

Best valid proposal ‘DistRLE’ cannot take into account hint(s)‘sparse’. Starting contrastive divergence estimation via CD-MCMLE: Iteration 1 of at most 60: Convergence test P-value:0e+00 The log-likelihood improved by 0.9066. Iteration 2 of at most 60: Convergence test P-value:1.5e-248 The log-likelihood improved by 0.4615. Iteration 3 of at most 60: Convergence test P-value:1.5e-176 The log-likelihood improved by 0.6357. Iteration 4 of at most 60: Convergence test P-value:6.1e-147 The log-likelihood improved by 0.4542. Iteration 5 of at most 60: Convergence test P-value:1.1e-92 The log-likelihood improved by 0.2675. Iteration 6 of at most 60: Convergence test P-value:2.4e-18 The log-likelihood improved by 0.04343. Iteration 7 of at most 60: Convergence test P-value:1.9e-05 The log-likelihood improved by 0.01517. Iteration 8 of at most 60: Convergence test P-value:4e-04 The log-likelihood improved by 0.01191. Iteration 9 of at most 60: Convergence test P-value:6e-03 The log-likelihood improved by 0.009803. Iteration 10 of at most 60: Convergence test P-value:1.7e-03 The log-likelihood improved by 0.01049. Iteration 11 of at most 60: Convergence test P-value:6.2e-02 The log-likelihood improved by 0.006888. Iteration 12 of at most 60: Convergence test P-value:6e-02 The log-likelihood improved by 0.007122. Iteration 13 of at most 60: Convergence test P-value:2.9e-01 The log-likelihood improved by 0.004809. Iteration 14 of at most 60: Convergence test P-value:1.5e-02 The log-likelihood improved by 0.00833. Iteration 15 of at most 60: Convergence test P-value:1.6e-01 The log-likelihood improved by 0.005631. Iteration 16 of at most 60: Convergence test P-value:4.3e-01 The log-likelihood improved by 0.004001. Iteration 17 of at most 60: Convergence test P-value:5.8e-01 Convergence detected. Stopping. The log-likelihood improved by 0.003348. Finished CD. Starting Monte Carlo maximum likelihood estimation (MCMLE): Iteration 1 of at most 200: Optimizing with step length 0.0863421108600662. The log-likelihood improved by 2.006. Estimating equations are not within tolerance region. Iteration 2 of at most 200: Optimizing with step length 0.0875588561417017. The log-likelihood improved by 2.751. Estimating equations are not within tolerance region. Iteration 3 of at most 200:

This runs fast but do not progress anymore.

score 0 · Answer 1 · answered Mar 31 '21 at 22:46

0

A few suggestions:

The development version uses adaptive MCMC by default and automates a lot of other things. I recommend removing all control parameters before (but not including) seed=. (Maybe seed as well, since reproducibility under parallel processing isn't in a great shape.)
To see get a better idea of what the function is doing, pass the verbose=TRUE or verbose=2 or higher. (Higher numbers mean more output, but they can also slow things down.)
To see what the sampler is doing, pass the MCMC.runtime.traceplot=TRUE control parameter.
Unfortunately, this sampler is not very efficient, so you may, in fact, need to run things for a long time. (There is a ticket for this: https://github.com/statnet/ergm.count/issues/5 .)

answered Mar 31 '21 at 22:46

Pavel N. Krivitsky

66
3

@ Pavel N. Krivitsky Thanks a lot for your suggestions, I have now run the ergm function: 1) removing all the control arguments before seed 2) added a "verbose = 4" argument (below 4 I had no additional output at the step where my code get slow) 3) The "MCMC.runtime.traceplot=TRUE" control parameter makes my R session systematically crash when starting the Monte Carlo maximum likelihood estimation so I had to remove it but the trace plots for the CD look ok – Louise Apr 06 '21 at 07:46
4) it is now running continuously for 5 days, after few minutes it reached iteration 25 of the Monte Carlo maximum likelihood estimation but only progressed of two iterations since then – Louise Apr 06 '21 at 07:47
Do you have an idea why the "blockdiag" term itself would slow down the process? If I am not wrong the same sampler is used when I run the ergm.count without constraints right? – Louise Apr 06 '21 at 07:49

Why introducing a 'blockdiag' constraints argument in exponential random graph models (ergm) drastically slow down the computation?

1 Answers1