Pivoted Cholesky

ldlt() in Eigen produces an object that has the following methods:

• ldlt().matrixL(): returns lower triangular factor
• ldlt().matrixU(): upper
• ldlt().vectorD(): diagonal elements (see below)
• ldlt().transpositionsP(): vector that implements permutation matrix P

$$ \Sigma = P^\top L D U P $$

C++ generation: should we translate tmp <- chol(matrix, pivot = TRUE) to have tmp be the Eigen result of ldlt() or tmp be a nimList that contains the Eigen result of ldlt() as well as $L, $U, $D, $P? If the former, then when a user does something like tmp$U we need to translate that to tmp.matrixU() as opposed to simple list operations. Or we could have the ldlt object be a nimbleFunction with methods such as backsolve, solve, %*%, etc (see below). Pulling out $U, $D, etc. might involve active bindings when done on the R side.

Difficulty is we don't currently allow result of a DSL call to be an Eigen type. So we could extend the DSL or take the ad hoc approach of returning a nimbleList with U,P,D elements and then when the pivoted chol is used in one of the limited ways it can be used, we map the elements back to Eigen types.

Cholesky methods: The natural things to do with the result of LDLT are as follows where U and L here are taken to be the result of the unpivoted Cholesky:

• backsolve(tmp, B): $$ U^{-1}B $$
• forwardsolve(tmp, B): $$ L^{-1}B $$
• solve(tmp, B): $$ {LU}^{-1}B $$
• tmp %*% B or crossprod(tmp, B): $$ LB $$

However, things are tricky because of the permutations. We can represent our decomposition as:

$$ \Sigma = P^\top L D^{1/2} P P^\top D^{1/2} U P $$ such that $$ P^\top L D^{1/2} P $$ is the matrix factor we are working with. Note that that is NOT triangular because it takes the triangular L and permutes rows and columns, but it is a legit factorization.

So the translation of our 4 core operations above are:

• backsolve(tmp, B): $$ P^\top U^{-1} P B $$
• forwardsolve(tmp, B): $$ P^\top L^{-1} P B $$
• solve(tmp, B): this is just ldlt().solve(B) so translation is simple
• tmp %*% B or crossprod(tmp, B): $$ P^\top L P B $$

I haven't fully figured out the Eigen syntax for these, but am a good chunk of the way there. In particular, the following implements forwardsolve, but is missing the diagonal matrix piece (see below):

myLDLT = A.ldlt();
myLDLT.transpositionsP().transpose() * myLDLT.matrixL().solve(myLDLT.transpositionsP() * B);

I'm still trying to figure out how to take vectorD (a vector) and do an efficient 'matrix' multiplication with D^{1/2} or D^{-1/2}. This presumably involves element-wise multiplication and therefore something like use of array() or cwiseProduct() but need to work on this more.