Lower/upper triangular matrix data types


#1

I have two variables (one data, one parameter) which should be upper or lower triangular matrices, but they aren’t Cholesky factors of anything. Is it ok to use the cholesky_factor_cov data type for this?

edit if so, can a prior be set using to_vector, e.g.

cholesky_factor_cov[5] tri;
to_vector(tri) ~ beta(a, b);

or would one need to loop through each element?

I was worried about declaring a matrix parameter and the other triangle being unconstrained for estimation (and generally slowing things down).

I searched the forum and found only hints, but I wanted to confirm from an implementation perspective this works.


#2

To partly answer my own question, to_vector(tri) ~ normal(0, 1) acts only on the lower triangle & diagonal. The upper triangle appears in output but always zeroed.


#3

That is fine, but a beta prior is not appropriate since cholesky_factor_cov is not constrained to (0,1).

Instead, you should start with a vector that is K + {K \choose 2} with whatever constraints are appropriate and use that to fill a lower triangular matrix.


#4

Ah, indeed the model does not initialize because of this. I had wanted to avoid the index computation of moving between triangular and vector shapes.


#5

That computation is negligible in all respects.


#6

That will make sure that L * L' is positive definite—it won’t give you general triangular matrices.

It’s hard for people who aren’t used to writing nested loops with tricky boundary conditions. For the sake of argument, assuming you want a diagonal and have a square matrix,

functions {
  matrix to_triangular(vector x, int K) {
    // could check rows(y) = K * (K + 1) / 2
    matrix[K, K] y;    
    int pos = 1;
    for (int i = 1; i < K; ++i) {
      for (int j = 1; j <= i; ++j) {
        y[i, j] = y_basis[pos];
        pos += 1;
      }
    }
    return y;
  }
}

And then to use it,

vector[K * (K + 1) / 2] y_basis;
...
matrix[K, K] y = to_triangular(y_basis, K);

Then hopefully you can get away using triangular-specific matrix functiosn to avoid all the unnecessary multiplies by zero and adds.


#7

Thanks! I had had the closed form formula in mind for index i, j to linear triangular, but eventually settled on exactly what you suggested, since I don’t need the random access.


#8

Yup, there’s a formula—it’s unlikely the integer computation is saturated, so either way should work just fine. I just find it easier not to have to think so hard.

Also, I should’ve mentioned I didn’t include the 0s. You should add those if you want to use it with non-triangular matrix functions.