Cholesky decomposition for sparse matrices

[sagetrac-r-gaia-cs]

The Cholesky decomposition are implemented in the file sage/matrix/matrix2.pyx for some subfield of the algebraic numbers.

In this implementation the base ring must be exact or, for numerical work, a
matrix with a base ring of RDF or CDF must be used.

For the numerical work it's used the Cholesky decomposition implemented in sage/matrix/matrix_double_dense.pyx and because of this a error raised when try to compute the numerical Cholesky decomposition of a sparse matrix.

sage: A = matrix(QQ, [[1, 1], [1, 2]]) 
sage: A.cholesky()                    
[1 0]
[1 1]
sage: A = matrix(QQ, [[1, 1], [1, 2]], sparse=True)
sage: A.cholesky()                                 
[1 0]
[1 1]
sage: A = matrix(RDF, [[1, 1], [1, 2]], sparse=True)  
sage: A = matrix(RDF, [[1, 1], [1, 2]])             
sage: A.cholesky()                                  
[1.0 0.0]
[1.0 1.0]
sage: A = matrix(RDF, [[1, 1], [1, 2]], sparse=True)
sage: A.cholesky()                                  
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/raniere/opt/sage/devel/sage-rcm/sage/matrix/<ipython console> in <module>()

/home/raniere/opt/sage/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.cholesky (sage/matrix/matrix2.c:47738)()
   9867             if not self.base_ring().is_exact():
   9868                 msg = 'base ring of the matrix must be exact, not {0}'
-> 9869                 raise TypeError(msg.format(self.base_ring()))
   9870             try:
   9871                 posdef = self.is_positive_definite()

TypeError: base ring of the matrix must be exact, not Real Double Field

For this solve this ticket the numerical sparce Cholesky decompostion need to be implemented.

For more information about this topic see https://groups.google.com/forum/?fromgroups=#!topic/sage-support/do55Fayur6U.

CC: @orlitzky @collares

Component: linear algebra

Keywords: matrix, decomposition, cholesky, sparse

Author: Siddharth Bhat, Michael Orlitzky

Branch: 954b9ba

Reviewer: Dima Pasechnik

Issue created by migration from https://trac.sagemath.org/ticket/13674

[bollu]

comment:5

I have a solution implemented at https://github.com/bollu/sage/tree/u/gh-bollu/jan-24-cholesky-for-sparse-numerical-matrices.

[bollu]

Branch: u/gh-bollu/jan-24-cholesky-for-sparse-numerical-matrices

[bollu]

Author: Siddharth Bhat

[bollu]

Commit: 4351684

[bollu]

New commits:

4351684

Trac #13674: implement Chokesly factorization for sparse numerical matrices

[bollu]

comment:8

Things I don't understand very well:

