Model does no contain samples when using vb with init function

[haansn08]

I'm trying to estimate a large correlation matrix (dimensions about 100x100), which makes Stan fail to find initial values, I think because filling all entries with random values from (-2, 2) doesn't give a valid correlation matrix, even if you try 100 times. My idea for a workaround was to choose the initial values using the R package randcorr. However it seems that vb doesn't like having a custom init function and returns a model without any samples. Here is a minimal example to reproduce the problem:

library(rstan)
library(randcorr)

model_code <- '
data {
 array[100] vector[100] y;
}
parameters {
 cholesky_factor_corr[100] y_corr_chol;
}
model {
 y_corr_chol ~ lkj_corr_cholesky(4);

 vector[100] zero = rep_vector(0, 100);
 y ~ multi_normal_cholesky(zero, y_corr_chol);
}
'

mod <- stan_model(model_code = model_code)

initf <- function() {
	list(
		y_corr_chol = t(chol(randcorr(100)))
	)
}

fit <- vb(
	mod,
	cores=1, chains=1,
	init = initf,
	data = list(
		y = matrix(rnorm(100 * 100), 100, 100)
	)
)

> fit
Stan model 'anon_model' does not contain samples.

> sessionInfo()
R version 4.5.1 (2025-06-13)
Platform: x86_64-pc-linux-gnu
Running under: Arch Linux

Matrix products: default
BLAS:   /usr/lib/libblas.so.3.12.0
LAPACK: /usr/lib/liblapack.so.3.12.0  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

time zone: Europe/Vienna
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] randcorr_1.0        rstan_2.32.7        StanHeaders_2.32.10

loaded via a namespace (and not attached):
 [1] vctrs_0.6.5         cli_3.6.5           rlang_1.1.6
 [4] processx_3.8.6      RcppParallel_5.1.10 glue_1.8.0
 [7] gridExtra_2.3       ps_1.9.1            pkgbuild_1.4.7
[10] stats4_4.5.1        scales_1.4.0        grid_4.5.1
[13] tibble_3.2.1        lifecycle_1.0.4     QuickJSR_1.7.0
[16] inline_0.3.21       compiler_4.5.1      codetools_0.2-20
[19] RColorBrewer_1.1-3  Rcpp_1.0.14         pkgconfig_2.0.3
[22] farver_2.1.2        R6_2.6.1            pillar_1.10.2
[25] parallel_4.5.1      callr_3.7.6         magrittr_2.0.3
[28] loo_2.8.0           tools_4.5.1         gtable_0.3.6
[31] matrixStats_1.5.0   ggplot2_3.5.2

[bob-carpenter]

Thanks for filing a reproducible bug---we really appreciate it.

I'm trying to estimate a large correlation matrix (dimensions about 100x100), which makes Stan fail to find initial values, I think because filling all entries with random values from (-2, 2) doesn't give a valid correlation matrix, even if you try 100 times.

This isn't the problem. We parameterize an N x N correlation matrix with (N choose 2) unconstrained parameters. Here's the specific transform we use for correlation matrices:

https://mc-stan.org/docs/reference-manual/transforms.html#cholesky-factors-of-correlation-matrices

For something this size, you should be able to sample with NUTS. When things get bigger, the numerics get tricky and we should have better transforms. I'm guessing it's just another issue with the way autodiff variational inference (ADVI) handles initialization.

And in case it wasn't clear from the doc, ADVI is only giving you a multivariate normal approximation to your posterior on the unconstrained scale. That is, we fit the transformed parameters, then sample there, then transform back to sample from the posterior. This is probably not going to be a great approximation for a correlation matrix because there is a lot of correlation among the parameters making up the correlation matrix due to the two constraints (positive define, rows are unit Euclidean length).

[bob-carpenter]

I forgot to say that it's unlikely anyone is going to work on patching ADVI. We're focused on other VI algorithms these days like normalizing flows.

[haansn08]

Thank you for the explanation. The above example also doesn't work with sampling, unless init=initf is passed:

fit <- sampling(
        mod,
        cores=1, chains=1,
        data = list(
                y = matrix(rnorm(100 * 100), 100, 100)
        )
)

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
...

So something is wrong with the way Stan chooses the initial values for the Cholesky factors per default...

Is there a way to try out the another VI algorithm? I can't seem to find information on "normalizing flows" in the documentation...

[bob-carpenter]

This isn't a bug so much as a numerical instability we don't know how to fix. It's really not intended to fit large correlation matrices like dimension 100. We should have a warning in the functions reference.

All of our code initializes the same way, which is uniform(-2, 2) on the unconstrained scale. You can reduce those bounds, which often leads to more stability on initialization. I tried it with argument inits=0.1 and it initialized and fit OK, though it's still slow as it's (100 x 100 choose 2), or roughly 5K, highly constrained parameters for a 100 x 100 Cholesky factor. The fit looks OK, but the first row is suspiciously not moving at all.

Here's the code I used that generates the test data directly in the transformed data block.

transformed data {
  vector[100] zero = rep_vector(0, 100);
  array[500] vector[100] y;
  for (n in 1:500) {
    for (m in 1:100) {
      y[n, m] = normal_rng(0, 1);
    }
  }
}
parameters {
  cholesky_factor_corr[100] Omega;
}
model {
  Omega ~ lkj_corr_cholesky(10);
  y ~ multi_normal_cholesky(zero, Omega);
}

And here's the calls I used in CmdStanPy to fit:

import cmdstanpy as csp
m = csp.CmdStanModel(stan_file='cholesky.stan')
f = f = m.sample(inits=0.1, iter_warmup=500, iter_sampling=500)
f.summary()

It took a couple minutes to run and the output is:

>>> f.summary(sig_figs=2)
                      Mean     MCSE  StdDev     MAD         5%         50%         95%  ESS_bulk  ESS_tail  R_hat
lp__           -25000.0000  1.70000  46.000  45.000 -25000.000 -25000.0000 -25000.0000     720.0    1200.0    1.0
Omega[1,1]          1.0000      NaN   0.000   0.000      1.000      1.0000      1.0000       NaN       NaN    NaN
Omega[1,2]          0.0000      NaN   0.000   0.000      0.000      0.0000      0.0000       NaN       NaN    NaN
Omega[1,3]          0.0000      NaN   0.000   0.000      0.000      0.0000      0.0000       NaN       NaN    NaN
Omega[1,4]          0.0000      NaN   0.000   0.000      0.000      0.0000      0.0000       NaN       NaN    NaN
...                    ...      ...     ...     ...        ...         ...         ...       ...       ...    ...
Omega[100,96]      -0.0650  0.00055   0.032   0.032     -0.120     -0.0640     -0.0100    3600.0    1400.0    1.0
Omega[100,97]       0.0065  0.00064   0.036   0.036     -0.052      0.0071      0.0670    3200.0    1500.0    1.0
Omega[100,98]       0.0460  0.00063   0.035   0.034     -0.014      0.0460      0.1100    3400.0    1400.0    1.0
Omega[100,99]      -0.0600  0.00055   0.034   0.034     -0.120     -0.0610     -0.0028    4100.0    1400.0    1.0
Omega[100,100]      0.8700  0.00032   0.014   0.014      0.850      0.8700      0.9000    1900.0    1400.0    1.0

I'm not sure what's going on with the first row as it's pinned to the correct estimate, but isn't mixing at all.

What we need to do to fix all of this is to use a different transformation for Cholesky factors than we use now.