Cholmod

[dpo]

I'm very sorry for adding this now. I know we've discussed it and I thought we wouldn't need it. I didn't realize that Scipy no longer supplies an interface to Cholmod. There is scikits.sparse but I had to hack it to make it compile. I pushed my changes to my fork.

Since the plan is to interface SPQR and that it uses the Cholmod data structure, let's also have a complete Cholmod interface. I happen to need a sparse Cholesky factorization in the class I'm currently teaching.

[ghost]

Not a problem. Will be done by end of June.

[dpo]

I wasn't aware of this, but Cholmod has an option to compute the LDL' factorization of a symmetric matrix A that isn't positive definite when possible. We make heavy use of that factorization in Sylvain's code an elsewhere. I became aware of this because Julia appears to have that feature, and uses Cholmod.

Fortunately, scikits.sparse also supports this feature:

In [91]: A = np.random.random((3,3)); A = np.dot(A.T, A);
In [92]: B = np.random.random((2,3))
In [93]: C = np.random.random((2,2)); C = -np.dot(C.T, C)
In [94]: K = np.zeros((5,5))
In [95]: K[:3,:3] = A; K[3:,:3] = B; K[:3,3:] = B.T; K[3:,3:] = C;  # K is SQD!
In [96]: Ksp = sp.sparse.csc_matrix(K)
In [97]: factor = cholesky(Ksp)
In [98]: L = factor.L_D()[0].toarray()
In [99]: D = factor.L_D()[1].toarray()
In [101]: np.linalg.norm(np.dot(L, np.dot(D, L.T)) - K)
Out[101]: 1.6319239955030595e-15

In [102]: factor(Ksp * np.ones(5))
Out[102]: array([ 1.,  1.,  1.,  1.,  1.])

[dpo]

Thought/question: Should we package Suite Sparse inside CySparse? We could have an option to build against the user's Suite Sparse, but build against the one we ship by default. @counterclocker @syarra

[ghost]

Yes, why not? Most of the packages if not all are under lgpl-2.1 or higher.

https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html

As you wish.

[ghost]

Will be implemented by end July...

[ghost]

End of September... ?

[ghost]

Nope, done before September! But needs testing and I didn't add the complex case yet...