Im(x)

Re(x)

Figure A.3: Hermite & sinc functions: straight lines, Im(x) = constant

x ∈ [’ ∞, ∞ ]

Im(x)

TBn(x)

Re(x)

Figure A.4: T Bn (x): pair of circles, the contours of “bipolar” coordinates in the complex

x-plane

A.10. GRAPHS OF CONVERGENCE DOMAINS IN THE COMPLEX PLANE 513

x ∈ [0, ∞ ]

LAGUERRE

Im(x)

Re(x)

Figure A.5: Laguerre functions: parabola with focus at x = 0

x ∈ [0, ∞ ]

TLn(x)

Im(x)

Re(x)

Figure A.6: T Ln (x): quartic curve with no standard name

Appendix B

Direct Matrix-Solvers

“When iterating, three things may happen and two of them are bad.”

” J. P. Boyd (paraphrase of football proverb)

In football, the three fates of a forward pass are “completion” (good), “incomplete” (bad) and “interception” (very

bad). Similarly, an iteration may converge, may converge too slowly to be useful, or diverge.

B.1 Matrix Factorizations

The most ef¬cient “direct” way to solve a matrix equation is to ¬rst simplify the square

matrix to a special factored form and then perform a “backsolve”. A “direct” method is a

non-iterative algorithm; iterative matrix-solvers are discussed in Chapter 15 above.

The three most useful factorizations are (i) the product of two triangular matrices L

and U , usually called the “LU” factorization (or “Crout reduction”) (ii) the product of a

“lower triangular” matrix L and its transpose [“Choleski factorization”], possible only for a

symmetric matrix and (iii) the product of an “orthogonal” matrix Q with an upper triangular

matrix R, commonly called the “QR” decomposition. The relative costs and merits of these

three methods are given Table B.1.

Very robust library software is available for all these algorithms, so the user need never

write his own. However, a little knowledge of matrix properties is important in choosing

the most ef¬cient subroutine.

Triangular matrices and the LU and Cholesky factorizations are de¬ned by following.

Table B.1: Cost of Matrix-Solving Methods for Dense (Full) Matrices

Note: In the cost estimates, N is the size of the (dense) matrix which is factored. The algo-

rithms all require equal numbers of additions/subtractions and multiplications/division;

the numbers are the sums of all operations.

Name Factorization Cost Backsolve Cost Comments

2

N3 2N 2

LU Default method

3

1

N3 2N 2

Cholesky Cheaper than LU, but only

3

applies to symmetric matrices

4

N3 4 N2

QR Good for ill-conditioned matrices, but

3

twice as costly as LU

514

B.1. MATRIX FACTORIZATIONS 515

De¬nition 45 (TRIANGULAR MATRIX) A matrix A is said to be in LOWER TRIANGULAR

or UPPER TRIANGULAR form, respectively, if

whenever j > i [LOWER TRIANGULAR] (B.1)

Aij = 0

whenever i < j [UPPER TRIANGULAR] (B.2)

Aij = 0

In other words, on one side, either the upper or lower, of the diagonal row of elements Aii , a trian-

gular matrix A has nothing but elements whose values are zero.

Theorem 40 (LU DECOMPOSITION OF A MATRIX) If its determinant is non-zero, then an

otherwise arbitary matrix A can always be factored into the product of a lower triangular matrix L

and an upper triangular matrix U. For 4—4 matrices, for example,

l11 0 0 0 1 u12 u13 u14 a11 a12 a13 a14

l21 l22 0 0 0 1 u23 u24 a a22 a23 a24

= 21 (B.3)

l31 l32 l33 0 0 0 1 u34 a31 a32 a33 a34

l41 l42 l43 l44 0 0 0 1 a41 a42 a43 a44

If the matrix A is REAL, and POSITIVE DEFINITE, then it can factored into

SYMMETRIC,

T

[Cholesky decomposition] (B.4)

LL = A

where superscript T denotes the transpose. The Cholesky decomposition needs only half the number

of operations and half the storage of the standard LU decomposition.

PROOF: The elements of L and U [uii ≡ 1] are given explicitly by

j’1

lij = aij ’ i≥j (B.5a)

lik ukj

k=1

i’1

1

aij ’ (B.5b)

uij = lik ukj i<j

lii

k=1

which is nothing more than a restatement of the fact that the elements of A are given by

the matrix product of L and U. The elements are computed in the order

li1 , u1j ; li2 , u2j ; . . . , li,N ’1 , uN ’1,j ; lN N ,

that is, we compute one column of L, then the corresponding row of U, and return to

evaluate the next column of L and so on. Q. E. D.

The proof ignores pivoting, which is rarely necessary for matrices arising from spectral

or ¬nite element discretizations. Good reviews of Gaussian elimination and factorization

