Estimate of the condition number

The condition number of square matrix A is defined as follows:

k(A) = ||A||·||A -1||

The condition number has the following meaning: if the machine precision equals ε, when solving a system of linear equations AX = b the order of relative error will be ε·k(A). Although the matrix condition number depends on the selected norm, if the matrix is well-conditioned, the condition number will be small in all norms, and if the matrix is ill-conditioned, the condition number will be large in all norms. Thus, the most convenient norm is usually selected. 1-norm, 2-norm and ∞-norm are used most frequently. They are as follows: Algorithm description

If A is given directly (i.e. we have matrix itself, not its triangular factorization), then algorithm has the following form:

1. 1-norm of A (or ∞-norm) is calculated using closed form expression with O(N 2) complexity
2. A is factorized using one of the triangular factorizations, O(N 3) complexity
3. norm of A -1 is estimated using iterative algorithm with O(N 2) complexity
4. k(A)=norm(A)·norm(inv(A)). Overall complexity: O(N 3)

If A is given by its truangular factorization (LU decomposition or Cholesky decomposition), we don't need to perform step (2). However, when performing step (1), we can't use closed form expression for matrix norm because it needs A itself, not its factorization. So following algorithm is used:

1. norm of A is estimated using iterative algorithm with O(N 2) complexity
2. norm of A -1 is estimated using iterative algorithm with O(N 2) complexity
3. k(A)=norm(A)·norm(inv(A)). Overall complexity: O(N 2)

A disadvantage of the iterative algorithm lies in the fact that it only obtains an estimate of the condition number, and a rather rough lower bound for that. An estimate is usually undersized by 5-10%, but sometimes the error is much bigger (during the numerical experiments, the estimate was lower than the condition number by 8% on average, and the maximum error was 87%). When the estimate is 10 times lower than a value itself, it is usually considered as an inaccurate estimate. Nevertheless, even such a precision may be sufficient.

Let's consider an example. Let the condition number of matrix A be equal to 10 4. This means that if the machine precision is 10 -15, the system of linear equations will be solved with precision 10 -11. Suppose that our estimate can vary by 90% and equals 10 3. This means that we can expect a higher precision: 10 -12. The difference is not very big because it is not the value of error that's important but its order. On average, 10% of the difference could be ignored since it is less than one decimal digit.

Subroutines

ALGLIB package contains a number of subroutines for condition number estimation:

• rmatrixrcond1
• rmatrixrcondinf
• cmatrixrcond1
• cmatrixrcondinf
• spdmatrixrcond
• hpdmatrixrcond
• rmatrixlurcond1
• rmatrixlurcondinf
• cmatrixlurcond1
• cmatrixlurcondinf
• spdmatrixcholeskyrcond
• hpdmatrixcholeskyrcond

Subroutine name depends on problem it solves and has the following form: MatrixTypeMatrix[DecompositionType]RCond[NormType].

• MatrixType determines type of matrix being processed and may be R (general real), C (general complex), SPD (real symmetric positive definite) or HPD (complex Hermitian positive definite)
• if subroutine accepts matrix itself, then DecompositionType is omitted. But if subroutine accepts triangular factorization of A, then factorization name as added to subroutine name (LU for general real/complex matrices, Cholesky - for SPD/HPD matrices)
• subroutine name ends with description of norm used to estimate condition number: 1 - for 1-norm, Inf - for ∞-norm. However, these norms are equal for SPD/HPD matrices, so SPD/HPD subroutines have no norm suffix in their names.

Example

Let's compare condition numbers of two matrices related to the polynomial interpolation: Vandermonde matrix and matrix of Chebyshev basis functions. We'll use equidistant grid on (-1,+1). For different N's we have:

```                 CONDITION NUMBERS
OF VANDERMONDE AND CHEBYSHEV INTERPOLATION MATRICES

VANDERMONDE   CHEBYSHEV
N      1-norm      1-norm
2         2,0         2,0
3         3,7         3,0
4        18,0         4,0
5        50,0         5,0
6       159,4         6,0
7       455,0         7,0
8      1460,5         8,0
9      4250,0         9,0
10     13625,2        10,0
11     40145,3        11,0
12    128360,3        12,2
13    381433,0        22,6
14   1216286,9       122,8
```

You may see that condition number of Vandermonde matrix grows very fast, but condition number of Chebyshev matrix grows significantly slower.

ALGLIB Project offers you two editions of ALGLIB:

ALGLIB Free Edition: delivered for free offers full set of numerical functionality extensive algorithmic optimizations no low level optimizations non-commercial license

ALGLIB Commercial Edition: flexible pricing offers full set of numerical functionality extensive algorithmic optimizations high performance (SMP, SIMD) commercial license with support plan

ALGLIB 3.15.0 for C++ C++ library.
Delivered with sources.
Monolithic design.
Extreme portability.
Editions:   FREE   COMMERCIAL

ALGLIB 3.15.0 for C# C# library with native kernels.
Delivered with sources.
VB.NET and IronPython wrappers.
Extreme portability.
Editions:   FREE   COMMERCIAL

ALGLIB 3.15.0 for Delphi Delphi wrapper around C core.
Delivered as precompiled binary.
Compatible with FreePascal.
Editions:   FREE   COMMERCIAL

ALGLIB 3.15.0 for CPython CPython wrapper around C core.
Delivered as precompiled binary.
Editions:   FREE   COMMERCIAL