Numerical Linear Algebra

Probably the simplest kind of problem.

Occurs in many contexts, often as part of larger problem.

Symbolic manipulation packages can do linear algebra "analytically" (e.g. Mathematica, Maple).

Numerical methods needed when:

Number of equations very large

Coefficients all numerical

Linear Systems

Write linear system as:

This system has n unknowns and m equations.

If n = m, system is closed.

If any equation is a linear combination of any others, equations are degenerate and system is singular.*

*see Singular Value Decomposition (SVD), NRiC 2.6.

Numerical Constraints

Numerical methods also have problems when:

Equations are degenerate "within round-off error".

Accumulated round-off errors swamp solution (magnitude of a's and x's varies wildly).

For n,m < 50, single precision usually OK.

For n,m < 200, double precision usually OK.

For 200 < n,m < few thousand, solutions possible only for sparse systems (lots of a's zero).

Matrix Form

Write system in matrix form:

where:

Columns

Rows

Matrix Data Representation

Recall, C stores data in row-major form:

a11, a12, ..., a1n; a21, a22, ..., a2n; ...; am1, am2, ..., amn

If using "pointer to array of pointers to rows" scheme in C, can reference entire rows by first index, e.g. 3rd row = a[2].

Recall in C array indices start at zero!!

FORTRAN stores data in column-major form:

a11, a21, ..., am1; a12, a22, ..., am2; ...; a1n, a2n, ..., amn

Note on Numerical Recipes in C

The canned routines in NRiC make use of special functions defined in nrutil.c (header nrutil.h).

In particular, arrays and matrices are allocated dynamically with indices starting at 1, not 0.

If you want to interface with the NRiC routines, but prefer the C array index convention, pass arrays by subtracting 1 from the pointer address (i.e. pass p-1 instead of p) and pass matrices by using the functions convert_matrix() and free_convert_matrix() in nrutil.c (see NRiC 1.2 for more information).

Tasks of Linear Algebra

We will consider the following tasks:

Solve Ax = b, given A and b.

Solve Axi = bi for multiple bi's.

Calculate A-1, where A-1A = I, the identity matrix.

Calculate determinant of A, det(A).

Large packages of routines available for these tasks, e.g. LINPACK, LAPACK (public domain); IMSL, NAG libraries (commercial).

We will look at methods assuming n = m.

The Augmented Matrix

The equation Ax = b can be generalized to a form better suited to efficient manipulation:

The system can be solved by performing operations on the augmented matrix.

The xi's are placeholders that can be omitted until the end of the computation.

Elementary Row Operations

The following row operations can be performed on an augmented matrix without changing the solution of the underlying system of equations:

Interchange two rows.

Multiply a row by a nonzero real number.

Add a multiple of one row to another row.

The idea is to apply these operations in sequence until the system of equations is trivially solved.

The Generalized Matrix Equation

Consider the generalized linear matrix equation:

Its solution simultaneously solves the linear sets:

Ax1 = b1, Ax2 = b2, Ax3 = b3, and AY = I,

where the xi's and bi's are column vectors.

Gauss-Jordan Elimination

GJE uses one or more elementary row operations to reduce matrix A to the identity matrix.

The RHS of the generalized equation becomes the solution set and Y becomes A-1.

Disadvantages:

Requires all bi's to be stored and manipulated at same time  memory hog.

Don't always need A-1.

Other methods more efficient, but good backup.

Gauss-Jordan Elimination: Procedure

Start with simple augmented matrix as example:

Divide first row (a1|b1) by first element a11.

Subtract ai1 (a1|b1) from all other rows:

Continue process for 2nd row, etc.

Row a1|b1

Pivot row

First column of identity matrix

GJE Procedure, Cont'd

Problem occurs if leading diagonal element ever becomes zero.

Also, procedure is numerically unstable!

Solution: use "pivoting" - rearrange remaining rows (partial pivoting) or rows & columns (full pivoting - requires permutation!) so largest coefficient is in diagonal position.

Best to "normalize" equations (implicit pivoting).

Gaussian Elimination with Backsubstitution

If, during GJE, only subtract rows below pivot, will be left with a triangular matrix:

Solution for x3 is then trivial: x3 = b3'/a33'.

Substitute into 2nd row to get x2.

Substitute x3 & x2 into 1st row to get x1.

Faster than GJE, but still memory hog.

"Gaussian Elimination"

LU Decomposition

Suppose we can write A as a product of two matrices: A = LU, where L is lower triangular and U is upper triangular:

Then Ax = (LU)x = L(Ux) = b, i.e. must solve,

(1) Ly = b; (2) Ux = y

Can reuse L & U for subsequent calculations.

LU Decomposition, Cont'd

Why is this better?

Solving triangular matrices is easy: just use forward substitution for (1), backsubstitution for (2).

Problem is, how to decompose A into L and U?

Expand matrix multiplication LU to get n2 equations for n2 + n unknowns (elements of L and U plus n extras because diagonal elements counted twice).

Get an extra n equations by choosing Lii = 1 (i = 1,n).

Then use Crout's algorithm for finding solution to these n2 + n equations "trivially" (NRiC 2.3).

LU Decomposition in NRiC

The routines ludcmp() and lubksb() perform LU decomposition and backsubstitution respectively.

Can easily compute A-1 (solve for the identity matrix column by column) and det(A) (find the product of the diagonal elements of the LU decomposed matrix) - see NRiC 2.3.

WARNING: for large matrices, computing det(A) can overflow or underflow the computer's floating-point dynamic range.

Iterative Improvement

For large sets of linear equations Ax = b, roundoff error may become a problem.

We want to know x but we only have x + x, which is an exact solution to A(x + x) = b + b.

Subtract the exact solution and eliminate b:

Ax = A(x + x) - b

The RHS is known, hence can solve for x. Subtract this from the wrong solution to get an improved solution (make sure to use doubles!).

Tridiagonal Matrices

Many systems can be written as (or reduced to):

ai xi-1 + bi xi + ci xi+1 = di	i = 1,n

i.e. a tridiagonal matrix:

Here a1 and cn are associated with "boundary conditions" (i.e. x0 and xn+1).

Sparse Matrices

LU decomposition and backsubstitution is very efficient for tri-di systems: O(n) operations as opposed to O(n3) in general case.

Operations on sparse systems can be optimized.

e.g. Tridiagonal

Band diagonal with bandwidth M

Block diagonal

Banded

See NRiC 2.7 for various systems & techniques.

Iterative Methods

For very large systems, direct solution methods (e.g. LU decomposition) are slow and RE prone.

Often iterative methods much more efficient:

Guess a trial solution x0

Compute a correction x1 = x0 + x

Iterate procedure until convergence, i.e. |x| < D

e.g. Congugate gradient method for sparse systems (NRiC 2.7).

Singular Value Decomposition

Can diagnose or (nearly) solve singular or near-singular systems.

Used for solving linear least-squares problems.

Theorem: any m n matrix A can be written:

A = UWV T

where U (m n) & V (n n) are orthogonal and W (n n) is a diagonal matrix.

Proof: buy a good linear algebra textbook.

SVD, Cont'd

The values Wi are zero or positive and are called the "singular values".

The NRiC routine svdcmp() returns U, V, & W given A. You have to trust it (or test it yourself!).

Uses Householder reduction, QR diagonalization, etc.

If A is square then we know:

A-1 = V [diag(1/Wi)] U T

This is fine so long as no Wi is too small (or 0).

Definitions

Condition number cond(A) = (max Wi)/(min Wi).

If cond(A) = , A is singular.

If cond(A) very large (106, 1012), A is ill-conditioned.

Consider Ax = b. If A is singular, there is some subspace of x (the nullspace) such that Ax = 0.

The nullity is the dimension of the nullspace.

The subspace of b such that Ax = b is the range.

The rank of A is the dimension of the range.

The Homogeneous Equation

SVD constructs orthonormal bases for the nullspace and range of a matrix.

Columns of U with corresponding non-zero Wi are an orthonormal basis for the range.

Columns of V with corresponding zero Wi are an orthonormal basis for the nullspace.

Hence immediately have solution for Ax = 0, i.e. the columns of V with corresponding zero Wi.

Residuals

If b ( 0) lies in the range of A, then the singular equations do in fact have a solution.

Even if b is outside the range of A, can get solution which minimizes residual r = |Ax - b|.

Trick: replace 1/Wi by 0 if Wi = 0 and compute

x = V [diag (1/Wi)] (U T b)

Similarly, can set 1/Wi = 0 if Wi very small.

Approximation of Matrices

Can write A = UWV T as:

If most of the singular values Wk are small, then A is well-approximated by only a few terms in the sum (strategy: sort Wk's in descending order).

For large memory savings, just store the columns of U and V corresponding to non-negligible Wk's.

Useful technique for digital image processing.

2001-02-15T09:13:37

Derek Richardson