CmdStan & Stan 2.30 release candidate

Developers
1

I am happy to announce that the latest release candidates of Cmdstan and Stan are now available on Github!

This release cycle brings complex vectors and matrices, 6 new functions and 3 new distributions.

You can find the release candidate for cmdstan here. Instructions for installing are given at the bottom of this post.

Please test the release candidate with your models and report if you experience any problems. We also kindly invite you to test the new features and provide feedback. If you feel some of the new features could be improved or should be changed before the release, please do not hesitate to comment.

The Stan development team appreciates your time and help in making Stan more efficient while maintaining a high level of reliability.

If everything goes according to plan, the 2.30 version will be released in about ten days.

Below are some of the highlights of the new release.

Complex vectors, row_vector and matrices

The Stan language gains three new types: complex_vector, complex_row_vector and complex_matrix with a number of functions that support them:

  • abs,
  • columns_dot_self,
  • cumulative_sum,
  • diag_matrix,
  • diagonal,
  • dot_self,
  • eigenvalues_sym,
  • eigenvectors_sym,
  • fft,
  • fft2,
  • inv_fft,
  • inv_ftt2,
  • prod,
  • reverse,
  • rows_dot_self,
  • singular_values,
  • sum,
  • svd_U,
  • svd_V,
  • symmetrize_from_lower_tri,
  • transpose,
  • trace

Other utility functions:

  • get_imag, get_real,
  • cols, rows, num_elements, size,
  • to_array_1d, to_matrix, to_row_vector, to_vector

fft() returns the Fourier transform of a vector, while fft2() can be used for the 2-dimensional case of the Fast Fourier transform. The inv_fft() and inv_fft2() provide the respective inverse transforms. More details on the FFT functions are temporarily available here. The rest of the listed functions are all existing functions that were extended to support the new complex types.

In addition to these, two new functions were added that return complex vectors/matrices: eigenvectors and eigenvalues. See the next section for details.

Note that real vectors/matrices automatically promote complex vector/matrices, which means that supplying a vector to fft(complex_vector) will work seamlessly, with no need to wrap the vector in a to_complex().

New functions

  • eigenvectors and eigenvalues:

eigenvalues returns the complex-valued vector of eigenvalues of the matrix A. The eigenvalues are repeated according to their algebraic multiplicity, so there are as many
eigenvalues as rows in the matrix. The eigenvalues are not sorted in any particular order. eigenvectors returns the matrix with the complex-valued (column) eigenvectors of
the matrix A in the same order as returned by eigenvalues.

Temporary documentation link.

  • inv_inc_beta

The inverse of the regularized incomplete beta function. The return value x is the value that solves p = inc_beta(alpha, beta, x).

Temporary documentation link.

  • log_determinant_spd

The log of the absolute value of the determinant of the symmetric, positive-definite matrix A. This function is faster when the input matrix is known to be symmetric and positive definite.

Temporary documentation link.

  • norm1, norm2

These functions return the L1 and L2 norms.

Temporary documentation link.

New distributions

  • Multivariate Student-t distribution, Cholesky parameterization: multi_student_t_cholesky_lpdf and multi_student_t_cholesky_rng

Temporary documentation link.

  • Wishart distribution, Cholesky Parameterization: wishart_cholesky_lpdf and wishart_cholesky_rng

Temporary documentation link.

  • Inverse Wishart distribution, Cholesky Parameterization: inv_wishart_cholesky_lpdf and inv_wishart_cholesky_rng

Temporary documentation link.

Deprecations

The fabs() function is now deprecated and will be removed in the 2.33.0 release (estimated to be released in June 2023). All usages of fabs() can be replaced with abs().

Discontinued support for Rtools 3.5 on Windows

The Stan 2.30 release can not be used with the Rtools 3.5 toolchain on Windows. The recommended toolchain on Windows is Rtools40, unless you are using Stan with R 4.2+, in which
case Rtools42 should be used.

Toolchain support on Mac and Linux systems should remain unchanged.

