Condition numbers of sparse matrices

Contributor & Development Discussion scipy
1

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.)

2

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.

3

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?

4

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.

5

Just for reference, here is the implementation that I am using to compute sparse condition estimate (based on splu and onenormest). But I agree that using the native SuperLU API is a better solution for SciPy.

import time
import numpy as np
import numpy.linalg as lna
import scipy.sparse as sps
import scipy.sparse.linalg as sla


def condest(A, splu_opt={}, onenormest_opt={}):
    """
    Compute an estimate of the 1-norm condition number of a sparse matrix.

    Parameters
    ----------
    A : (M, M) sparse matrix
        square matrix to be inverted
    splu_opt : dict, optional
        Additional named arguments to `splu`.
    onenormest_opt : dict, optional
        Additional named arguments to `onenormest`.

    Returns
    -------
    c : {float, inf}
        The condition number of the matrix. May be infinite.

    References
    ----------
    .. [1] Nicholas J. Higham and Francoise Tisseur (2000),
           "A Block Algorithm for Matrix 1-Norm Estimation,
           with an Application to 1-Norm Pseudospectra."
           SIAM J. Matrix Anal. Appl. Vol. 21, No. 4, pp. 1185-1201.

    .. [2] William W. Hager (1984), "Condition Estimates."
           SIAM J. Sci. Stat. Comput. Vol. 5, No. 2, pp. 311-316.

    Examples
    --------
    >>> from numpy.linalg import cond
    >>> from scipy.sparse import csc_matrix
    >>> A = csc_matrix([[1., 0., 0.], [5., 8., 2.], [0., -1., 0.]], dtype=float)
    >>> A.toarray()
    array([[ 1.,  0.,  0.],
           [ 5.,  8.,  2.],
           [ 0., -1.,  0.]])
    >>> condest(A)
    45.0
    >>> cond(A.toarray(), p=1)
    45.0
    """

    # Check the input sparse matrix.
    if A.ndim != 2:
        raise ValueError('expected the matrix to have to dimensions')
    if A.shape[0] != A.shape[1]:
        raise ValueError('expected the matrix to be square')
    if A.shape[0] == 0:
        raise ValueError('cond is not defined on empty arrays')

    # Get LU decomposition of the matrix.
    try:
        decomposition = sla.splu(A, **splu_opt)
    except RuntimeError:
        return np.inf

    # Function for solving the equation system (original matrix).
    def matvec(rhs): return decomposition.solve(rhs, trans="N")

    # Function for solving the equation system (Hermitian matrix).
    def rmatvec(rhs): return decomposition.solve(rhs, trans="H")

    # Create a linear operator for the matrix inverse.
    op = sla.LinearOperator(A.shape, matvec=matvec, rmatvec=rmatvec)

    # Compute the 1-norm of the matrix inverse (estimate).
    nrm_inv = sla.onenormest(op, **onenormest_opt)

    # Compute the 1-norm of the matrix (estimate).
    nrm_ori = sla.onenormest(A, **onenormest_opt)

    # Compute an estimate of the condition number.
    c = nrm_ori*nrm_inv

    return c


if __name__ == "__main__":
    # test matrix size
    n_mat = 5000

    # number of nonzero diagonals
    n_diag = 5

    # create a random test matrix
    mat = sps.csc_matrix((n_mat,  n_mat))
    for idx in range(n_diag):
        vec = np.random.rand(n_mat-idx)
        mat += sps.diags_array(vec, offsets=+idx)
        mat += sps.diags_array(vec, offsets=-idx)

    # evaluate the sparse condition estimate
    timestamp = time.time()
    c = condest(mat.tocsc())
    print("sparse / cond = %f / time = %f seconds" % (c, time.time()-timestamp))

    # evaluate the dense condition number
    c = lna.cond(mat.todense(), p=1)
    print("dense / cond = %f / time = %f seconds" % (c, time.time()-timestamp))