# Trouble setting up QR

**URL:** <https://discourse.mc-stan.org/t/trouble-setting-up-qr/1612>\
**Category:** General\
**Tags:** stanc\
**Created:** [August 19, 2017, 4:52am UTC](https://discourse.mc-stan.org/t/trouble-setting-up-qr/1612 "2017-08-19T04:52:41Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![bhomass](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bhomass/32/407_2.png) [@bhomass](https://discourse.mc-stan.org/u/bhomass)\
**Post date:** [August 19, 2017, 4:52am UTC](https://discourse.mc-stan.org/t/trouble-setting-up-qr/1612/1 "2017-08-19T04:52:41Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![bhomass](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bhomass/32/407_2.png) [@bhomass](https://discourse.mc-stan.org/u/bhomass)\
**Post date:** [August 23, 2017, 4:37pm UTC](https://discourse.mc-stan.org/t/trouble-setting-up-qr/1612/3 "2017-08-23T16:37:49Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Charles\_Driver](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charles_driver/32/6524_2.png) [@Charles\_Driver](https://discourse.mc-stan.org/u/Charles_Driver)\
**Post date:** [August 23, 2017, 7:51pm UTC](https://discourse.mc-stan.org/t/trouble-setting-up-qr/1612/4 "2017-08-23T19:51:45Z")

</div>

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]