Lower/upper triangular matrix data types

Algorithms
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.

1 Like
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.

9

Just came across this. Is there any way to create a triangular matrix from a vector of its values using a pointer?
Using map reduce I need to pass a covariance matrix to the function as a vector, but I don’t want so much copying to slow me down.

Also, to stay in the example:

cholesky_factor_cov[5] tri;
to_vector(tri)

It seems that to_vector(tri) includes '0’s for the 0 elements…

10

Only if you decompose and recompose it yourself.

We really need to add triangular data type support.

2 Likes
11

that would be awesome!

12

Hi,
I am wondering if you have an idea how to declare an upper triangular matrix with simplex rows please?
Maybe this was already solved. I saw the discussion dates back to 2020😊