• What's the correct solution for dense matrices? while the above code works, there should be a fast way to supply an entire dense matrix to cvxopt? The best idea I had was to use cvxopt.matrix(sagemat.to_numpy()); perhaps there are faster methods I am unaware of?
• Similarly, when building the sparse matrix, I coerce the matrix entries into a float. What is the correct way to hand the matrix entries off to cvxopt?
• On the google groups thread(https://groups.google.com/g/cvxopt/c/xQ-lR9ESijg/discussion), there was some discussion about also getting a good permutation for cholesky. I'd like to submit this patch next once I figure the API out. Any pointers to this would be appreciated!

[dimpase]

Reviewer: Dima Pasechnik

[orlitzky]

comment:10

Ohhkay. So, the title of this ticket is a bit misleading. It's not a cholesky for sparse matrices that's missing, but rather a cholesky decomposition for matrices over inexact rings. For example,

sage: A = matrix(RR, [[1, 1], [1, 2]])                                          
sage: A.cholesky()
...
TypeError: base ring of the matrix must be exact, not Real Field with 53 bits of precision

The use of RealField(100) or any similar ring produces the same result. The problem only appears related to sparse matrices because there is a special case implemented for dense matrices over RDF, but not sparse matrices over RDF. Thus, over RDF, the sparse implementation is indeed just "missing" (it could use the dense implementation without hurting anything).

But regardless, the real problem to be solved here is that we need a numerically-stable cholesky factorization that works for most inexact rings. That implementation can then be used (in the meantime) for both the dense and sparse methods, until someone decides to come along and speed it up in the sparse case. I've essentially already done this on #10332, which adds a numerically-stable block_ldlt() method for all Hermitian matrices. When your matrix is positive-definite, the block-LDLT factorization can be turned into a cholesky factorization after taking the square root of D.

So my suggestion for this ticket is to use block_ldlt() from the other ticket instead of adding a special case that relies on cvxopt. By using block_ldlt(), you allow the cholesky() method to work on fields like RR that cvxopt doesn't know about. It should also simplify the code a little bit, at the price of some administrative overhead: you'll have to,

1. Rebase your git branch onto u/mjo/ticket/10332
2. Make this ticket depend on isPositiveSemiDefinite not accessible #10332
3. Update the patch to use block_ldlt() instead of calling cvxopt.

I think (3) is pretty self-explanatory after reading the docs for block_ldlt(), but if not, I'm happy to help. E.g.

sage: A = matrix(QQ, [[9, 0, 3], [0, 1, 0], [3, 0, 25/4]])                      
sage: A.is_positive_definite()                                                  
True
sage: P,L,D = A.block_ldlt()                                                    
sage: L*D.apply_map(sqrt)                                                       
[           3            0            0]
[           0            1            0]
[           1            0 1/2*sqrt(21)]
sage: A.cholesky()                                                              
[                 3                  0                  0]
[                 0                  1                  0]
[                 1                  0 2.291287847477920?]

You might be able to save a few microseconds by using the internal _block_ldlt() method instead of the nice public one as well.

[bollu]

comment:11

Thanks a lot for the link! My only concern is with that of sparseness; the reason I choose to jump hoops with cvxopt is for the performance wins for large sparse matrices, that I get from running discrete differential geometry algorithms such as geodesics in heat: https://arxiv.org/abs/1204.6216

I'll be glad to benchmark the code, and see which method is faster; I think the route you propose is performant for dense inexact matrices for sure. I'm unsure about the sparse case. Could you please shed some light on the performance characteristics of the linked implementation?

[orlitzky]

comment:12

block_ldlt() will be pretty slow on large sparse matrices since the implementation scans the matrix (without regard for sparsity) looking for pivots. It's main benefit (with respect to cholmod) is that it would work for other fields, like sparse matrices over RR.

If you need the performance for sparse matrices over RDF then something like cholmod is indeed the way to go. I'm surprised we don't have a special matrix subclass for sparse RDF matrices already; my first attempt would be to override cholesky() there, similar to how cholesky() is overridden in matrix_double_dense.pyx for dense matrices.

[orlitzky]

comment:15

Replying to @dimpase:

In general, I'd be looking at wrapping up arb's cholesky and LDLs,
as something more robust than a naive attempt to implement this stuff.

For things like cholesky over RR this would likely be better than a naive implementation based on block_ldlt(), but the big problem I was solving by reimplementing block-LDLT was to gain a factorization that works on indefinite matrices.

[mkoeppe]

comment:16

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[orlitzky]

comment:17

Ticket #31619 allows cholesky() and is_positive_definite() to work over inexact rings as promised, but it will be slower than necessary for sparse RDF matrices. A matrix subclass that delegates to cholmod() in that case is the last missing piece.

[mkoeppe]

comment:18

Setting a new milestone for this ticket based on a cursory review.

[orlitzky]

comment:19

cvxopt is now a pseudo-optional package, but we should still be able to use the approach in comment:14. We can super() if cvxopt isn't available.

[orlitzky]

Changed author from Siddharth Bhat to Siddharth Bhat, Michael Orlitzky

[orlitzky]

Changed commit from 4351684 to 954b9ba

[orlitzky]

comment:20

This should work for both RDF and CDF, but it would be nice to have some more serious examples for test cases. The cvxopt interface is a little sketchy so I'd like to be sure that we're using it "correctly," insofar as is possible for an undocumented interface.

[orlitzky]

Changed branch from u/gh-bollu/jan-24-cholesky-for-sparse-numerical-matrices to u/mjo/ticket/13674

[dimpase]

comment:21

lgtm

[vbraun]

Changed branch from u/mjo/ticket/13674 to 954b9ba

[collares]

comment:23

I am seeing failures such as the one below on aarch64:

sage -t --long --random-seed=138452687149883420730489596915102319785 /nix/store/1jyscb1slmz6134mlsfs9gfjs4kv8w8i-sage-src-9.5.beta8/src/sage/matrix/matrix_double_sparse.pyx
**********************************************************************
File "/nix/store/1jyscb1slmz6134mlsfs9gfjs4kv8w8i-sage-src-9.5.beta8/src/sage/matrix/matrix_double_sparse.pyx", line 95, in sage.matrix.matrix_double_sparse.Matrix_double_sparse.cholesky
Failed example:
    L = A.cholesky()
Exception raised:
    Traceback (most recent call last):
      File "/nix/store/vwd2z6p52kzhidwwvwavgw9jxp1165qh-python3-3.9.6-env/lib/python3.9/site-packages/sage/doctest/forker.py", line 694, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/nix/store/vwd2z6p52kzhidwwvwavgw9jxp1165qh-python3-3.9.6-env/lib/python3.9/site-packages/sage/doctest/forker.py", line 1096, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.matrix.matrix_double_sparse.Matrix_double_sparse.cholesky[20]>", line 1, in <module>
        L = A.cholesky()
      File "sage/matrix/matrix_double_sparse.pyx", line 110, in sage.matrix.matrix_double_sparse.Matrix_double_sparse.cholesky (build/cythonized/sage/matrix/matrix_double_sparse.c:2820)
        raise ValueError("matrix is not Hermitian")
    ValueError: matrix is not Hermitian

[collares]

Changed commit from 954b9ba to none

[dimpase]

comment:24

@collares - please open a new ticket for this.

[collares]

comment:25

Opened #33023 for the test failure, thanks for letting me know about the correct procedure. I didn't know the "Commit" field would be cleared when I posted my first comment, sorry about that.

[orlitzky]

comment:26

Replying to @collares:

Opened #33023 for the test failure, thanks for letting me know about the correct procedure. I didn't know the "Commit" field would be cleared when I posted my first comment, sorry about that.

Thanks, the commit field thing is no big deal, that always happens. I'll see if I can reproduce the problem. The matrix is Hermitian by construction (unless I've made some typo I can't see) so it should be interesting.