This chapter provides methods for the solution of systems of simultaneous linear equations, and associated computations. It provides methods for

– matrix factorizations;
– solution of linear equations;
– estimating matrix condition numbers;
– computing error bounds for the solution of linear equations;
– matrix inversion;
– computing scaling factors to equilibrate a matrix.

Methods are provided for both real and complex data.

The methods in this chapter (F07 class) handle only dense and band matrices (not matrices with more specialized structures, or general sparse matrices).

The methods in this chapter have all been derived from the LAPACK project (see Anderson et al. (1999)). They have been designed to be efficient on a wide range of high-performance computers, without compromising efficiency on conventional serial machines.

Syntax

C#
public static class F07

Visual Basic (Declaration)
Public NotInheritable Class F07

Visual C++
public ref class F07 abstract sealed

F#
[<AbstractClassAttribute>] [<SealedAttribute>] type F07 = class end

Background to the Problems

This section is only a brief introduction to the numerical solution of systems of linear equations. Consult a standard textbook, for example Golub and Van Loan (1996) for a more thorough discussion.

Notation

We use the standard notation for a system of simultaneous linear equations:

Ax=b

(1)

where A is the coefficient matrix, b is the right-hand side, and x is the solution. A is assumed to be a square matrix of order n.

If there are several right-hand sides, we write

AX=B

(2)

where the columns of B are the individual right-hand sides, and the columns of X are the corresponding solutions.

We also use the following notation, both here and in the method documents:

x^	a computed solution to Ax=b, (which usually differs from the exact solution x because of round-off error)
r = b-Ax^	the residual corresponding to the computed solution x^
x∞= maxi xi	the ∞-norm of the vector x
x1 = ∑ j=1 n xj	the 1-norm of the vector x
A∞ = maxi ∑j aij	the ∞-norm of the matrix A
A1 = maxj ∑ i=1 n aij	the 1-norm of the matrix A
x	the vector with elements xi
A	the matrix with elements aij

Inequalities of the form A≤B are interpreted component-wise, that is aij≤bij for all i,j.

Matrix Factorizations

If A is upper or lower triangular, Ax=b can be solved by a straightforward process of backward or forward substitution.

Otherwise, the solution is obtained after first factorizing A, as follows.

General matrices (LU factorization with partial pivoting)

A=PLU

where P is a permutation matrix, L is lower-triangular with diagonal elements equal to 1, and U is upper-triangular; the permutation matrix P (which represents row interchanges) is needed to ensure numerical stability.

Symmetric positive-definite matrices (Cholesky factorization)

A=UTU or A=LLT

where U is upper triangular and L is lower triangular.

Symmetric indefinite matrices (Bunch–Kaufman factorization)

A = PUD UT PT or A = PLD LT PT

where P is a permutation matrix, U is upper triangular, L is lower triangular, and D is a block diagonal matrix with diagonal blocks of order 1 or 2; U and L have diagonal elements equal to 1, and have 2 by 2 unit matrices on the diagonal corresponding to the 2 by 2 blocks of D. The permutation matrix P (which represents symmetric row-and-column interchanges) and the 2 by 2 blocks in D are needed to ensure numerical stability. If A is in fact positive-definite, no interchanges are needed and the factorization reduces to A=UDUT or A=LDLT with diagonal D, which is simply a variant form of the Cholesky factorization.

Solution of Systems of Equations

Given one of the above matrix factorizations, it is straightforward to compute a solution to Ax=b by solving two subproblems, as shown below, first for y and then for x. Each subproblem consists essentially of solving a triangular system of equations by forward or backward substitution; the permutation matrix P and the block diagonal matrix D introduce only a little extra complication:

General matrices ( LU factorization)

Ly=PTb Ux=y

Symmetric positive-definite matrices (Cholesky factorization)

UTy=b Ux=y or Ly=b LTx=y

Symmetric indefinite matrices (Bunch–Kaufman factorization)

PUDy=b UTPTx=y or PLDy=b LTPTx=y

Sensitivity and Error Analysis

Normwise error bounds

Frequently, in practical problems the data A and b are not known exactly, and it is then important to understand how uncertainties or perturbations in the data can affect the solution.

