# Fastest algorithm to compute the condition number of a large matrix in Matlab/Octave

From the definition of condition number it seems that a matrix inversion is needed to compute it, I'm wondering if for a generic square matrix (or better if symmetric positive definite) is possible to exploit some matrix decomposition to compute the condition number in a faster way.

Computing the condition number (even approximating it within a factor of 2) seems to have the same complexity as computing a factorization, though there are no theorems in this direction.

From a sparse Cholesky factor $R$ of a symmetric positive definite matrix, or from a sparse $QR$ factorization (with implicit $Q$) of a general square matrix, one can obtain the condition number in the Frobenius norm by computing the sparse inverse subset of $(R^TR)^{-1}$, which is much faster than computing the full inverse. (Related to this is my paper: Hybrid norms and bounds for overdetermined linear systems, Linear Algebra Appl. 216 (1995), 257-266. http://www.mat.univie.ac.at/~neum/scan/74.pdf)

Edit: If $A=QR$ then with respect to any unitarily invariant norn, $$cond(A)=cond(R)=\sqrt{cond(R^TR)}.$$ For the computation of sparse QR factorizations see, e.g.,
http://dl.acm.org/citation.cfm?id=174408.
For the computation of the sparse inverse, see, e.g., my paper: Restricted maximum likelihood estimation of covariances in sparse linear models, Genetics Selection Evolution 30 (1998), 1-24.
https://www.mat.univie.ac.at/~neum/ms/reml.pdf The cost is about 3 times the cost for the factorization.

• So you are suggesting the following: Given a matrix $\mathbf{A}$ compute its QR of the form $\mathbf{A}=\mathbf{Q}\mathbf{R}$ where $\mathbf{R}$ is an upper triangular matrix and $\mathbf{Q}$ is an orthogonal matrix and then the condition number is given by $\textrm{cond}(\mathbf{A})=||A|| ||A^{-1}|| (\mathbf{R}^T\mathbf{R})^{-1}$ The point here is how to find a fast method to compute a QR factorization. Am I right? – linello Oct 22 '12 at 19:21
• @linello: not quite; see my edit. – Arnold Neumaier Oct 22 '12 at 19:58
• Thanks! I'm going to check it, btw what is the cost of this step? – linello Oct 22 '12 at 20:02
• @linello: For a full matrix, $O(n^3)$; for a sparse matrix, it depends a lot on the sparsity structure. – Arnold Neumaier Oct 22 '12 at 20:23

It's certainly easy to use the eigenvalue/eigenvector decomposition of a symmetric matrix or the SVD of a general matrix to compute the condition number, but these aren't particularly fast ways to proceed.

There are iterative algorithms that can compute an estimate of the condition number that is useful for most purposes without going to all of the work of computing $A^{-1}$. See for example the condest function in MATLAB.

• But the estimate is sometimes significantly too small. Computing the condition number (even approximating it within a factor of 2) seems to have the same complexity as computing a factorization, though there are no theorems in this direction. – Arnold Neumaier Oct 22 '12 at 14:27

For sparse Hermitian matrices $H$, you can use Lanczos algorithm to compute its eigenvalues. If $H$ is not Hermitian, you can compute its singular values by computing the eigenvalues of $H^TH$.

Since the largest and smallest eigenvalues/singularvalues can be found very fast (long before the tridiagonalization is complete), the Lanczos method is particular useful to compute the condition number.

• I've always wondered where to find a readable matlab code for lanczos iteration which clarifies how to get the smallest or largest eigenvalue. Can you suggest me one? – linello Oct 23 '12 at 8:16
• I don't have the MATLAB codes for Lanczos algorithm. – chaohuang Oct 23 '12 at 8:50