Naive question on handling high-dimensional cov/precision matrices

[Alp Kucukelbir]

Hi folks,

I am not accustomed to having to worry about numerical accuracy in my coding. However, I'm quickly realizing that I really should be.

I've learned about log_sum_exp trick, and I've been trying to compute things in log-space to avoid numerical underflow. (Any other general tips would be appreciated!)

My main question is: what do I need to know about high-dimensional symmetric matrices? I'd like to compute log determinants and the such. I'd also like to compute quadratic forms. (Presumably, there's a smart way to compute these without inverting the matrix?)

I'm seeing that Stan seems so use the LDLT Factor. Is this what I should read up on?

Thanks very much in advance,
Alp

P.S. This whole discussion began as I realized my algorithms were only stable until dimension = 5. I think this is because Eigen has closed-form implementations of things like determinants up to that dimension. (I may be wrong.) Because when I go to dimension = 6, I seem to lose a lot of numerical stability.

[Marcus Brubaker]

Hi Alp,

Yes, read up on Cholesky factorizations in particular, LDLT is mild generalization of Cholesky which sometimes has better stability.  Internally for things like multi_normal and multi_gp Stan is generally doing what is numerical best practice for log determinants, quadratic forms, etc so they're a good place to look for how to deal with things if you're writing C++ code.  If you have more specific questions definitely ask as, much like the log_sum_exp trick, some of these issues are not well explained or documented anywhere.

The Stan model API that is exposed to users for implementing these things themselves is somewhat limited and may be a bit less optimal both numerically and in terms of efficiency.  Part of that is simply due to a lack of time, although in most cases it's simply not trivial to expose certain types of functionality.  E.G., the LDLT factorization would be very hard to expose to users in an (autodiff) efficient way.

Cheers,

Marcus

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

If you are using a covariance matrix (or precision matrix) as a parameter in a model, then it is much better to reparameterize the model in terms of a Cholesky factor of a covariance matrix (or precision matrix). There are several Cholesky factorizations. The one I am referring to here is the lower-triangular matrix with strictly positive diagonal entries, denoted L, such that L * transpose(L) is equal to the covariance matrix (or precision matrix) of interest. This has as much numerical stability as anything in Stan. However, in most cases, substantive discussion of the model will involve the covariance matrix (or precision matrix) and thus in order to figure out your prior for L, you have to take your prior for the covariance matrix (or precision matrix) and multiply by a fairly simple Jacobian determinant.

If you are talking about a known covariance matrix or something internal that you tune in your algorithm, then it is slightly better to use the LDL' form of the Cholesky factorization, which expresses a covariance matrix in terms of L --- which is lower-triangular and has ones on its diagonal --- and D which is a diagonal matrix with strictly diagonal entries. This form is better for matrices that are singular or nearly singular because Eigen does pivoting to avoid dividing by (nearly) zero. So, really it is PLDL'P' where P is a permutation matrix and you may have to keep track of the P matrix. It isn't necessary to worry about the P matrix in some cases, like solving least squares problems.

The other matrix decompositions, such as QR, eigenvalues / vectors, etc. are either not as fast or not as stable as the Cholesky factorizations, so you wouldn't ordinarily need to use them.

Ben

[Alp Kucukelbir]

Thanks Marcus and Ben!

Both of your responses were very helpful. I will switch to LDL' form then.

In the meantime, if anyone else has any tricks like the "log_sum_exp", please feel free to share!

Cheers,
Alp

[Bob Carpenter]

There are lots of such "tricks" and I don't think we have
a list anywhere. Everything from using
log1p(x) instead of log(1 + x) for values of x near 0 to
making sure exp() doesn't overflow anywhere by pushing the
logs through with algebra, to preferring A \ B to inverse(A) * B.

- Bob

[Michael Betancourt]

Alp,

I think the best answer is to write up the equations that you’d like to implement and
then we can help you develop and implementation that’s most efficient and accurate.
As Bob notes there are lots of tricks so trying to address general strategies can be tricky.

[Alp Kucukelbir]

Thanks Bob and Michael,

I might take you up on your offer, Michael. But already, I'd say that this discussion has increased my awareness. Which is a good thing.

Cheers,
Alp

> You received this message because you are subscribed to a topic in the Google Groups "stan development mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-dev/zoLMKrrk4vU/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-dev+u...@googlegroups.com.

[Alp Kucukelbir]

Hi folks,

Thanks for the all the help.

Would you guys mind looking over the attached PDF? I think I've figured out how to do log abs determinants and quadratic forms. (Please do let me know if I got something wrong.) But I'm a bit confused by a trace that shows up in the inverse Wishart distribution.

Cheers,
Alp

[Alp Kucukelbir]

To follow up, I've taken a shot at how that trace should be calculated.

Still not %100 confident though...

Cheers,
Alp

[Marcus Brubaker]

Everything you have in there looks pretty good from a numerics perspective.  I should note that in some cases there are (slightly) faster ways to compute those expressions if you don't have to worry about autodiff but I'm not entirely clear what the context you're working in is.  If you're implementing another distribution in Stan, then following those patterns is pretty close to as good as you can do.

Cheers,

Marcus

[Michael Betancourt]

The quadratic form can be cleaned up a bit,

LLT<MatrixXd> L = H.selfadjointView<Eigen::Lower>().ldlt();

Eigen::VectorXd y = L.solve(x);
quad = x.squaredNorm();

The trace is a tricker one and I’m surprised more people don’t write it out more efficiently. Calling trace will be O(N^{3}) because of the intermediate matrix multiplication, but you can actually do it in O(N^{2}). In components,

Tr[A . B] = \sum_{ij} A_{ij} B_{ji} = A_{i-} B_{-i}

In other words, the trace is the inner product of the ith row of A with the ith column of B (or the other way around, it doesn’t matter) summed over i.

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <latex-image-1.pdf>

[Marcus Brubaker]

On Wed, Sep 3, 2014 at 1:58 PM, Michael Betancourt <betan...@gmail.com> wrote:

The quadratic form can be cleaned up a bit,

LLT<MatrixXd> L = H.selfadjointView<Eigen::Lower>().ldlt();

Eigen::VectorXd y = L.solve(x);
quad = x.squaredNorm();

This won't work if you're using ldlt(), as you need to divide by the square root of the diagonals of D. I'm guessing you meant to call llt() given the type of L but in that case L.solve() isn't computing what you mean to compute.  You want L.matrixL().solve().  We don't have LLT implemented for autodiff right now, so that's not really an option if autodiff is needed.

If using LDLT and dividing by the square root of the diagonals, it will be faster for simple double/float based arithmetic, but will be slower for autodiff because the square roots and divisions will add extra variables to the autodiff stack.  For autodiff with LDLT the quadratic form as it's in Alps notes is about as good as it gets with our code right now.

 

The trace is a tricker one and I'm surprised more people don't write it out more efficiently.  Calling trace will be O(N^{3}) because of the intermediate matrix multiplication, but you can actually do it in O(N^{2}).  In components,

Tr[A . B] = \sum_{ij} A_{ij} B_{ji} = A_{i-} B_{-i}

In other words, the trace is the inner product of the ith row of A with the ith column of B (or the other way around, it doesn't matter) summed over i.

It is written explicitly in the case of the inverse wishart because A is the inverse of Sigma so the trace can't easily be computed as a dot product without computing A explicitly.  Numerically it's more stable to compute A . B and then take its trace and it's basically as efficient.

Cheers,

Marcus

[Alp Kucukelbir]

Interesting. The subtleties are deep...

Another related question:

Let's say I have an LDLT decomposition of some matrix A.

I want to compute scalar * A.

Eigen doesn't appear to allow setting/adjusting the L and D matrices without re-computing everything.

Is there a trick for this that I'm not seeing?

Cheers,
Alp

> You received this message because you are subscribed to a topic in the Google Groups "stan development mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-dev/zoLMKrrk4vU/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-dev+u...@googlegroups.com.

[Marcus Brubaker]

The best solution there is probably to recognize that basically everything you want to compute using scalar * A can be computed using A and then rescaling it appropriately.  E.G., (scalar * A) \ b = (1.0/scalar) * (A \ b) or det(scalar * A) = scalar^N * det(A) where N is the size of A.

Alternatively, if you *really* need to update the LDLT, you simply need to rescale the entries of D but I don't know if Eigen's API allows this (matrixD() may be const) and I know for sure that the Stan wrapping of LDLT doesn't.

Cheers,

Marcus

[Bob Carpenter]

Can't you just pull out the L and D matrices and do it
externally? Or do you need it in the LDLT structure for
some reason?

- Bob

[Marcus Brubaker]

Yes, you could pull out L and D and update and use them externally, although I'd argue that that is less convenient and more error prone than simply using LDLT directly and rescaling the results appropriately.

[Alp Kucukelbir]

Marcus: I'd like to take a gradient step. So I think I need to update the LDLT, but perhaps there's some obvious workaround that I'm missing...

Bob: If I pull out L and D, Eigen doesn't seem to let me put the updated versions back in. I'd have to use them externally, as Marcus suggests, which is not ideal.

Cheers,
Alp

[Marcus Brubaker]

Take a gradient step in what variable?  A?  That is much more than just rescaling A in that case and you'll almost certainly have to recompute the LDLT factorization.

Cheers,

Marcus

Cheers,
Alp

[Alp Kucukelbir]

Sorry. I should be clearer. (BTW: I appreciate all the help guys -- this is immensely useful.)

I have A_old at iteration t-1. It is in LDLT form.

I have computed it's gradient. Call it G. I have a closed-form expression for it.

I want to compute A_new at iteration t in LDLT form.

Currently, it seems like I need to do something like

A_new = LDLT<MatrixXd>(A_old.reconstructedMatrix() + G);

I guess what I'm asking is: there's no "trick" to avoid that LDLT, right?

Cheers,
Alp

> You received this message because you are subscribed to a topic in the Google Groups "stan development mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-dev/zoLMKrrk4vU/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-dev+u...@googlegroups.com.

[Marcus Brubaker]

On Wed, Sep 3, 2014 at 3:12 PM, Alp Kucukelbir <akucu...@gmail.com> wrote:

Sorry. I should be clearer. (BTW: I appreciate all the help guys -- this is immensely useful.)

I have A_old at iteration t-1. It is in LDLT form.

I have computed it's gradient. Call it G. I have a closed-form expression for it.

I want to compute A_new at iteration t in LDLT form.

Currently, it seems like I need to do something like

A_new = LDLT<MatrixXd>(A_old.reconstructedMatrix() + G);

I guess what I'm asking is: there's no "trick" to avoid that LDLT, right?

That's right.  The only possible trick is if you somehow know that G is a rank one update, e.g., G = u u^T for some vector u.  In which case the LDLT factorization can be (relatively) cheaply updated.  However, for a general gradient this is unlikely to be the case.

Cheers,

Marcus

[Ben Goodrich]

Or a sequence of rank-one updates, which is still unlikely to be applicable. You don't want to call the reconstructedMatrix() method because that will literally do the multiplication. It is presumably worth the trouble to retain the matrix whose factorization is A_old and add G to that before factoring again.

Also, remember to keep track of the pivoting.

Ben

[Alp Kucukelbir]

Marcus: that makes sense. I came to a similar conclusion by looking at what Eigen offered for the LDLT class.

Ben: that's a really good point. However, I missed the subtlety about pivoting. Do I need to keep track of pivoting if I'm using "reconstrucedMatrix()" too?

Cheers,
Alp

[Ben Goodrich]

On Wednesday, September 3, 2014 9:12:22 PM UTC-4, Alp Kucukelbir wrote:

Marcus: that makes sense. I came to a similar conclusion by looking at what Eigen offered for the LDLT class.

Ben: that's a really good point. However, I missed the subtlety about pivoting. Do I need to keep track of pivoting if I'm using "reconstrucedMatrix()" too?

The reconstructedMatrix method reverses the pivoting so that it reconstructs the original matrix. But you already have the original matrix, so there is no reason to use that method, which mostly exists for unit testing.

Ben

[Alp Kucukelbir]

Ben: OK got it. So I guess I only need to worry about pivoting if I'm somehow keeping a copy of L outside of Eigen's LDLT class.

Great.

Thanks!
Alp