About the Project
3 Numerical MethodsAreas

§3.2 Linear Algebra

Contents
  1. §3.2(i) Gaussian Elimination
  2. §3.2(ii) Gaussian Elimination for a Tridiagonal Matrix
  3. §3.2(iii) Condition of Linear Systems
  4. §3.2(iv) Eigenvalues and Eigenvectors
  5. §3.2(v) Condition of Eigenvalues
  6. §3.2(vi) Lanczos Tridiagonalization of a Symmetric Matrix
  7. §3.2(vii) Computation of Eigenvalues

§3.2(i) Gaussian Elimination

To solve the system

3.2.1 𝐀𝐱=𝐛,

with Gaussian elimination, where 𝐀 is a nonsingular n×n matrix and 𝐛 is an n×1 vector, we start with the augmented matrix

3.2.2 [a11a1nb1an1annbn].

By repeatedly subtracting multiples of each row from the subsequent rows we obtain a matrix of the form

3.2.3 [u11u12u1ny10u22u2ny200unnyn].

During this reduction process we store the multipliers jk that are used in each column to eliminate other elements in that column. This yields a lower triangular matrix of the form

3.2.4 𝐋=[1002110n1n,n11].

If we denote by 𝐔 the upper triangular matrix comprising the elements ujk in (3.2.3), then we have the factorization, or triangular decomposition,

3.2.5 𝐀=𝐋𝐔.

With 𝐲=[y1,y2,,yn]T the process of solution can then be regarded as first solving the equation 𝐋𝐲=𝐛 for 𝐲 (forward elimination), followed by the solution of 𝐔𝐱=𝐲 for 𝐱 (back substitution).

For more details see Golub and Van Loan (1996, pp. 87–100).

Example

3.2.6 [123231312]=[100210351][1230150018].

In solving 𝐀𝐱=[1,1,1]T, we obtain by forward elimination 𝐲=[1,1,3]T, and by back substitution 𝐱=[16,16,16]T.

In practice, if any of the multipliers jk are unduly large in magnitude compared with unity, then Gaussian elimination is unstable. To avoid instability the rows are interchanged at each elimination step in such a way that the absolute value of the element that is used as a divisor, the pivot element, is not less than that of the other available elements in its column. Then |jk|1 in all cases. This modification is called Gaussian elimination with partial pivoting.

For more information on pivoting see Golub and Van Loan (1996, pp. 109–123).

Iterative Refinement

When the factorization (3.2.5) is available, the accuracy of the computed solution 𝐱 can be improved with little extra computation. Because of rounding errors, the residual vector 𝐫=𝐛𝐀𝐱 is nonzero as a rule. We solve the system 𝐀𝛿𝐱=𝐫 for 𝛿𝐱, taking advantage of the existing triangular decomposition of 𝐀 to obtain an improved solution 𝐱+𝛿𝐱.

§3.2(ii) Gaussian Elimination for a Tridiagonal Matrix

Tridiagonal matrices are ones in which the only nonzero elements occur on the main diagonal and two adjacent diagonals. Thus

3.2.7 𝐀=[b1c10a2b2c2an1bn1cn10anbn].

Assume that 𝐀 can be factored as in (3.2.5), but without partial pivoting. Then

3.2.8 𝐋=[100210n1100n1],
3.2.9 𝐔=[d1u100d2u20dn1un100dn],

where uj=cj, j=1,2,,n1, d1=b1, and

3.2.10 j =aj/dj1,
dj =bjjcj1,
j=2,,n.

Forward elimination for solving 𝐀𝐱=𝐟 then becomes y1=f1,

3.2.11 yj=fjjyj1,
j=2,,n,

and back substitution is xn=yn/dn, followed by

3.2.12 xj=(yjujxj+1)/dj,
j=n1,,1.

For more information on solving tridiagonal systems see Golub and Van Loan (1996, pp. 152–160).

§3.2(iii) Condition of Linear Systems

The p-norm of a vector 𝐱=[x1,,xn]T is given by

3.2.13 𝐱p =(j=1n|xj|p)1/p,
p=1,2,,
𝐱 =max1jn|xj|.

The Euclidean norm is the case p=2.

The p-norm of a matrix 𝐀=[ajk] is

3.2.14 𝐀p=max𝐱𝟎𝐀𝐱p𝐱p.

The cases p=1,2, and are the most important:

3.2.15 𝐀1 =max1knj=1n|ajk|,
𝐀 =max1jnk=1n|ajk|,
𝐀2 =ρ(𝐀𝐀T),

where ρ(𝐀𝐀T) is the largest of the absolute values of the eigenvalues of the matrix 𝐀𝐀T; see §3.2(iv). (We are assuming that the matrix 𝐀 is real; if not 𝐀T is replaced by 𝐀H, the transpose of the complex conjugate of 𝐀.)

The sensitivity of the solution vector 𝐱 in (3.2.1) to small perturbations in the matrix 𝐀 and the vector 𝐛 is measured by the condition number

3.2.16 κ(𝐀)=𝐀p𝐀1p,

