Generating orthogonal latents

Modeling
18

@spinkney Tried to use the Gram-Schmidt with a probabilistic PCA approach. Chains did not mix well and the largest Rhat was above 1.55. (some effective samples were really low). However, it did seem to recover the covariance matrix well in my case (though the data I simulated were dominated by a single factor). I think I might stick with other approaches over this.

functions {
  // Computes the Gram-Schmidt process to orthogonalize a matrix
  matrix gram_schmidt(matrix A) {
    int J = rows(A);
    int K = cols(A);
    matrix[J, K] Q;
    if(J < K) reject("number of rows must be >= the number of columns!");
    for (k in 1:K) {
      vector[J] v = A[:, k];
      for (j in 1:k - 1) {
        vector[J] u = Q[:, j];
        v -= dot_product(v, u) * u;
      }
      Q[:, k] = v / sqrt(dot_self(v));
    }
    return Q;
  }
}
data {
    int<lower=1> T;
    int<lower=2> N;
    vector[N] X[T];
    int<lower=1> N_K;
}
transformed data {
    int N_nonzero_loadings = (N_K * (N_K - 1)) / 2 + N_K * (N - N_K);
    vector[N] zeros = rep_vector(0.0, N);
}
parameters {
    vector<lower=0>[N_K] diag_ele;
    vector[N_nonzero_loadings] off_diag;
    real<lower=0> sd_e;
    vector<lower=0>[N_K] sd_F;
}
transformed parameters {
    matrix[N, N_K] Q; // orthogonal matrix
    cov_matrix[N] Sigma;
    {
        matrix[N, N_K] L;
        L[1, 1] = diag_ele[1];
        {
            int count = 0;
            for (j in 2:N_K) {
                L[j, j] = diag_ele[j];
                for (k in 1:j - 1) {
                    count += 1;
                    L[j, k] = off_diag[count];
                    L[k, j] = 0;
                }
            }
            for (j in (N_K + 1):N) {
                for (k in 1:N_K) {
                    count += 1;
                    L[j, k] = off_diag[count];
                }
            }
        }
        Q = gram_schmidt(L);
    }
    Sigma = tcrossprod(diag_post_multiply(Q, sd_F));
    for (i in 1:N) {
        Sigma[i, i] += square(sd_e);
    }
}
model {
    // prior on L
    off_diag ~ std_normal();
    diag_ele ~ lognormal(0, 1);
    // other priors
    sd_e ~ cauchy(0, 1);
    sd_F ~ cauchy(0, 1);
    // model
    X ~ multi_normal(zeros, Sigma);
}
2 Likes