defining a covariance matrix from a correlation matrix

[terrance savitsky]

Hi,  I finally got one of my monster multivariate models lofted in stan.  (see attached script).   I constructed a set of covariance matrices from an associated set of correlation matrices, the latter under uniform priors.   The manner in which I was able to get this to run was to declare the covariance matrices as unstructured matrix objects, rather than cov_matrix objects.  I had been receiving error messages under the cov_matrix declaration of the first-used matrix that the sampler failed for violation of symmetry.   That said, the error occurred so rapidly that I wonder if there it never got to the sampler.   These covariance matrices are, in turn, used to generate multivariate gaussian samples.   The cholesky decomposition of each covariance matrix is rendered for that purpose.   So the chains may yet crash under an SPD error, but I wanted to let you know of this event.  

Lastly, I think previous run attempts crashed R because I input a NULL data element.

--
Thank you, Terrance Savitsky

[Ben Goodrich]

Can you give us the script / data that produces the symmetry error? We need to get that fixed and it may be the case that the check is too strict.

Ben

[terrance savitsky]

i will send, tomorrow.  thanks for asking.  the only difference is that i would replace the "matrix" declarations for those matrices (all labeled with capital "S") formed from correlation matrices.  something that suggested to me that changing from "cov_matrix" to "matrix" might be a good idea is the example in the stan users manual that builds a covariance matrix from a mixture of known data elements and unknown parameters.  that matrix was declared as an unrestricted "matrix".  something odd, though.  before changing the declaration i had used matt's suggestion of replacing covariance S = 0.5 * ( S + S') to enforce symmetry and that didn't work.

--
 
 

[terrance savitsky]

Ben,  I attach a simplified model that renders the error.  In particular, if the cov_matrix[D] S_t is changed to matrix[D,D] S_t, all is well.

Terrance.

On Sat, Oct 27, 2012 at 8:06 PM, Ben Goodrich <goodri...@gmail.com> wrote:

--
 
 

[terrance savitsky]

The attached monster model ran under stan - much faster than i'd thought, though I wasn't able to check for convergence because I received the following error in (64-bit) R: Error: cannot allocate vector of size 16 Kb.    There are some things I may try in R, but I write to see if anyone has experience managing memory allocations to arrays in R.   The model is generally a memory hog (what with so many covariance matrices) and consumes most of the available 16GB of RAM on my machine during the run.  I know I may improve the speed from the attached script by directly augmenting lp__ in a fashion that composes the quadratic product without performing inverses of covariance matrices.   Since I still need the cholesky decomposition and the inverse of that, I'm not clear that this approach would conserve RAM.   Also, even if I could conserve RAM, I bet I'd still hit the vector allocation size limit.

I also received a few warning statements (during 2000 iterations) that said this:

Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal (really) issue (and please ignore that what comes next might say 'error'): Error in function stan::prob::ordered_logistic(N4stan5agrad3varE): Cut points parameter is 1.#INF:0, but must be finite!

If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.

 If the problem does not persist, the resulting samples will still be drawn from the psoterior.

Cutpoints are generally notorious highly intracorrelated under ordered logistic models, but I do something that I want to check if you think it is a problem.   I fix the first element of a vector of ordered cutpoints to 0 to achieve hard identifiability.   the other elements are random.

Thanks, Terrance.

[Ben Goodrich]

On Tuesday, October 30, 2012 3:16:46 PM UTC-4, tds151 wrote:

The model is generally a memory hog (what with so many covariance matrices) and consumes most of the available 16GB of RAM on my machine during the run.  I know I may improve the speed from the attached script by directly augmenting lp__ in a fashion that composes the quadratic product without performing inverses of covariance matrices.   Since I still need the cholesky decomposition and the inverse of that, I'm not clear that this approach would conserve RAM.  

In R, you can pass a pars argument to stan() that names the posterior quantities you would like to keep. Using that would give R more headroom to calculate afterward.

But the main place where I think RAM is consumed unnecessarily is with the covariance matrices in the transformed parameters block:

        /* covariance matrices */
        matrix[D,D] S_t;
        matrix[Ds,Ds] S_s;
        matrix[Do,Do] S_v;
        matrix[Do,Do] S_r;
        matrix[D1,D1] S_z1;
        matrix[D2,D2] S_z2;
        matrix[D3,D3] S_z3;
        matrix[D,D] L_t;
        matrix[Ds,Ds] L_s;
        matrix[Do,Do] L_v;
        matrix[Do,Do] L_r;
        matrix[D1,D1] L_z1;
        matrix[D2,D2] L_z2;
        matrix[D3,D3] L_z3;