If x is the exact solution to Ax=b, and x+δx is the exact solution to a perturbed problem A+δAx+δx=b+δb, then

δx x ≤κA δA A +δb b +⋯second-order terms

where κA is the condition number of A defined by

κA = A.A-1 .

(3)

In other words, relative errors in A or b may be amplified in x by a factor κA. [Estimating condition numbers] discusses how to compute or estimate κA.

Similar considerations apply when we study the effects of rounding errors introduced by computation in finite precision. The effects of rounding errors can be shown to be equivalent to perturbations in the original data, such that δA A and δb b are usually at most pnε, where ε is the machine precision and pn is an increasing function of n which is seldom larger than 10n (although in theory it can be as large as 2n-1).

In other words, the computed solution x^ is the exact solution of a linear system A+δAx^=b+δb which is close to the original system in a normwise sense.

Estimating condition numbers

The previous section has emphasised the usefulness of the quantity κA in understanding the sensitivity of the solution of Ax=b. To compute the value of κA from equation (3) is more expensive than solving Ax=b in the first place. Hence it is standard practice to estimate κA, in either the 1-norm or the ∞-norm, by a method which only requires On2 additional operations, assuming that a suitable factorization of A is available.

The method used in this chapter is Higham's modification of Hager's method (see Higham (1988)). It yields an estimate which is never larger than the true value, but which seldom falls short by more than a factor of 3 (although artificial examples can be constructed where it is much smaller). This is acceptable since it is the order of magnitude of κA which is important rather than its precise value.

Because κA is infinite if A is singular, the methods in this chapter actually return the reciprocal of κA.

Scaling and Equilibration

The condition of a matrix and hence the accuracy of the computed solution, may be improved by scaling; thus if D1 and D2 are diagonal matrices with positive diagonal elements, then

B = D1 A D2

is the scaled matrix. A general matrix is said to be equilibrated if it is scaled so that the lengths of its rows and columns have approximately equal magnitude. Similarly a general matrix is said to be row-equilibrated (column-equilibrated) if it is scaled so that the lengths of its rows (columns) have approximately equal magnitude. Note that row scaling can affect the choice of pivot when partial pivoting is used in the factorization of A.

A symmetric or Hermitian positive-definite matrix is said to be equilibrated if the diagonal elements are all approximately equal to unity.

For further information on scaling and equilibration see Section 3.5.2 of Golub and Van Loan (1996), Section 7.2, 7.3 and 9.8 of Higham (1988) and Section 5 of Chapter 4 of Wilkinson (1965).

Methods are provided to return the scaling factors that equilibrate a matrix for general, general band, symmetric and Hermitian positive-definite and symmetric and Hermitian positive-definite band matrices.

Componentwise error bounds

A disadvantage of normwise error bounds is that they do not reflect any special structure in the data A and b – that is, a pattern of elements which are known to be zero – and the bounds are dominated by the largest elements in the data.

Componentwise error bounds overcome these limitations. Instead of the normwise relative error, we can bound the relative error in each component of A and b:

maxijk δaij aij , δbk bk ≤ω

where the component-wise backward error bound ω is given by

ω= maxi ri A. x^+bi .

Methods are provided in this chapter which compute ω, and also compute a forward error bound which is sometimes much sharper than the normwise bound given earlier:

x- x^ ∞ x∞ ≤ A-1 . r ∞ x∞ .

Care is taken when computing this bound to allow for rounding errors in computing r. The norm A-1.r∞ is estimated cheaply (without computing A-1) by a modification of the method used to estimate κA.

Iterative refinement of the solution

If x^ is an approximate computed solution to Ax=b, and r is the corresponding residual, then a procedure for iterative refinement of x^ can be defined as follows, starting with x0=x^:

for i=0,1,…, until convergence

compute ri=b-Axi
solve Adi=ri
compute xi+1=xi+di

In F04 (not in this release), methods are provided which perform this procedure using additional precision to compute r, and are thus able to reduce the forward error to the level of machine precision.

The methods in this chapter do not use additional precision to compute r, and cannot guarantee a small forward error, but can guarantee a small backward error (except in rare cases when A is very ill-conditioned, or when A and x are sparse in such a way that A.x has a zero or very small component). The iterations continue until the backward error has been reduced as much as possible; usually only one iteration is needed.

