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