Miscellaneous

  • Downstream libraries were updated: Sundials to 6.1.1 and Boost to 1.78.0
  • ./ and .*, as well as /, and * now support the same input types. The support was not consistent in previous versions.
  • Exposed an overload for cumulative_sum which takes and returns an array of integers.

How to install?

Download the tar.gz file from the link above, extract it and use it the way you use any Cmdstan release. We also have an online Cmdstan guide available at Redirecting…

If you are using cmdstanpy you can install the release candidate using

cmdstanpy.install_cmdstan(version='2.30.0-rc1')

With CmdStanR you can install the release candidate using

cmdstanr::install_cmdstan(version = "2.30.0-rc1", cores = 4)

And then select the RC with

cmdstanr::set_cmdstan_path(file.path(Sys.getenv("HOME"), ".cmdstan", "cmdstan-2.30.0-rc1"))
3 Likes
2

Edit: The multi_student_t has a call to trace_inv_quad_form_ldlt whereas the cholesky version has calls to dot_self and mdivide_left_tri_low. I’m just surprised at how much slower this is. I take it with larger K the ldlt decomposition is going to start impacting the results. Anyway, I think all this just means I need to get the derivatives in!

To my surprise the multi_student_t_cholesky is slower (on this one model). There’s gotta be something simple that’s causing this as the cpp implementation is not all that different between the two.

# multi_student_t_cholesky
> mod_out$profiles()
[[1]]
                      name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t_cholesky 0x1010a4580    6.86995      2.92227      3.94769   122552664       52245960          33384                 6

[[2]]
                      name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t_cholesky 0x104590580      6.653       2.8124       3.8406   117549091       50112865          32021                 6

[[3]]
                      name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t_cholesky 0x10447c580    6.21732      2.66844      3.54888   111587387       47571305          30397                 4

[[4]]
                      name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t_cholesky 0x10254c580     7.0145      2.99608      4.01842   125504148       53504220          34188                 5

# multi_student_t
> mod_out_reg$profiles()
[[1]]
             name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t 0x1026c8580     3.0147      1.36047      1.65422   104223834       17359596          33129                 1

[[2]]
             name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t 0x1026d0580    2.72917      1.23019      1.49898    94424044       15727336          30014                 1

[[3]]
             name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t 0x102eb8580    3.13896      1.41395      1.72501   108798118       18121492          34583                 1

[[4]]
             name   thread_id total_time forward_time reverse_time chain_stack no_chain_stack autodiff_calls no_autodiff_calls
1 multi_student_t 0x1027b8580    2.98872      1.34738      1.64134   103402728       17222832          32868                 1

The Stan model is from Simon Jackman’s Bayesian Model Examples in Stan. Here’s the code I’m testing.

data {
  // multivariate outcome
  int N;
  int K;
  array[N] vector[K] y;
  // covariates
  int P;
  array[N] vector[P] X;
  // prior
  vector[K] alpha_loc;
  vector[K] alpha_scale;
  array[K] vector[P] beta_loc;
  array[K] vector[P] beta_scale;
  real Sigma_corr_shape;
  real Sigma_scale_scale;
  int cholesky_lik;
}
parameters {
  // regression intercept
  vector[K] alpha;
  // regression coefficients
  array[K] vector[P] beta;
  // Cholesky factor of the correlation matrix
  cholesky_factor_corr[K] Sigma_corr_L;
  vector[K] Sigma_scale;
  // student-T degrees of freedom
  real<lower=0> nu;
}
transformed parameters {
  array[N] vector[K] mu;
 
  for (i in 1 : N) {
    for (k in 1 : K) {
      mu[i, k] = alpha[k] + dot_product(X[i], beta[k]);
    }
  }
}
model {

  for (k in 1:K) {
    alpha[k] ~ normal(alpha_loc[k], alpha_scale[k]);
    beta[k] ~ normal(beta_loc[k], beta_scale[k]);
  }
  
  nu ~ gamma(2, 0.1);
  Sigma_scale ~ cauchy(0., Sigma_scale_scale);
  Sigma_corr_L ~ lkj_corr_cholesky(Sigma_corr_shape);
  if (cholesky_lik == 1) {
    matrix[K, K] L = diag_pre_multiply(Sigma_scale, Sigma_corr_L);
    profile("multi_student_t_cholesky") {
    y ~ multi_student_t_cholesky(nu, mu, L);
    }
  } else {
    matrix[K, K] Sigma = multiply_lower_tri_self_transpose(diag_pre_multiply(Sigma_scale, Sigma_corr_L));
    profile("multi_student_t") {
    y ~ multi_student_t(nu, mu, Sigma);
    }
  }
}
3

