Eigen-decomposition-based covariance-matrix transforms?
Bob Carpenter
as specified in terms of Cholesky factors by Ben.
I'm still blown away that the CPCs generate a
correlation matrix when all the dust settles.
Having just boned up on positive definiteness,
it appears another way to specify a covariance
matrix is through an Eigen decomposition in
terms of a set of K orthonormal vectors and
K positive ordered Eigenvalues. We already have
the data type and transform for the Eigenvalues.
Conceptually, we should be able to generate a
set of orthonormal vectors by
choose normal x[1] from R^K
choose normal x[2] from (K-1)-dim null space of { x[1] }
choose normal x[3] from (K-2)-dim null space of { x[1], x[2] }
...
Is there an efficient way to do that without
running solvers? For instance, in 3D, generating
the first vector is easy, with
x[1,1]^2 in +/- (0,1)
x[1,2]^2 in +/- (0,1 - x[1,1]^2)
x[1,3]^2 = 1 - x[1,1]^2 - x[1,2]^2
x[2,] is then constrained so that
1 = x[2,]' * x[2,]
and
0 = x[1,]' * x[2,]
So we have to choose an x[2,1] that ensures the above
is solvable (basically, the absolute value of
x[1,1] * x[2,1] can't be too large so that the remaining
coordinates can't bring the solution back to 0 while
retaining unit length).
Then once we choose x[2,1], the values x[2,2] and x[2,3]
will be solvable for because of the unit-length constraint.
Then x[3,] is fully determined by the unit-length constraint.
That leaves (K choose 2) free parameters in total,
plus K eigenvalue parameters.
- Bob
P.S. In the current setup, we can use the cov_matrix
transform to define all sorts of priors on covariance
matrices. For example, we could put a prior on the
positive ordered sequence of Eigenvalues. One way to
do this is choose the smallest Eigenvalue, the condition
[ratio of largest to smallest], then fill in the ones
in between choosing uniformly on the remaining interval.
Ben Goodrich
I finally got my head around the LKJ transform
as specified in terms of Cholesky factors by Ben.
I'm still blown away that the CPCs generate a
correlation matrix when all the dust settles.
If you like LKJ, you will love Lie algebra.
Having just boned up on positive definiteness,
it appears another way to specify a covariance
matrix is through an Eigen decomposition in
terms of a set of K orthonormal vectors and
K positive ordered Eigenvalues. We already have
the data type and transform for the Eigenvalues.Conceptually, we should be able to generate a
set of orthonormal vectors bychoose normal x[1] from R^K
choose normal x[2] from (K-1)-dim null space of { x[1] }
choose normal x[3] from (K-2)-dim null space of { x[1], x[2] }
...Is there an efficient way to do that without
running solvers?
http://www.mit.edu/~wingated/introductions/stiefel-mfld.pdf
which is sufficient and perhaps necessary for unconstrained algorithms involving eigenvector matrices with more than 2 or 3 columns.
In short, this is a job for a genuine math person, although Stan is probably the right tool for a genuine math person with a problem that requires something like this.
Ben
Marcus Brubaker
I don't know anything about this LKJ transform, could you send along a reference? Google seems to have failed me on that one.
At any rate, a couple things came to mind for representing orthonormal vectors for covariance matrices.
1. Grahm-Schmidt orthogonalization on a set of non-orthonormal vectors
2. Ignore orthogonality and just take use a set of vectors. The odds of n randomly chosen vectors actually being linearly dependent is small. (This approach is taken in, e.g., collaborative filtering problems where enforcing the orthonormal constraint is expensive and/or hard.)
3. Use the skew-symmetric representation of the tangent space or orthogonal matrices. The mapping to an orthogonal matrix is then the matrix exponential of a skew-symmetric matrix. See http://en.wikipedia.org/wiki/Skew-symmetric_matrix#Infinitesimal_rotations
4. (Shameless self plug) Implement non-linearly constrained HMC. See http://www.cs.toronto.edu/~mbrubake/projects/cmcmc/ If I ever get some more free time, I might implement this myself and would be happy to help someone else implement it. There are even some clever specializations one could do for orthogonality constraints. It would be a decidedly unique feature for stan.
Those ideas aside, is there a reason to replace the Cholesky factor based representation? It seems easy to implement and has some advantages when computing elliptical distributions which would need the Cholesky factors anyway.
Cheers,
Marcus
Ben Goodrich
Hey Bob,
I don't know anything about this LKJ transform, could you send along a reference? Google seems to have failed me on that one.
That is because Bob and I coined LKJ because we couldn't remember how to spell L or K's name
http://scholar.google.com/scholar?cluster=5449838284935626139&hl=en&as_sdt=0,33
But LKJ relies on scalar algebra and didn't actually work out the Cholesky-based form in Stan. I did.
Those ideas aside, is there a reason to replace the Cholesky factor based representation? It seems easy to implement and has some advantages when computing elliptical distributions which would need the Cholesky factors anyway.
Agreed.
Ben
Bob Carpenter
On 4/2/12 4:43 PM, Marcus Brubaker wrote:
> Hey Bob,
>
> I don't know anything about this LKJ transform, could you send along a
> reference? Google seems to have failed me on that one.
:-) See Ben's response. Lewandowski, Kurowicka and Joe. I'm a good
speller, but it's a mouthful.
As Ben says, we used his version of the transform based on
Cholesky factors. Ben also has a note he sent me earlier with
more details. That should probably at least be made into
an Arxiv paper to cite in our doc.
> At any rate, a couple things came to mind for representing orthonormal
> vectors for covariance matrices.
>
> 1. Grahm-Schmidt orthogonalization on a set of non-orthonormal vectors
I thought of that, but isn't there a problem because
the orthogonalization is many to one? I've never done
non-bijective transforms, but hey, I had never done any
before last year.
> 2. Ignore orthogonality and just take use a set of vectors. The odds of
> n randomly chosen vectors actually being linearly dependent is small.
> (This approach is taken in, e.g., collaborative filtering problems where
> enforcing the orthonormal constraint is expensive and/or hard.)
Approach 1. would rely on this feature, too.
Don't the linearly dependent vectors have measure zero among
all the sets of vectors?
> 3. Use the skew-symmetric representation of the tangent space or
> orthogonal matrices. The mapping to an orthogonal matrix is then the
> matrix exponential of a skew-symmetric matrix. See
> http://en.wikipedia.org/wiki/Skew-symmetric_matrix#Infinitesimal_rotations
Thanks -- hadn't seen that. We'll see if I can
understand it.
> 4. (Shameless self plug) Implement non-linearly constrained HMC. See
> http://www.cs.toronto.edu/~mbrubake/projects/cmcmc/ If I ever get some
> more free time, I might implement this myself and would be happy to help
> someone else implement it. There are even some clever specializations
> one could do for orthogonality constraints. It would be a decidedly
> unique feature for stan.
I'll take a look. Michael Betancourt's also been urging
us to think about implementing at least simple constraints.
I bcc-ed him on this because I don't think he's on the stan-dev
list.
> Those ideas aside, is there a reason to replace the Cholesky factor
> based representation?
No. I should've made it clear that this was simply
an exercise I set for myself to help me understand
all the linear algebra winging about.
> It seems easy to implement and has some
> advantages when computing elliptical distributions which would need the
> Cholesky factors anyway.
I didn't even understand this point until I
finally understood the Cholesky factors and their
utility for specifying positive definite matrices
and multivariate distributions. This is the point Ben's
been trying to get through my thick skull for around
a year now. Thanks again for that patience, Ben!
- Bob
Bob Carpenter
I was looking for! I got the feeling from Strang's
linear algebra book that Gram-Schmidt would be the
efficient way to do the solving.
I've always wanted to be enough of a mathematician
to have a use for algebraic topology. Algebra and
topology were two of my fave undergrad subjects, and
though I took linear algebra and ring theory, I never
really did matrices in an applied math setting, like
say, physics or statistics.
- Bob
Marcus Brubaker
:-) See Ben's response. Lewandowski, Kurowicka and Joe. I'm a good
On 4/2/12 4:43 PM, Marcus Brubaker wrote:
Hey Bob,
I don't know anything about this LKJ transform, could you send along a
reference? Google seems to have failed me on that one.
speller, but it's a mouthful.
As Ben says, we used his version of the transform based on
Cholesky factors. Ben also has a note he sent me earlier with
more details. That should probably at least be made into
an Arxiv paper to cite in our doc.I thought of that, but isn't there a problem because
At any rate, a couple things came to mind for representing orthonormal
vectors for covariance matrices.
1. Grahm-Schmidt orthogonalization on a set of non-orthonormal vectors
the orthogonalization is many to one? I've never done
non-bijective transforms, but hey, I had never done any
before last year.
Yes, it would make computing the log-determinant of the jacobian impossible as the Jacobian would be rectangular.
Approach 1. would rely on this feature, too.
2. Ignore orthogonality and just take use a set of vectors. The odds of
n randomly chosen vectors actually being linearly dependent is small.
(This approach is taken in, e.g., collaborative filtering problems where
enforcing the orthonormal constraint is expensive and/or hard.)
Don't the linearly dependent vectors have measure zero among
all the sets of vectors?
Almost, the set of nxn singular matrices is of measure zero. A set of n n-dimensional vectors is almost certain to a) be somewhat linearly dependent and b) still span the full space. b) is just a restatement of the measure zero argument. a) follows by considering that, with n-1 orthonormal vectors, there are exactly 2 choices of the last vector, which must have measure zero.
Thanks -- hadn't seen that. We'll see if I can
3. Use the skew-symmetric representation of the tangent space or
orthogonal matrices. The mapping to an orthogonal matrix is then the
matrix exponential of a skew-symmetric matrix. See
http://en.wikipedia.org/wiki/Skew-symmetric_matrix#Infinitesimal_rotations
understand it.
In 3D this is called the exponential map representation of rotations which might help with intuition.
Basically, you represent an n dimensional orthonormal matrix as an n(n − 1)/2-dimensional vector whose entries comprise the off diagonal entries of the skew-symmetric matrix. The magnitude of the vector is the angle of the rotation and the direction of the vector is the axis of rotation. This is a minimal representation of O(n), but has singularities when the magnitude is zero or a multiple of 2pi. These are almost unavoidable (see hairy ball theorem) and rarely a problem in practice, but can come up.
Cheers,
Marcus
Ben Goodrich
On Monday, April 2, 2012 5:19:37 PM UTC-4, Bob Carpenter wrote:
Neat. The Stiefel manifold's exactly what
I was looking for! I got the feeling from Strang's
linear algebra book that Gram-Schmidt would be the
efficient way to do the solving.
I remembered that paper getting into the Lie algebra stuff, but actually I must have been thinking of even more beautiful work by Mark Plumbley:
http://scholar.google.com/scholar?hl=en&q=author%3AM-Plumbley+exponential+orthogonal+matrix&btnG=Search&as_sdt=0%2C33&as_ylo=&as_vis=0
Markus' #3
3. Use the skew-symmetric representation of the tangent space or orthogonal matrices. The mapping to an orthogonal matrix is then the matrix exponential of a skew-symmetric matrix. See http://en.wikipedia.org/wiki/Skew-symmetric_matrix#Infinitesimal_rotations
Ben
Ben Goodrich
As Ben says, we used his version of the transform based onCholesky factors. Ben also has a note he sent me earlier with
more details. That should probably at least be made into
an Arxiv paper to cite in our doc.
Ben
Ben Goodrich
If Sigma is an unknown KxK covariance matrix, then it has a Cholesky decomposition
Sigma = L * L'
where L is lower-triangular with positive diagonal entries. Conversely, for any such L, L * L' is a valid covariance matrix.
So, Stan could simply parameterize a covariance matrix in terms of the K choose 2 strictly lower-triangular elements of L --- which are unconstrained --- and the K diagonal elements of L, which are unconstrained aside from needing to be positive. At the moment, Stan does *not* do that, but we could, particularly if the user wants to place a {scaled, inverse} Wishart prior on Sigma. Conversely, it would be very difficult to specify a prior on the unbounded parameters and understand what that implies for Sigma.
What Stan actually does at the moment is to write
Sigma = D * L * L' * D
where D is a diagonal matrix of standard deviations and L is parameterized in terms of so-called canonical partial correlations (CPCs) such that L * L' is a correlation matrix. Furthermore, we use exp() to form a positive standard deviation from an unbounded parameter and tanh() to bind a CPC to the (-1,1) interval from an unbounded parameter. In other words, that entails three more levels of transformation that are not strictly necessary and are in fact wasteful if the user wants to place a {scaled, inverse} Wishart prior on Sigma.
The motivation for the extra three levels of transformation --- and the ultimate contribution of the LKJ paper --- is that you can write the joint density of a correlation matrix in a simple, useful way:
f(L(eta) * L(eta)' | eta) = constant * det[ L(eta) * L(eta)' ]^(eta - 1)
where eta > 0 is a hyperprior that affects the shape of the density for each CPC. This gives you a firm basis for specifying a prior on a correlation matrix: Iff eta = 1, then the joint density is constant and all correlation matrices are equally likely, and as eta -> oo, L(eta) * L(eta)' has to converge to the identity matrix. And it is pretty easy to supplement this with a prior on the standard deviations to imply (with an easy change-of-variables) a prior on a covariance matrix.
I think that is a substantively better way to specify a prior on a covariance matrix, but certainly the {scaled, inverse} Wishart priors are more popular because they were the only practical choices until a few years ago. So, if a user is going to use a {scaled, inverse} Wishart prior on Sigma, is it worth our coding effort to provide a simpler parameterization in terms of the lower-triangular elements of L? I'd say probably. And to get back to Bob's question, is there any other decomposition of a covariance matrix that is both computationally efficient and lends itself to specifying reasonable priors? Not that I know of, which is why I'd say to wait on that until a situation presents itself where an alternate decomposition is necessary.
Ben
Bob Carpenter
after fussing around with builds and writing doc (though
I love writing doc more than just about anything -- it
feels like polishing a nice piece of furniture we just
made).
On 4/3/12 1:03 PM, Ben Goodrich wrote:
> Bob's post about parameterizing a covariance matrix in terms of a
> decomposition raises another point that I had been meaning to bring up.
>
> If Sigma is an unknown KxK covariance matrix, then it has a Cholesky
> decomposition
>
> Sigma = L * L'
>
> where L is lower-triangular with positive diagonal entries. Conversely,
> for any such L, L * L' is a valid covariance matrix.
Awesome -- that's exactly what I was looking for.
We should definitely use that as our basic parameterization
for covar_matrix!
I didn't realize there weren't any constraints on
the resulting Cholesky factor other than positive diagonals.
I'm still working on my merit badge in matrices and linear
algebra.
Is there an easy way to constrain L so that
L * L' is a correlation matrix? Basically each row
of L needs to be unit length so that
(L * L')[i,i] = L[i] dot L[i]
= 1.
I know we could just divide, but we want to do this in such
a way that the transform is bijective so we can easily
compute the change of variables.
How should we build data types into Stan to support
LKJ? We can now separate the issue from that of specifying
a basic covariance matrix transform.
In the worst case, we'd rount trip through the
factorizations every leapfrog step, but that would be
very very slow as Ben's made clear in previous messages.
It's not that hard to add basic data types to Stan. All
we need is:
mapping to C++ type and getters/setters (trivial if matrix-like,
but we should think about triangular and symmetric matrix
types)
unconstrained -> constrained transform & log abs Jacobian det
constrained -> unconstrained transform
So I'm thinking we need to do this for Cholesky factors
that lead to positive-definitite covariance matrices (and/or
correlation matrices).
Any suggestions for a name? Then we could use multivariate
densities defined directly in terms of Cholesky factors,
which should be much speedier as Ben's pointed out before.
I spent last night playing around with Bingham-von Mises-Fisher
distributions on the hypersphere and was thinking adding a
unit-K-hypersphere data type might be useful. (Especially if
we want to study paleomagnetism, which is the main app I saw
sited.)
> So, Stan could simply parameterize a covariance matrix in terms of the K
> choose 2 strictly lower-triangular elements of L --- which are
> unconstrained --- and the K diagonal elements of L, which are
> unconstrained aside from needing to be positive.
That makes sense. And now I see the problem -- the K diagonal
elements don't represent the scale any more. See question (2) above.
> At the moment, Stan
> does *not* do that, but we could, particularly if the user wants to
> place a {scaled, inverse} Wishart prior on Sigma.
That's one motivation, but the real motiviation is just to
guarantee the positive definiteness constraint (up to arithmetic
singularities, of course).
I also want it in case they want to do something like put
a prior on the eigenvalues. But that'd be expensive as it'd
still require an Eigen-decomposition of the L * L' positive
definite matrix.
> Conversely, it would
> be very difficult to specify a prior on the unbounded parameters and
> understand what that implies for Sigma.
>
> What Stan actually does at the moment is to write
>
> Sigma = D * L * L' * D
>
> where D is a diagonal matrix of standard deviations and L is
> parameterized in terms of so-called canonical partial correlations
> (CPCs) such that L * L' is a correlation matrix. Furthermore, we use
> exp() to form a positive standard deviation from an unbounded parameter
> and tanh() to bind a CPC to the (-1,1) interval from an unbounded
> parameter. In other words, that entails three more levels of
> transformation that are not strictly necessary and are in fact wasteful
> if the user wants to place a {scaled, inverse} Wishart prior on Sigma.
>
> The motivation for the extra three levels of transformation --- and the
> ultimate contribution of the LKJ paper --- is that you can write the
> joint density of a correlation matrix in a simple, useful way:
>
> f(L(eta) * L(eta)' | eta) = constant * det[ L(eta) * L(eta)' ]^(eta - 1)
>
> where eta > 0 is a hyperprior that affects the shape of the density for
> each CPC.
What is L(eta) here?
> This gives you a firm basis for specifying a prior on a
> correlation matrix: Iff eta = 1, then the joint density is constant and
> all correlation matrices are equally likely, and as eta -> oo, L(eta) *
> L(eta)' has to converge to the identity matrix.
I agree that this is much easier to understand and
more natural conceptually than the (scaled) (inverse)
Wishart business.
> And it is pretty easy to
> supplement this with a prior on the standard deviations to imply (with
> an easy change-of-variables) a prior on a covariance matrix.
This is what we'd like to make easy to specify in
Stan. So we want something like a positive vector
of scales and a correlation matrix.
> I think that is a substantively better way to specify a prior on a
> covariance matrix, but certainly the {scaled, inverse} Wishart priors
> are more popular because they were the only practical choices until a
> few years ago. So, if a user is going to use a {scaled, inverse} Wishart
> prior on Sigma, is it worth our coding effort to provide a simpler
> parameterization in terms of the lower-triangular elements of L? I'd say
> probably.
Especially if we want to compete with BUGS/JAGS in speed
on simpler models.
> And to get back to Bob's question, is there any other
> decomposition of a covariance matrix that is both computationally
> efficient and lends itself to specifying reasonable priors? Not that I
> know of,which is why I'd say to wait on that until a situation presents
> itself where an alternate decomposition is necessary.
Agreed.
- Bob
Ben Goodrich
Is there an easy way to constrain L so that
L * L' is a correlation matrix? Basically each row
of L needs to be unit length so that(L * L')[i,i] = L[i] dot L[i]
= 1.
I know we could just divide, but we want to do this in such
a way that the transform is bijective so we can easily
compute the change of variables.
How should we build data types into Stan to support
LKJ? We can now separate the issue from that of specifying
a basic covariance matrix transform.
If we don't go whole hog to a CovMatrix class, then we somehow need to figure out from the .stan file if the user wants a LKJ prior and in that case, we need to do the more complicated LKJ transform. Ergo, I don't like this route as much.
Any suggestions for a name? Then we could use multivariate
densities defined directly in terms of Cholesky factors,
which should be much speedier as Ben's pointed out before.
> joint density of a correlation matrix in a simple, useful way:
>
> f(L(eta) * L(eta)' | eta) = constant * det[ L(eta) * L(eta)' ]^(eta - 1)
>
> where eta > 0 is a hyperprior that affects the shape of the density for
> each CPC.What is L(eta) here?
In the simpler parameterization where
Sigma = L * L'
Ben
Marcus Brubaker
So I think this discussion is conflating two different issues/use cases:
1. Someone wants to define a covariance matrix for use in a model in arbitrary or yet to be determined ways.
2. Someone wants to specify a prior in terms of the LKJ representation, the cholesky factorization, the eigenvalue decomposition or some other arbitrary representation of the covariance matrix.
1 means the covariance matrix is primal and how it's used is arbitrary or unpredictable. Hence we want the simplest, most compact and efficient way to represent a (positive definite) covariance matrix. Here it seems the L * L' form with a log transformation on the diagonal is the most sensible.
2 on the other hand means that from a covariance matrix we'll need to compute things which may be inconvenient or expensive and we may just want to represent the covariance matrix directly in those terms.
For me it seems that there should be different types for these cases in the Stan modelling language. A basic definition of a cov_matrix corresponds to use case 1. The LKJ form could be handled with a corr_matrix type and vector which is constrained to have positive elements and the user would define a transformed parameter for the full covariance. The Cholesky case, the eigenvalue case and anything else could be handled similarly.
An example showing three different representations to make things clear:
parameters {
cov_matrix<N> Sigma_basic;
corr_matrix<N> Sigma_LKJ_correlations;
vector<N> Sigma_LKJ_logstddevs; // or add some way to make the entries positive
chol_factor<N> Sigma_cholesky_L;
}
transformed parameters {
matrix<N,N> Sigma_LKJ;
matrix<N,N> Sigma_cholesky;
Sigma_LKJ <- lkj_covariance_matrix(Sigma_LKJ_correlations,exp(Sigma_LKJ_logstddevs));
Sigma_cholesky <- Sigma_cholesky_L * Sigma_cholesky_L';
}
model {
x1 ~ multi_normal(mu,Sigma_basic);
lp__ <- lp__ + arbitrary_sigma_prior(Sigma_basic,...); // or whatever
x2 ~ multi_normal(mu,Sigma_LKJ);
Sigma_LKJ_correlations ~ lkj_corr_prior(...);
x3 ~ multi_normal(mu,Sigma_cholesky); // or multi_normal_chol(mu,Sigma_cholesky_L);
Sigma_cholesky_L ~ inverse_wishart_chol(...);
Etc.
}
In other words, if someone wants an LKJ prior, then they can define their model appropriately.
It places a somewhat larger burden on the user but avoids having to add heuristics to determine which representation to use based on how the model is specified which just seems ugly. However, if the user knows they want an LKJ prior, then presumably they know what the LKJ representation means, at least in terms of a correlation matrix and standard deviations.
Adding a CovMatrix class makes it possible to speed some things up automatically under the hood, but would still require deciding how best to represent the covariance matrix for a given model. That would mean the same kind of heuristics or some kind of hinting in the .stan file provided by the user.
Cheers,
Marcus
Ben Goodrich
That will teach me not to read the code first, I had assumed that Stan was using the L * L' representation of a covariance matrix.
That isn't so incorrect. At the moment, Stan does represent a covariance matrix as
Sigma = (D * L) * (L' * D)
where D is a diagonal matrix of standard deviations and L * L' is a correlation matrix but circumvents LKJ's recursion formula (thus D * L is still lower triangular). But we should unbundle the LKJ transform from the LKJ prior for the reasons you have given, which I basically agree with.
Adding a CovMatrix class makes it possible to speed some things up automatically under the hood, but would still require deciding how best to represent the covariance matrix for a given model. That would mean the same kind of heuristics or some kind of hinting in the .stan file provided by the user.
The way I had been envisioning it, the CovMatrix class would be kind of like the S4 Matrix class in R or the LLT class in Eigen in the sense that it "holds" redundant information about Sigma in various different ways. So, the public attributes of the CovMatrix class might be
Sigma
Sigma_is_initialized
L
L_is_initialized
Sigma_inverse
Sigma_inverse_is_initialized
L_inverse
L_inverse_is_initialized
eigenvalues
eigenvalues_are_initialized
eigenvectors
eigenvectors_are_initialized
determinant
determinant_is_calculated
transformation_scheme
...
The density function (or whatever function is downstream) could first check whether what it wants has been initialized. If not, initialize it. Then, calculate the return value.
So, for example, the multi_log function would test L_inverse_is_initialized, and in many cases, it wouldn't be. Then, test whether L_is_initialized, and in many cases, it would be. (If not, Sigma should have been initialized, so Cholesky decompose it and update L and L_is_initialized). Forward solve L to get L_inverse, cache it, and change L_inverse_is_initialized. Test determinant_is_calculated and if that fails, calculate and cache the determinant of Sigma from L. Now you have all the pieces to evaluate the density. And, moreover, everything is ready in case this CovMatrix needs to be used again (i.e., you need to evaluate the density for the next observation).
In short, we could have fewer verbosely named functions and allow each density function to figure out what they need from the CovMatrix class.
Ben
Ben Goodrich
In the simpler parameterization where
Sigma = L * L'
L is not a function of anything (except for the diagonal elements which are the antilog of something), so there is less for the computer to do. Although, we do need to code up the Jacobian for the transformation from the lower triangle of L to the lower triangle of Sigma.
Note to self mostly: This handout has a bunch of useful Jacobians
http://web.mit.edu/18.325/www/handouts/handout2.pdf
including the Jacobian determinant of vech(L) to vech(Sigma) at the end of section 4.1.
Ben
Marcus Brubaker
My main point was that if we want the model code to do the right thing we either a) need to figure out from the model specification which representation to use when parameterizing the model (I think this is going to be hard and error prone for everything but really simple cases) or have the user tell us, either by hinting when defining a cov_matrix parameter, or by requiring to the user to build the covariance matrix themselves in the model (what I had proposed).
Cheers,
Marcus
Bob Carpenter
On 4/3/12 3:33 PM, Marcus Brubaker wrote:
> That will teach me not to read the code first, I had assumed that Stan was using the L * L' representation of a
> covariance matrix.
>
> So I think this discussion is conflating two different issues/use cases:
We're obviously not communicating very well. These are exactly
the issues we've been talking about:
> 1. Someone wants to define a covariance matrix for use in a model in arbitrary or yet to be determined ways.
> 2. Someone wants to specify a prior in terms of the LKJ representation, the cholesky factorization, the eigenvalue
> decomposition or some other arbitrary representation of the covariance matrix.
The confusion was my lack of linear algebra skill that lead
to my thinking the LKJ representation was the only feasible
way to get from an arbitrary parameterization to a covariance
matrix.
> 1 means the covariance matrix is primal and how it's used is arbitrary or unpredictable. Hence we want the simplest,
> most compact and efficient way to represent a (positive definite) covariance matrix. Here it seems the L * L' form with
> a log transformation on the diagonal is the most sensible.
Exactly -- that's why I'm so keen to switch to the L*L' rep.
I think we may provide a type for that, too. But
we won't need to tie them together.
It still looks like LKJ might be the only simple way
to parameterize a correlation matrix, so maybe we can
keep the LKJ correlation matrix, which I think lets Ben
do everything he wants, or close to it, and then use
the L*L' version of covariance matrices as the basic
data type.
For instance, we could build in what you're doing in the model
above.
> It places a somewhat larger burden on the user but avoids having to add heuristics to determine which representation to
> use based on how the model is specified which just seems ugly. However, if the user knows they want an LKJ prior, then
> presumably they know what the LKJ representation means, at least in terms of a correlation matrix and standard deviations.
Right. I don't mind adding new, little-used types.
> Adding a CovMatrix class makes it possible to speed some things up automatically under the hood, but would still require
> deciding how best to represent the covariance matrix for a given model. That would mean the same kind of heuristics or
> some kind of hinting in the .stan file provided by the user.
Right.
- Bob
Ben Goodrich
I actually think a CovMatrix class is a good idea since they are such an integral part of almost any model.
My main point was that if we want the model code to do the right thing we either a) need to figure out from the model specification which representation to use when parameterizing the model (I think this is going to be hard and error prone for everything but really simple cases) or have the user tell us, either by hinting when defining a cov_matrix parameter, or by requiring to the user to build the covariance matrix themselves in the model (what I had proposed).
I agree. The only difference is that I was thinking something along these lines
Case 1:
parameters {
cov_matrix(N) Sigma;
}
model {
Sigma ~ wishart();
// something with Sigma
}
In this situation, we should use a plain Cholesky parameterization behind the scenes without the user knowing it or having to specify it explicitly. The C++ constructor for the CovMatrix would only initialize the L attribute.
In other words, the following more explicit version should be legal, but we shouldn't force the user to do
parameters {
chol_factor(N) L;
}
transformed parameters {
cov_matrix(N) Sigma;
Sigma <- L * L'
}
model {
Sigma ~ wishart();
// something with Sigma
}
But if the user wants to do LKJ, then along the lines of what you said
Case 2:
parameters {
vector CPCs(N * (N - 1));
vector sds(N);
}
transformed parameters {
cov_matrix(N) Sigma;
Sigma <- LKJ_transformation(CPCs, sds);
}
model {
double eta;
eta <- 1.0;
CPCs ~ LKJ_prior(eta);
for(i in 1:N) {
sds[i] ~ something();
}
// something with Sigma
}
Again, the Cholesky factorization is hidden from the user but integral to the C++.
Ben
Ben Goodrich
parameters {
double rho(-1,1);
double variance(0,);
}
transformed parameters {
cov_matrix(N) Sigma;
Sigma <- AR1(rho, variance);
}
model {
rho ~ uniform(-1,1)
variance ~ something();
// something with Sigma
}
In the AR1 case, it is usually better *not* to calculate a Cholesky factor because the covariance matrix, its inverse (which is tridiagonal), and its determinant are known in closed form as a function of rho and variance. So, the CovMatrix constructor in C++ might initialize these three things, but not L or L_inverse. But whatever density function is using Sigma doesn't need to care much about about that. It just needs to check whether the desired attributed is initialized before it tries to use it in a calculation.
Ben
Andrew Gelman
A
Marcus Brubaker
The CovMatrix class would allow Stan to be smarter about how covariance matrices are used when computing likelihoods.
So in summary, this is proposing to add a chol_factor and corr_matrix type to the modeling language. Beyond that maybe also a way to specify a vector of positive reals? The CovMatrix class can be added somewhat independently of this.
Cheers,
Marcus
Ben Goodrich
Speaking as a user, my first goal is to be able to fit multilevel models with the scaled-inverse-Wishart class of priors with dimension up to, say, 10 or so. But I'm sure I'll have other goals once that gets working!
The dimensionality is not a major issue for the proposed code refactoring as far as I can tell. But the dimensionality is relevant to whether a chosen prior is reasonable (which Andrew knows well because we co-wrote a paper with a couple of other people about how scaled-inverse-Wishart with dimensionality > 2 and low degrees of freedom puts more mass in the corners of the parameter space where covariance matrices are nearly singular). The {scaled, inverse} Wishart priors are overly popular due to their computational convenience in Gibbs samplers, but I suspect they are computationally inconvenient for HMC due to the induced high dependence among parameters. Maybe NUTS 3 will make this feasible, but it should already be feasible in NUTS 1 with an LKJ prior and eta >= 1.
Ben
Bob Carpenter
On 4/3/12 4:49 PM, Marcus Brubaker wrote:
> I think we're in agreement: a cov_matrix defaults to a Cholesky factor parameterization, and other forms can be
> explicitly instantiated by the user.
>
> The CovMatrix class would allow Stan to be smarter about how covariance matrices are used when computing likelihoods.
>
> So in summary, this is proposing to add a chol_factor and corr_matrix type to the modeling language. Beyond that maybe
> also a way to specify a vector of positive reals? The CovMatrix class can be added somewhat independently of this.
The currently supported types (arranged by ISA/assignability)
are:
matrix(M,N)
cov_matrix(K)
corr_matrix(K)
row_vector(N)
vector(M)
simplex(K)
pos_ordered(K)
real
real(lb,)
real(,ub)
real(lb,ub)
int
int(lb,)
int(,ub)
int(lb,ub)
plus arrays of arbitrary dimensions with entries
declared to be any of these types. Stan doesn't
currently forbid
cov_matrix(4) Sigma;
matrix(4,4) M;
...
Sigma <- M;
but every use case is checked (though not at
intermediate stages of processing).
There's not a positive vector type, but there is
a positive array of reals declared with lb = 0.
We could add bounds to the basic matrix, row vector
and vector types.
Three new types of square matrices are under discussion:
symmetric_matrix(K)
triangular_matrix(K)
cov_matrix_cholesky(K)
and a new vector for points on the unit
hypersphere:
hypersphere(K)
I'm about to reparameterize cov_matrix using the
Cholesky factors.
I've also thought some typedefs would be useful:
typedef real(0,) positive;
typedef real(0,1) probability;
typedef real(-1,1) correlation;
typedef int(0,) natural;
typedef int(0,1) binary;
I'm pretty sure we can handle correlation matrices
by using first a bound and then a hyperspherical
transform on the rows so that the sum of the square
entries is 1. That'll get us positive-definiteness
and a unit diagonal. All other entries must have
absolute values below 1 because they're all products
of unit vectors and only the diagonals are from
ones going in the same direction. I think I got
my geometrical reasoning right: if x and y are
unit-length vectors, then
x' * x >= x' * y >= x' * -x
- Bob
Bob Carpenter
On 4/3/12 4:31 PM, Ben Goodrich wrote:
> On Tuesday, April 3, 2012 4:11:34 PM UTC-4, Marcus Brubaker wrote:
>
> I actually think a CovMatrix class is a good idea since they are such an integral part of almost any model.
>
> My main point was that if we want the model code to do the right thing we either a) need to figure out from the
> model specification which representation to use when parameterizing the model (I think this is going to be hard and
> error prone for everything but really simple cases) or have the user tell us, either by hinting when defining a
> cov_matrix parameter, or by requiring to the user to build the covariance matrix themselves in the model (what I had
> proposed).
>
>
> I agree. The only difference is that I was thinking something along these lines
These 3 examples are in line with how we've been doing
things up to now (up to parameterizations and some syntax
in variable declarations for vectors).
> Case 1:
>
> parameters {
> cov_matrix(N) Sigma;
> }
> model {
> Sigma ~ wishart();
> // something with Sigma
> }
>
> In this situation, we should use a plain Cholesky parameterization behind the scenes without the user knowing it or
> having to specify it explicitly. The C++ constructor for the CovMatrix would only initialize the L attribute.
>
Right -- the key is that we allow a user to
declare a cov_matrix param Sigma and then do
whatever they want with it. All this transformation
stuff's in the background (which is why it'll be
easy to swap out LKJ for Cholesky).
> In other words, the following more explicit version should be legal, but we shouldn't force the user to do
>
> parameters {
> chol_factor(N) L;
> }
> transformed parameters {
> cov_matrix(N) Sigma;
> Sigma <- L * L'
> }
> model {
> Sigma ~ wishart();
> // something with Sigma
> }
The nice thing here is that they could also
do something with the entries of L for a prior
then use the resulting covariance matrix as
a parameter, as in:
L ~ some_funky_distro_on_cholesky_factors();
for (i in 1:I)
y[i] ~ multi_normal(mu,Sigma);
A funky distro could be designed for L along the
lines of:
for (n in 1:N)
for (n2 in 1:n)
L[n,n2] ~ log_normal(...);
As Ben illustrates in the next case:
> But if the user wants to do LKJ, then along the lines of what you said
>
> Case 2:
>
> parameters {
> vector CPCs(N * (N - 1));
> vector sds(N);
> }
> transformed parameters {
> cov_matrix(N) Sigma;
> Sigma <- LKJ_transformation(CPCs, sds);
> }
> model {
> double eta;
> eta <- 1.0;
> CPCs ~ LKJ_prior(eta);
> for(i in 1:N) {
> sds[i] ~ something();
> }
> // something with Sigma
> }
>
> Again, the Cholesky factorization is hidden from the user but integral to the C++.
This is the key for letting users write models flexibly
without themselves having to worry too much about the
constraints.
- Bob
Andrew Gelman
Ben Goodrich
Regarding scaled-inv-wishart, I agree that it has convenient conditional conjugacy, but the scaling (thus, the redundant parameterization) appeals to me because it allows the correlations to be controlled (i.e., smoothed, i.e., partially pooled) while allowing a reasonable prior distribution for the scale parameters. Plain old unscaled inverse-Wishart doesn't do that; the degrees of freedom simultaneously determine the precisions for the corrs and the scales.
It looks like Google ate my post with attachments again. Resending in a minute.
Ben
Ben Goodrich
The way I had been envisioning it, the CovMatrix class would be kind of like the S4 Matrix class in R or the LLT class in Eigen in the sense that it "holds" redundant information about Sigma in various different ways. So, the public attributes of the CovMatrix class might be
Sigma
Sigma_is_initialized
L
L_is_initialized
Sigma_inverse
Sigma_inverse_is_initialized
L_inverse
L_inverse_is_initialized
eigenvalues
eigenvalues_are_initialized
eigenvectors
eigenvectors_are_initialized
determinant
determinant_is_calculated
transformation_scheme
...
A corollary of the idea to parameterize
Sigma = L * L'
The user may want to store the posterior distribution of Sigma, in which case we would have to calculate [L^(-1)]^(-1) = L and then L * L'. But, we could do that after casting L^(-1) from var to double, so we wouldn't have to calculate a zillion auto-diffs.
Ben
Ben Goodrich
BenThat's interesting, thanks for sending. The pictures are pretty. Usually I don't get much from 3d scatterplots but these really work!The inv-Wishart with d+1 df gives uniform marginals on the corrs, but I can see how, to get the uniform marginal [sic, AG meant joint], you have to pile up mass on the corners. One solution would be to increase the df, right?
Yes. The expectation of the inverse Wishart is S / (d - k - 1) iff d > k + 1. So, I picked the popular but extreme case that the expectation doesn't exist. Calculating the expectation of an implied correlation is a bit involved, but as d increases, the standard deviation of the correlations shrinks.
library(MCMCpack)
k <- 3
I <- diag(k)
mark <- lower.tri(I, diag = FALSE)
REPS <- 10000
dof <- (k+1):20
results <- sapply(dof, FUN = function(d) replicate(REPS, cov2cor(riwish(d, I))[mark]) )
colnames(results) <- dof
plot(dof, apply(results, 2, sd), type = "l", xlab = "Degrees of freedom",
ylab = "Standard deviation of correlations")
However, I continue to think that increasing the shape parameter of a LKJ distribution (starting at eta = 1 -> joint uniform) is more intuitive than increasing the degrees of freedom of an inverse Wishart (starting at d = k + 1 -> no moments) in part because the "degrees of freedom" is kind of a nondescript number, while the LKJ correlation matrix is a composition of independent symmetric betas on (-1,1) whose shape parameters are a function of eta.
More fundamentally, though, we don't have a great way of thinking about these multivariate models. Once d>2, I have very little intuition about what sorts of models would make sense. Maybe the joint uniform is fine, I just can't get my head around what would be the best parameterization.
I think what you said is generally true: Most people don't have a good understanding of joint distributions, correlation matrices, etc. and even if they understand them, they don't have strong priors about correlation matrices, particularly if it is a correlation matrix among parameters rather than data. So, in that situation, a joint uniform distribution on a correlation matrix accurately reflects their prior beliefs. And insofar as they have informative prior beliefs, they probably slightly favor the identity matrix over the nearly singular correlation matrices, which can be accomplished with an eta slightly greater than 1.0.
A related issue is that in a multilevel model, if you shift the x-variables (for example, by centering them), this changes all the correlations.
Yes, if the correlation matrix in question is among parameters.
Based on what you're saying here, the uniform on the corrs does seem to make the most sense. Is there a parameterized version of this (just as the inv-wishart has that df parameterization)?
Yes, what Bob and I call the LKJ distribution over a correlation matrix (with shape parameter eta > ) in stan is based on this paper
http://scholar.google.com/scholar?cluster=5449838284935626139&hl=en&as_sdt=0,33
which gives conditions under which the PDF of a correlation matrix is
f(Sigma | eta) = 1 / c(k,eta) * det(Sigma)^(eta - 1)
where Sigma is a correlation matrix (and a function of eta) and c(k,eta) only depends on k and eta. And so eta = 1 -> f(Sigma | eta) = 1 / c(k,eta) is constant.
Also, getting back to the uniform model (which I think was also discussed by Barnard, Meng, and McCulloch), let's think briefly about high dimensions. I know that earlier I said I'm only interested in d up to 10, but it would be good to understand the limit of very high dimensions, if only to get a sense of what the model is doing. If d is high, then the joint uniform distribution on the corrs implies a marginal dist on each corr that is concentrated near 0. There's a conceptual problem with such a model: imagine you start with two variables, y1 and y2, with a corr rho_12 which has a uniform prior. Now we add more and more variables: y3, y4, y5, etc. If we go with the inverse-Wishart, we're fine, all these models are coherent (although they have other problems, as you've noted!). But if we go with the joint uniform, then adding y3, y4, y5, etc. has the effect of changing the prior on rho_12. Which doesn't seem right to me, somehow, that the addition of more variables into the analysis changes the model for an existing variable. This sort of paradox shoudn't stop us from possibly using the model, but I'd like to have a good way of thinking about it.
I don't think it is a problem with the prior; it is a problem with high dimensions. An earlier paper by Joe
http://scholar.google.com/scholar?cluster=8483085447913922429&hl=en&as_sdt=0,33
proves that the c(k,eta) term in the PDF when eta = 1 is
c(k) = c(k - 1) * 2^(k^2 - 2 * k + 1) * Beta(k/2,k/2)^(k - 1)
where Beta() is the Beta function. The hypervolume of correlation matrices is c(k), which eventually starts to decrease as Beta(k/2,k/2)^(k - 1) shrinks faster than 2^(k^2 - 2 * k + 1) increases:
k c(k)
2 2.0
3 4.934802
4 11.69731
5 22.53256
6 31.11388
7 27.85823
8 14.87740
9 4.411544
10 0.682269
So, given that correlation matrices must be positive definite, as k increases, a Bayesian robot should update its beliefs about the correlations to make them closer to zero. In other words, for large k, the set of correlation matrices is small relative to the set of merely symmetric matrices with ones on the diagonal, and the correlation matrices in that set tend to have small off-diagonals. But people are not Bayesian robots and can fall into the trap of thinking that their model stays "the same" as k increases, instead of conditioning their prior beliefs on the known k.
Ben
<covariance.png>
Here is a demonstration as to why I think scaled inverse Wishart priors are going to be problematic with HMC, relative to scaled LKJ priors. This is similar to some of the plots in that paper Andrew and I did, except here I am only doing 3 dimensions. Higher dimensions have the same intuition but are harder to visualize. In every case, I draw 10,000 times from a distribution, and plot the sampling distribution of some quantity. For inverse Wishart, I always use identity scale and 4 degrees of freedom. For LKJ, I always set the shape parameter equal to 1.0, so that the correlation matrix has a joint uniform distribution.
In the top left, I rescale the inverse Wishart covariance matrix into a correlation matrix and plot the three correlations. The main point here is that rarely would anyone have a prior belief that places so much mass in the corners of this constrained space. In contrast (top right), in many situations, people do not have strong prior beliefs about a correlation matrix, so a joint uniform accurately reflects their ignorance.
[sic, BG messed up the indices: It should be L_{21}, L_{31}, and L_{32}]
In the bottom left, I take the Cholesky decomposition of an inverse Wishart covariance matrix and plot L_{12}, L_{13}, and L_{23}. This is the unconstrained space (ignoring the diagonal) that stan would be drawing from, and its tails are scary. I had to blow up all the axes just to hold the occasional extreme points.
In the LKJ case, we do not actually draw from an unconstrained L, but in the bottom right, I plot L_{12}, L_{13}, and L_{23} anyway where L * L' is a jointly uniform correlation matrix. If we were drawing from this unconstrained space, it doesn't look like it would be a big problem, and this intuition carries over to the unconstrained space we are actually drawing from in the LKJ case.
To summarize, the popular scaled inverse Wishart parameterization
Sigma = S * W * S
where S is a diagonal matrix and W has an inverse Wishart distribution is similar in form to the unpopular scaled LKJ parameterization
Sigma = D * (L * L') * D
where D is a diagonal matrix of standard deviations and L * L' is a correlation matrix (often with a uniform distribution).
The main differences are:
1) D has well-defined standard deviations on its diagonal, whereas S has weird positive numbers that are not standard deviations and harder to form a marginal prior over.
2) Few understand the inverse Wishart distribution and fewer genuinely believe it reflects their priors over a covariance matrix, whereas many people can understand what it means to draw uniformally from the set of all correlation matrices and in many cases that reflects their well-founded prior ignorance over a correlation matrix.
3) The unconstrained space that corresponds to the constrained space for an inverse Wishart is computationally scary, whereas the unconstrained space --- the atanh() of symmetric betas on (-1,1) with shape parameter >= 1.0 --- that corresponds to the constrained space for a jointly uniform correlation matrix is pretty tame.
In short, with a scaled inverse Wishart prior, you have to get bailed out by the data or fiddle with the degrees of freedom parameter. But again, it is more straightforward to fiddle with the shape parameter for the LKJ prior: Starting from eta = 1.0, increasing eta puts more mass on the identity matrix and less mass in the corners of the constrained space.
Ben
<covariance.R>
Andrew Gelman
(replying inline to stan-dev so that everyone can follow)On Wed, Apr 4, 2012 at 11:17 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
BenThat's interesting, thanks for sending. The pictures are pretty. Usually I don't get much from 3d scatterplots but these really work!The inv-Wishart with d+1 df gives uniform marginals on the corrs, but I can see how, to get the uniform marginal [sic, AG meant joint], you have to pile up mass on the corners. One solution would be to increase the df, right?
Ben Goodrich
On Thursday, April 5, 2012 4:21:45 PM UTC-4, Andrew Gelman wrote:
See below, I think there is some confusion.
On Apr 5, 2012, at 2:01 AM, Ben Goodrich wrote:
(replying inline to stan-dev so that everyone can follow)
On Wed, Apr 4, 2012 at 11:17 PM, Andrew Gelman wrote:
BenThat's interesting, thanks for sending. The pictures are pretty. Usually I don't get much from 3d scatterplots but these really work!The inv-Wishart with d+1 df gives uniform marginals on the corrs, but I can see how, to get the uniform marginal [sic, AG meant joint], you have to pile up mass on the corners. One solution would be to increase the df, right?
Andrew Gelman
Yes. The expectation of the inverse Wishart is S / (d - k - 1) iff d > k + 1. So, I picked the popular but extreme case that the expectation doesn't exist. Calculating the expectation of an implied correlation is a bit involved, but as d increases, the standard deviation of the correlations shrinks.
library(MCMCpack)
k <- 3
I <- diag(k)
mark <- lower.tri(I, diag = FALSE)
REPS <- 10000
dof <- (k+1):20
results <- sapply(dof, FUN = function(d) replicate(REPS, cov2cor(riwish(d, I))[mark]) )
colnames(results) <- dof
plot(dof, apply(results, 2, sd), type = "l", xlab = "Degrees of freedom",
ylab = "Standard deviation of correlations")
However, I continue to think that increasing the shape parameter of a LKJ distribution (starting at eta = 1 -> joint uniform) is more intuitive than increasing the degrees of freedom of an inverse Wishart (starting at d = k + 1 -> no moments) in part because the "degrees of freedom" is kind of a nondescript number, while the LKJ correlation matrix is a composition of independent symmetric betas on (-1,1) whose shape parameters are a function of eta.
Most people don't have a good understanding of joint distributions, correlation matrices, etc. and even if they understand them, they don't have strong priors about correlation matrices, particularly if it is a correlation matrix among parameters rather than data. So, in that situation, a joint uniform distribution on a correlation matrix accurately reflects their prior beliefs.
And insofar as they have informative prior beliefs, they probably slightly favor the identity matrix over the nearly singular correlation matrices, which can be accomplished with an eta slightly greater than 1.0.
Ben Goodrich
I don't think this is quite right, but I'll have to think further about how to express my argument about increasing dimensionality.
Maybe another way of thinking about it is in terms of Stein's paradox. Suppose you have two standardized variables x and y. You could estimate the correlation between them with the MLE. Then suppose you get K - 2 more standardized variables, so you could form a K by K correlation matrix. Should you change your estimate of the correlation between x and y? Yes, no matter how arbitrarily chosen the other K - 2 variables are. Specifically, the MLE is dominated in expected aggregate loss by many estimators, most of which have the flavor of shrinking toward the mean / median / modal / common correlation estimate. Assuming what you shrink toward is zero, then you should favor something smaller in magnitude than the MLE correlation between x and y. That doesn't mean that it is a better estimate of the population correlation between x and y, even in expectation, but it is better in expected loss aggregated over all K choose 2 correlations. Now all that is for estimates, but if you know it is going to happen with estimates, then your prior should evolve in a similar fashion as K increases.
It seems that having a prior belief about the marginal correlation between any two variables isn't that helpful for forming a joint prior over K choose 2 correlations, even though the margin is the only thing that is easy to grasp. But for K > 3, the union of someone's marginal priors is usually inadmissible for a correlation matrix, so what good are the marginal priors? I suppose you could have a prior about the marginal correlation between x and y and condition on that when building up the joint prior. That can be handled in the LKJ framework but conditioning on more than one margin gets tedious and introduces the possibility that the joint is impossible.
Ben