Condition numbers of sparse matrices

Hi all,

I’d like to add a feature to estimate the reciprocal condition number of sparse matrices

RCOND = 1 / ( norm(A) * norm(inv(A)) )

as proposed by @loganbvh in #18969. As @loganbvh suggested, we can use SuperLU’s dgscon function for this, which takes a LU factorization of the sparse matrix as input arguments. It works with the 1- and infinity-norm and implements the algorithm by Hager and Higham also used in the Matlab CONDEST function, according to the SuperLU docs.

(Continuing below because new users can only put five links into one post.)

The SuperLU function gstrf to compute the LU factorization is already wrapped as scipy.sparse.linalg.splu and returns a scipy.sparse.linalg.SuperLU object (see _superluobject.h for the C definition). Currently SuperLU objects only expose a single member function solve to the Python side, implemented in _superluobject.c. Regarding the API for the new sparse condition number feature, @loganbvh suggested to implement another member function

scipy.sparse.linalg.SuperLU.estimate_rcond(A_norm: float, norm=np.inf [or 1]) -> float

which takes the norm of A (to be computed separately by the caller, e.g. using scipy.sparse.linalg.norm) and the information whether 1- or infinity-norms should be used for the condition number.

I would like to raise at least the following points to discuss:

  • Do you think a SuperLU member function is the right place? (I think yes.)
  • It’s called estimate_rcond because dgscon approximates the reciprocal condition number. We can think about computing the inverse and calling the function estimate_cond instead. I don’t know the traditions/conventions and the relevant numerical or floating-point aspects though.
  • I think it is a bit weird to ask users to supply the matrix norm A_norm themselves. In the code of dgscon, it is only used as an additional factor for the condition number. We might think about setting it to 1.0, renaming the function into something like estimate_rinvnorm (“reciprocal norm of the inverse”), and writing into the docs that the (reciprocal) condition number can be computed by multiplication/division with the norm of A.

What do you think?

Thanks for following up @maxaehle !

numpy computes the condition number in linalg.cond. For simplicity, I think we should to do the same for sparse arrays. The inverse condition number is the important one to detect singularity, I know, but API consistency is more important here in my opinion.