Correlation matrix with positively constrained off-diagonals

Modeling
1

Hi, everyone.

I would like to estimate a correlation matrix with it off-diagonal values are positively constrained. So far, I tried

// snipet 
parameter{
  corr_matrix<lower=0>[N] Omega; 
}
model{
  Omega ~ lkj_corr(1);
}

but this returned a compilation error. Then, I tried

// snipet 
parameter{
  matrix<lower=0>[N,N] Omega; 
}
model{
  Omega ~ lkj_corr(1);
}

which compiled but returned an initialization error. To remove an initialization error, I provided good initial values, but I got lots of sampling errors.

By searching this forum, I found Dealing with constraints on correlations. Is this the only way to set up a positively constrained correlation matrix?

Thank you for your help and kindness. Any help is greatly appreciated.

2 Likes
2

Hi,
I think this specific case should be possible, but to do this, you will need to implement the constraining transform + Jacobian adjustment yourself. I’ll take the default transformation as a starting point 10.9 Correlation matrices | Stan Reference Manual

Conjecture: the cholesky decomposition of a correlation matrix with all positive elements has all elements positive. Trivially all positive elements of cholesky decomposition means all positive matrix. Not sure about the converse, can’t quickly prove it. If the converse does not hold, the construction below will be overly restrictive.

Note that in the construction of w (the cholesky decomposition) in the Stan manual link, the diagonals are ensured to be positive (to make the decomposition valid) and the sign of off-diagonal w_{i,j} equals the sign of z_{i,j}. So we need to ensure z_{i,j} = \tanh y_k > 0 which is equivalent to y_k > 0. So if you can reimplement the constraint and Jacobian yourself, but start from parameters constrained to be positive, you will get a correlation matrix with all entries positive.

EDIT: The C++ code for the constraints Stan does could serve as a useful starting point,
the code for the cholesky_factor_corr type is a bit more readable: math/cholesky_corr_constrain.hpp at develop · stan-dev/math · GitHub while the version for corr_matrix decomposes into multiple function calls, see: math/corr_matrix_constrain.hpp at develop · stan-dev/math · GitHub

In both cases lp accumulates the required log jacobian adjustment to be added to target.

Hope that helps at least a bit.

2 Likes
3

I had most of this code already, just constraining tanh(x) to positive.

functions {
  vector corr_constrain_lp(vector x) {
    int N = num_elements(x);
    vector[N] tanh_x = tanh(x);
    
    target += sum(log1m(square(tanh_x)));
    
    return tanh_x;
  }

  matrix cholesky_corr_constrain_lp (vector y, int K) {
    int k_choose_2 = (K * (K - 1)) / 2;
    vector[k_choose_2] z = corr_constrain_lp(y);
    matrix[K, K] x = rep_matrix(0, K, K);
    int iter = 0;
    
    x[1, 1] = 1;

    for (i in 2:K) {
      real sum_sqs;
      iter += 1;
      x[i, 1] = z[iter];
      sum_sqs = square(x[i, 1]);
     if (i > 2) {
      for (j in 2:i - 1) {
        iter += 1;
        target += 0.5 * log1m(sum_sqs);
        x[i, j] = z[iter] * sqrt(1.0 - sum_sqs);
        sum_sqs += square(x[i, j]);
      }
     }
    x[i, i] = sqrt(1.0 - sum_sqs);
  }
  return x;
  }
}
data {
  int K;
  real eta;
}
transformed data {
  int k_choose_2 = (K * (K - 1)) / 2;
}
parameters {
  vector<lower=0>[k_choose_2] y_raw;
}
transformed parameters {
  matrix[K, K] y = cholesky_corr_constrain_lp(y_raw, K);
}
model {
  y ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(y);
}

Run with

library(cmdstanr)

fp <- file.path("../positive_corr.stan")
mod <- cmdstan_model(fp)

K <- 10
mod_out <- mod$sample(
  data = list(K = K,
              eta = 100),
  chains = 2,
  adapt_delta = 0.8,
  parallel_chains = 2,
  iter_warmup = 200,
  iter_sampling = 200
)

matrix(mod_out$summary("Omega")$mean, K, K)
3 Likes
4

@martinmodrak and @spinkney
Thank you for your detailed explanations on this issue.

@martinmodrak I will read the suggested chapter to understand the math behind the correlation matrices.
@spinkney Thank you very much for providing your code. It is greatly appreciated and helped me a lot.