where p is one of the matrix norms. For any norm (3.2.14) we have κ(𝐀)1. The larger the value κ(𝐀), the more ill-conditioned the system.

Let 𝐱* denote a computed solution of the system (3.2.1), with 𝐫=𝐛𝐀𝐱* again denoting the residual. Then we have the a posteriori error bound

3.2.17 𝐱*𝐱p𝐱pκ(𝐀)𝐫p𝐛p.

For further information see Brezinski (1999) and Trefethen and Bau (1997, Chapter 3).

§3.2(iv) Eigenvalues and Eigenvectors

If 𝐀 is an n×n matrix, then a real or complex number λ is called an eigenvalue of 𝐀, and a nonzero vector 𝐱 a corresponding (right) eigenvector, if

3.2.18 𝐀𝐱=λ𝐱.

A nonzero vector 𝐲 is called a left eigenvector of 𝐀 corresponding to the eigenvalue λ if 𝐲T𝐀=λ𝐲T or, equivalently, 𝐀T𝐲=λ𝐲. A normalized eigenvector has Euclidean norm 1; compare (3.2.13) with p=2.

The polynomial

3.2.19 pn(λ)=det[λ𝐈𝐀]

is called the characteristic polynomial of 𝐀 and its zeros are the eigenvalues of 𝐀. The multiplicity of an eigenvalue is its multiplicity as a zero of the characteristic polynomial (§3.8(i)). To an eigenvalue of multiplicity m, there correspond m linearly independent eigenvectors provided that 𝐀 is nondefective, that is, 𝐀 has a complete set of n linearly independent eigenvectors.

§3.2(v) Condition of Eigenvalues

If 𝐀 is nondefective and λ is a simple zero of pn(λ), then the sensitivity of λ to small perturbations in the matrix 𝐀 is measured by the condition number

3.2.20 κ(λ)=1|𝐲T𝐱|,

where 𝐱 and 𝐲 are the normalized right and left eigenvectors of 𝐀 corresponding to the eigenvalue λ. Because |𝐲T𝐱|=|cosθ|, where θ is the angle between 𝐲T and 𝐱 we always have κ(λ)1. When 𝐀 is a symmetric matrix, the left and right eigenvectors coincide, yielding κ(λ)=1, and the calculation of its eigenvalues is a well-conditioned problem.

§3.2(vi) Lanczos Tridiagonalization of a Symmetric Matrix

Let 𝐀 be an n×n symmetric matrix. Define the Lanczos vectors 𝐯j and coefficients αj and βj by 𝐯0=𝟎, a normalized vector 𝐯1 (perhaps chosen randomly), α1=𝐯1T𝐀𝐯1, β1=0, and for j=1,2,,n1 by the recursive scheme

3.2.21 𝐮 =𝐀𝐯jαj𝐯jβj𝐯j1,
βj+1 =𝐮2,
𝐯j+1 =𝐮/βj+1,
αj+1 =𝐯j+1T𝐀𝐯j+1.

Then 𝐯jT𝐯k=δj,k, j,k=1,2,,n. The tridiagonal matrix

3.2.22 𝐁=[α1β20β2α2β3βn1αn1βn0βnαn]

has the same eigenvalues as 𝐀. Its characteristic polynomial can be obtained from the recursion

3.2.23 pk+1(λ)=(λαk+1)pk(λ)βk+12pk1(λ),
k=0,1,,n1,

with p1(λ)=0, p0(λ)=1.

In the case that the orthogonality condition is replaced by 𝐒-orthogonality, that is, 𝐯jT𝐒𝐯k=δj,k, j,k=1,2,,n, for some positive definite matrix 𝐒 with Cholesky decomposition 𝐒=𝐋T𝐋, then the details change as follows. Start with 𝐯0=𝟎, vector 𝐯1 such that 𝐯1T𝐒𝐯1=1, α1=𝐯1T𝐀𝐯1, β1=0. Then for j=1,2,,n1

3.2.24 𝐮 =(𝐀αj𝐒)𝐯jβj𝐒𝐯j1,
βj+1 =(𝐋1)T𝐮2,
𝐯j+1 =𝐋1(𝐋1)T𝐮/βj+1,
αj+1 =𝐯j+1T𝐀𝐯j+1.

For more details see Guan et al. (2007).

For numerical information see Stewart (2001, pp. 347–368).

Lanczos’ method is related to Gauss quadrature considered in §3.5(v). When the matrix 𝐀 is replaced by a scalar x, the recurrence relation in the first line of (3.2.21) with 𝐮=βj+1𝐯j+1 is similar to the one in (3.5.30_5). Also, the recurrence relations in (3.2.23) and (3.5.30) are similar, as well as the matrix 𝐁 in (3.2.22) and the Jacobi matrix 𝐉n in (3.5.31).

§3.2(vii) Computation of Eigenvalues

Many methods are available for computing eigenvalues; see Golub and Van Loan (1996, Chapters 7, 8), Trefethen and Bau (1997, Chapter 5), and Wilkinson (1988, Chapters 8, 9).