# Initialization of QR reparameterized bounded variables

Hello,

I am revisiting a multiple linear regression problem that I attempted to solve using Stan. I can get Stan to run but it is horrendously slow. To resolve this I have tried to use the QR reparameterization trick as described here.

``````data {
int<lower=1> N;
int<lower=1> K;
vector[N] d;
vector[N] err;
matrix[N,K] W;
}

transformed data {
matrix[N,K] A;
matrix[N,N] H;
vector[N] b;
real sig_b;

matrix[N,K] Q_ast;
matrix[K,K] R_ast;
matrix[K,K] R_ast_inv;

//Define centering matrix
H = diag_matrix(rep_vector(1.0,N)) - rep_matrix(1.0/N,N,N);

//Standardize data
b = H*(d ./ err);
sig_b = sd(b);
b = b/sig_b;

//Center Matrix
for (i in 1:K){
A[:,i] = W[:,i] ./ err;
}
A = (1.0/sig_b)*(H*A);

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

}

parameters {
vector[K] X_ast;
}

transformed parameters {
vector<lower=0>[K] X = R_ast_inv*X_ast;
}

model {
print(R_ast_inv*(X_ast));
X ~ lognormal(30,5);
b ~ normal(Q_ast*X_ast , 1/sig_b);
}

generated quantities {
real F_tot;
F_tot = sum(X);
}
``````

However I run into problems with the initialization since the transformed parameter, `X`, is lower-bounded to be positive. I have tried setting the initial value of `X_ast` to `R_ast*rep_vector(3.5e15,500)` to guaranteed that `X` is positive but for some reason the values do not work or the initialization is ignored and I get the following error

``````.
.
.
data
file = log_tomography_1.data.R
init = log_tomography_1.init.R
random
seed = 3361586922
output
file = output.csv (Default)
diagnostic_file =  (Default)
refresh = 100 (Default)

Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: validate transformed params: X is -1.12091e+18, but must be greater than or equal to 0  (in '/home/lstagner/orbit_tom/tmp/log_tomography.stan' at line 44)

Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: validate transformed params: X is -1.12091e+18, but must be greater than or equal to 0  (in '/home/lstagner/orbit_tom/tmp/log_tomography.stan' at line 44)
.
.
.
``````

Sorry this didn’t get answered before. I don’t understand the model, but why are you trying to assume `X` is positive? In general, that won’t be true for matrix products and putting that constraint in turns Stan into a rejection sampler.

I’m not sure why our normal has `1 / sig_b`—Stan parameterizes on the standard deviation scale for the second parameter and usually `sigma` denotes a standard deviation.

I’m not sure, but you might find more stable performance with `R_ast \ X_ast` as opposed to `inverse(R_ast) * X_ast`. My only concern is that it’ll be slower as you’ll have to factor that matrix each call in the transformed parameters block.

Also, you may want to constrain the scale a bit, beause `lognormal(30, 5)` is has an enormous range of `exp(20)` or so—think `exp(normal(30, 5))`, so the scale gets exponentiated.

If the model doesn’t fit the data well, it can also be slow.

You can define `X` as a local variable in the model block to avoid saving it.