are given in Conte and de Boor (1980) and Forsythe and Moler (1967).

The matrices which are the Chebyshev or Legendre pseudospectral discretizations of

derivatives may be somewhat ill-conditioned, that is, generate a lot of round-off error in

the LU method, when N is large (> 100) and/or the order of the derivatives is high (third

derivatives or higher). One remedy is to use a different factorization

APPENDIX B. DIRECT MATRIX-SOLVERS

516

Theorem 41 (QR Factorization)

• An arbitrary matrix A, even if rectangular or singular, can always be factored as the product

of an orthogonal matrix Q multiplied on the right by a triangular matrix R. (A matrix is

“orthogonal” or “unitary” if its transpose is also its inverse.)

• The matrix Q is (implicitly) the product of a series of rank-one orthogonal matrices (“House-

holder transformations”). Because all transformations are unitary, QR factorization is much

more resistant to accumulation of roundoff error than the LU factorization.

• It is expensive to explicitly compute Q. To solve Ax = b, one alternative is to apply the

Householder transformations simultaneously to both A and b to transform the matrix into

T

R and the vector into Q b. The solution x is then computed by backsolving the triangular

system Rx = b.

• Alternatively, one can use the “semi-normal equations”: one multiplication of a vector by a

square matrix followed by two triangular backsolves:

T T

(B.6)

R Rx =A b

This has less numerical stability (to roundoff) than the “Q-and-R” or “simultaneous trans-

formation” approach, so the semi-normal backsolve is usually repeated with b replaced by

b ’ Ax, the residual of the linear system, to provide one step of “iterative re¬nement”.

T

• R is a Cholesky factor of A A, that is

T T

(B.7)

R R=A A

This implies that the QR factorization doubles the bandwidth of a banded matrix; the ratio of

(QR cost)/(LU cost) is even bigger for sparse matrices than the ratio of two for dense matrices.

• If A contains even a single full row, R will be a dense matrix.

T T

• For overdetermined systems, the “normal equations” A Ax = A b and the QR factoriza-

tion will both give bounded solutions if the columns of A are all linearly independent (“full

column rank”). However, accumulation of roundoff error is often disastrous in the normal

equations; QR factorization is much more stable.

• If the matrix is intrinsically singular, and not merely apparently singular because of accu-

mulated roundoff, R will have one or more zero diagonal elements and the QR algorithm will

fail.

QR may be a useful fallback algorithm when LU software crashes with a message like

“Matrix is nearly singular”. However, my personal experience is that 95% of the occur-

rences of an ill-conditioned pseudospectral matrix arise because of a coding error. Dennis

and Schnabel (1983) give a brief and readable treatment of QR.

De¬nition 46 (SPARSE & DENSE) A matrix A is said to be SPARSE if most of its matrix el-

ements are zero. A matrix whose elements are all (or mostly) non-zero is said to be DENSE or

FULL.

B.1. MATRIX FACTORIZATIONS 517

Figure B.1: Schematic of sparse matrix with disks representing nonzero matrix elements.

The heavy curves are the “skyline” which bounds the non-zero elements; it is (by de¬-

nition) symmetric with respect to the diagonal of the matrix. The right half of the ¬gure

shows that the lower and upper triangular factors of a matrix have the same “skyline”.

Spectral methods usually generate dense matrices (alas!). The costs of the factorizations

for such dense matrices are given in Table B.1. However, domain decomposition generates

sparse matrices (with a few dense blocks). Galerkin discretizations of special differential

equations, (usually constant coef¬cient and linear) , also generate sparse matrices (Chapter

15, Sec. 10). The LU and Cholesky factorizations can then be modi¬ed to omit operations

on the zero elements, thereby greatly reducing cost.

Because matrix sparsity occurs in many patterns and can generate enormous savings, a

huge literature on “sparse matrix” methods has arisen. We cannot possibly do it justice in

an appendix, but must instead refer the reader to a good book on numerical linear algebra.

However, some patterns are so common that it is worthwhile to brie¬‚y describe them here.

The most general sparse matrix software which can be explained in simple terms is a

“skyline solver”. Kikuchi (1986) gives both FORTRAN and BASIC codes and also a good

explanation of the rather complicated “skylines” that can arise from ¬nite element and

spectral element discretizations. The concept of the “skyline”, which is de¬ned formally

below , is important because it is preserved by LU factorization as also illustrated in the

graph.

De¬nition 47 (SKYLINE of a MATRIX) Let all the nonzero elements of a matrix A satisfy the

inequality

Aij , Aji = 0 ∀|i ’ j| ≥ σ(i)

where σ(i) is chosen to be as small as possible for each i. Then σ(i) is the “SKYLINE” of the matrix

A.

The crucial point is that if the matrix is sparse, many of the operations above involve

only elements known to be zero. The cost of the factorization can be greatly reduced by