Spatial distributions?

[John]

Has anyone tried to implement in Stan the kind of spatial distributions found in Geobugs, such as proper CAR?  Would Stan likely be much faster for these models?

Thanks, and good luck with the project,

John

[Ben Goodrich]

 On Wednesday, September 19, 2012 6:21:32 PM UTC-4, John wrote:

Has anyone tried to implement in Stan the kind of spatial distributions found in Geobugs, such as proper CAR?  Would Stan likely be much faster for these models?

A coauthor and I have been using one such model. The important thing to keep in mind is that Stan doesn't do conditional distributions, inclusive of the conditional autoregressive (CAR) distribution. Stan does joint distributions. Some flavors of BUGS implement CAR because it admits Gibbs sampling. Stan does variants of Hamiltonian Monte Carlo. The implications of this difference are that 1) you can use whatever priors and distributional assumptions reflect your beliefs rather than the priors and distributional assumptions that lead to known full-conditional distributions and 2) you need to express your priors in terms of the joint distribution rather than the conditional distributions. In our case, it turned out that the normal CAR thing that we would have used if we were using BUGS corresponded to an improper multivariate normal joint distribution (improper because the covariance matrix was singular). And that eventually ran reasonably fast. Your case might well imply some other multivariate normal (or maybe t) with a nonsingular but structured covariance matrix. At the moment, the C++ code in Stan is a bit thin in this area. The multivariate normal and t distributions exist, but are not optimized for structured covariance (or precision) matrix. So, in order to get optimal speed, you might have to do some parts within the .stan file with

lp__ <- lp__ + whatever;

syntax. See the chapter called "Assignments" in the reference manual or this thread

https://groups.google.com/d/topic/stan-users/HcvYaDu71_Y/discussion

Ben

[terrance savitsky]

To your general comment on structured covariance matrices of a multivariate Gaussian, is your point that if there are an underlying set of fewer parameters generates some covariance matrix, then the gradient may be more efficiently processed by stating the multivate Gaussian prior in terms of those parameters, rather than the resulting covariance matrix?   For example, what if I had a factor analytic structure where Sigma = Lambda * Lambda' + diag(sigma)?  Same idea, again, if I specified a proper CAR covariance structure?

--
Thank you, Terrance Savitsky

[terrance savitsky]

Response from Ben Goodrich:

On Wednesday, September 19, 2012 9:15:10 PM UTC-4, tds151 wrote:

To your general comment on structured covariance matrices of a multivariate Gaussian, is your point that if there are an underlying set of fewer parameters generates some covariance matrix, then the gradient may be more efficiently processed by stating the multivate Gaussian prior in terms of those parameters, rather than the resulting covariance matrix?   For example, what if I had a factor analytic structure where Sigma = Lambda * Lambda' + diag(sigma)?  Same idea, again, if I specified a proper CAR covariance structure?

In short, yes. TL;DR version as follows:

In general, I would say that if you are using a multivariate normal or a multivariate t distribution and your covariance matrix has a particular structure, then there are computational gains to be had if you feel comfortable with handling that part inside the .stan file with

lp__ <- lp__ + whatever;

syntax as in this thread:

https://groups.google.com/d/topic/stan-users/HcvYaDu71_Y/discussion

The reason is that the multivariate normal or a multivariate t distributions in Stan are currently only parameterized in terms of the covariance matrix. So, the C++ code that implements these two distributions inverts the covariance matrix (it actually takes a Cholesky decomposition and inverts the Cholesky factor) which is costly for large matrices and generates a huge expression tree to keep track of the auto-differentiation.

So, if you can do something clever to invert the covariance matrix, you can go faster with less RAM used. In the link above, we had the case where the covariance matrix was a scalar variance parameter times a known matrix. For a factor structure like

Sigma <- Lambda * Phi * Lambda' + Theta;

with a diagonal Theta, you could use the Woodbury Identity

https://en.wikipedia.org/wiki/Woodbury_matrix_identity

to get something like

