5

scipy.linalg.solve, in its newer versions, has a parameter assume_a that can be used to specify that the matrix $A$ is symmetric or positive definite; in these cases, LDL or Cholesky are used rather than LU (Lapack's sysv and posv rather than gesv).

Is there a similar interface for sparse solvers? As far as I understand, scipy.sparse.linalg.spsolve does not support assume_a and always uses LU. What is the recommended way to use a symmetric sparse direct solver in Scipy, then (if there is any at all)?

I have seen that there is also sksparse.cholmod, but it is a separate package with a different interface, and from the documentation it looks like it does not handle indefinite matrices at all.

Federico Poloni
  • 11,344
  • 1
  • 31
  • 59
  • For Cholesky, $A$ must of course be nonnegative definite, eigenvalues $\ge 0$. E.g. sksparse.cholmod.cholesky( A, beta=1e-6 ) does $A + 10^{-6} I \to L L'$. Is that what you want ? – denis Jan 12 '20 at 18:03
  • @denis I was hoping to find also something for truly indefinite matrices like LDL (something that can solve a system with $\begin{bmatrix}0 & 1 \ 1 & 0\end{bmatrix}$, for instance). And I was hoping to find it in scipy rather than in another package, but that may be asking too much. :) – Federico Poloni Jan 12 '20 at 18:10
  • Would you know of test cases for such problems, with matrices, on the web ? Thanks – denis Jan 15 '20 at 14:26
  • @denis You can find a lot of them if you select "numerical symmetry property: symmetric indefinite" in the Matrix market search tool. For instance these ones. – Federico Poloni Jan 15 '20 at 15:23
  • Would you know of problems for generalized eigenvalues, pairs $A, M$ ? Matrixmarket is from 1999 / 2004, a long time ago (Moore's laws of cpus, memory, software, maybe eigenvalue algorithms too). Thanks – denis Jan 24 '20 at 15:10
  • @denis I am only familiar with Matrix market and with this collection, so I cannot help you much more with that. – Federico Poloni Jan 24 '20 at 15:52
  • @denis Anyway, your question looks like it should be a separate question here on this site, not a comment, if it does not refer specifically to my question. – Federico Poloni Jan 24 '20 at 15:53

0 Answers0