Adaptive MCMC, especially with complex models, spends a lot of time calculating the ACF.

[krivit]

Here's an example of Rprof results for a recent run of a 60-parameter multilayer model:

It was run on a cluster with parallel=8. The runtime was:

> proc.time()
     user    system   elapsed 
29897.472   771.917 83982.011 

so about a third of the run was spent inside the master process, and of that, about 2/3 was spent calculating the ACFs for the ar() function for adaptive burnin and precision estimation:

> lapply(fit.prof, head, 30)
$by.self
                      self.time self.pct total.time total.pct
"acf"                  20649.44    67.55   20786.08     67.99
"tcrossprod"            5886.48    19.26    6120.02     20.02
"solve.yw"               897.60     2.94    1606.74      5.26
"serialize"              697.84     2.28     697.98      2.28
"cov"                    342.04     1.12     342.50      1.12
"%*%"                    305.80     1.00     305.80      1.00
"NextMethod"             300.52     0.98     300.68      0.98
".Fortran"               231.42     0.76     231.42      0.76
"aperm.default"          197.32     0.65     198.64      0.65
"t"                      105.52     0.35     455.78      1.49
"cal.resid"               89.06     0.29    6267.92     20.50
"rbind"                   82.26     0.27      93.00      0.30
"unclass"                 81.28     0.27      81.28      0.27
"t.default"               58.74     0.19      58.74      0.19
"array"                   50.22     0.16      58.04      0.19
"apply"                   49.48     0.16      97.10      0.32
"ar.yw.default"           45.86     0.15   28920.72     94.60
".Call"                   45.22     0.15      45.22      0.15
"as.matrix.mcmc.list"     40.52     0.13      48.36      0.16
"<Anonymous>"             36.26     0.12      43.40      0.14
"FUN"                     32.00     0.10   21967.70     71.86
"colnames<-"              22.00     0.07      22.06      0.07
"colMeans"                19.72     0.06      28.26      0.09
"ts"                      19.36     0.06      24.70      0.08
"matrix"                  17.44     0.06      18.14      0.06
"unserialize"             16.86     0.06      16.90      0.06
"is.na"                   14.96     0.05      14.96      0.05
"var"                     13.20     0.04      19.62      0.06
"determinant.matrix"      13.16     0.04      13.62      0.04
"any"                     11.58     0.04      11.58      0.04

$by.total
                             total.time total.pct self.time self.pct
"ergm"                         30569.90    100.00      0.00     0.00
"ergm.MCMLE"                   30552.06     99.94      0.04     0.00
"tryCatch"                     30088.82     98.43      0.56     0.00
"tryCatchList"                 30088.76     98.43      0.40     0.00
"tryCatchOne"                  30088.68     98.43      0.16     0.00
"doTryCatch"                   30088.46     98.42      1.14     0.00
"spectrum0.mvar"               29384.16     96.12      9.46     0.03
"withCallingHandlers"          29366.80     96.06      0.90     0.00
".catchToList"                 28921.60     94.61      0.06     0.00
"ar"                           28921.38     94.61      0.10     0.00
"ar.yw"                        28920.92     94.61      0.20     0.00
"ar.yw.default"                28920.72     94.60     45.86     0.15
"ergm_MCMC_sample"             27211.52     89.01      1.98     0.01
"eval"                         22496.04     73.59      1.18     0.00
"try"                          22452.30     73.45      0.02     0.00
"FUN"                          21967.70     71.86     32.00     0.10
"lapply"                       21965.92     71.85      6.40     0.02
"approx.hotelling.diff.test"   21736.26     71.10      0.40     0.00
"within.list"                  21734.30     71.10      0.00     0.00
"within"                       21734.30     71.10      0.00     0.00
"ERRVL"                        21730.14     71.08      0.08     0.00
"acf"                          20786.08     67.99  20649.44    67.55
"sapply"                       20674.72     67.63      0.00     0.00
".find_OK_burnin"              20673.46     67.63      0.00     0.00
"suppressWarnings"             20651.72     67.56      0.00     0.00
"geweke.diag.mv"               20651.32     67.55      0.00     0.00
"cal.resid"                     6267.92     20.50     89.06     0.29
"tcrossprod"                    6120.02     20.02   5886.48    19.26
"solve.yw"                      1606.74      5.26    897.60     2.94
"ergm.estimate"                 1492.76      4.88      1.02     0.00

$sample.interval
[1] 0.02

$sampling.time
[1] 30570.04

That means that about This makes sense, since acf runtime is proportional to the sample size, the order, and the square of the number of parameters. What optimisations could be made?

[handcock]

I'd suggest replacing the acf() call with a Fast Fourier Transform based computation. You could, in addition, subsample the sample with, e.g., lag 10 as well and only use those in the computation. There are other ideas too as it is called repeatedly. Alos, it is a guide only so accuracy is not required.

[krivit]

Thanks, Mark. One limitation is that acf() is being called from within ar(), though I guess we could copy it out (again).

[krivit]

For now, perhaps we should use first-order AR only for adaptive MCMC. Although that would produce biased ESS estimates, they are meant to be largely heuristic anyway.