for(k in 1:K) {
  for(j in 1:K) {
    Theta_inverse[j,k] = 0.0;
  }
  Theta_inverse[k,k] <- 1/Theta[k,k];
}
Sigma_inverse <- Theta_inverse - Theta_inverse * Lambda *
                 inverse(inverse(Phi) + Lambda' * Theta_inverse * Lambda) *
                 Lambda' * Theta_inverse;

(and in your case Phi is the identity matrix so you don't actually need to do inverse(Phi))

In which case, you could do some of this in the .stan file

lp__ <- lp__ - 0.5 * (x - mu)' * Sigma_inverse * (x - mu);

But if that is all you need Sigma_inverse for, then you could consolidate without explicitly forming Sigma_inverse

model {
  vector[K] x_mu_ThetaInv;
  real sum_of_squares;
  for(k in 1:K) {
    x_mu_ThetaInv[k] <- (x - mu) / Theta[k,k];
  }
  // (x - mu)' * Sigma_inverse * (x - mu)
  sum_of_squares <- dot_product(x_mu_ThetaInv,x_mu_ThetaInv) -       
                    x_mu_ThetaInv' * Lambda *

                    inverse(inverse(Phi) + Lambda' * Theta_inverse * Lambda) *
                    Lambda' * x_mu_ThetaInv;
  lp__ <- lp__ - 0.5 * sum_of_squares;

But, you need to remember the (log) determinant part of the multivariate normal (or t). You could do log(determinant(Sigma)) or log(1/determinant(Sigma_inverse)) but the determinant() function is a computational hazard because it involves a QR decomposition of a general parameterized matrix and thus a giant expression tree of auto-diff stuff. So, instead it is a bit better to calculate the determinant yourself using a Cholesky decomposition to exploit the positive definiteness of Sigma, as in

L <- cholesky_decompose(Sigma);
// det(Sigma) == det(L)^2 so 0.5 * log(det(Sigma)) == sum(log(diagonal(L)))
lp__ <- lp__ - sum(log(diagonal(L)));

That would probably help a good bit, but the inverse() function also does not exploit positive definiteness. So, you could gain a little more speed by doing the inversion via a Cholesky factorization. Except it looks like we don't actually expose the triangular solver to the Stan parser for some reason, so never mind.

For a covariance matrix with AR1 structure, equation 5.2.15 in Hamilton's textbook

http://books.google.com/books?id=B8_1UBmqVUoC&lpg=PP1&ots=2oJHZoc1xP&dq=Hamilton%20Time&pg=PA120#v=onepage&q=%22could%20be%20written%20down%20immediately%22&f=false

implies something like

parameters {
  real phi<lower=-1,upper=1>;
  real sigma_squared<lower=0>;
  vector[K] mu;
}
transformed parameters {
  real denominator;
  cov_matrix[K] Omega;
  denominator <- 1 - square(phi);
  for(i in 1:(K - 1)) {
    Omega[i,i] <- sigma_squared / denominator;
    for(j in (i+1):K) {
      Omega[i,j] <- sigma_squared * pow(rho, abs(j - i)) / denominator;
      Omega[j,i] <- Omega[i,j];
    }
  }
  Omega[K,K] <- sigma_squared / denominator;
}
model {
  x ~ multi_normal(mu, V); // entails slow inversion of V
}

But Equation 5.2.18 in Hamilton also implies that a faster way to get (x - mu)' * inverse(V) * (x - mu) is to use the Cholesky decomposition where

V^-1 = L' * L
Omega^-1 = 1/sigma_squared * V^-1

and L is lower triangular. So, we get

parameters {
  real phi<lower=-1,upper=1>;
  real sigma_squared<lower=0>;
  vector[K] mu;
}
model {
  vector[K] x_mu;
  real sum_of_squares;
  x_mu <- x - mu;
  // post-multiply L by x_mu without explicitly forming L, square, and sum
  sum_of_squares <- (1 - square(phi)) * square(x_mu[1]); // (row(L,1) * x_mu)^2
  for(k in 2:K) {
    sum_of_squares <- sum_of_squares + square(-phi * x_mu[k - 1] + x_mu[k]);
  }
  sum_of_squares <- sum_of_squares / sigma_squared;
  lp__ <- lp__ - 0.5 * sum_of_squares; // kernel of multivariate normal
  // add 0.5 * log(det(1/sigma_squared * V^-1)) in the form of:
  lp__ <- lp__ + 0.5 * log( 1 - square(phi) ) + 0.5 * K * log(1 / sigma_squared);
}

Hopefully, that is sufficiently clear without too many typos. This will get easier for the user once we get around to implementing a covariance_matrix class in C++.

Ben

[Ben Goodrich]

Thanks Terrance! (Google Groups wouldn't let me post to the list last night). I found (at least) one mistake in what I said:

But if that is all you need Sigma_inverse for, then you could consolidate without explicitly forming Sigma_inverse

model {
  vector[K] x_mu_ThetaInv;
  real sum_of_squares;
  for(k in 1:K) {
    x_mu_ThetaInv[k] <- (x - mu) / Theta[k,k];
  }
  // (x - mu)' * Sigma_inverse * (x - mu)
  sum_of_squares <- dot_product(x_mu_ThetaInv,x_mu_ThetaInv) -       
                    x_mu_ThetaInv' * Lambda *

                    inverse(inverse(Phi) + Lambda' * Theta_inverse * Lambda) *
                    Lambda' * x_mu_ThetaInv;
  lp__ <- lp__ - 0.5 * sum_of_squares;

That was wrong because dot_product(x_mu_ThetaInv,x_mu_ThetaInv) != (x-mu)' * ThetaInv * (x-mu). And there is a typo because (x - mu) is a vector. So, better would be

model {
  vector[K] x_mu_ThetaInv;
  real sum_of_squares;

  sum_of_squares <- 0.0;
  for(k in 1:K) {
    x_mu_ThetaInv[k] <- (x[k] - mu[k]) / Theta[k,k];
    sum_of_squares <- sum_of_squares + x_mu_ThetaInv[k] * (x[k] - mu[k]);
  }

  // (x - mu)' * Sigma_inverse * (x - mu)

  sum_of_squares <- sum_of_squares -       
                    x_mu_ThetaInv' * Lambda *

                    inverse(inverse(Phi) + Lambda' * Theta_inverse * Lambda) *
                    (Lambda' * x_mu_ThetaInv);

  lp__ <- lp__ - 0.5 * sum_of_squares;

Ben

[terrance savitsky]

Ben, your note and tutorial is *unbelievably* helpful.   i have a much better insight into how i may leverage the flexibility of Stan to mitigate RAM or speed issues with the auto-differentiation.  my particular example employs a covariance matrix composed from the kronecker product of 2 others, one of which employs a factor analytic construction.   after reading your tutorial, i realize that I may simplify my work by employing a matrix variate Gaussian that allows me to separately work with each component (of the kronecker product) covariance matrix to formulate my sum_of_squares update to lp_.   this is a nice procedure.

Thanks much for consultation.

Terrance.

[Andrew Gelman]

If it's incredibly helpful, I guess it should go into the Stan manual (if it's not already there…)

[terrance savitsky]

Funny, I thought the same - just didn't want to be presumptuous to say so.  Ben lays out a nice tutorial.   

[Ben Goodrich]

On Thursday, September 20, 2012 5:05:56 PM UTC-4, tds151 wrote:

Funny, I thought the same - just didn't want to be presumptuous to say so.  Ben lays out a nice tutorial.   

On Thu, Sep 20, 2012 at 2:03 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

If it's incredibly helpful, I guess it should go into the Stan manual (if it's not already there…)

On Sep 20, 2012, at 5:02 PM, terrance savitsky wrote:

Ben, your note and tutorial is *unbelievably* helpful.   i have a much better insight into how i may leverage the flexibility of Stan to mitigate RAM or speed issues with the auto-differentiation.  my particular example employs a covariance matrix composed from the kronecker product of 2 others, one of which employs a factor analytic construction.   after reading your tutorial, i realize that I may simplify my work by employing a matrix variate Gaussian that allows me to separately work with each component (of the kronecker product) covariance matrix to formulate my sum_of_squares update to lp_.   this is a nice procedure.

Thanks much for consultation.

Terrance.

You are the second or third person to mention Kronecker products, so hopefully Stan will have a kronecker() function soon. I guess you are doing the Kronecker product by hand, in which case, yes it makes sense wherever possible to manipulate the smaller blocks. Basically, anytime you have to do an O(K^3) matrix decomposition, stop and think if there might be a shortcut.

As regards to the manual, hopefully everything I said will be obsolete soon. Once we have a covariance matrix class in C++, then it will just be a matter of specifying inherited classes with various special structures.

Ben

[terrance savitsky]

Interesting.  So eventually we won't lose computational efficiency by directly employing a covariance matrix built from a structured configuration in a multivariate normal distribution once a covariance matrix class is defined in C++?  The woodbury matrix identity wouldn't be employed, tho??

In my case re: the kronecker product, I use your approach to jujitsu my way out of direct computation of the kronecker product.   in particular, I directly augment lp__ with the matrix variate distribution of Dawid, where for a D x M parameter matrix, Y,  the matrix variate normal is notated with:

Y - C ~ N(P,Q) where P is a D x D covariance matrix and Q is an M x M covariance matrix.   C is a D x M mean.

The density is proportional to:  det(P)^{-M/2} * det(Q)^{-D/2} * exp(-0.5*trace[P^{-1} * Y * Q^{-1} * Y^{'}]).

In my case, P is specified under a sparse, factor analytic construction.  

[Ben Goodrich]

On Thursday, September 20, 2012 8:52:50 PM UTC-4, tds151 wrote:

Interesting.  So eventually we won't lose computational efficiency by directly employing a covariance matrix built from a structured configuration in a multivariate normal distribution once a covariance matrix class is defined in C++?  The woodbury matrix identity wouldn't be employed, tho??

#include<forward_looking_statements_disclaimer>

Eventually, there will be an abstract covariance matrix class and a bunch of subclasses with specialized methods for common operations. The public attributes of the abstract covariance matrix class would be things like the correlation matrix, the standard deviations, the Cholesky factor, the precision matrix, the Cholesky factor of the precision matrix, etc. plus a bunch of flags indicating whether each attribute has been allocated yet. Then when some function needs one of these attributes to do some calculation, it would first have to check whether it has been allocated and if not, do it in an efficient way and then cache the result so that subsequent functions do not need to recreate it.

So, in your case, we might have an EFA structure covariance matrix class that inherits from the abstract covariance matrix class and has additional attributes for the loadings and unique variances. In which case, we would specialize the inverse() function to do Woodbury-based calculations (as opposed to Cholesky-based calculations for a generic covariance matrix or LU-based calculations for a generic square matrix). Similarly, we would create a function called quadratic_form() that takes a vector (x), a covariance matrix (Sigma), and a bool indicating whether to calculate x' * Sigma * x or x' * Sigma^-1 * x and it would spit out that sum of squares in the best way we know how for that class of covariance matrix.

So, for an AR1 structured covariance matrix class, there would be additional attributes for phi and sigma_squared. The inverse() function would be implemented analytically to exploit the tridiagonality and ideally would produce a sparse matrix. And quadratic_form() would exploit the bidiagonality of the Cholesky factor of the precision matrix, and so on.

In my case re: the kronecker product, I use your approach to jujitsu my way out of direct computation of the kronecker product.   in particular, I directly augment lp__ with the matrix variate distribution of Dawid, where for a D x M parameter matrix, Y,  the matrix variate normal is notated with:

Y - C ~ N(P,Q) where P is a D x D covariance matrix and Q is an M x M covariance matrix.   C is a D x M mean.

The density is proportional to:  det(P)^{-M/2} * det(Q)^{-D/2} * exp(-0.5*trace[P^{-1} * Y * Q^{-1} * Y^{'}]).

In my case, P is specified under a sparse, factor analytic construction.

You might look into the "bag of trace tricks" also. Taken literally, Stan will construct the whole matrix, P^{-1} * Y * Q^{-1} * Y^{'} , even though it only needs the diagonal to calculate the trace. If you define A' = P^{-1} * Y and B = Q^{-1} * Y^{'}, then

https://en.wikipedia.org/wiki/Trace_%28linear_algebra%29#Trace_of_a_product
 
says that trace(A' * B) is equal to the sum of the Hadamaard product. So, in Stan syntax it would be

sum(A .* B)

We use that trick in src/stan/prob/distributions/multivariate/continuous/{inv_}wishart.hpp .

Ben