. 105
( 136 .)



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

x ∈ [’ ∞, ∞ ]



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

x ∈ [0, ∞ ]


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

x ∈ [0, ∞ ]


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
N3 2N 2
LU Default method
N3 2N 2
Cholesky Cheaper than LU, but only
applies to symmetric matrices
N3 4 N2
QR Good for ill-conditioned matrices, but
twice as costly as LU


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
[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
lij = aij ’ i≥j (B.5a)
lik ukj

aij ’ (B.5b)
uij = lik ukj i<j

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

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
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:
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”.
• R is a Cholesky factor of A A, that is
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.
• 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
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

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

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

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

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


. 105
( 136 .)