...
        /* Generate covariance matrices for multivariate random effects from associated correlation matrices and scale parameters */
        for( d in 1:D ) for( e in 1:D ) S_t[d,e]        <- R_t[d,e] * sqrt(s2_t[d]) * sqrt(s2_t[e]);
        for( d in 1:Ds ) for( e in 1:Ds ) S_s[d,e]      <- R_s[d,e] * sqrt(s2_s[d]) * sqrt(s2_s[e]);
        for( d in 1:Do ) for( e in 1:Do ) S_v[d,e]      <- R_v[d,e] * sqrt(s2_v[d]) * sqrt(s2_v[e]);
        for( d in 1:Do ) for( e in 1:Do ) S_r[d,e]      <- R_r[d,e] * sqrt(s2_r[d]) * sqrt(s2_r[e]);
        for( d in 1:D1 ) for( e in 1:D1 ) S_z1[d,e]     <- R_z1[d,e] * sqrt(s2_z1[d]) * sqrt(s2_z1[e]);
        for( d in 1:D2 ) for( e in 1:D2 ) S_z2[d,e]     <- R_z2[d,e] * sqrt(s2_z2[d]) * sqrt(s2_z2[e]);
        for( d in 1:D3 ) for( e in 1:D3 ) S_z3[d,e]     <- R_z3[d,e] * sqrt(s2_z3[d]) * sqrt(s2_z3[e]);

        /* render cholesky decompositions of covariance matrices */
        L_t             <- cholesky_decompose(S_t);
        L_s             <- cholesky_decompose(S_s);
        L_v             <- cholesky_decompose(S_v);
        L_r             <- cholesky_decompose(S_r);
        L_z1            <- cholesky_decompose(S_z1);
        L_z2            <- cholesky_decompose(S_z2);
        L_z3            <- cholesky_decompose(S_z3);

If you think of a covariance matrix as a correlation matrix multiplied from both sides by a diagonal matrix of standard deviations, then it is not necessary to first calculate the covariance matrix in order to then calculate its Cholesky factor. You could just make the covariance matrix temporary by doing

L_ <- cholesky_decompose(diag_matrix(s_) * R_ * diag_matrix(s_));

However, doing so entails that you define each s_ as a vector (rather than real) of standard deviations (rather than variances) and thus you need to be using v1.0.3 in order to put the non-negativity constraint on a vector.

But there is an even better way. The Cholesky factor of a covariance matrix is the Cholesky factor of the correlation matrix pre-multiplied by the diagonal matrix of standard deviations. So, you could do something like

L_ <- diag_matrix(s_) * cholesky_decompose(R_);

This way entails less burden on the auto-differentiator.
 