Matrix Inversion

It is seldom necessary to compute an explicit inverse of a matrix. In particular, do not attempt to solve Ax=b by first computing A-1 and then forming the matrix-vector product x=A-1b; the procedure described in [Solution of Systems of Equations] is more efficient and more accurate.

However, methods are provided for the rare occasions when an inverse is needed, using one of the factorizations described in [Matrix Factorizations].

Packed Storage

Methods which handle symmetric matrices are usually designed so that they use either the upper or lower triangle of the matrix; it is not necessary to store the whole matrix. If the upper or lower triangle is stored conventionally in the upper or lower triangle of a two-dimensional array, the remaining elements of the array can be used to store other useful data. However, that is not always convenient, and if it is important to economize on storage, the upper or lower triangle can be stored in a one-dimensional array of length nn+1/2; in other words, the storage is almost halved.

This storage format is referred to as packed storage; it is described in [Packed storage]. It may also be used for triangular matrices.

Methods designed for packed storage perform the same number of arithmetic operations as methods which use conventional storage, but they are usually less efficient, especially on high-performance computers, so there is then a trade-off between storage and efficiency.

Band and Tridiagonal Matrices

A band matrix is one whose nonzero elements are confined to a relatively small number of subdiagonals or superdiagonals on either side of the main diagonal. A tridiagonal matrix is a special case of a band matrix with just one subdiagonal and one superdiagonal. Algorithms can take advantage of bandedness to reduce the amount of work and storage required. The storage scheme used for band matrices is described in [Band storage].

The LU factorization for general matrices, and the Cholesky factorization for symmetric and Hermitian positive-definite matrices both preserve bandedness. Hence methods are provided which take advantage of the band structure when solving systems of linear equations.

The Cholesky factorization preserves bandedness in a very precise sense: the factor U or L has the same number of superdiagonals or subdiagonals as the original matrix. In the LU factorization, the row-interchanges modify the band structure: if A has kl subdiagonals and ku superdiagonals, then L is not a band matrix but still has at most kl nonzero elements below the diagonal in each column; and U has at most kl+ku superdiagonals.

The Bunch–Kaufman factorization does not preserve bandedness, because of the need for symmetric row-and-column permutations; hence no methods are provided for symmetric indefinite band matrices.

The inverse of a band matrix does not in general have a band structure, so no methods are provided for computing inverses of band matrices.

Block Partitioned Algorithms

Many of the methods in this chapter use what is termed a block partitioned algorithm. This means that at each major step of the algorithm a block of rows or columns is updated, and most of the computation is performed by matrix-matrix operations on these blocks. The matrix-matrix operations are performed by calls to the Level 3 BLAS (see F06 class), which are the key to achieving high performance on many modern computers. See Golub and Van Loan (1996) or Anderson et al. (1999) for more about block partitioned algorithms.

The performance of a block partitioned algorithm varies to some extent with the blocksize – that is, the number of rows or columns per block. This is a machine-dependent parameter, which is set to a suitable value when the library is implemented on each range of machines. You do not normally need to be aware of what value is being used. Different block sizes may be used for different methods. Values in the range 16 to 64 are typical.

On some machines there may be no advantage from using a block partitioned algorithm, and then the methods use an unblocked algorithm (effectively a blocksize of 1), relying solely on calls to the Level 2 BLAS (see F06 class again).

Mixed Precision LAPACK Routines

Some LAPACK routines use mixed precision arithmetic in an effort to solve problems more efficiently on modern hardware. They work by converting a double precision problem into an equivalent single precision problem, solving it and then using iterative refinement in double precision to find a full precision solution to the original problem. The method may fail if the problem is too ill-conditioned to allow the initial single precision solution, in which case the routines fall back to solve the original problem entirely in double precision. The vast majority of problems are not so ill-conditioned, and in those cases the technique can lead to significant gains in speed without loss of accuracy. This is particularly true on machines where double precision arithmetic is significantly slower than single precision.

References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Higham N J (1988) Algorithm 674: Fortran codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation ACM Trans. Math. Software 14 381–396

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

Inheritance Hierarchy

System..::.Object
NagLibrary..::.F07

compute	ri=b-Axi
solve	Adi=ri
compute	xi+1=xi+di