Cholesky for sparse matrix

[Raniere Gaia Silva]

Hi,
I need to performe a (numerical) cholesky factorization of a sparse
matrix (20x20) but I'm getting a error.
I look at this thread
https://groups.google.com/forum/?fromgroups=#!searchin/sage-devel/cholesky/sage-devel/AW4pmKx49H4/7iuet3rWYQgJ
but it didn't help very much (and is a little old).

Bellow a minimal example of the problem.

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

/home/raniere/documents/cnpq_126874_2012-3/src_sage/<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:47860)()

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

Thanks,
Raniere

[Dima Pasechnik]

On 2012-10-29, Raniere Gaia Silva <r.ga...@gmail.com> wrote:
> Hi,
> I need to performe a (numerical) cholesky factorization of a sparse
> matrix (20x20) but I'm getting a error.

the error message is clear: there is no implementation for RDF matrices,
only for exact matrices. So you can do this for QQ, for instance.

This is of course a deficiency in Sage matrices. There are Cholesky
implementations available for RDF matrices from numpy, for instance,
it's just absense of proper hooks in Sage.
You can do the following as a workaround:

sage: import numpy as np

sage: D = matrix(RDF, 2, 2, [[1, 1],[1,2]])

sage: np.linalg.cholesky(D.numpy())
array([[ 1., 0.],
[ 1., 1.]])
sage:

you can get Sage matrix out of it, too:
sage: L=matrix(np.linalg.cholesky(D.numpy())); L
[1.0 0.0]
[1.0 1.0]

HTH,
Dmitrii

[Jason Grout]

I think the point was that cholesky was not implemented for *sparse*
matrices. Dense RDF matrices in Sage do have cholesky implemented (and
that calls scipy).

Of course, converting sparse to dense RDF and using cholesky on the
dense matrix works if the matrix is small enough. Scipy had some
interest in adding sparse cholesky [1], but it seems that nothing came
of it.

On the other hand, it would be cool if we wrapped something like CHOLMOD
[2] to compute sparse cholesky decompositions. Even cooler would be
supporting the SuiteSparse package [3] and its various algorithms.
Sparse matrices in Sage need some attention.

Thanks,

Jason

[1] http://projects.scipy.org/scipy/ticket/261
[2] http://www.cise.ufl.edu/research/sparse/cholmod/
[3] http://www.cise.ufl.edu/research/sparse/SuiteSparse/

[Jason Grout]

On 10/30/12 6:54 AM, Jason Grout wrote:
> Scipy had some interest in adding sparse cholesky [1], but it seems that
> nothing came of it.
>
> On the other hand, it would be cool if we wrapped something like CHOLMOD
> [2] to compute sparse cholesky decompositions.

It looks like the scikits-sparse package wrapped CHOLMOD to provide a
sparse Cholesky decomposition: https://github.com/njsmith/scikits-sparse

This scikit isn't in Sage. However, you can install it into Sage using
these instructions:
http://packages.python.org/scikits.sparse/overview.html#introduction.

You could then use it like this:
http://packages.python.org/scikits.sparse/cholmod.html

Thanks,

Jason

[Raniere Gaia Silva]

Thanks all reply.

I really like to make a little improve to Sage can handle the sparse
matrix. Should I open a ticket?

Raniere

[Dima Pasechnik]

ok, I might have misread the post in question.

>
> Of course, converting sparse to dense RDF and using cholesky on the
> dense matrix works if the matrix is small enough. Scipy had some
> interest in adding sparse cholesky [1], but it seems that nothing came
> of it.
>
> On the other hand, it would be cool if we wrapped something like CHOLMOD
> [2] to compute sparse cholesky decompositions. Even cooler would be
> supporting the SuiteSparse package [3] and its various algorithms.
> Sparse matrices in Sage need some attention.

cvxopt bundles and wraps SuiteSparce.
More precisely: cholmod:
http://www.cise.ufl.edu/research/sparse/cholmod/
see
http://abel.ee.ucla.edu/cvxopt/userguide/spsolvers.html?highlight=cholmod#cvxopt.cholmod.numeric

Let me see if I can get this to work in Sage...

Dima

[Dima Pasechnik]

On 2012-10-29, Raniere Gaia Silva <r.ga...@gmail.com> wrote:

> Hi,
> I need to performe a (numerical) cholesky factorization of a sparse
> matrix (20x20) but I'm getting a error.

OK, my previous reply didn't make much sense, so here is one that is
meaningful. In cvxopt (a standard Sage spkg) there is an interface to sparse Cholesky
factorization.
There aren't too many hooks for cvxopt in Sage, and, moreover, cvxopt
does not provide a straightforward interface, you need to solve matrix
linear equations to get L, see
http://abel.ee.ucla.edu/cvxopt/userguide/spsolvers.html#positive-definite-linear-equations

so it's a bit of a walk:

sage: preparser(False) # somehow, preparser breaks cvxopt :(
sage: import cvxopt
sage: from cvxopt import cholmod
sage: A=cvxopt.sparse(cvxopt.matrix([[1,1],[1,2]])); A
<2x2 sparse matrix, tc='d', nnz=4>
sage: print A
[ 1.00e+00 1.00e+00]
[ 1.00e+00 2.00e+00]
sage: F = cholmod.symbolic(A)
sage: cholmod.numeric(A,F)
sage: Lt=cvxopt.cholmod.spsolve(F,A,sys=4); Lt
<2x2 sparse matrix, tc='d', nnz=3>
sage: print Lt
[ 1.00e+00 1.00e+00]
[ 0 1.00e+00]

Achtung! This only works in the case the decomposition PAP^T=LL^T
computed by cholmod (and stored in F), where P is a permutation
matrix, has P=identity. More generally, you need to be more careful:
either cvxopt.cholmod.symbolic() must be called with the parameter
p specified to be the identity (permutation) matrix, or it can (perhaps)
be retrieved using cvxopt.cholmod.spsolve() with appropriate B and sys=...

Well, perhaps you actually do not need to know L explicitly, just to
solve these matrix equations...

By the way, you can extract entries of Lt using syntax like this:
sage: Lt[1,0]
0.0
sage: Lt[1,1]
1.0

note that
sage: type(Lt[1,0])
<type 'float'>

so they can be converted back to Sage types.

HTH,
Dmitrii

PS. cvxopt has its own google group:
https://groups.google.com/forum/?fromgroups=#!forum/cvxopt

[Dima Pasechnik]

please see this https://groups.google.com/d/topic/cvxopt/xQ-lR9ESijg/discussion

far a way to make it more sane. It turns out that cvxopt has un undocumented function to 

extract the Cholesky factor directly.

[Dima Pasechnik]

Hi,

after a discussion here: https://groups.google.com/d/topic/cvxopt/xQ-lR9ESijg/discussion

it turns out that there is enough CHOLMOD interface already in Sage, it's just not documented properly.

Should I volunteer to get sparse  Cholesky implemented for non-exact sparse Sage matrices, via this route?

Dima

[Jason Grout]

On 10/31/12 11:53 PM, Dima Pasechnik wrote:

> Should I volunteer

+1. Always +1 to that question. :)

Jason

[Raniere Gaia Silva]

Dima,
I never work with CHOLMOD cvxopt directily, but if there is any thing that I can do to help you to get sparse Cholesky implemented for non-exact sparse Sage matrices just say.

Raniere

Jason


--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To post to this group, send email to sage-s...@googlegroups.com.
To unsubscribe from this group, send email to sage-support+unsubscribe@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support?hl=en.