But you can also make the Cholesky factors temporary. Rather than looping as in

        /* Re-scaled and located multivariate random effects */
        for( t in 1:T ) Delta[t]        <- (L_t * Delta_norm[t]')';
        for( s in 1:S ) Phi[s]          <- (L_s * Phi_norm[s]')';
        for( v in 1:V ) Theta[v]        <- (L_v * Theta_norm[v]')';
        for( r in 1:Rr ) Beta[r]        <- (L_r * Beta_norm[r]' + mean_Beta)';
        for( z in 1:VR1 ) Zeta1[z]      <- (L_z1 * Zeta1_norm[z]')';
        for( z in 1:VR2 ) Zeta2[z]      <- (L_z2 * Zeta2_norm[z]')';
        for( z in 1:VR3 ) Zeta3[z]      <- (L_z3 * Zeta3_norm[z]')';

if you define everything as matrices, you could do something like

Delta <- Delta_norm * transpose(cholesky_decompose(R_t)) * diag_matrix(s_t);

and so on. For Beta, it is a little annoying because Stan currently does not support assigning to a row or column of a matrix. So, you could deal with Beta separately or make another matrix whose columns are all mean_Beta and tack that onto Beta.

Ben

[Ben Goodrich]

On Tuesday, October 30, 2012 4:17:35 PM UTC-4, Ben Goodrich wrote:

On Tuesday, October 30, 2012 3:16:46 PM UTC-4, tds151 wrote:

The model is generally a memory hog (what with so many covariance matrices) and consumes most of the available 16GB of RAM on my machine during the run.  I know I may improve the speed from the attached script by directly augmenting lp__ in a fashion that composes the quadratic product without performing inverses of covariance matrices.   Since I still need the cholesky decomposition and the inverse of that, I'm not clear that this approach would conserve RAM.  

In R, you can pass a pars argument to stan() that names the posterior quantities you would like to keep. Using that would give R more headroom to calculate afterward.

But the main place where I think RAM is consumed unnecessarily is with the covariance matrices in the transformed parameters block:

A smaller waste of RAM is with the cutpoints thing in the transformed parameters {} block

/* truncation points for ordered_logistic on y1[i,d], y2[i,d] and y3[i,d]  */
// construct ordered cut points, gamma[d,c], where gamma[d,1] = 0 for identifiability,
// else gamma[d,c] = gamma_u[d,(c-1)] for c >= 2

You should do that stuff inside the model {} block so that gamma is not stored, only gamma_u.

Ben

[Ben Goodrich]

Is there a particular seed that produces the error? I've run it for 40 -- 60 iterations a few times without hitting it yet.

Ben

[Bob Carpenter]

On 10/30/12 5:08 PM, Ben Goodrich wrote:
> On Tuesday, October 30, 2012 4:17:35 PM UTC-4, Ben Goodrich wrote:
>
> On Tuesday, October 30, 2012 3:16:46 PM UTC-4, tds151 wrote:
>
> The model is generally a memory hog (what with so many covariance matrices) and consumes most of the available
> 16GB of RAM on my machine during the run. I know I may improve the speed from the attached script by directly
> augmenting lp__ in a fashion that composes the quadratic product without performing inverses of covariance
> matrices. Since I still need the cholesky decomposition and the inverse of that, I'm not clear that this
> approach would conserve RAM.
>
>
> In R, you can pass a pars argument to stan() that names the posterior quantities you would like to keep. Using that
> would give R more headroom to calculate afterward.
>
> But the main place where I think RAM is consumed unnecessarily is with the covariance matrices in the transformed
> parameters block:

These are two different kinds of memory requirements.

The transformed params (and anything else involved in the
model) is a transient memory demand. As soon as you get
the gradient, the memory's recovered.

The parameters that are saved in memory (in RStan, not
in Stan itself, which streams the samples to file) are
a constant, growing memory demand.

I actually don't think covar matrix transforms cause too
much of a hit, but it could be improved.

Anything that involves solving is going to be way worse,
because we build up the full expression tree.

The best thing we can do throughout Stan to improve
efficiency is to partially evaluate auto-dif. We could
do more of this in all of our matrix ops. Right now,
there are some efficient products and dot products, and
a few functions like sum(), and quite few auto-dif efficient
vectorized distributions.

- Bob

[Bob Carpenter]

On 10/30/12 3:16 PM, terrance savitsky wrote:
> The attached monster model ran under stan - much faster than i'd thought, though I wasn't able to check for convergence
> because I received the following error in (64-bit) R: Error: cannot allocate vector of size 16 Kb.

Sounds like R's out of memory.

By dumping your data, you can run Stan from
the command line, which will stream out samples
rather than trying to keep them in memory.

I'd be curious as to hwo much of the memory is R
and how much is Stan's model execution.

You can try more thinning within R.

The profile should be gradually increasing memory
due to saving samples, plus spikes as the model
is evaluated every leapfrog step.

That should use as much memory the first time
being executed as any other time, so you can
probably watch whatever OS your using's system
monitor to watch memory.

> There are some
> things I may try in R, but I write to see if anyone has experience managing memory allocations to arrays in R. The
> model is generally a memory hog (what with so many covariance matrices) and consumes most of the available 16GB of RAM
> on my machine during the run. I know I may improve the speed from the attached script by directly augmenting lp__ in a
> fashion that composes the quadratic product without performing inverses of covariance matrices. Since I still need the
> cholesky decomposition and the inverse of that, I'm not clear that this approach would conserve RAM. Also, even if I
> could conserve RAM, I bet I'd still hit the vector allocation size limit.
>
> I also received a few warning statements (during 2000 iterations) that said this:
>
> Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal
> (really) issue (and please ignore that what comes next might say 'error'): Error in function
> stan::prob::ordered_logistic(N4stan5agrad3varE): Cut points parameter is 1.#INF:0, but must be finite!
> If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.
> If the problem does not persist, the resulting samples will still be drawn from the psoterior.

I was hoping that message would be less scary than
"(Non-fatal Warning)" followed by "Error in function..."

I'll fix the spelling of "posterior".

Any other suggestions for how to say this given that the "Error in function..
must be finite!" part is fixed.

> Cutpoints are generally notorious highly intracorrelated under ordered logistic models, but I do something that I want
> to check if you think it is a problem. I fix the first element of a vector of ordered cutpoints to 0 to achieve hard
> identifiability. the other elements are random.

It's not so much an identifiability problem as a problem
with the discretization in the Hamiltonian simulation.

For example, if you have a small value for a scale parameter
sigma bounded below at zero, and a gradient toward zero and
a large step size, it's possible to overshoot zero in one
leapfrog step inside the Hamiltonian simulation. Lower step
sizes can always be found (up to arithmetic issues), but we
can't get around the fact that we're using a first order
(derivative only) discretization of the Hamiltonian.

- Bob

[terrance savitsky]

it seems like R is holding on to the memory allocated to the temporary objects from the auto diff as it proceeds to write out the sampled objects.  i'm guessing garbage collection / unprotecting objects in rcpp doesn't take place until after the c++ objects are returned to the R layer??  all to say, returning similarly-sized objects from JAGS doesn't hit the RAM ceiling within R.

something that would be useful is a user option to pre-select the parameters that will be returned / written to file.   such a feature may be useful for managing memory for running stan within R.   

- Bob

--

[Jiqiang Guo]

Just to highlight that Ben has mentioned the following.

[terrance savitsky]

whoops.  i see the 'pars' option for stan function.  thanks for reiterating.

--
 
 

[terrance savitsky]

a related question. i don't see a stan function option to discard the burnin.  is that done, by default?

On Thu, Nov 1, 2012 at 6:04 PM, Jiqiang Guo <guo...@gmail.com> wrote:

--
 
 

[Jiqiang Guo]

There is no such option now.  So all warmup samples are saved.  

--

Jiqiang

--
 
 

[terrance savitsky]

Hi, I ran my model, again, after some help from Ben and Bob such that less memory is drawn by the autodiff and I employed the 'pars' option to restrict returned parameters.   I still ran out of memory (using R) as stan attempted to write out objects to file.  I'm able to store a lot more parameters for the same number of stored iterations in JAGS without a problem.  Since I assume both JAGS and stan store the sampled parameters in similar (S4) objects, I'm wondering if the storing the last state of the sampler requires also storing the last gradient computation and if this may be an issue.  Also, this time when the sampler finished, I had a message indicating the iterations were completed as well as the out of memory message.  Yet, the object existed and when I attempted to extract samples, it told me that there were no samples.   I did not receive any indication that the sampler encountered issues, however.

- Bob

--

[Bob Carpenter]

On 11/12/12 4:29 PM, terrance savitsky wrote:
> Hi, I ran my model, again, after some help from Ben and Bob such that less memory is drawn by the autodiff and I
> employed the 'pars' option to restrict returned parameters. I still ran out of memory (using R) as stan attempted to
> write out objects to file. I'm able to store a lot more parameters for the same number of stored iterations in JAGS
> without a problem. Since I assume both JAGS and stan store the sampled parameters in similar (S4) objects, I'm
> wondering if the storing the last state of the sampler requires also storing the last gradient computation and if this
> may be an issue. Also, this time when the sampler finished, I had a message indicating the iterations were completed as
> well as the out of memory message. Yet, the object existed and when I attempted to extract samples, it told me that
> there were no samples. I did not receive any indication that the sampler encountered issues, however.

Thanks for the report.

You can always write the models to file rather than
storing them in R using the command-line version of Stan.
That just streams the draws to a file, so there's no overhead
at all for number of samples.

I don't know if there's an easy way to read the results
back into a Stan object.

I don't know how we're storing data in R (but I'm guessing
we're about to find out). The original plan was
to store the data in a compact C++ object and access it
through Rcpp, but that proved to be a huge hassle for saving
and restoring across sessions. I still have some hope we
could write custom serialization but I'm not really enough of
an R programmer to do it myself and it didn't sound easy after
a query to the Rcpp list.

Because users really wanted to be able to save sessions,
we switched the code to store the data in R rather than in
C (hopefully it's not being duplicated -- if so, we could
probably cut out a lot of fat there).

- Bob

[Andrew Gelman]

We really shouldn't be running out of memory in R for models that can be fit in Jags!
A

> --
>

[Jiqiang Guo]

The gradient is not stored in the fitted object.  We do store some other information such as the tree_depth, stepsize, lp__.  I do not think it will be an issue in general.  My guess about that there is no sample in the fitted object for your case if because of the error; in the whole sampling process, whenever error is encountered, an empty object is returned instead. 

Another option is that set pars to be as smallest as possible (currently, it does not allow to be empty) and specify sample_file argument for stan(), so that the samples are written to text files as using Stan from the command line. 

--

Jiqiang

--
 
 

[Bob Carpenter]

On 11/12/12 5:14 PM, Jiqiang Guo wrote:
> The gradient is not stored in the fitted object.

The only issue with the gradient is the space needed
for the expression tree gets stored in memory in order
to compute gradients. It's transient memory that's
allocated and used per sample and then freed.

> We do store some other information such as the tree_depth, stepsize,
> lp__. I do not think it will be an issue in general. My guess about that there is no sample in the fitted object for
> your case if because of the error; in the whole sampling process, whenever error is encountered, an empty object is
> returned instead.

I think this makes sense to keep consistency.

If there's a way to keep the "samples so far" in
a consistent state that might help, but I don't think it
should be a high priority.

> Another option is that set pars to be as smallest as possible (currently, it does not allow to be empty) and specify
> sample_file argument for stan(), so that the samples are written to text files as using Stan from the command line.

Is there an easy way to read the results back into
a stanfit object? I looked at the rstan doc index
and didn't see anything. If we could do that, then we
could use sflist2stanfit to put multiple samples together.

Does anyone recall how JAGS reads data? I could always
download r2jags and look at the source.

- Bob

[terrance savitsky]

Thanks, Jiqiang.  I'll give this a try (writing to a file instead of storing in an R session) as I'm also wondering if there's something I'm doing on my end that is a problem.   Getting back a set of reasonable samples will let me know.  I will read the saved-to-file objects back into R.   If there is anything I need to know or do (for example using the extract function), I'd appreciate your letting me know.

Thanks, Terrance

[Bob Carpenter]

On 11/12/12 5:37 PM, Jiqiang Guo wrote:
> [to stan-dev]
>
> 1. reading csv files of samples in RStan
>
> Currently, no functionality in RStan to do that. I guess I should work on it so maybe in the release after 1.0.4, we
> can have it.

This would be really useful for scalability in R or even
just for running from Stan command line and then
reading into R for later analysis.

> 2. keeping `samples so far'
>
> I do not know now. I have never encounter an error in R of `out of memory`.

Lucky you.

If you use 32-bit R and try to run JAGS models,
you'll run into it fast enough. I rewrote some models in
Stan that took up way less memory than their JAGS counterparts
in R.

Any model that calls iterative solvers on matrices may
take up a lot of memory during sampling the way things are
implemented now. The issue is that we build up the entire
expression tree over the iterations before propagating the
chain rule backwards for calculating gradients. We plan to
make this more efficient by extracting the double values,
calling the iterative algorithsm on those, then partially
evaluating the gradients using all the rules in The Matrix
Cookbook and Magnus and Neudecker's book.

> I do not know how R or OS deals with this
> case. And in the case if it is R calling a C++ function that is running out of memory, I am not sure we can do
> something in R then.

It's going to be handled by process. When R calls the
C function, it's all in the same process, and so sharing
all the memory from the OSes point of view. Thus if you run
R in 32-bit mode, the OS is going to put a low ceiling on
how much memory can be used. A C function called from R
might not be well behaved with respect to R -- it can call its
own malloc() and play havoc with R's memory model.

I don't know how well Rcpp deals with these issues and how
much is up to the user.

You can try something like object.size() in R. memory.profile()
also looks useful.

There may also be extra objects hanging around in the R environment
taking up space that have nothing to do with Stan. R uses
garbage collection, but doesn't appear to do it very well. So
one probably shouldn't use long-running R sessions even if you get
rid of local variables and call gc() to try to clean things up.

> In short, it is difficult. Unless this happens a lot, I would not put the investigation on top
> priority.

We should probably test the scalability to make sure we're
not doing something sub-optimally. If it's a deep issue with
R, I agree that it's not a priority for us to try to fix it!

- Bob