Trouble setting up QR

General
1

I am implementing the QR reparametrization of section 9.2. However, my model calls for an array of matrices.

data {
    int<lower=0> N;
 	int T; // number of days in campaign
	int K; // number of predictors
	
    matrix[N, K] x[T];
    vector[N] y[T];
}

So for the transformed data block I want to set up T Q_ast matrices with the modification that each QR parameter is a array of dimension T.

	for (t in 1:T){
		Q_ast[t] = qr_Q(z[t])[, 1:K] * sqrt(N - 1);
	R_ast[t] = qr_R(z[t])[1:K, ] / sqrt(N - 1);
	R_ast_inverse[t] = inverse(R_ast[t]);
	}

However, I am not sure how I would recover beta using the formula from the manual.

beta = R_ast_inverse * theta

since in my case both R_ast_inverse and theta are arrays of dimension T, I will end up with T sets of betas. Do I take the average?

1 Like
3

I am trying the QR technique out. I ran into non-invertible matrix problem

 Q_ast = qr_Q(x)[, 1:K] * sqrt(N - 1);
  R_ast = qr_R(x)[1:K, ] / sqrt(N - 1);
  R_ast_inverse = inverse(R_ast);

But because my data contains most 0’s, the inverse calculation fails

Q_ast = [[0.154673,0,-0.0379848,-0.121529,0.00411132],[-0,3,-0,0,0],[-0,0,3.53977e-06,0.0282982,-0.324708],[0.12117,0,-0.0297572,2.99507,0.0032208],[-0,0,0.324708,0.00306288,2.96485],[0.719883,0,2.89347,-0.000357266,-0.313177],[-0,0,1.13273e-06,1.06847e-08,-1.22602e-07],[5.82244e-05,0,0.10423,0.000980936,-0.0112814],[2.45722,0,-0.603449,-0.100074,0.0653149],[1.55087,0,-0.380866,-0.0631614,0.0412234]]

R_ast = [[0.648872,-0,0.0799872,-0,-0],[0,0,0,0,9.26125e-36],[-0,-0,0.325705,-0,-0],[0,0,0,0,0],[0,0,0,0,0]]

R_ast_inverse = [[nan,nan,nan,nan,nan],[nan,nan,nan,nan,nan],[nan,nan,nan,nan,nan],[nan,nan,nan,nan,nan],[nan,nan,nan,nan,inf]]

Does any one have tricks to get around this problem in order to still leverage QR?

4

I wrote a lower triangular inversion function recently, perhaps that helps? edit: no probably not, but will leave it here anyway!

matrix inverse_lowtri(matrix l, int unitdiag){ //inverse of lower triangular matrix
  matrix[rows(l),rows(l)] b ;
  matrix[rows(l),rows(l)] x;
  x = diag_matrix(rep_vector(1,rows(l)));
  b = x;

    for(j in 1:cols(x)){
      if(unitdiag != 1) x[j,j] = b[j,j] / l[1,1];
      
     if(rows(x)>1){
      for(m in 2:rows(x)){
      x[m,j]  =  (b[m,j] - (l[m,1:(m-1)] * x[1:(m-1),j]) ) / l[m,m];
      }
     }
    }
  return x;
}

[edit: added triple back ticks before and after code to format it]