Initialization of QR reparameterized bounded variables

Modeling
1

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[2] 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[2] 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)
.
.
.
2

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.