sparse cholesky design

[paciorek]

Eigen computes sparse Cholesky based on a reordering to preserve as much sparsity as possible. This means it doesn't make sense to try to give the user the sparse Cholesky factor, either in compiled or uncompiled code.

Here is a draft design:

1. Provide a sparseCholFactor operator that returns a shared pointer to a simplicialLLT object by initializing the simplicialLLT (e.g., called llt) and then calling llt->compute(inputMat). Note that we will no overload chol(A) because the output is different when A is sparse versus dense.
2. In R, we will instead call Matrix::Cholesky which does a similar reordering (but I don't know it is identical).
3. Provide the following operators than can use the simplicialLLT:
   i. sparseCholeskySolve(b): this will call llt->solve(b) in Eigen code, solving the system of equations Ax=b.
   ii. sparseCholeskyForwardsolve(b): this will call llt->matrixL().triangularView<Eigen::Lower>().solve(b) to do $L^{-1}b$.
   iii. sparseCholeskyDiag(): this will call code that will extract the diag from matrixL and inverse permute. This will presumably mostly be used for determinants.

I have a draft of items 1 and 2 working on branch sparse_chol. And I have the Eigen code to do (3) from Claude.

@perrydv happy to have your inputs, but also happy to proceed on this path.

LDLT for sparse and dense will be future work (possibly post-initial release). I will want to experiment to better understand cases where LDLT works and LLT does not. Need to look at some GP cases where matrix is nearly non-p.d., such as Gaussian correlation or Matern with larger nu values.

[paciorek]

Oh, I guess we need L %*% x as well (e.g., if $A$ is a covariance), presumably via an operator sparseCholeskyLmult or some such.

[paciorek]

Ok, after more puzzling over what methods make sense given the permutation, I think they are:

• sparseCholLogdet
• sparseCholMult (e.g., for generating MVN with sparse covariance)
• sparseCholBacksolve (e.g., for generating MVN with sparse precision)
• sparseCholSolve (e.g., for quadratic forms with sparse covariance and more general systems of equations)

The result of sparseCholMult and sparseCholBacksolve will be the same in distribution as if using a dense matrix, but not the same exact values.

I have a draft of sparseCholLogdet working. I think others will be straightforward given that as a starting point.

[paciorek]

The current usage is like this:

ch <- sparseCholFactor(Q)
out <- sparseCholLogdet(ch)

Not sure if we want to provide ch$sparseCholLogdet(), which would require creating our own class (wrapping around the Eigen SimplicialLLT) for the sparse factor.

We'll also want to think about whether we want to overload operators to allow, e.g., solve(ch, x).

[perrydv]

A few quick thoughts:

• All the shared pointers we use in C++ are part of the nClass scheme. Potentially the natural way to do this for a special type like a sparse Cholesky result is to make a small nClass (probably predefined) that holds the custom type needed as a member variable. Then operators using it can have handlers that insert access to that member variable. These steps would allow such a result to be held in R (even if only to pass back into another compiled function), included in serialization, included in other types (nClass members or nList elements), and then used as if they are the underlying variable (with some handler magic).
• There is (coming, when I merge the big chunk I'm working on) some overloading resolution based on first argument type, which could be useful here.
• Compile-time arguments are another option for choosing syntax, like "chol(Q, sparse = TRUE)" where sparse is marked as a compile-time argument. This allows chol(Q, sparse = my_sparse_setting), where my_sparse_setting is evaluated in R at compile-time and so provides an easy switch for experimenting with (say) sparse vs non-sparse methods. I'm not sure how relevant or general that could get.
• My initial thought would have been to keep the the same operator names and let them default to making sparse results from sparse inputs unless one over-rides that. So ch <- sparse(Q), defaulting to sparse or dense according to Q, unless manually over-ridden by a compile-time argument.

Editing on some additional notes after looking at branch sparse_chol:

• I think for now (as you make code to get the math working) it would be better to pass and return the actual objects (not shared_ptrs). These are the same level of object as e.g. Eigen::Tensors and we don't move those around by shared_ptr.
• We may have a general need to allow first class types in C++ that we don't do anything with in R except hold pointers to. This could be used for the variety of Eigen math-related types. I can start to see how that could be fit in to our system somewhat cleanly.
• I guess we don't handle even a dense chol result as its native Eigen object, and that comes up naturally as a need as well.

[paciorek]

Ok, thanks for the extensive thoughts.

Yes, I see the appeal of being able to pass the sparse Cholesky result around easily.

Some additional considerations/concerns:

1. I tried to pass the SimplicialLLT objects around directly but got C++ compilation errors about
   "‘Eigen::SimplicialLLT<Eigen::SparseMatrix<double, 0, int> >::SimplicialLLT(Eigen::SimplicialLLT<Eigen::SparseMatrix<double, 0, int> >&&)’ is implicitly deleted because the default definition would be ill-formed:
   338 | class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering>"
   that Claude told me had to do with " SimplicialLLT cannot be moved (or copied) — its move constructor is implicitly deleted. This almost always means your code is trying to return a SimplicialLLT by value from a function". I'm assuming that if we had the nClass wrapper that the embedded simplicialLLT would need to be a reference.

2. My concern about using the same operator names is that chol(Q) in the dense case returns a matrix but a more complicated object in the sparse case. And that things like diag(ch) wouldn't then make sense in the sparse case, and that ch %*% x in the sparse case might be difficult to get to work in uncompiled operation (e.g., in R I think the way we would do it is to use Matrix operations with this syntax: solve(ch, expand1(ch, "L") %*% x, system = "Pt").

3. Regarding chol(Q) returning a matrix in the dense case, that's how Josh set it up, but given that further use of that matrix should use Eigen's special triangular matrix handling, it may not be what we want. So maybe for both dense and sparse we want to restrict the user to the limited subsequent operations I list in my previous comments. So that is consistent with your idea to use the same operator names and with my thinking that we need specific new operators that restrict users to only do sensible things with the Cholesky, disallowing things like diag(ch) and ch %*% x.

Thoughts? Feel free to defer until you are back; talking in person may ultimately be what we need to do here...

[paciorek]

Ok, a bit of tweaking of my comments. I guess that objects of the sparse chol wrapper class would be shared pointers.

I'm a bit confused about the simplicialLLT member data in the class. That object needs to actually live somewhere; I guess it would be a pointer to a simplicialLLT object that the C++ initializer for the class would initialize.

[paciorek]

Ok, I think the following may be about what we need.

# `nClass` code

sparseCholFactor <- nClass(
    classname = 'sparseCholFactor',
    predefined = quote(...),
    Cpublic = list(
        llt = 'simplicialLLT'
    )
)

// generated C++ class:
class sparseCholFactor {
    public:
      sparseCholFactor();
      Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> llt;
};

// permanent C++ code

std::shared_ptr<sparseCholFactor> chol(const Eigen::SparseMatrix<double>& x) {
    auto ch = std::make_shared<sparseCholFactor>();  // or use `new`?
    ch->llt.compute(x);
    if (ch->llt.info() != Eigen::Success)
        throw std::runtime_error("Cholesky factorization failed");
    return ch;
}

double sparseCholLogdet(const sparseCholFactor& ch) {
    Eigen::SparseMatrix<double> L = ch.llt.matrixL();
    return L.diagonal().array().log().sum();
}

// add other methods for solve,mult, etc.

# example nFunction code:

   fun = function(A = 'nSparseMatrix') {
      ch <- chol(A) 
      ld <- sparseCholLogdet(ch)
    }

// C++ code generated from nFunction:
             
std::shared_ptr<sparseCholFactor> ch = chol(A);
double ld = sparseCholLogdet(*ch);