I think so yes.

4

Simulating some data and I’m finding that sampling is extremely slow with dimensions > 25 on the non-cholesky version due to numerical issues with LDLt. So, if you have relatively high dimensional student t, the Cholesky version is the lpdf for you!

3 Likes
5

I study matrix models to study population dynamics of species that are structured in age/stage classes. The growth characteristics in these models is governed by the dominant eigenvalue and corresponding vectors. Since the eigenvalues are not sorted in any order, is there a way to identify the dominant eigenvalue and vector

6

There’s the max function and there’s 5.5 Sorting functions | Stan Functions Reference.

1 Like
7

The behavior of our function is the same as the underlying Eigen library’s implementation, which doesn’t order the results. For complex numbers there isn’t really a natural order you can impose

8

Just compiled cmdstan-2.30-rc1 and noticed 2 small issues:

  1. For a very long time I think I could get away with below erroneous model:
bernoulli_model = "
data {
  int<lower=1> N;
  real y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

This now reports:

Ill-typed arguments supplied to function 'bernoulli':
(array[] real, real)
Available signatures:
(int, row_vector) => real
  The first argument must be int but got array[] real
(int, vector) => real
  The first argument must be int but got array[] real
(int, array[] real) => real
  The first argument must be int but got array[] real
(int, real) => real
  The first argument must be int but got array[] real
(array[] int, row_vector) => real
  The first argument must be array[] int but got array[] real
(Additional signatures omitted)

which is correct but just wanted to confirm this is an intended change. My JSON input is {"N":10,"y":[0,1,0,1,0,0,0,0,0,1]}.

Replacing real by int resolves this.

As usually, the .csv file reports cmdstan v2.29.2. But that is usually updated pretty late.

All other tests seem to run fine.

Best, Rob

Edit: New keyboard …

1 Like
9

Thanks @Rob_J_Goedman. I quickly tried a few recent versions (back to 2.27) and the model in 1. failed to compile in all of them. It’s quite possible that it compiled in some previous version, but I would consider that previous behaviour a bug, as that would be doing truncation from reals to integers.

Yeah, we update that one prior to the actual release.

1 Like
10

Thanks @rok_cesnovar , you’re right. For quite some time I’ve had array[N] int<lower=0, upper=1> y; in the model and somehow dropped the int<lower=0, upper=1> part when testing 2.30.1 and must have reentered as real y[N]. Should have gone back to the Github history.

11

I updated httpstan to use the release candidates. Everything seems to be working.

2 Likes
12

Thanks! If things keep trending as they are, the release will happen on Thursday.

1 Like
13

have the docs been frozen as well?
I have a some updates to the installation instructions for CmdStan

14

Docs dont fall under the freeze, so fire those updates away.

15

The path should have .cmdstan not .cmdstanr (see Change default directory to .cmdstan by jgabry · Pull Request #520 · stan-dev/cmdstanr · GitHub)

2 Likes
16

Fixed, thanks! I copied that section from an old release candidate announcement before the PR you linked.

Also an update, the official release should be out today.

2 Likes
17

Thanks for this tip. It works great. I really want to thank the development team. With complex vector and eigen function support we can explore a really wide range of spatial models that were not possible (alleast not in computationally scalable way)

2 Likes