Linear Equations and Eigensystems

CHAPTER

LINEAR EQUATIONS AND EIGENSYSTEMS

2

Abstract When physical systems are modeled mathematically, they are sometimes described by linear equation systems or alternatively by eigensystems and in this chapter we examine how such equation systems are solved. M ATLAB is an ideal environment for studying linear algebra, including linear equation systems and eigenvalue problems, because M ATLAB functions and operators can work directly on vectors and matrices. M ATLAB is rich in functions and operators which facilitate the manipulation of matrices.

2.1 INTRODUCTION We now discuss linear equation systems and defer discussion of eigensystems until Section 2.15. To illustrate how linear equation systems arise in the modeling of certain physical problems, we will consider how current flows are calculated in a simple electrical network. The necessary equations can be developed using one of several techniques; here we use the loop-current method together with Ohm’s law and Kirchhoff’s voltage law. A loop current is assumed to circulate around each loop in the network. Thus, in the network given in Fig. 2.1, the loop current I1 circulates around the closed loop abcd. Thus the current I1 − I2 flows in the link connecting b to c. Ohm’s law states that the voltage across an ideal resistor is proportional to the current flow through the resistor. For example, for the link connecting b to c Vbc = R2 (I1 − I2 ) where R2 is the value of the resistor in the link connecting b to c. Kirchhoff’s voltage law states that the algebraic sum of the voltages around a loop is zero. Applying these laws to the circuit abcd of Fig. 2.1 we have Vab + Vbc + Vcd = V Substituting the product of current and resistance for voltage gives R1 I1 + R2 (I1 − I2 ) + R4 I1 = V Numerical Methods. https://doi.org/10.1016/B978-0-12-812256-3.00011-7 Copyright © 2019 Elsevier Inc. All rights reserved.

73

74

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

FIGURE 2.1 Electrical network.

We can repeat this process for each loop to obtain the following four equations: (R1 + R2 + R4 ) I1 − R2 I2 = V (R1 + 2R2 + R4 ) I2 − R2 I1 − R2 I3 = 0 (R1 + 2R2 + R4 ) I3 − R2 I2 − R2 I4 = 0 (R1 + R2 + R3 + R4 ) I4 − R2 I3 = 0

(2.1)

Letting R1 = R4 = 1 , R2 = 2 , R3 = 4  and V = 5 volts, (2.1) becomes 4I1 − 2I2 = 5 −2I1 + 6I2 − 2I3 = 0 −2I2 + 6I3 − 2I4 = 0 −2I3 + 8I4 = 0 This is a system of linear equations in four variables, I1 , ..., I4 . In matrix notation it becomes ⎡ ⎤⎡ ⎤ ⎡ ⎤ 5 4 −2 0 0 I1 ⎢ −2 ⎥ ⎢ I 2 ⎥ ⎢ 0⎥ 6 −2 0 ⎢ ⎥⎢ ⎥ = ⎢ ⎥ ⎣ 0 −2 6 −2 ⎦ ⎣I3 ⎦ ⎣0⎦ I4 0 0 0 −2 8

(2.2)

This equation has the form Ax = b where A is a square matrix of known coefficients, in this case relating to the values of the resistors in the circuit. The vector b is a vector of known coefficients, in this case the voltage applied to each current loop. The vector x is the vector of unknown currents. Although this set of equations can be solved by hand, the process is time consuming and error prone. Using M ATLAB we simply enter matrix A and vector b and use the command A\b as follows: >> A = [4 -2 0 0;-2 6 -2 0;0 -2 6 -2;0 0 -2 8]; >> b = [5 0 0 0].’; >> x = A\b

2.1 INTRODUCTION

75

x = 1.5426 0.5851 0.2128 0.0532 >> A*x ans = 5.0000 -0.0000 0.0000 -0.0000

The sequence of operations that are invoked by this apparently simple command is examined in Section 2.3. In many electrical networks the ideal resistors of Fig. 2.1 are more accurately represented by electrical impedances. When a harmonic alternating current (AC) supply is connected to the network, electrical engineers represent the impedances by complex quantities. This is to account for the effect of capacitance and/or inductance. To illustrate this we will replace the 5 volt DC supply to the network of Fig. 2.1 by a 5 volt AC supply and replace the ideal resistors R1 , ..., R4 by impedances Z1 , ..., Z4 . Thus, (2.1) becomes (Z1 + Z2 + Z4 )I1 − Z2 I2 = V (Z1 + 2Z2 + Z4 )I2 − Z2 I1 − Z2 I3 = 0 (Z1 + 2Z2 + Z4 )I3 − Z2 I2 − Z2 I4 = 0

(2.3)

(Z1 + Z2 + Z3 + Z4 )I4 − Z2 I3 = 0 At the frequency of the 5 volt AC supply √ we will assume that Z1 = Z4 = (1 + 0.5j ) , Z2 =√(2 + 0.5j ) , Z3 = (4 + 1j )  where j = −1. Electrical engineers prefer to use j rather than ı for −1. This avoids any possible confusion with i (or I ) which is normally used to denote the current in a circuit. Thus, (2.3) becomes (4 + 1.5j )I1 − (2 + 0.5j )I2 = 5 −(2 + 0.5j )I1 + (6 + 2.0j )I2 − (2 + 0.5j )I3 = 0 −(2 + 0.5j )I2 + (6 + 2.0j )I3 − (2 + 0.5j )I4 = 0 −(2 + 0.5j )I3 + (8 + 2.5j )I4 = 0 This system of linear equations becomes, in matrix notation, ⎡

⎤⎡ (4 + 1.5j ) −(2 + 0.5j ) 0 0 I1 ⎢ −(2 + 0.5j ) ⎥ ⎢ I2 (6 + 2.0j ) −(2 + 0.5j ) 0 ⎢ ⎥⎢ ⎣ 0 −(2 + 0.5j ) (6 + 2.0j ) −(2 + 0.5j ) ⎦ ⎣ I3 I4 0 0 −(2 + 0.5j ) (8 + 2.5j )

⎤ 5 ⎥ ⎢ 0 ⎥ ⎥=⎢ ⎥ ⎦ ⎣ 0 ⎦ 0 ⎤

⎡

(2.4)

76

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

Note that the coefficient matrix is now complex. This does not present any difficulty for M ATLAB because the operation A\b works directly with both real and complex numbers. Thus: >> p = 4+1.5i; q = -2-0.5i; >> r = 6+2i; s = 8+2.5i; >> A = [p q 0 0;q r q 0;0 q r q;0 0 q s]; >> b = [5 0 0 0].’; >> A\b ans = 1.3008 - 0.5560i 0.4560 - 0.2504i 0.1530 - 0.1026i 0.0361 - 0.0274i

Note that strictly we have no need to re-enter the values in vector b, assuming that we have not cleared the memory, reassigned the vector b or quit M ATLAB. The answer shows that currents flowing in the network are complex. This means that there is a phase difference between the applied harmonic voltage and the currents flowing. We will now begin a more detailed examination of linear equation systems.

2.2 LINEAR EQUATION SYSTEMS In general, a linear equation system can be written in matrix form as Ax = b

(2.5)

where A is an n × n matrix of known coefficients, b is a column vector of n known coefficients and x is the column vector of n unknowns. We have already seen an example of this type of equation system in Section 2.1 where the matrix equation (2.2) is the matrix equivalent of the set of linear equations (2.1). The equation system (2.5) is called homogeneous if b = 0 and inhomogeneous if b = 0. Before attempting to solve an equation system it is reasonable to ask if it has a solution and if so is it unique? A linear inhomogeneous equation system may be consistent and have one or an infinity of solutions or be inconsistent and have no solution. This is illustrated in Fig. 2.2 for a system of three equations in three variables x1 , x2 , and x3 . Each equation represents a plane surface in the x1 , x2 , x3 space. In Fig. 2.2A the three planes have a common point of intersection. The coordinates of the point of intersection give the unique solution for the three equations. In Fig. 2.2B the three planes intersect in a line. Any point on the line of intersection represents a solution so there is no unique solution but an infinite number of solutions satisfying the three equations. In Fig. 2.2C two of the surfaces are parallel to each other and therefore they never intersect while in Fig. 2.2D the line of intersection of each pair of surfaces is different. In both of these cases there is no solution and the equations these surfaces represent are inconsistent.

2.2 LINEAR EQUATION SYSTEMS

77

FIGURE 2.2 Three intersecting planes representing three equations in three variables. (A) Three plane surfaces intersecting in a point. (B) Three plane surfaces intersecting in a line. (C) Three plane surfaces, two of which do not intersect. (D) Three plane surfaces intersecting in three lines.

To obtain an algebraic solution to the inhomogeneous equation system (2.5) we multiply both sides of (2.5) by a matrix called the inverse of A, denoted by A−1 , thus: A−1 Ax = A−1 b

(2.6)

A−1 A = AA−1 = I

(2.7)

where A−1 is defined by

and I is the identity matrix. Thus, we obtain x = A−1 b

(2.8)

The standard algebraic formula for the inverse of A is A−1 = adj(A)/|A|

(2.9)

78

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

where |A| is the determinant of A and adj(A) is the adjoint of A. The determinant and the adjoint of a matrix are defined in Appendix A. Eqs. (2.8) and (2.9) are algebraic statements allowing us to determine x but they do not provide an efficient means of solving the system because computing A−1 using (2.9) is extremely inefficient, involving an order (n + 1)! multiplications where n is the number of equations. However, (2.9) is theoretically important because it shows that if |A| = 0 then A does not have an inverse. The matrix A is then said to be singular and a unique solution for x does not exist. Thus, establishing that |A| is non-zero is one way of showing that an inhomogeneous equation system is a consistent system with a unique solution. It is shown in Sections 2.6 and 2.7 that (2.5) can be solved without formally determining the inverse of A. An important concept in linear algebra is the rank of a matrix. For a square matrix, the rank is the number of independent rows or columns of the matrix. Independence can be explained as follows. The rows (or columns) of a matrix can clearly be viewed as a set of vectors. A set of vectors is said to be linearly independent if none of them can be expressed as a linear combination of any of the others. By linear combination we mean a sum of scalar multiples of the vectors. For example, the matrix ⎡  ⎤ ⎡ ⎡ ⎡ ⎤ ⎤ ⎡ ⎤ ⎡ ⎤ ⎤ 1 2 3 1 2 3 1 2 3 ⎥ ⎢  ⎣ −2 1 4 ⎦ or ⎢ −2 1 4 ⎥ or ⎣ ⎣ −2 ⎦ ⎣ 1 ⎦ ⎣ 4 ⎦ ⎦ ⎣  ⎦ −1 3 7 7 −1 3 −1 3 7 has linearly dependent rows and columns. This is because row 3 − row 1 − row 2 = 0 and column 3 − 2(column 2) + column 1 = 0. There is only one equation relating the rows and thus there are two independent rows. Similarly with the columns. Hence, this matrix has a rank of 2. Now consider ⎡ ⎤ 1 2 3 ⎣ 2 4 6 ⎦ 3 6 9 Here row 2 = 2(row 1) and row 3 = 3(row 1). There are two equations relating the rows and hence only one row is independent so that the matrix has a rank of 1. Note that the number of independent rows and columns in a square matrix is identical; that is, its row rank and column rank are equal. In general, matrices may be non-square and the rank of an m × n matrix A is written rank(A). Matrix A is said to be of full rank if rank(A) = min(m, n); otherwise rank(A) < min(m, n) and A is said to be rank deficient. M ATLAB provides the function rank which works with both square and non-square matrices to find the rank of a matrix. In practice, M ATLAB determines the rank of a matrix from its singular values; see Section 2.10. For example, consider the following M ATLAB statements: >> D = [1 2 3;3 4 7;4 -3 1;-2 5 3;1 -7 6] D = 1 3 4 -2 1

2 4 -3 5 -7

3 7 1 3 6

2.2 LINEAR EQUATION SYSTEMS

79

>> rank(D) ans = 3

Thus D is of full rank, since its rank is equal to the minimum number of rows or column of D. A useful operation in linear algebra is the conversion of a matrix to its reduced row echelon form (RREF). The RREF is defined in Appendix A. In M ATLAB we can use the rref function to compute the RREF of a matrix thus: >> rref(D) ans = 1

0

0

0

1

0

0

0

1

0

0

0

0

0

0

It is a property of the RREF of a matrix that the number of rows with at least one non-zero element equals the rank of the matrix. In this example we see that there are three rows in the RREF of the matrix containing a non-zero element, confirming that the matrix rank is 3. The RREF also allows us to determine whether a system has a unique solution or not. We have discussed a number of important concepts relating to the nature of linear equations and their solutions. We now summarize the equivalencies between these concepts. Let A be an n × n matrix. If Ax = b is consistent and has a unique solution, then all of the following statements are true: Ax = 0 has only the trivial solution x = 0. A is non-singular and det(A) = 0. The RREF of A is the identity matrix. A has n linearly independent rows and columns. A has full rank, that is rank(A) = n. In contrast, if Ax = b is either inconsistent or consistent but with more than one solution, then all of the following statements are true: Ax = 0 has more than one solution. A is singular and det(A) = 0. The RREF of A contains at least one zero row. A has linearly dependent rows and columns. A is rank deficient, that is rank(A) < n. So far we have only considered the case where there are as many equations as unknowns. Now we consider the cases where there are fewer or more equations than the number of unknown variables. If there are fewer equations than unknowns, then the system is said to be under-determined. The equation system does not have a unique solution; it is either consistent with an infinity of solutions, or inconsistent with no solution. These conditions are illustrated by Fig. 2.3. The diagram shows two plane surfaces in three-dimensional space, representing two equations in three variables. It is seen that

80

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

FIGURE 2.3 Planes representing an under-determined system of equations. (A) Two plane surfaces intersecting in a line. (B) Two plane surfaces which do not intersect.

the planes either intersect in a line so that the equations are consistent with an infinity of solutions represented by the line of intersection, or the surfaces do not intersect and the equations they represent are inconsistent. Consider the following system of equations: ⎡ ⎤

 x1

 ⎥ 1 2 3 4 ⎢ ⎢ x2 ⎥ = 1 −4 2 −3 7 ⎣ x3 ⎦ 3 x4 This under-determined system can be rearranged thus: 





  1 2 x1 3 4 1 x3 + = x2 x4 −4 2 −3 7 3 or



1 2 −4 2



x1 x2



=

1 3



−

3 4 −3 7



x3 x4



Thus, we have reduced this to a system of two equations in two unknowns, provided values are assumed for x3 and x4 . Since x3 and x4 can take an infinity of values, the problem has an infinity of solutions. If a system has more equations than unknowns, then the system is said to be over-determined. Fig. 2.4 shows four plane surfaces in three-dimensional space, representing four equations in three variables. Fig. 2.4A shows all four planes intersecting in a single point so that the system of equations is consistent with a unique solution. Fig. 2.4B shows all the planes intersecting in a line and this represents a consistent system with an infinity of solutions. Fig. 2.4D shows planes that represent an inconsistent system of equations with no solution. In Fig. 2.4C the planes do not intersect in a single point and so the system of equations is inconsistent. However, in this example the points of intersection of groups of three planes, (S1, S2, S3), (S1, S2, S4), (S1, S3, S4), and (S2, S3, S4), are close to each other and a mean point of intersection could be determined and used as an approximate solution. This example of marginal inconsistency often arises because the coefficients in the equations are determined experimentally; if the coefficients were known exactly, it is likely that the equations would be consistent with a unique solution. Rather than accepting that the system is inconsistent, we

2.3 OPERATORS \ AND / FOR SOLVING Ax = b

81

FIGURE 2.4 Planes representing an over-determined system of equations. (A) Four plane surfaces intersecting in a point. (B) Four plane surfaces intersecting in a line. (C) Four plane surfaces not intersecting at a single point. Point of intersection of (S1, S2, S3) and (S1, S2, S4) are visible. (D) Four plane surfaces representing inconsistent equations.

may ask which solution best satisfies the equations approximately. In Sections 2.11 and 2.12 we deal with the problem of over-determined and under-determined systems in more detail.

2.3 OPERATORS \ AND / FOR SOLVING Ax = b The purpose of this section is to introduce the reader to the M ATLAB operator \. A detailed discussion of the algorithms behind its operation will be given in later sections. This operator is a very powerful one which provides a unified approach to the solution of many categories of linear equation systems. The operators / and \ perform matrix “division” and have identical effects. Thus, to solve Ax = b we may write either x=A\b or x’=b’/A’. In the latter case the solution x is expressed as a row rather than a column vector. The algorithm chosen to solve a system of linear equations is dependent on the form of the matrix A. These cases are outlined as follows: • if A is a triangular matrix, the system is solved by back or forward substitution alone, described in Section 2.6.

82

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

• elseif A is a positive definite, square symmetric or Hermitian matrix (see Appendix A.7 for definitions of these matrices) Cholesky decomposition (described in Section 2.8) is applied. When A is sparse, Cholesky decomposition is preceded by a symmetric minimum degree preordering (described in Section 2.14). • elseif A is a square matrix, general LU decomposition (described in Section 2.7) is applied. If A is sparse, this is preceded by a non-symmetric minimum degree preordering (described in Section 2.14). • elseif A is a full non-square matrix, QR decomposition (described in Section 2.9) is applied. • elseif A is a sparse non-square matrix, it is augmented and then a minimum degree preordering is applied, followed by sparse Gaussian elimination (described in Section 2.14). The M ATLAB \ operator can also be used to solve AX = B where B and the unknown X are m × n matrices. This could provide a simple method of finding the inverse of A. If we make B the identity matrix I then we have AX = I and X must be the inverse of A since AA−1 = I. Thus, in M ATLAB we could determine the inverse of A by using the statement A\eye(size(A)). Alternatively, M ATLAB provides the function inv(A) to find the inverse of a matrix. It is important to stress that the inverse of a matrix should only be determined if it is specifically required. If we require the solution of a set of linear equations it is more efficient to use the operators \ or /. We now give some example which show how the \ operator works, beginning with the solution of a system where the system matrix is triangular. The experiment in this case examines the time taken by the operator \ to solve a system when it is full and then when the same system is converted to triangular form by zeroing appropriate elements to produce a triangular matrix. The script e4s201, used for this experiment, is % e4s201.m disp(’

n

full-time

full-time/n^3

tri-time

tri-time/n^2’);

A = [ ]; b = [ ]; for n = 2000:500:6000 A = 100*rand(n); b = [1:n].’; tic, x = A\b; t1 = toc; t1n = 5e9*t1/n^3; for i = 1:n for j = i+1:n A(i,j) = 0; end end tic, x = A\b; t2 = toc; t2n = 1e9*t2/n^2; fprintf(’%6.0f %9.4f %12.4f %12.4f %11.4f\n’,n,t1,t1n,t2,t2n) end

The results for a series of randomly generated n × n matrices are as follows:

2.3 OPERATORS \ AND / FOR SOLVING Ax = b

n 2000 2500 3000 3500 4000 4500 5000 5500 6000

full-time 0.3845 0.7833 0.9172 1.3808 2.0325 2.4460 3.5037 4.7222 5.9742

full-time/n^3 0.2403 0.2507 0.1698 0.1610 0.1588 0.1342 0.1401 0.1419 0.1383

tri-time 0.0080 0.0089 0.0122 0.0178 0.0243 0.0277 0.0344 0.0394 0.0556

83

tri-time/n^2 2.0041 1.4182 1.3560 1.4504 1.5215 1.3696 1.3741 1.3016 1.5444

Column 1 of this table gives n, the size of the square matrix. Columns 2 and 4 give the times taken to solve the full system and the triangular system respectively. Note that the triangular times are much smaller. Column 3 gives the time divided by n3 and shows that the time take is closely proportional to n3 , whereas Column 5 is multiplied by 109 and divided by n2 . We see that the time required to solve a triangular system is approximately proportional to n2 . We now perform experiments to examine the effects of using the operator \ with positive definite symmetric systems. This is a more complex problem than those previously discussed and the script e4s202.m implements this test. It is based on comparing the application of the \ operator to a positive definite system and a non-positive definite system of equations. We can create a positive define matrix by letting A = M*M’ where M is any matrix, but in this case M will be a matrix of random numbers, A will then be a positive define matrix. To generate a non-positive definite system we add a random matrix to the positive define matrix and we compare the time required to solve the two forms of matrix. The script for e4s202.m is of the form % e4s202.m disp(’ n time-pos time-pos/n^3 time-npos time-b/n^3’); for n = 500:500:5000 A = [ ]; M = 100*randn(n,n); A = M*M’; b = [1:n].’; tic, x = A\b; t1 = toc*1000; t1d = t1/n^3; A = A+rand(size(A)); tic, x = A\b; t2 = toc*1000; t2d = t2/n^3; fprintf(’%4.0f %10.4f %14.4e %11.4f %13.4e\n’,n,t1,t1d,t2,t2d) end

The result of running this script is n 500 1000 1500 2000 2500

time-pos 9.9106 46.7403 122.6701 191.2798 374.3233

time-pos/n^3 7.9285e-08 4.6740e-08 3.6347e-08 2.3910e-08 2.3957e-08

time-npos 10.6241 60.0436 164.9225 334.4048 736.9797

time-b/n^3 8.4993e-08 6.0044e-08 4.8866e-08 4.1801e-08 4.7167e-08

84

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

3000 3500 4000 4500 5000

619.7430 769.8043 1040.6866 1474.6506 1934.6456

2.2953e-08 1.7955e-08 1.6261e-08 1.6183e-08 1.5477e-08

974.7610 1354.9658 1976.1037 2518.8017 3437.0971

3.6102e-08 3.1603e-08 3.0877e-08 2.7641e-08 2.7497e-08

Column 1 of this table gives n, the size of the matrix. Column 2 gives the time taken multiplied by 1000 for the positive definite matrix and column 4 gives the time taken, again multiplied by 1000, for the non-positive definite matrix. These results show that the time taken to determine the solution for the system is somewhat faster for the positive definite system. This is because the operator \ checks to see if the matrix is positive definite and if so uses the more efficient Cholesky decomposition. Columns 3 and 5 give the times divided by the size of the matrix cubed to illustrate that the processing times are approximately proportional to n3 . The next test we perform examines how the operator \ succeeds with the very badly conditioned Hilbert matrix. The test gives the time taken to solve the system and the accuracy of the solution given by the Euclidean norm of the residuals, that is norm(Ax − b). For the definition of the norm see Appendix A, Section A.10. In addition the test compares these results for the \ operator with the results obtained using the inverse, that is x = A−1 b. The script e4s203.m performs this test. % e4s203.m disp(’ n time-slash acc-slash for n = 4:2:20 A = hilb(n); b = [1:n].’; tic, x = A\b; t1 = toc; t1 = nm1 = norm(b-A*x); tic, x = inv(A)*b; t2 = toc; nm2 = norm(b-A*x); c = cond(A); fprintf(’%2.0f %10.4f %10.2e end

time-inv

acc-inv

condition’);

t1*10000; t2 = t2*10000;

%8.4f %11.2e %11.2e \n’,n,t1,nm1,t2,nm2,c)

Script e4s203.m produces the following table of results: n 4 6 8 10 12 14 16 18 20

time-slash acc-slash 0.8389 4.11e-14 0.4730 2.89e-12 0.4462 1.89e-10 0.4641 9.81e-09 10.2279 5.66e-07 9.0498 1.33e-05 8.7687 4.57e-05 8.7776 2.25e-05 9.5094 6.91e-05

time-inv 0.6961 0.6649 0.7229 0.8479 9.4157 8.5991 8.6035 8.6660 8.9784

acc-inv 8.89e-13 3.74e-10 5.60e-07 2.97e-04 3.42e-01 5.33e+01 1.01e+03 1.92e+02 7.04e+01

condition 1.55e+04 1.50e+07 1.53e+10 1.60e+13 1.62e+16 2.55e+17 4.89e+17 1.35e+18 2.11e+18

This output has been edited to remove warnings about the ill-conditioning of the matrix for n >= 10. Column 1 gives the size of the matrix. Columns 2 and 3 give the time taken multiplied by 10,000 and

2.4 ACCURACY OF SOLUTIONS AND ILL-CONDITIONING

85

the accuracy when using the \ operator, respectively. Columns 4 and 5 give the same information when using the inv function. Column 6 gives the condition number of the system matrix. When the condition number is large, the matrix is nearly singular and the equations are ill-conditioned. Ill-conditioning is fully described in Section 2.4. The results in the preceding table demonstrate convincingly the superiority of the \ operator over the inv function for solving a system of linear equations. It is considerably more accurate than using matrix inversion. However, it should be noted that the accuracy falls off as the matrix becomes increasingly ill-conditioned. The M ATLAB operator \ can also be used to solve under- and over-determined systems. In this case the operator \ uses a least squares approximation, discussed in detail in Section 2.12.

2.4 ACCURACY OF SOLUTIONS AND ILL-CONDITIONING We now consider factors which affect the accuracy of the solution of Ax = b and how any inaccuracies can be detected. A further discussion on the accuracy of the solution of this equation system is given in Appendix B, Section B.3. We begin with the following examples. Example 2.1. Consider the following M ATLAB statements: >> A = [3.021 2.714 6.913;1.031 -4.273 1.121;5.084 -5.832 9.155] A = 3.0210

2.7140

6.9130

1.0310

-4.2730

1.1210

5.0840

-5.8320

9.1550

>> b = [12.648 -2.121 8.407].’ b = 12.6480 -2.1210 8.4070 >> A\b ans = 1.0000 1.0000 1.0000

This result is correct and easily verified by substitution into the original equations.

86

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

Example 2.2. This example is obtained from Example 2.1 with A(2,2) changed from −4.2730 to −4.2750 thus: >> A(2,2) = -4.2750 A = 3.0210

2.7140

6.9130

1.0310

-4.2750

1.1210

5.0840

-5.8320

9.1550

>> A\b ans = -1.7403 0.6851 2.3212

Here we have a solution which is very different from that of Example 2.1, even though the only change in the equation system is less than 0.1% in coefficient a(2,2). The two examples just shown have dramatically different solutions because the coefficient matrix A is ill-conditioned. Ill-conditioning can be interpreted graphically by representing each of the equation systems by three plane surfaces, in the manner shown in Fig. 2.2. In an ill-conditioned system at least two of the surfaces will be almost parallel so that the point of intersection of the surfaces will be very sensitive to small changes in slope, caused by small changes in coefficient values. A system of equations is said to be ill-conditioned if a relatively small change in the elements of the coefficient matrix A causes a relatively large change in the solution. Conversely, a system of equations is said to be well-conditioned if a relatively small change in the elements of the coefficient matrix A causes a relatively small change in the solution. Clearly, we require a measure of the condition of a system of equations. We know that a system of equations without a solution – the very worst condition possible – has a coefficient matrix with a determinant of zero. It is therefore tempting to think that the size of the determinant of A can be used as a measure of condition. However, if Ax = b and A is an n × n diagonal matrix with each element on the leading diagonal equal to s, then A is perfectly conditioned, regardless of the value of s. But the determinant of A in this case is s n . Thus, the size of the determinant of A is not a suitable measure of condition because in this example it changes with s even though the condition of the system is constant. Two of the functions M ATLAB provides to estimate the condition of a matrix are cond and rcond. The function cond is a sophisticated function and is based on singular value decomposition, discussed in Section 2.10. For a perfect condition cond is unity but gives a large value for a matrix which is illconditioned. The function rcond is less reliable but usually faster. This function gives a value between zero and one. The smaller the value, the worse the conditioning. The reciprocal of rcond is usually of the same order of magnitude as cond. We now illustrate these points with two examples.

2.4 ACCURACY OF SOLUTIONS AND ILL-CONDITIONING

87

Example 2.3. Illustration of a perfectly conditioned system: >> A = diag([20 20 20]) A = 20 0 0 0 20 0 0 0 20 >> [det(A) rcond(A) cond(A)] ans = 8000

1

1

Example 2.4. Illustration of a badly conditioned system: >> A = [1 2 3;4 5 6;7 8 9.000001]; >> format short e >> [det(A) rcond(A) 1/rcond(A) cond(A)] ans = -3.0000e-06

6.9444e-09

1.4400e+08

1.0109e+08

Note that the reciprocal of the rcond value is close to the value of cond. Using the M ATLAB functions cond and rcond we now investigate the condition number of the Hilbert matrix (defined in Problem 2.1), using the script e4s204.m: % e4s204.m Hilbert matrix test. disp(’ n cond rcond log10(cond)’) for n = 4:2:20 A = hilb(n); fprintf(’%5.0f %16.4e’,n,cond(A)); fprintf(’%16.4e %10.2f\n’,rcond(A),log10(cond(A))); end

Running script e4s204.m gives n n 4 6 8 10 12 14 16 18 20

cond cond 1.5514e+04 1.4951e+07 1.5258e+10 1.6025e+13 1.6212e+16 2.5515e+17 4.8875e+17 1.3500e+18 2.1065e+18

rcond rcond 3.5242e-05 3.4399e-08 2.9522e-11 2.8285e-14 2.5334e-17 6.9286e-19 1.3275e-18 2.5084e-19 1.2761e-19

log10(cond) log10(cond) 4.19 7.17 10.18 13.20 16.21 17.41 17.69 18.13 18.32

88

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

This shows that the Hilbert matrix is ill-conditioned even for relatively small values of n, the size of the matrix. The last column of the output of e4s204.m gives the value of log10 of the condition number of the appropriate Hilbert matrix. This gives a rule of thumb estimate of the number of significant figures lost in solving an equation system with this matrix or inverting the matrix. The Hilbert matrix of order n was generated in the preceding script using the M ATLAB function hilb(n). Other important matrices with interesting structure and properties, such as the Hadamard matrix and the Wilkinson matrix can be obtained using, in these cases the M ATLAB function hadamard(n) and wilkinson(n) where n is the required size of the matrix. In addition, many other interesting matrices can be accessed using the gallery function. In almost every case, we can choose the size of the matrix and in many cases we can also choose other parameters within the matrix. Example calls are gallery(’hanowa’,6,4) gallery(’cauchy’,6) gallery(’forsythe’,6,8)

The next section begins the detailed examination of one of the algorithms used by the \ operator.

2.5 ELEMENTARY ROW OPERATIONS We now examine the operations that can usefully be carried out on each equation of a system of equations. Such a system will have the form a11 x1 + a12 x2 + . . . + a1n xn = b1 a21 x2 + a22 x2 + . . . + a2n xn = b2 ............... an1 xn + an2 x2 + . . . + ann xn = bn or in matrix notation Ax = b where ⎡

a11 ⎢a21 ⎢ A=⎢ . ⎣ .. an1

a12 . . . a22 . . . .. . an2 . . .

⎤ ⎡ ⎤ ⎡ ⎤ b1 a1n x1 ⎢ b2 ⎥ ⎢ x2 ⎥ a2n ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ b=⎢ . ⎥ x=⎢ . ⎥ .. ⎥ . ⎦ ⎣ ⎦ ⎣ .. ⎦ . . ann xn bn

A is called the coefficient matrix. Any operation performed on an equation must be applied to both its left and right sides. With this in mind it is helpful to combine the coefficient matrix A with the right

2.6 SOLUTION OF Ax = b BY GAUSSIAN ELIMINATION

side vector b thus:

⎡

a11 ⎢a21 ⎢ A=⎢ . ⎣ .. an1

a12 . . . a1n a22 . . . a2n .. .. . . an2 . . . ann

89

⎤ b1 b2 ⎥ ⎥ .. ⎥ .⎦ bn

This new matrix is called the augmented matrix and we will write it as [A b]. We have chosen to adopt this notation because it is consistent with M ATLAB notation for combining A and b. Note that if A is an n × n matrix, then the augmented matrix is an n × (n + 1) matrix. Each row of the augmented matrix holds all the coefficients of an equation and any operation must be applied to every element in the row. The three elementary row operations described in 1, 2, and 3 can be applied to individual equations in a system without altering the solution of the equation system. They are: 1. Interchange the position of any two rows (that is any two equations). 2. Multiply a row (that is an equation) by a non-zero scalar. 3. Replace a row by the sum of the row and a scalar multiple of another row. These elementary row operations can be used to solve some important problems in linear algebra and we now discuss an application of them.

2.6 SOLUTION OF Ax = b BY GAUSSIAN ELIMINATION Gaussian elimination is an efficient way to solve equation systems, particularly those with a nonsymmetric coefficient matrix having a relatively small number of zero elements. The method depends entirely on using the three elementary row operations, described in Section 2.5. Essentially the procedure is to form the augmented matrix for the system and then reduce the coefficient matrix part to an upper triangular form. To illustrate the systematic use of the elementary row operations we consider the application of Gaussian elimination to solve the following equation system: ⎡

3 6 ⎣ 2 (4 + p) −3 −4

⎤⎡ ⎤ ⎡ ⎤ 9 3 x1 2 ⎦ ⎣x2 ⎦ = ⎣ 4 ⎦ x3 −5 −11

(2.10)

where the value of p is known. Table 2.1 shows the sequence of operations, beginning at Stage 1 with the augmented matrix. In Stage 1 the element in the first column of the first row (enclosed in a box in the table) is designated the pivot. We wish to make the elements of column 1 in rows 2 and 3 zero. To achieve this, we divide row 1 by the pivot and then add or subtract a suitable multiple of the modified row 1 to or from rows 2 and 3. The result of this is shown in Stage 2 of the table. We then select the next pivot. This is the element in the second column of the new second row, which in Stage 2 is equal to p. If p is large, this does not present a problem, but if p is small, then numerical problems may arise because we will be dividing all the elements of the new row 2 by this small quantity p. If p is small then we are in effect multiplying the pivot by a large number. This not only multiplies the pivot element by this value but also amplifies any rounding error in the element by the same amount. If p

90

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

Table 2.1 Gaussian elimination used to transform an augmented matrix to upper triangular form A1 3 6 9 3 A2 2 (4 + p) 2 4 −3 −4 −11 −5 A3 A1 3 6 9 3 B2 = A2 − 2(A1)/3 0 p −4 2 B3 = A3 + 3(A1)/3 0 2 −2 −2 A1 3 6 9 3 0 2 −2 −2 B3 0 p −4 2 B2 A1 3 6 9 3 0 2 −2 −2 B3 C3 = B2 − p(B3)/2 0 0 (−4 + p) (2 + p)

Stage 1: Initial matrix Stage 2: Reduce col 1 of row 2 & 3 to zero Stage 3: Interchange rows 2 & 3 Stage 4: Reduce col 2 of row 3 to zero

is zero, then we have an impossible situation because we cannot divide by zero. This difficulty is not related to ill-conditioning; indeed this particular equation system is quite well conditioned when p is zero. To circumvent these problems the usual procedure is to interchange the row in question with the row containing the element of largest modulus in the column below the pivot. In this way we provide a new and larger pivot. This procedure is called partial pivoting. If we assume in this case that p < 2, then we interchange rows 2 and 3 as shown in Stage 3 of the table to replace p by 2 as the pivot. From row 3 we now subtract row 2 divided by the pivot and multiplied by a coefficient in order to make the element of column 2, row 3, zero. Thus, in Stage 4 of the table it can be seen that the original coefficient matrix has been reduced to an upper triangular matrix. If, for example, p = 0, we obtain 3x1 + 6x2 + 9x3 = 3 2x2 − 2x3 = −2 −4x3 = 2

(2.11) (2.12) (2.13)

We can obtain the values of the unknowns x1 , x2 , and x3 by a process called back substitution. We solve the equations in reversed order. Thus, from (2.13), x3 = −0.5. From (2.12), knowing x3 , we have x2 = −1.5. Finally from (2.11), knowing x2 and x3 , we have x1 = 5.5. It can be shown that the determinant of a matrix can be evaluated from the product of the elements on the main diagonal provided at Stage 3 in Table 2.1. This product must be multiplied by (−1)m where m is the number of row interchanges used. For example, in the preceding problem, with p = 0, one row interchange is used so that m = 1 and the determinant of the coefficient matrix is given by 3 × 2 × (−4) × (−1)1 = 24. A method for solving a linear equation system that is closely related to Gaussian elimination is Gauss–Jordan elimination. The method uses the same elementary row operations but differs from Gaussian elimination because elements both below and above the leading diagonal are reduced to zero. This means that back substitution is avoided. For example, solving system (2.10) with p = 0 leads to the

2.7 LU DECOMPOSITION

following augmented matrix:

⎡

3 ⎣ 0 0

91

⎤ 0 0 16.5 2 0 −3.0 ⎦ 0 −4 2.0

Thus x1 = 16.5/3 = 5.5, x2 = −3/2 = −1.5, and x3 = 2/ − 4 = −0.5. Gaussian elimination requires order n3 /3 multiplications followed by back substitution requiring order n2 multiplications. Gauss–Jordan elimination requires order n3 /2 multiplications. Thus, for large systems of equations (say n > 10), Gauss–Jordan elimination requires approximately 50% more operations than Gaussian elimination.

2.7 LU DECOMPOSITION LU decomposition (or factorization) is a similar process to Gaussian elimination and is equivalent in terms of elementary row operations. The matrix A can be decomposed so that A = LU

(2.14)

where L is a lower triangular matrix with a leading diagonal of ones and U is an upper triangular matrix. Matrix A may be real or complex. Compared with Gaussian elimination, LU decomposition has a particular advantage when the equation system we wish to solve, Ax = b, has more than one right side or when the right sides are not known in advance. This is because the factors L and U are obtained explicitly and they can be used for any right sides as they arise without recalculating L and U. Gaussian elimination does not determine L explicitly but rather forms L−1 b so that all right sides must be known when the equation is solved. The major steps required to solve an equation system by LU decomposition are as follows. Since A = LU, then Ax = b becomes LUx = b where b is not restricted to a single column. Letting y = Ux leads to Ly = b Because L is a lower triangular matrix this equation is solved efficiently by forward substitution. To find x we then solve Ux = y Because U is an upper triangular matrix, this equation can also be solved efficiently by back substitution. We now illustrate the LU decomposition process by solving (2.10) with p = 1. We are not concerned with b and we do not form an augmented matrix. We proceed exactly as with Gaussian elimination, see Table 2.1, except that we keep a record of the elementary row operations performed at the ith stage in T(i) and place the results of these operations in a matrix U(i) rather than over-writing A.

92

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

We begin with the matrix

⎡

3 ⎣ A= 2 −3

⎤ 6 9 5 2 ⎦ −4 −11

(2.15)

Following the same operations as used in Table 2.1, we will create a matrix U(1) with zeros below the leading diagonal in the first column using the following elementary row operations: row 2 of U(1) = row 2 of A − 2(row 1 of A)/3

(2.16)

row 3 of U(1) = row 3 of A + 3(row 1 of A)/3

(2.17)

and T(1)

U(1)

as follows: Now A can be expressed as the product ⎡ ⎤ ⎡ ⎤⎡ 3 6 9 1 0 0 3 6 ⎣ 2 5 2 ⎦ = ⎣ 2/3 1 0 ⎦ ⎣ 0 1 −3 −4 −11 −1 0 1 0 2

⎤ 9 −4 ⎦ −2

Note that row 1 of A and row 1 of U(1) are identical. Thus row 1 of T(1) has a unit entry in column 1 and zero elsewhere. The remaining rows of T(1) are determined from (2.16) and (2.17). For example, row 2 of T(1) is derived by rearranging (2.16); thus: row 2 of A = row 2 of U(1) − 2(row 1 of A)/3

(2.18)

row 2 of A = 2(row 1 of U(1) )/3 + row 2 of U(1)

(2.19)

or is identical to row 1 of A. Hence row 2 of is [2/3 1 0]. since row 1 of We now move to the next stage of the decomposition process. In order to bring the largest element of column 2 in U(1) onto the leading diagonal we must interchange rows 2 and 3. Thus U(1) becomes the product T(2) U(2) as follows: ⎤⎡ ⎤ ⎡ ⎤ ⎡ 1 0 0 3 6 9 3 6 9 ⎣ 0 1 −4 ⎦ = ⎣ 0 0 1 ⎦ ⎣ 0 2 −2 ⎦ 0 1 0 0 1 −4 0 2 −2 U(1)

T(1)

Finally, to complete the process of obtaining an upper triangular matrix we make row 3 of U = row 3 of U(2) − (row 2 of U(2) )/2 Hence, U(2) becomes the product T(3) U as follows: ⎡ ⎤ ⎡ ⎤⎡ ⎤ 3 6 9 1 0 0 3 6 9 ⎣ 0 2 −2 ⎦ = ⎣ 0 1 0 ⎦ ⎣ 0 2 −2 ⎦ 0 1 −4 0 1/2 1 0 0 −3

2.7 LU DECOMPOSITION

93

Thus A = T(1) T(2) T(3) U, implying that L = T(1) T(2) T(3) as follows: ⎡ ⎤⎡ ⎤⎡ ⎤ ⎡ ⎤ 1 0 0 1 0 0 1 0 0 1 0 0 ⎣ 2/3 1 0 ⎦ ⎣ 0 0 1 ⎦ ⎣ 0 1 0 ⎦ = ⎣ 2/3 1/2 1 ⎦ −1 0 1 0 1 0 0 1/2 1 −1 1 0 Note that owing to the row interchanges L is not strictly a lower triangular matrix but it can be made so by interchanging rows. M ATLAB implements LU factorization by using the function lu and may produce a matrix that is not strictly a lower triangular matrix. However, a permutation matrix P may be produced, if required, such that LU = PA with L lower triangular. We now show how the M ATLAB function lu solves the example based on the matrix given in (2.15): >> A = [3 6 9;2 5 2;-3 -4 -11] A = 3

6

2

5

9 2

-3

-4

-11

To obtain the L and U matrices, we must use that M ATLAB facility of assigning two parameters simultaneously as follows: >> [L1 U] = lu(A) L1 = 1.0000

0

0

0.6667

0.5000

1.0000

-1.0000

1.0000

0

U = 3.0000

6.0000

9.0000

0

2.0000

-2.0000

0

0

-3.0000

Note that the L1 matrix is not in lower triangular form, although its true form can easily be deduced by interchanging rows 2 and 3 to form a triangle. To obtain a true lower triangular matrix we must assign three parameters as follows: >> [L U P] = lu(A) L = 1.0000

0

0

-1.0000

1.0000

0

0.6667

0.5000

1.0000

94

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

U = 3.0000

6.0000

9.0000

0

2.0000

-2.0000

0

0

-3.0000

P = 1

0

0

0

0

1

0

1

0

In the preceding output, P is the permutation matrix such that L*U = P*A or P’*L*U = A. Thus P’*L is equal to L1. The M ATLAB operator \ determines the solution of Ax = b using LU factorization. As an example of an equation system with multiple right sides we solve AX = B where ⎤ ⎤ ⎡ 3 4 −5 1 3 A = ⎣ 6 −3 4 ⎦ and B = ⎣ 9 5 ⎦ 8 9 −2 9 4 ⎡

Performing LU decomposition such that LU = A gives ⎡

⎤ ⎡ ⎤ 0.375 −0.064 1 8 9 −2 L = ⎣ 0.750 1 0 ⎦ and U = ⎣ 0 −9.75 5.5 ⎦ 1 0 0 0 0 −3.897 Thus LY = B is given by ⎡

⎤⎡ 0.375 −0.064 1 y11 ⎣ 0.750 1 0 ⎦ ⎣ y21 y31 1 0 0

⎤ ⎡ ⎤ 1 3 y12 y22 ⎦ = ⎣ 9 5 ⎦ y32 9 4

We note that implicitly we have two systems of equations which when separated can be written: ⎡

⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 3 y11 y12 L ⎣ y21 ⎦ = ⎣ 9 ⎦ and L ⎣ y22 ⎦ = ⎣ 5 ⎦ y31 y32 9 4 In this example, L is not strictly a lower triangular matrix owing to the reordering of the rows. However, the solution of this equation is still found by forward substitution. For example, 1y11 = b31 = 9, so that y11 = 9. Then 0.75y11 + 1y21 = b21 = 9. Hence y21 = 2.25, etc. The complete Y matrix is ⎡

⎤ 9.000 4.000 Y = ⎣ 2.250 2.000 ⎦ −2.231 1.628

2.8 CHOLESKY DECOMPOSITION

95

Finally solving UX = Y by back substitution gives ⎡ ⎤ 1.165 0.891 X = ⎣ 0.092 −0.441 ⎦ 0.572 −0.418 The M ATLAB function det determines the determinant of a matrix using LU factorization as follows. Since A = LU then |A| = |L| |U|. The elements of the leading diagonal of L are all ones so that |L| = 1. Since U is upper triangular, its determinant is the product of the elements of its leading diagonal. Thus, taking account of row interchanges the appropriately signed product of the diagonal elements of U gives the determinant.

2.8 CHOLESKY DECOMPOSITION Cholesky decomposition or factorization is a form of triangular decomposition that can only be applied to either a positive definite symmetric matrix or a positive definite Hermitian matrix. A symmetric matrix A is said to be positive definite if xT Ax > 0 for any non-zero x. Similarly, if A is Hermitian, then xH Ax > 0. A more useful definition of a positive definite matrix is one that has all eigenvalues greater than zero. The eigenvalue problem is discussed in Section 2.15. If A is symmetric or Hermitian, we can write A = PT P (or A = PH P when A is Hermitian)

(2.20)

where P is an upper triangular matrix. The algorithm computes P row by row by equating coefficients of each side of (2.20). Thus p11 , p12 , p13 , ..., p22 , p23 , ... are determined in sequence, ending with pnn . Coefficients on the leading diagonal of P are computed from expressions that involve determining a square root. For example, 2 p22 = a22 − p12 A property of positive definite matrices is that the term under the square root is always positive and so the square root will be real. Furthermore, row interchanges are not required because the dominant coefficients will always be on the main diagonal. The whole process requires only about half as many multiplications as LU decomposition. Cholesky factorization is implemented for positive definite symmetric matrices in M ATLAB by the function chol. For example, consider the Cholesky factorization of the following positive definite Hermitian matrix: >> A = [2 -i 0;i 2 0;0 0 3] A = 2.0000 + 0.0000i

0.0000 - 1.0000i

0.0000 + 0.0000i

0.0000 + 1.0000i

2.0000 + 0.0000i

0.0000 + 0.0000i

0.0000 + 0.0000i

0.0000 + 0.0000i

3.0000 + 0.0000i

>> P = chol(A)

96

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

P = 1.4142 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i

0.0000 - 0.7071i 1.2247 + 0.0000i 0.0000 + 0.0000i

0.0000 + 0.0000i 0.0000 + 0.0000i 1.7321 + 0.0000i

When the operator \ detects a symmetric positive definite or Hermitian positive definite system matrix, it solves Ax = b using the following sequence of operations. A is factorized into PT P, y is set to Px; then PT y = b. The algorithm solves for y by forward substitution since PT is a lower triangular matrix. Then x can be determined from y by backward substitution since P is an upper triangular matrix. We can illustrate the steps in this process by the following example: ⎤ ⎤ ⎡ ⎡ 2 2 3 4 A=⎣ 3 6 7 ⎦ and b = ⎣ 4 ⎦ 8 4 7 10 Then by Cholesky factorization ⎡

⎤ 1.414 2.121 2.828 P=⎣ 0 1.225 0.817 ⎦ 0 0 1.155 Now since PT y = b solving for y by forward substitution gives ⎡ ⎤ 1.414 y = ⎣ 0.817 ⎦ 2.887 Finally solving Px = y by back substitution gives ⎡

⎤ −2.5 x = ⎣ −1.0 ⎦ 2.5

We now compare the performance of the operator \ with the function chol. Clearly, their performance should be similar in the case of a positive definite matrix. To generate a symmetric positive definite matrix in script e4s205.m, we multiply a matrix by its transpose: % e4s205.m disp(’ n time-backslash time-chol’); for n = 1000:100:2000 A = [ ]; M = 100*randn(n,n); A = M*M’; b = [1:n].’; tic, x = A\b; t1 = toc; tic, R = chol(A); v = R.’\b; x = R\b;

2.9 QR DECOMPOSITION

97

t2 = toc; fprintf(’%4.0f %14.4f %13.4f \n’,n,t1,t2) end

Running script e4s205.m gives n

time-backslash

time-chol

1000

0.0320

0.0232

1100

0.0388

0.0246

1200

0.0447

0.0271

1300

0.0603

0.0341

1400

0.0650

0.0475

1500

0.0723

0.0420

1600

0.0723

0.0521

1700

0.0965

0.0552

1800

0.1142

0.0679

1900

0.1312

0.0825

2000

0.1538

0.1107

The similarity in performance of the function chol and the operator \ is borne out in the preceding table. In this table, column 1 gives the size of the matrix and column 2 gives the time taken using the \ operator. Column 3 gives the time taken using Cholesky decomposition to solve the same problem. Cholesky factorization can be applied to a symmetric matrix which is not positive definite but the process does not possess the numerical stability of the positive definite case. Furthermore, one or more rows in P may be purely imaginary. For example, ⎡

⎤ ⎡ 1 2 3 1 If A = ⎣ 2 −5 9 ⎦ then P = ⎣ 0 3 9 4 0

⎤ 2 3 3ı −ı ⎦ 0 2ı

This is not implemented in M ATLAB.

2.9 QR DECOMPOSITION We have seen how a square matrix can be decomposed or factorized into the product of a lower and an upper triangular matrix by the use of elementary row operations. An alternative decomposition of A is into an upper triangular matrix and an orthogonal matrix if A is real, or into an upper triangular matrix and a unitary matrix if A is complex. This is called QR decomposition. Thus A = QR where R is the upper triangular matrix and Q is the orthogonal, or the unitary matrix. If Q is orthogonal, Q−1 = QT , and if Q is unitary, Q−1 = QH . The preceding properties are very useful.

98

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

There are several procedures which provide QR decomposition; here we present Householder’s method. To decompose a real matrix, Householder’s method begins by defining a matrix P thus: P = I − 2wwT

(2.21)

where w is a column vector and P is symmetrical matrix. Provided wT w = 1, P is also orthogonal. The orthogonality can easily be verified by expanding the product PT P = PP as follows:

  PP = I − 2wwT I − 2wwT

 = I − 4wwT + 4wwT wwT

 = I − 4wwT + 4w wT w wT = I since wT w = 1. To decompose A into QR, we begin by forming the vector w1 from the coefficients of the first column of A as follows: wT1 = μ1 [(a11 − s1 ) a21 a31 . . . an1 ] where ⎛ ⎞1/2 n  1 aj21 ⎠ and s1 = ± ⎝ μ1 = √ 2s1 (s1 − a11 ) j =1 By substituting for μ1 and s1 in w1 it can be verified that the necessary orthogonality condition, wT1 w1 = 1, is satisfied. Substituting w1 into (2.21) we generate an orthogonal matrix P(1) . The matrix A(1) is now created from the product P(1) A. It can easily be verified that all elements in the first column of A(1) are zero except for the element on the leading diagonal which is equal to s1 . Thus, ⎡ ⎤ s1 + · · · + ⎢ 0 + ··· + ⎥ ⎢ ⎥ ⎢ (1) (1) .. .. ⎥ A = P A = ⎢ ... . . ⎥ ⎢ ⎥ ⎣ 0 + ··· + ⎦ 0 + ··· + In the matrix A(1) , + indicates a non-zero element. We now begin the second stage of the orthogonalization process by forming w2 from the coefficients of the second column of A(1) thus:    (1) (1) (1) (1) wT2 = μ2 0 a22 − s2 a32 a42 · · · an2

2.9 QR DECOMPOSITION

99

where aij are the coefficients of A and 1 μ2 = (1) 2s2 (s2 − a22 )

⎛ ⎞1/2 n  (1) and s2 = ± ⎝ (aj 2 )2 ⎠ j =2

Then the orthogonal matrix P(2) is generated from P(2) = I − 2w2 wT2 The matrix A(2) is then created from the product P(2) A(1) as follows: ⎡ ⎢ ⎢ ⎢ A(2) = P(2) A(1) = P(2) P(1) A = ⎢ ⎢ ⎣

s1 0 .. . 0 0

+ s2 .. . 0 0

⎤ ··· + ··· + ⎥ ⎥ .. ⎥ . ⎥ ⎥ ··· + ⎦ ··· +

Note that A(2) has zero elements in its first two columns except for the elements on and above the leading diagonal. We can continue this process n − 1 times until we obtain an upper triangular matrix R. Thus, R = P(n−1) . . . P(2) P(1) A

(2.22)

Note that since P(i) is orthogonal, the product P(n−1) ... P(2) P(1) is also orthogonal. We wish to determine the orthogonal matrix Q such that A = QR. Thus R = Q−1 A or R = QT A. Hence, from (2.22), QT = P(n−1) . . . P(2) P(1) Apart from the signs associated with the columns of Q and the rows of R, the decomposition is unique. These signs are dependent on whether the positive or negative square root is taken in determining s1 , s2 , etc. Complete decomposition of the matrix requires 2n3 /3 multiplications and n square roots. To illustrate this procedure consider the decomposition of the matrix ⎡ ⎤ 4 −2 7 A=⎣ 6 2 −3 ⎦ 3 4 4 Thus, s1 = μ1 =

√

√ 42 + 62 + 32 = 7.8102 1 2×7.8102×(7.8102−4)

= 0.1296

wT1 = 0.1296[(4 − 7.8102) 6 3] = [−0.4939 0.7777 0.3889]

100

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

Using (2.21) we generate P(1) and hence A(1) thus: ⎡ ⎤ 0.5121 0.7682 0.3841 P(1) = ⎣ 0.7682 −0.2097 −0.6049 ⎦ 0.3841 −0.6049 0.6976 ⎡ ⎤ 7.8102 2.0486 2.8168 A(1) = P(1) A = ⎣ 0 −4.3753 3.5873 ⎦ 0 0.8123 7.2936 Note that we have reduced the elements of the first column of A(1) below the leading diagonal to zero. We continue with the second stage thus: s2 = μ2 =

√ √

−4.37532 + 0.81232 = 4.4501

1 2×4.4501×(4.4501+4.3753)

= 0.1128

wT2 = 0.1128 [0 (−4.3753 − 4.4501) 0.8123] = [0 −0.9958 0.0917] ⎡ ⎤ 1 0 0 P(2) = ⎣ 0 −0.9832 0.1825 ⎦ 0 0.1825 0.9832 ⎡ ⎤ 7.8102 2.0486 2.8168 R = A(2) = P(2) A(1) = ⎣ 0 4.4501 −2.1956 ⎦ 0 0 7.8259 Note that we have now reduced the first two columns of A(2) below the leading diagonal to zero. This completes the process to determine the upper triangular matrix R. Finally, we determine the orthogonal matrix Q as follows: ⎡ ⎤ 0.5121 −0.6852 0.5179

T Q = P(2) P(1) = ⎣ 0.7682 0.0958 −0.6330 ⎦ 0.3841 0.7220 0.5754 It is not necessary for the reader to carry out the preceding calculations since M ATLAB provides the function qr to carry out this decomposition. For example, >> A = [4 -2 7;6 2 -3;3 4 4] A = 4

-2

7

6

2

-3

3

4

4

>> [Q R] = qr(A)

2.10 SINGULAR VALUE DECOMPOSITION

101

Q = -0.5121

0.6852

0.5179

-0.7682

-0.0958

-0.6330

-0.3841

-0.7220

0.5754

-7.8102

-2.0486

-2.8168

0

-4.4501

2.1956

0

0

7.8259

R =

Note that the matrices Q and R in the M ATLAB output are the negative of the hand calculations of Q and R above. This is not significant since their product is equal to A, and in the multiplication, the signs cancel. One advantage of QR decomposition is that it can be applied to non-square matrices, decomposing an m × n matrix into an m × m orthogonal matrix and an m × n upper triangular matrix. Note that if m > n, the decomposition is not unique.

2.10 SINGULAR VALUE DECOMPOSITION The singular value decomposition (SVD) of an m × n matrix A is given by A = USVT (or A = USVH if A is complex)

(2.23)

where U is an orthogonal m × m matrix and V is an orthogonal n × n matrix. If A is complex then U and V are unitary matrices. In all cases S is a real diagonal m × n matrix. The elements of the leading diagonal of this matrix are called the singular values of A. Normally they are arranged in decreasing value so that s1 > s2 > ... > sn . Thus, ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ S=⎢ ⎢ ⎢ ⎢ ⎢ ⎣

s1 0 .. . 0 0 .. . 0

0 s2 .. . 0 0 .. . 0

... ... ... ... ...

0 0 .. . sn 0 .. . 0

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

The singular values are the non-negative square roots of the eigenvalues of AT A. Because AT A is symmetric or Hermitian these eigenvalues are real and non-negative so that the singular values are

102

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

also real and non-negative. Algorithms for computing the SVD of a matrix are given by Golub and Van Loan (1989). The SVD of a matrix has several important applications. In Section 2.4, we introduced the reduced row echelon form of a matrix and explained how the M ATLAB function rref gave information from which the rank of a matrix can be deduced. However, rank can be more effectively determined from the SVD of a matrix since its rank is equal to the number of its non-zero singular values. Thus, for a 5 × 5 matrix of rank 3, s4 and s5 would be zero. In practice, rather than counting the non-zero singular values, M ATLAB determines rank from the SVD by counting the number of singular values greater than some tolerance value. This is a more realistic approach to determining rank than counting any non-zero value, however small. To illustrate how singular value decomposition helps us to examine the properties of a matrix we will use the M ATLAB function svd to carry out a singular value decomposition and compare it with the function rref. Consider the following example in which a Vandermonde matrix is created using the M ATLAB function vander. The Vandermonde matrix is known to be ill-conditioned. SVD allows us to examine the nature of this ill-conditioning. In particular a zero or a very small singular value indicates rank deficiency and this example shows that the singular values are becoming relatively close to this condition. In addition SVD allows us to compute the condition number of the matrix. In fact, the M ATLAB function cond uses SVD to compute the condition number and this gives the same values as obtained by dividing the largest singular value by the smallest singular value. Additionally, the Euclidean norm of the matrix is supplied by the first singular value. Comparing the SVD with the RREF process in the following script, we see that the result of using the M ATLAB functions rref and rank give the rank of this special Vandermonde matrix as 5 but tells us nothing else. There is no warning that the matrix is badly conditioned. >> c = [1 1.01 1.02 1.03 1.04]; >> V = vander(c) V = 1.0000 1.0406 1.0824 1.1255 1.1699

1.0000 1.0303 1.0612 1.0927 1.1249

>> format long >> s = svd(V) s = 5.210367051037899 0.101918335876689 0.000699698839445 0.000002352380295 0.000000003294983 >> norm(V)

1.0000 1.0201 1.0404 1.0609 1.0816

1.0000 1.0100 1.0200 1.0300 1.0400

1.0000 1.0000 1.0000 1.0000 1.0000

2.10 SINGULAR VALUE DECOMPOSITION

103

ans = 5.210367051037898 >> cond(V) ans = 1.581303244929480e+09 >> s(1)/s(5) ans = 1.581303244929480e+09 >> rank(V) ans = 5 >> rref(V) ans = 1 0 0 0 0

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

The following example is very similar to the preceding one but the Vandermonde matrix has now been generated to be rank deficient. The smallest singular value, although not zero, is zero to machine precision and rank returns the value of 4. >> c = [1 1.01 1.02 1.03 1.03]; >> V = vander(c) V = 1.0000 1.0406 1.0824 1.1255 1.1255

1.0000 1.0303 1.0612 1.0927 1.0927

>> format long e >> s = svd(V)

1.0000 1.0201 1.0404 1.0609 1.0609

1.0000 1.0100 1.0200 1.0300 1.0300

1.0000 1.0000 1.0000 1.0000 1.0000

104

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

s = 5.187797954424026e+00 8.336322098941423e-02 3.997349250041025e-04 8.462129966394409e-07 1.261104176071465e-23 >> format short >> rank(V) ans = 4 >> rref(V) ans = 1.0000 0 0 0 0

0 1.0000 0 0 0

0 0 1.0000 0 0

0 0 0 1.0000 0

-0.9424 3.8262 -5.8251 3.9414 0

>> cond(V)

ans = 4.1137e+23

The rank function does allow the user to vary the tolerance using a second parameter, that is rank(A,tol) where tol gives the tolerance for acceptance of the singular values of A to be taken as non-zero. However, tolerance should be used with care since the rank function counts the number of singular values greater than tolerance and this gives the rank of the matrix. If tolerance is very small, that is, smaller than the machine precision, the rank may be miscounted.

2.11 THE PSEUDO-INVERSE Here we discuss the pseudo-inverse and in Section 2.12 we apply it to solve over- and under-determined systems. If A is an m × n rectangular matrix, then the system Ax = b

(2.24)

cannot be solved by inverting A, since A is not a square matrix. Assuming an equation system with more equations than variables, that is m > n, then by pre-multiplying (2.24) by AT we can convert the

2.11 THE PSEUDO-INVERSE

105

system matrix to a square matrix as follows: AT Ax = AT b The product AT A is square and, provided it is non-singular, it can be inverted to give the solution to (2.24) thus:

−1 x = AT A AT b (2.25) Let

−1 A+ = A T A AT

(2.26)

The matrix A+ is called the Moore–Penrose pseudo-inverse of A or just the pseudo inverse. Thus, the solution of (2.24) is   (2.27) x = A+ b This definition of the pseudo-inverse, A+ , requires A to have full rank. If A is full rank and m > n, then rank(A) = n. Now rank(AT A) = rank(A) and hence rank(AT A) = n. Since AT A is an n × n array, AT A is automatically full rank and A+ is then a unique m × n array. If A is rank deficient, then AT A is rank deficient and cannot be inverted. If A is square and non-singular, then A+ = A−1 . If A is complex then

−1 AH A+ = A H A

(2.28)

where AH is the conjugate transpose, described in Appendix A, Section A.6. The product AT A has a condition number which is the square of the condition number of A. This has implications for the computations involved in A+ . It can be shown that the pseudo-inverse A+ has the following properties: 1. A(A+ )A = A, 2. (A+ )A(A+ ) = A+ , 3. (A+ )A and A(A+ ) are symmetrical matrices. We must now consider the situation that occurs when A of (2.24) is an m × n array with m < n, that is an equation system with more variables than equations. If A is a full rank, then rank(A) = m. Now rank(AT A) = rank(A) and hence rank(AT A) = m. Since AT A is an n × n matrix, AT A is rank deficient and cannot be inverted, even though A is of full rank. We can avoid this problem by recasting (2.24) as follows:

 −1 b Ax = AAT AAT and hence, multiplying by A−1

−1 x = AT AAT b

106

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

Thus,

  x = A+ b

where A+ = AT (AAT )−1 and is the pseudo-inverse. Note that AAT is an m × m array with rank m and can thus be inverted. Example 2.5. Consider the following matrix: ⎡

1 ⎢ 4 A=⎢ ⎣ 5 −2

2 5 6 3

⎤ 3 9 ⎥ ⎥ 7 ⎦ 1

Computing the pseudo-inverse of A using a M ATLAB implementation of (2.26) we have >> A = [1 2 3;4 5 9;5 6 7;-2 3 1]; >> rank(A) ans = 3

We note that A is full rank. Thus, >> A_plus = inv(A.’*A)*A.’ A_plus = -0.0747 -0.0378 0.0858

-0.1467 -0.2039 0.2795

0.2500 0.2500 -0.2500

-0.2057 0.1983 -0.0231

The M ATLAB function pinv provides this result directly and with greater accuracy. A*A_plus*A ans = 1.0000 4.0000 5.0000 -2.0000

2.0000 5.0000 6.0000 3.0000

3.0000 9.0000 7.0000 1.0000

0.2841 0.9096 0.0000 -0.0387

0.0000 0.0000 1.0000 -0.0000

>> A*A_plus ans = 0.1070 0.2841 0.0000 0.1218

0.1218 -0.0387 -0.0000 0.9834

2.11 THE PSEUDO-INVERSE

107

>> A_plus*A ans = 1.0000

0.0000

0.0000

1.0000

0.0000 0.0000

-0.0000

-0.0000

1.0000

Note that these calculations verify that A*A_plus*A equals A and that both A*A_plus and A_plus*A are symmetrical. It has been shown that if A is rank deficient, then (2.26) cannot be used to determine the pseudoinverse of A. This doesn’t mean that the pseudo-inverse does not exist; it always exists but we must use a different method to evaluate it. When A is rank deficient, or close to rank deficient, A+ is best calculated from the singular value decomposition (SVD) of A. If A is real, the SVD of A is USVT where U is an orthogonal m × m matrix and V is an orthogonal n × n matrix and S is an n × m matrix of singular values. Thus, the SVD of AT is VST UT so that AT A = (VST UT )(USVT ) = VST SVT since UT U = I Hence, A+ = (VST SVT )−1 VST UT = V−T (ST S)−1 V−1 VST UT = V(ST S)−1 ST UT

(2.29)

We note that because of orthogonality VVT = I and hence VVT V−T = V−T . Since V is an m × m matrix, U is an n × n matrix and S is an n × m matrix, then (2.29) is conformable, that is matrix multiplication is possible; see Appendix A, Section A.5. In this situation, however, ST S cannot be inverted because of the very small or zero singular values. To deal with this problem we take only the r non-zero singular values of the matrix so that S is an r × r matrix where r is the rank of A. To make the multiplications of (2.29) conformable (that is, possible) we take the first r columns of V and the first r rows of UT , that is the first r columns of U. This is illustrated in Example 2.6 where the pseudo-inverse of A is determined. Example 2.6. Consider the following rank deficient matrix: ⎡ ⎢ ⎢ A=⎢ ⎢ ⎣

Using M ATLAB, we have

⎤ 1 2 3 4 5 9 ⎥ ⎥ 7 11 18 ⎥ ⎥ −2 3 1 ⎦ 7 1 8

108

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

>> A = [1 2 3;4 5 9;7 11 18;-2 3 1;7 1 8] A = 1

2

4

5

3 9

7

11

18

-2

3

1

7

1

8

>> rank(A) ans = 2

Since A has a rank of 2 and is thus rank deficient. Hence we cannot use (2.26) to determine its pseudoinverse. We now find the SVD of A as follows >> [U S V] = svd(A) U = -0.1381

0.0839

-0.1198

0.4305

-0.8799

-0.4115

0.0215

0.9092

0.0489

-0.0333

-0.8258

0.2732

-0.3607

-0.3341

0.0413

-0.0524

0.5650

-0.0599

0.7077

0.4165

-0.3563

-0.7737

-0.1588

0.4469

0.2224

S = 26.8394

0

0

0

6.1358

0

0

0

0.0000

0

0

0

0

0

0

-0.3709

-0.7274

0.5774

-0.4445

0.6849

0.5774

-0.8154

-0.0425

-0.5774

V =

We now select the two significant singular values for use in the subsequent computation: >> SS = S(1:2,1:2) SS = 26.8394

0

0

6.1358

2.12 OVER- AND UNDER-DETERMINED SYSTEMS

109

To make the multiplication conformable we use only the first two columns of U and V thus: >> A_plus = V(:,1:2)*inv(SS.’*SS)*SS.’*U(:,1:2).’ A_plus = -0.0080 0.0117 0.0036

0.0031 0.0092 0.0124

-0.0210 0.0442 0.0232

-0.0663 0.0639 -0.0023

0.0966 -0.0805 0.0162

This result can be obtained directly using the pinv function which is based on the singular value decomposition of G. >> A*A_plus ans = 0.0261 0.0586 0.1369 0.0546 -0.0157

0.0586 0.1698 0.3457 0.0337 0.1300

0.1369 0.3457 0.7565 0.1977 0.0829

-0.3333 0.6667 0.3333

0.3333 0.3333 0.6667

0.0546 0.0337 0.1977 0.3220 -0.4185

-0.0157 0.1300 0.0829 -0.4185 0.7256

>> A_plus*A ans = 0.6667 -0.3333 0.3333

Note that A*A_plus and A_plus*A are symmetric. In Section 2.12, we will apply these methods to solve over- and under-determined systems and discuss the meaning of the solution.

2.12 OVER- AND UNDER-DETERMINED SYSTEMS We will begin by examining examples of over-determined systems, that is, systems of equations in which there are more equations than unknown variables. Although over-determined systems may have a unique solution, most often we are concerned with equation systems that are generated from experimental data which can lead to a relatively small degree of inconsistency between the equations. For example, consider the following over-determined system of linear equations: x1 + x2 = 1.98 2.05x1 − x2 = 0.95 3.06x1 + x2 = 3.98 −1.02x1 + 2x2 = 0.92 4.08x1 − x2 = 2.90

(2.30)

110

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

FIGURE 2.5 Plot of an inconsistent equation system (2.30).

Fig. 2.5 shows that (2.30) is such a system; the lines do not intersect in a point, although there is a point that nearly satisfies all the equations. We would like to choose the best of all possible solutions of (2.30) in the region defined by the intersections. One criterion for doing this is that the chosen solution should minimize the sum of squares of the residual errors (or residuals) of the equations. For example, consider the equation system (2.30). Letting r1 , ..., r5 be the residuals, then x1 + x2 − 1.98 = r1 2.05x1 − x2 − 0.95 = r2 3.06x1 + x2 − 3.98 = r3 −1.02x1 + 2x2 − 0.92 = r4 4.08x1 − x2 − 2.90 = r5 In this case the sum of the residuals squared is given by S=

5 

ri2

i=1

We wish to minimize S and we can do this by making ∂S = 0, k = 1, 2 ∂xk Now  ∂ri ∂S = 2ri , k = 1, 2 ∂xk ∂xk 5

i=1

(2.31)

2.12 OVER- AND UNDER-DETERMINED SYSTEMS

111

and thus 5  i=1

ri

∂ri = 0, k = 1, 2 ∂xk

(2.32)

It can be shown that minimizing the sum of the squares of the residuals using (2.32) gives an identical solution to that given by the pseudo-inverse method for solving the equation system. When solving a set of over-determined equations, determining the pseudo-inverse of the system matrix is only part of the process and normally we do not require this interim result. The M ATLAB operator \ solves over-determined systems automatically. Thus, the operator may be used to solve any linear equation system. In the following example, we compare the results obtained using the operator \ and using the pseudo-inverse for solving (2.30). The M ATLAB script e4s206.m follows: % e4s206.m A = [1 1;2.05 -1;3.06 1;-1.02 2;4.08 -1]; b = [1.98;0.95;3.98;0.92;2.90]; x = pinv(A)*b norm_pinv = norm(A*x-b) x = A\b norm_op = norm(A*x-b)

Running script e4s206.m gives the following numeric output: x = 0.9631 0.9885 norm_pinv = 0.1064 x = 0.9631 0.9885 norm_op = 0.1064

Here, both the M ATLAB operator \ and the function pinv have provided the same “best fit” solution for the inconsistent set of equations. Fig. 2.6 shows the region where these equations intersect in greater detail than Fig. 2.5. The symbol “+” indicates the M ATLAB solution which can be seen to lie in this region. The norm of Ax − b is the square root of the sum of the squares of the residuals and provides a measure of how well the equations are satisfied.

112

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

FIGURE 2.6 Plot of inconsistent equation system (2.30) showing the region of intersection of the equations, where + indicates “best” solution.

The M ATLAB operator \ does not solve an over-determined system by using the pseudo-inverse, as given in (2.26). Instead, it solves (2.24) directly by QR decomposition. QR decomposition can be applied to both square and rectangular matrices providing the number of rows is greater than the number of columns. For example, applying the M ATLAB function qr to solve the over-determined system (2.30) we have >> A = [1 1;2.05 -1;3.06 1;-1.02 2;4.08 -1]; >> b = [1.98 0.95 3.98 0.92 2.90].’; >> [Q R] = qr(A) Q = -0.1761

0.4123

-0.7157

-0.2339

-0.4818

-0.3610

-0.2702

0.0998

0.6751

-0.5753

-0.5388

0.5083

0.5991

-0.2780

-0.1230

0.1796

0.6839

-0.0615

0.6363

0.3021

-0.7184

-0.1756

-0.3394

0.0857

0.5749

R = -5.6792

0.7237

0

2.7343

0

0

0

0

0

0

In the equation Ax = b we have replaced A by QR so that QRx = b. Let Rx = y. Thus, we have y = Q−1 b = QT b since Q is orthogonal. Once y is determined we can efficiently determine x by back substitution since R is upper triangular. Thus, continuing the previous example,

2.12 OVER- AND UNDER-DETERMINED SYSTEMS

113

>> y = Q.’*b y = -4.7542 2.7029 0.0212 -0.0942 -0.0446

Using the second row of R and the second row of y we can determine x2 . From the first row of R and the first row of y we can determine x1 since x2 is known. Thus, −5.6792x1 + 0.7237x2 = −4.7542 2.7343x2 = 2.7029 gives x1 = 0.9631 and x2 = 0.9885, as before. The M ATLAB operator \ implements this sequence of operations. We now consider a case where the coefficient matrix of the over-determined system is rank deficient. The following example is rank deficient and represents a system of parallel lines. x1 + 2x2 = 1.00 x1 + 2x2 = 1.03 x1 + 2x2 = 0.97 x1 + 2x2 = 1.01 In M ATLAB this becomes >> A = [1 2;1 2;1 2;1 2] A = 1 1 1 1

2 2 2 2

>> b = [1 1.03 0.97 1.01].’ b = 1.0000 1.0300 0.9700 1.0100 >> y = A\b Warning: Rank deficient, rank = 1,

tol =

3.552714e-015.

114

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

y = 0 0.5012 >> norm(y) ans = 0.5012

The user is warned that this system is rank deficient. We have solved the system using the \ operator and now solve it using the pinv function as follows: >> x = pinv(A)*b x = 0.2005 0.4010 >> norm(x) ans = 0.4483

We see that when the pinv function and \ operator are applied to rank deficient systems the pinv function gives the solution with the smallest Euclidean norm; see Appendix A, Section A.10. Clearly there is no unique solution to this system since it represents a set of parallel lines. We now turn to the problem of under-determined systems. Here there is insufficient information to determine a unique solution. For example, consider the equation system x1 + 2x2 + 3x3 + 4x4 = 1 −5x1 + 3x2 + 2x3 + 7x4 = 2 Expressing these equations in M ATLAB we have >> A = [1 2 3 4;-5 3 2 7]; >> b = [1 2].’; >> x1 = A\b x1 = -0.0370 0 0 0.2593 >> x2 = pinv(A)*b

2.12 OVER- AND UNDER-DETERMINED SYSTEMS

115

x2 = -0.0780 0.0787 0.0729 0.1755

We calculate the norms: >> norm(x1) ans = 0.2619 >> norm(x2) ans = 0.2199

The first solution vector, x1, is a solution which satisfies the system; the second solution vector, x2, satisfies the system of equations but also gives the solution with the minimum norm. The definition of the Euclidean or 2-norm of a vector is the square root of the sum of the squares of the elements of the vector; see Appendix A, Section A.10. The shortest distance between a point in space and the origin is given by Pythagoras’s theorem as the square root of the sum of squares of the coordinates of the point. Thus, the Euclidean norm of a vector which is a point on a line, surface, or hypersurface may be interpreted geometrically as the distance between this point and the origin. The vector with the minimum norm must be the point on the line, surface, or hypersurface which is closest to the origin. It can be shown that the line joining this vector to the origin must be perpendicular to the line, surface, or hypersurface. Giving the minimum norm solution has the advantage that, whereas there are an infinite number of solutions to an under-determined problem, there is only one minimum norm solution. This provides a standard result. To complete the discussion of over- and under-determined systems we consider the use of the lsqnonneg function which solves the least squares problem with a non-negative solution. This solves the problem of finding a solution x to the linear equation system Ax = b, subject to x ≥ 0 where A and b must be real. This is equivalent to the problem of finding the vector x that minimizes norm(Ax − b) subject to x ≥ 0. We can call the M ATLAB function lsqnonneg for a specific problem using the statement x = lsqnonneg(A,b)

116

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

where A corresponds to A in our definition and b to b. The solution is given by x. Consider the example which follows. Solve ⎤ ⎡ x1

 ⎢ x2 ⎥

 ⎥ 7 1 1 1 1 0 ⎢ ⎢ x3 ⎥ = ⎥ 12 1 2 3 0 1 ⎢ ⎣ x4 ⎦ x5 subject to xi ≥ 0, where i = 1, 2, ..., 5. In M ATLAB this becomes >> A = [1 1 1 1 0;1 2 3 0 1]; >> b = [7 12].’;

Solving this system gives >> x = lsqnonneg(A,b) x = 0 0 4 3 0

We can also solve this using \ but this will not ensure non-negative values for x. >> x2 = A\b x2 = 0 0 4.0000 3.0000 0

In this case, we do obtain a non-negative solution but this is fortuitous and cannot be guaranteed. The following example illustrates how the lsqnonneg function forces a non-negative solution which best satisfies the equations: ⎤ ⎤ ⎡ ⎡ 2.5 3.0501 4.8913

 ⎢ 2.5 ⎥ ⎢ 3.2311 −3.2379 ⎥ x1 ⎥ ⎥ ⎢ (2.33) =⎢ ⎣ 0.5 ⎦ ⎣ 1.6068 7.4565 ⎦ x2 2.5 2.4860 −0.9815 >> A = [3.0501 4.8913;3.2311 -3.2379; 1.6068 7.4565;2.4860 -0.9815]; >> b = [2.5 2.5 0.5 2.5].’;

We can compute the solution using \ or lsqnonneg function thus:

2.13 ITERATIVE METHODS

117

>> x1 = A\b x1 = 0.8307 -0.0684 >> x2 = lsqnonneg(A,b) x2 = 0.7971 0 >> norm(A*x1-b) ans = 0.7040 >> norm(A*x2-b) ans = 0.9428

Thus, the best fit is given by using the operator \, but if we require all components of the solution to be non-negative, then we must use the lsqnonneg function.

2.13 ITERATIVE METHODS Except in special circumstances, it is unlikely that any function or script developed by the user will outperform a function or operator that is an integral part of M ATLAB. We cannot expect to develop a function that will determine the solution of Ax = b more efficiently than by using the M ATLAB operation A\b. However, we describe iterative methods here for the sake of completeness. These methods are attractive if A is diagonally dominant. Iterative methods of solution are developed as follows. We begin with a system of linear equations a11 x1 + a12 x2 + . . . a21 x1 + a22 x2 + . . . .. .. . . an1 x1 + an2 x2 + . . . These can be rearranged to give

+a1n xn +a2n xn .. . +ann xn

= b1 = b2 .. . = bn

118

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

x1 = (b1 − a12 x2 − a13 x3 − . . . − a1n xn ) /a11 x2 = (b2 − a21 x1 − a23 x3 − . . . − a2n xn ) /a22 .. .. .. .. . . . .   xn = bn − an1 x1 − an2 x2 − . . . − an,n−1 xn−1 /ann

(2.34)

If we assume initial values for xi , where i = 1, ..., n, and substitute these values into the right side of the preceding equations, we may determine new values for the xi from (2.34). The iterative process is continued by substituting these values of xi into the right side of the equations, etc. There are several variants of the process. For example, we can use old values of xi in the right side of the equations to determine all the new values of xi in the left side of the equation. This is called Jacobi or simultaneous iteration. Alternatively, we may use a new value of xi in the right side of the equation as soon as it is determined, to obtain the other values of xi in the right side. For example, once a new value of x1 is determined from the first equation of (2.34), it is used in the second equation, together with the old x3 , ..., xn to determine x2 . This is called Gauss–Seidel or cyclic iteration. The conditions for convergence for this type of iteration are |aii | >>

n 

  aij  for i = 1, 2, . . . , n

j =1, j =i

Thus, these iterative methods are only guaranteed to work when the coefficient matrix is diagonally dominant. An iterative method based on conjugate gradients for the solution of systems of linear equations is discussed in Chapter 9.

2.14 SPARSE MATRICES Sparse matrices arise in many problems of science and engineering – for example, in linear programming and the analysis of structures. Indeed, most large matrices that arise in the analysis of physical systems are sparse and the recognition of this fact and its efficient exploitation makes the solution of linear systems with millions of coefficients feasible. The aim of this section is to give a brief description of the extensive sparse matrix facilities available in M ATLAB and to give practical illustrations of their value through examples. For background information on how M ATLAB implements the concept of sparsity, see Gilbert et al. (1992). It is difficult to give a simple quantitative answer to the question: when is a matrix sparse? A matrix is sparse if it contains a high proportion of zero elements. However, this is significant only if the sparsity is of such an extent that we can utilize this feature to reduce the computation time and storage facilities required for operations used on such matrices. One major way in which time can be saved in dealing with sparse matrices is to avoid unnecessary operations on zero elements. M ATLAB does not automatically treat a matrix as sparse and the sparsity features of M ATLAB are not introduced until invoked. Thus, the user determines whether a matrix is in the sparse class or the full class. If the user considers a matrix to be sparse and wants to use this fact to advantage, the matrix must first be converted to sparse form. This is achieved by using the function sparse. Thus B = sparse(A) converts the matrix A to sparse form and assigns it to B and M ATLAB operations

2.14 SPARSE MATRICES

119

on B will take account of this sparsity. If we wish to return this matrix to full form, we simply use C = full(B). However, the sparse function can also be used directly to generate sparse matrices. It is important to note that binary operators *, +, -, /, and \ produce sparse results if both operands are sparse. Thus, the property of sparsity may survive a long sequence of matrix operations. In addition such functions as chol(A) and lu(A) produce sparse results if the matrix A is sparse. However, in mixed cases, where one operand is sparse and the other is full, the result is generally a full matrix. Thus, the property of sparsity may be inadvertently lost. Notice in particular that eye(n) is not in the sparse class of matrices in M ATLAB but a sparse identity matrix can be created using speye(n). Thus, the latter should be used in manipulations with sparse matrices. We will now introduce some of the key M ATLAB functions for dealing with sparse matrices, describe their use and, where appropriate, give examples of their application. The simplest M ATLAB function which helps in dealing with sparsity is the function nnz(a) which provides the number of non-zero elements in a given matrix a, regardless of whether it is sparse or full. A function which enables us to examine whether a given matrix has been defined or has been propagated as sparse is the function issparse(a) which returns the value 1 if the matrix a is sparse or 0 if it is not sparse. The function spy(a) allows the user to view the structure of a given matrix a by displaying symbolically only its non-zero elements; see Fig. 2.7 later in the chapter for examples. Before we can illustrate the action of these and other functions, it is useful to generate some sparse matrices. This is easily done using a different form of the sparse function. This time the function is supplied with the location of the non-zero entries in the matrix, the value of these entries, the size of the sparse matrix and the space allocated for the non-zero entries. This function call takes the form sparse(i, j, nzvals, m, n, nzmax). This generates an m × n matrix and allocates the non-zero values in the vector nzvals to the positions in the matrix given by the vectors i and j. The row position is given by i and the column position by j. Space is allocated for nzmax non-zeros. Since all but one parameter is optional, there are many forms of this function. We cannot give examples of all these forms but the following cases illustrate its use. >> colpos = [1 2 1 2 5 3 4 3 4 5]; >> rowpos = [1 1 2 2 2 4 4 5 5 5]; >> value = [12 -4 7 3 -8 -13 11 2 7 -4]; >> A = sparse(rowpos,colpos,value,5,5)

These statements give the following output: A = (1,1)

12

(2,1)

7

(1,2)

-4

(2,2)

3

(4,3)

-13

(5,3)

2

(4,4)

11

(5,4)

7

(2,5)

-8

(5,5)

-4

120

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

We see that a 5 × 5 sparse matrix with 10 non-zero elements has been generated with the required coefficient values in the required positions. This sparse matrix can be converted to a full matrix as follows: >> B = full(A) B = 12 7 0 0 0

-4 3 0 0 0

0 0 0 -13 2

0 0 0 11 7

0 -8 0 0 -4

This is the equivalent full matrix. Now the following statements test to see if the matrices A and B are in the sparse class and give the number of non-zeros they contain. >> [issparse(A) issparse(B) nnz(A) nnz(B)] ans = 1

0

10

10

Clearly, these functions give the expected results. Since A is a member of the class of sparse matrices, the value of issparse(A) is 1. However, although B looks sparse, it is not stored as a sparse matrix and hence is not in the class sparse within the M ATLAB environment. The next example shows how to generate a large 10,000 × 10,000 sparse matrix and compares the time required to solve a linear system of equations involving this sparse matrix with the time required for the equivalent full matrix. This is provided in script e4s207.m. % e4s207.m Generates a sparse triple diagonal matrix n = 10000; rowpos = 2:n; colpos = 1:n-1; values = 2*ones(1,n-1); Offdiag = sparse(rowpos,colpos,values,n,n); A = sparse(1:n,1:n,4*ones(1,n),n,n); A = A+Offdiag+Offdiag.’; %generate full matrix B = full(A); %generate arbitrary right side for system of equations rhs = [1:n].’; tic, x = A\rhs; f1 = toc; tic, x = B\rhs; f2 = toc; fprintf(’Time to solve sparse matrix = %8.5f\n’,f1); fprintf(’Time to solve full matrix = %8.5f\n’,f2);

Script e4s207.m provides the results Time to solve sparse matrix = Time to solve full matrix =

0.00078 9.66697

2.14 SPARSE MATRICES

121

In this example there is a massive reduction in the time taken to solve the system when using the sparse class of matrix. We now perform a similar exercise, this time to determine the lu decomposition of a 10,000 × 10,000 matrix, using script e4s208: % e4s208.m n = 10000; offdiag = sparse(2:n,1:n-1,2*ones(1,n-1),n,n); A = sparse(1:n,1:n,4*ones(1,n),n,n); A = A+offdiag+offdiag’; %generate full matrix B = full(A); %generate arbitrary right side for system of equations rhs = [1:n]’; tic, lu1 = lu(A); f1 = toc; tic, lu2 = lu(B); f2 = toc; fprintf(’Time for sparse LU = %8.4f\n’,f1); fprintf(’Time for

full

LU = %8.4f\n’,f2);

The time taken to solve the systems using script e4s208 are Time for sparse LU = Time for

full

LU =

0.0051 14.9198

Again this provides a considerable reduction in the time taken. An alternative way to generate sparse matrices is to use the functions sprandn and sprandsym. These provide random sparse matrices and random sparse symmetric matrices respectively. The call A = sprandn(m,n,d)

produces an m × n random matrix with normally distributed non-zero entries of density d. The density is the proportion of the non-zero entries to the total number of entries in the matrix. Thus d must be in the range 0 to 1. To produce a symmetric random matrix with normally distributed non-zero entries of density d, we use A = sprandsys(n,d)

Examples of calls of these functions are given by >> A = sprandn(5,5,0.25) A = (2,1)

-0.4326

(3,3)

-1.6656

(5,3)

-1.1465

(4,4)

0.1253

(5,4)

1.1909

(4,5)

0.2877

122

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

>> B = full(A) B = 0 -0.4326 0 0 0

0 0 0 0 0

0 0 -1.6656 0 -1.1465

0 0 0 0.1253 1.1909

0 0 0 0.2877 0

0 0 0 0 -0.0376

0 0 0.1746 -0.0376 1.1892

>> As = sprandsym(5,0.25) As = (3,1) (1,3) (5,3) (5,4) (3,5) (4,5) (5,5)

0.3273 0.3273 0.1746 -0.0376 0.1746 -0.0376 1.1892

>> Bs = full(As) Bs = 0 0 0.3273 0 0

0 0 0 0 0

0.3273 0 0 0 0.1746

An alternative call for sprandsym is given by A = sprandsym(n,density,r)

If r is a scalar, then this produces a random sparse symmetric matrix with a condition number equal to 1/r. Remarkably, if r is a vector of length n, a random sparse matrix with eigenvalues equal to the elements of r is produced. Eigenvalues are discussed in Section 2.15. A positive definite matrix has all its eigenvalues positive and consequently such a matrix can be generated by choosing each of the n elements of r to be positive. An example of this form of call is >> Apd = sprandsym(6,0.4,[1 2.5 6 9 2 4.3]) Apd = (1,1) (2,1) (4,1)

1.0058 -0.0294 -0.0879

2.14 SPARSE MATRICES

(1,2) (2,2) (4,2) (3,3) (5,3) (1,4) (2,4) (4,4) (3,5) (5,5) (6,6)

123

-0.0294 8.3477 -1.9540 5.4937 -1.3300 -0.0879 -1.9540 3.1465 -1.3300 2.5063 4.3000

>> Bpd = full(Apd) Bpd = 1.0058 -0.0294 0 -0.0879 0 0

-0.0294 8.3477 0 -1.9540 0 0

0 0 5.4937 0 -1.3300 0

-0.0879 -1.9540 0 3.1465 0 0

0 0 -1.3300 0 2.5063 0

0 0 0 0 0 4.3000

This provides an important method for generating test matrices with required properties since, by providing a list of eigenvalues with a range of values, we can produce positive definite matrices that are very badly conditioned. We now examine further the value of using sparsity. The reasons for the very high level of improvement in computing efficiency when using the \ operator, illustrated in the example at the beginning of this section, are complex. The process includes a special preordering of the columns of the matrix. This special preordering, called minimum degree ordering, is used in the case of the \ operator. This preordering takes different forms depending on whether the matrix is symmetric or non-symmetric. The aim of any preordering is to reduce the amount of fill-in from any subsequent matrix operations. Fill-in is the introduction of additional non-zero elements. We can examine this preordering process using the spy function and the function symamd which implements symmetric minimum degree ordering in M ATLAB. The function is automatically applied when working on matrices which belong to the class of sparse matrices for the standard functions and operators of M ATLAB. However, if we are required to use this preordering in non-standard applications, then we may use the symmmd function. The following examples illustrate the use of this function. We first consider the simple process of multiplication applied to a full and a sparse matrix. The sparse multiplication uses the minimum degree ordering. The following script e4s209.m generates a sparse matrix, obtains a minimum degree ordering for it and then examines the result of multiplying the matrix by itself transposed. This is compared with the same operations carried out on the full matrix, and the time required for each operation is compared. % e4s209.m % generate a sparse matrix

124

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

n = 5000; offdiag = sparse(2:n,1:n-1,2*ones(1,n-1),n,n); offdiag2 = sparse(4:n,1:n-3,3*ones(1,n-3),n,n); offdiag3 = sparse(n-5:n,1:6,7*ones(1,6),n,n); A = sparse(1:n,1:n,4*ones(1,n),n,n); A = A+offdiag+offdiag’+offdiag2+offdiag2’+offdiag3+offdiag3’; A = A*A.’; % generate full matrix B = full(A); m_order = symamd(A); tic spmult = A(m_order,m_order)*A(m_order,m_order).’; flsp = toc; tic, fulmult = B*B.’; flful = toc; fprintf(’Time for sparse mult = %6.4f\n’,flsp) fprintf(’Time for full mult = %6.4f\n’,flful)

Running this script e4s209.m results in the following output: Time for sparse mult = 0.0073 Time for full mult = 3.2559

We now perform a similar experiment to that of e4s209.m, but for a more complex numerical process than multiplication. In the script that follows we examine LU decomposition. We consider the result of using a minimum degree ordering on the LU decomposition process by comparing the performance of the lu function with and without preordering. This is illustrated in script e4s210.m. % e4s210.m % generate a sparse matrix n = 100; offdiag = sparse(2:n,1:n-1,2*ones(1,n-1),n,n); offdiag2 = sparse(4:n,1:n-3,3*ones(1,n-3),n,n); offdiag3 = sparse(n-5:n,1:6,7*ones(1,6),n,n); A = sparse(1:n,1:n,4*ones(1,n),n,n); A = A+offdiag+offdiag’+offdiag2+offdiag2’+offdiag3+offdiag3’; A = A*A.’; A1 = flipud(A); A = A+A1; n1 = nnz(A) B = full(A); %generate full matrix m_order = symamd(A); lud = lu(A(m_order,m_order)); n2 = nnz(lud) fullu = lu(B); n3 = nnz(fullu) subplot(2,2,1), spy(A,’k’);

2.14 SPARSE MATRICES

125

FIGURE 2.7 Effect of minimum degree ordering on LU decomposition. The spy function shows the matrix, the ordered matrix, and LU decomposition with and without preordering.

title(’Original matrix’) subplot(2,2,2), spy(A(m_order,m_order),’k’) title(’Ordered Matrix’) subplot(2,2,3), spy(fullu,’k’) title(’LU decomposition,unordered matrix’) subplot(2,2,4), spy(lud,’k’) title(’LU decomposition, ordered matrix’)

Running script e4s210.m gives n1 = 2096 n2 = 1307 n3 = 4723

As expected, by using a sparse operation we achieve a reduction in the time taken to determine the LU decomposition. Fig. 2.7 shows the original matrix with 2096 non-zero elements, the reordered matrix (which has the same number of non-zeros), and the LU decomposition structure both with and without

126

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

minimum degree ordering. Notice that the number of non-zeros in the LU matrices with preordering is 1307 and without is 4465. Thus, there is a large increase in the number of non-zero elements in the LU matrices without preordering. In contrast, LU decomposition of the preordered matrix has produced fewer non-zeros than the original matrix. The reduction of fill-in is an important feature of sparse numerical processes and may ultimately lead to great savings in computational effort. Note that if the size of the matrices are increased from 100 × 100 to 3000 × 3000 then the output from e4s210 is n1 = 65896 n2 = 34657 n3 = 510266 Time for sparse lu = 0.0131 Time for full lu = 0.8313

Here we obtain a more substantial reduction by using sparse operations. The M ATLAB function symamd provides a minimum degree ordering for symmetric matrices. For non-symmetric matrices M ATLAB provides the function colmmd which gives the column minimum degree ordering for non-symmetric matrices. An alternative ordering, which is used to reduce bandwidth, is the reverse Cuthill–MacGee ordering. This is implemented in M ATLAB by the function symrcm. The execution of the statement p = symrcm(A) provides the permutation vector p to produce the required preordering and A(p,p) is the reordered matrix. We have shown that in general taking account of sparsity will provide savings in floating-point operations. However, these savings fall off as the matrices on which we are operating become less sparse as script e4s211.m illustrates. % e4s211.m n = 1000; b = 1:n; disp(’ density time_sparse time_full’); for density = 0.004:0.003:0.039 A = sprandsym(n,density)+0.1*speye(n); density = density+1/n; tic, x = A\b’; f1 = toc; B = full(A); tic, y = B\b’; f2 = toc; fprintf(’%10.4f %12.4f %12.4f\n’,density,f1,f2); end

In script e4s211 a diagonal of elements has been added to the randomly generated sparse matrix. This is done to ensure that each row of the matrix contains a non-zero element; otherwise the matrix may be singular. Adding this diagonal modifies the density. If the original n × n matrix has a density of d, then, assuming that this matrix has only zeros on the diagonal, the modified density is d + 1/n.

2.14 SPARSE MATRICES

density 0.0050 0.0080 0.0110 0.0140 0.0170 0.0200 0.0230 0.0260 0.0290 0.0320 0.0350 0.0380

time_sparse 0.0204 0.0329 0.0508 0.0744 0.0892 0.1064 0.1179 0.1317 0.1444 0.1516 0.1789 0.1627

127

time_full 0.1907 0.1318 0.1332 0.1399 0.1351 0.1372 0.1348 0.1381 0.1372 0.1369 0.1404 0.1450

This output shows that the advantage of using a sparse class of matrix diminishes as the density increases. Another application where sparsity is important is in solving the least squares problem. This problem is known to be ill-conditioned and hence any saving in computational effort is particularly beneficial. This is directly implemented by using A\b where A is non-square and sparse. To illustrate the use of the \ operator with sparse matrices and compare its performance when no account is taken of sparsity we use script e4s212.m: % e4s212.m % generate a sparse triple diagonal matrix n = 2000; rowpos = 2:n; colpos = 1:n-1; values = ones(1,n-1); offdiag = sparse(rowpos,colpos,values,n,n); A = sparse(1:n,1:n,4*ones(1,n),n,n); A = A+offdiag+offdiag’; %Now generate a sparse least squares system Als = A(:,1:n/2); %generate full matrix Cfl = full(Als); rhs = 1:n; tic, x = Als\rhs’; f1 = toc; tic, x = Cfl\rhs’; f2 = toc; fprintf(’Time for sparse least squares solve = %8.4f\n’,f1) fprintf(’Time for full least squares solve = %8.4f\n’,f2)

Script e4s212.m provides the following results Time for sparse least squares solve = Time for full least squares solve =

Again we see the advantage of using sparsity.

0.0027 0.9444

128

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

FIGURE 2.8 Mass-spring system with three degrees of freedom.

We have not covered all aspects of sparsity or described all the related functions. However, we hope this section provides a helpful introduction to this difficult but important and valuable feature of M ATLAB.

2.15 THE EIGENVALUE PROBLEM Eigenvalue problems arise in many branches of science and engineering. For example, the vibration characteristics of structures are determined from the solution of an algebraic eigenvalue problem. Here, we consider a particular example of a system of masses and springs shown in Fig. 2.8. The equations of motion for this system are m1 q¨1 + (k1 + k2 + k4 ) q1 − k2 q2 − k4 q3 = 0 m2 q¨2 − k2 q1 + (k2 + k3 ) q2 − k3 q3 = 0 m3 q¨3 − k4 q1 − k3 q2 + (k3 + k4 ) q3 = 0

(2.35)

If we assume a where m1 , m2 , and m3 are the system masses and k1 , ..., k4 are the spring stiffnesses. √ harmonic solution for each coordinate, then qi (t) = ui exp(j ωt) where j = −1, for i = 1, 2, and 3. Hence d 2 qi /dt 2 = −ω2 ui exp(j ωt). Substituting in (2.35) and canceling the common factor exp(j ωt) gives −ω2 m1 u1 + (k1 + k2 + k4 ) u1 − k2 u2 − k4 u3 = 0 −ω2 m2 u2 − k2 u1 + (k2 + k3 ) u2 − k3 u3 = 0 −ω2 m3 u3

(2.36)

− k4 u1 − k3 u2 + (k3 + k4 ) u3 = 0

If m1 = 10 kg, m2 = 20 kg, m3 = 30 kg, k1 = 10 kN/m, k2 = 20 kN/m, k3 = 25 kN/m, and k4 = 15 kN/m, then (2.36) becomes −ω2 10u1 + 45,000u1 − 20,000u2 − 15,000u3 = 0 −ω2 20u1 − 20,000u1 + 45,000u2 − 25,000u3 = 0 −ω2 30u1 − 15,000u1 − 25,000u2 + 40,000u3 = 0

2.15 THE EIGENVALUE PROBLEM

129

This can be expressed in matrix notation as

where

−ω2 Mu + Ku = 0

(2.37)

⎤ ⎡ ⎤ 10 0 0 45 −20 −15 M = ⎣ 0 20 0 ⎦ kg and K = ⎣ −20 45 −25 ⎦ kN/m 0 0 30 −15 −25 40

(2.38)

⎡

Eq. (2.37) can be rearranged in a variety of ways. For example, it can be written Mu = λKu where λ =

1 ω2

(2.39)

This is an algebraic eigenvalue problem and solving it determines values for the eigenvectors u and the eigenvalues, λ. M ATLAB provides a function eig to solve the eigenvalue problem and to illustrate its use we apply it to the solution of (2.37). >> M = [10 0 0;0 20 0;0 0 30]; >> K = 1000*[45 -20 -15;-20 45 -25;-15 -25 40]; >> lambda = eig(M,K).’ lambda = 0.0002

0.0004

0.0073

>> omega = sqrt(1./lambda) omega = 72.2165

52.2551

11.7268

This result tells us that the system of Fig. 2.8 vibrates with natural frequencies 11.72, 52.25, and 72.21 rad/s. In this example we have chosen not to determine u, the eigenvectors. We will discuss further the use of the function eig in Section 2.17. Having provided an example of an eigenvalue problem, we consider the standard form of this problem thus: Ax = λx

(2.40)

This equation is an algebraic eigenvalue problem where A is a given n × n matrix of coefficients, x is an unknown column vector of n elements and λ is an unknown scalar. Eq. (2.40) can alternatively be written as (A − λI)x = 0

(2.41)

Our aim is to discover the values of x, called the characteristic vectors or eigenvectors, and the corresponding values of λ, called the characteristic values or eigenvalues. The values of λ that satisfy (2.41) are given by the roots of the equation |A − λI| = 0

(2.42)

130

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

These values of λ are such that (A − λI) is singular. Since (2.41) is a homogeneous equation system, non-trivial solutions exist for these values of λ. Evaluation of the determinant (2.42) leads to an nth degree polynomial in λ which is called the characteristic polynomial. This characteristic polynomial has n roots, some of which may be repeated, giving the n values of λ. In M ATLAB we can create the coefficients of the characteristic polynomial using the function poly, and the roots of the resulting characteristic polynomial can be found using the function roots. For example, ⎡ ⎤ 1 2 3 A=⎣ 4 5 −6 ⎦ 7 −8 9 then we have >> A = [1 2 3;4 5 -6;7 -8 9]; >> p = poly(A) p = 1.0000

-15.0000

-18.0000

360.0000

Hence the characteristic equation is λ3 − 15λ2 − 18λ + 360 = 0. To find the roots of this equation we use the statement >> roots(p).’ ans = 14.5343

-4.7494

5.2152

We can verify this result using the function eig, thus: >> eig(A).’ ans = -4.7494

5.2152

14.5343

Having obtained the eigenvalues we can substitute back into (2.41) to obtain linear equations for the characteristic vectors thus: (A − λi I)x = 0 for i = 1, 2, ..., n

(2.43)

This homogeneous system provides n non-trivial solutions for x. However, the use of (2.42) and (2.43) is not a practical means of solving eigenvalue problems. If A is a symmetric real matrix or a Hermitian matrix, (such that A = AH ), then its eigenvalues are real, but not necessarily positive, and the corresponding eigenvectors are real if A is real and symmetric, and complex if A is Hermitian. For such a system, if λi , xi and λj , xj satisfy the eigenvalue problem and λi and λj are distinct, then xH i xj = 0 i = j

(2.44)

2.15 THE EIGENVALUE PROBLEM

131

and xH i Axj = 0 i = j

(2.45)

Eqs. (2.44) and (2.45) are called the orthogonality relationships. If A is real and symmetric, then in (2.44) and (2.45) we take the transpose of the vectors denoted by (T ), rather than the complex conjugate H transpose denoted by (H ). If i = j , then in general xH i xi and xi Axi are non-zero. The vector xi includes an arbitrary scalar multiplier because the vector multiplies both sides of (2.40). Hence the product xH i xi must be arbitrary. However, if the arbitrary scalar multiplier is adjusted so that xH x = 1 then the vectors i i are said to be normalized. In fact, the M ATLAB function eig normalizes the eigenvectors, so that xH i xi = 1 for i = 1, 2, ..., n

(2.46)

When the vectors are normalized so that (2.46) is valid, then xH i Axi = λi for i = 1, 2, ..., n

(2.47)

Note that the eigenvectors are orthogonal, that is (2.44) and (2.45) are valid, then they are orthogonal irrespective of whether or not the eigenvectors are normalized according to (2.46). Sometimes the eigenvalues are not distinct and the eigenvectors associated with these equal or repeated eigenvalues are not necessarily orthogonal. If λi = λj and the other eigenvalues, λk , are distinct, then ⎫ ⎬ x = 0 xH i k k = 1, 2, ..., n, k = i, k = j (2.48) x H xk = 0 ⎭ j

For consistency, we can choose to make xH i xj = 0. When λi = λj , the eigenvectors xi and xj are not unique and a linear combination of them, (αxi + γ xj ) where α and γ are arbitrary constants, also satisfies the eigenvalue problem. We can combine all the eigensolutions into a single matrix equation by placing all the normalized eigenvectors into a single matrix, X and the eigenvalues in a diagonal matrix, D thus:  X = x 1 x2 . . . xn (2.49) ⎡ ⎢ ⎢ D=⎢ ⎣

λ1 0 .. . 0

0 ... λ2 . . . .. .. . . 0 ...

0 0 .. . λn

⎤ ⎥ ⎥ ⎥ ⎦

(2.50)

If A is an Hermitian matrix, then X is a unitary matrix, that is XH = X−1 . Its determinant is either +1 or −1 and its eigenvalues are complex but lie on a unit circle in the complex plane, that is their amplitudes all equal one but their phases differ. If A is an n × n real symmetric matrix, then X is an orthogonal matrix, that is XT = X−1 . The determinant of X is +1 and its eigenvalues are complex conjugate pairs if n is even, or complex conjugate pairs plus one real eigenvalue if n is odd. Again, all the eigenvalues have an absolute value of one and thus lie on a unit circle in the complex plane.

132

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

The n solutions (2.49) and (2.50) satisfy the eigenvalue problem (2.40). Assembling all the eigensolutions into one equation, we have: AX = XD

(2.51)

Using similar reasoning to that above, and substituting (2.49) and (2.50) in (2.46) gives XH X = I

(2.52)

Pre-multiplying (2.51) by XH , and noting the relationship given by (2.52) we have: XH AX = D

(2.53)

Example 2.7. Determine the eigenvalues of the Hermitian matrix

A=

1 + 0ı 2 − 2ı

2 + 2ı 3 + 0ı



Using M ATLAB function eig we can determine both the eigenvalues and eigenvectors of this matrix as follows: >> A = [1 2;2 3]+j*[0 2;-2 0]; >> [X,D] = eig(A) X = 0.5774 + 0.5774i

0.4082 + 0.4082i

-0.5774 + 0.0000i

0.8165 + 0.0000i

D = -1.0000

0

0

5.0000

where the diagonal elements of D are the eigenvalues and the corresponding eigenvector are given by X. We can verify, for this problem, the relationships given by (2.52) and (2.53) as follows: >> X’*A*X ans = -1.0000 + 0.0000i

-0.0000 - 0.0000i

-0.0000 + 0.0000i

5.0000 + 0.0000i

>> X’*X ans = 1.0000 + 0.0000i

-0.0000 - 0.0000i

-0.0000 + 0.0000i

1.0000 + 0.0000i

2.15 THE EIGENVALUE PROBLEM

133

In summary, if A is real, both its eigenvalues and eigenvectors are real, but the eigenvalues are not necessarily positive. If A is a Hermitian matrix, its eigenvalues are real, but not necessarily positive, but its eigenvectors are generally complex. Consider now the case when A is a general complex square matrix (that is, not Hermitian) with a set of eigenvalues xi and corresponding eigenvalues λi . A pair of related eigenvalue problems exist as follows: Ax = λx

(2.54)

AT y = βy

(2.55)

and

The second eigenvalue problem, (2.55), can be transposed to give yT A = βyT

(2.56)

The vectors x and y are called the right and left vectors of A, respectively. The equations |A − λI| = 0 and |AT − βI| = 0 must have the same solutions for λ and β because the determinant of a matrix and the determinant of its transpose are equal. Thus, the eigenvalues of A and AT are identical but the eigenvectors x and y will, in general, differ from each other. The eigenvalues and eigenvectors of a non-symmetric real matrix are either real or pairs of complex conjugates. If λi , xi , yi and λj , xj , yj are solutions that satisfy the eigenvalue problems of (2.54) and (2.55) and λi and λj are distinct, then yTi xj = 0 i = j

(2.57)

yTi Axj = 0 i = j

(2.58)

and

Eqs. (2.57) and (2.58) are called the bi-orthogonal relationships. As with (2.46) and (2.47) if, in these equations, i = j , then in general yTi xi and yTi Axi are not zero. The eigenvectors xi and yi include arbitrary scaling factors and so the product of these vectors will also be arbitrary. However, if the vectors are adjusted so that yTi xi = 1 for i = 1, 2, ..., n

(2.59)

yTi Axi = λi for i = 1, 2, ..., n

(2.60)

then

We cannot, in these circumstances, describe either xi or yi as normalized; the vectors still include an arbitrary scale factor; only their product can be chosen uniquely. The eigensolutions can be combined us: X= Y=

 

x1

x2

...

xn

y1

y2

...

yn



(2.61) (2.62)

134

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

⎡ ⎢ ⎢ D=⎢ ⎣

0 ... 0 λ1 0 λ2 . . . 0 .. .. .. .. . . . . 0 0 . . . λn

⎤ ⎥ ⎥ ⎥ ⎦

(2.63)

The n solutions (2.61) and (2.63) satisfy the eigenvalue problem (2.54). Assembling all the eigensolutions into one equation, we have: AX = XD

(2.64)

Similarly, noting that βi = λi , the n solutions (2.62) and (2.63) satisfy the eigenvalue problem (2.55). Assembling all the eigensolutions into one equation, we have: AT Y = YD or YT A = DYT

(2.65)

Using (2.49) and (2.50) we can write (2.59) in matrix form: YT X = I

(2.66)

Pre-multiplying (2.64) by YT and noting the relationship given by (2.66) we have YT AX = D

(2.67)

In summary, if A is real (but not symmetric) or complex (but not Hermitian) then both its eigenvalues and eigenvectors will be complex. Example 2.8. determine the eigenvalues of the following complex matrix,   1 + 2ı 2 + 3ı A= 4 − 1ı 3 + 1ı >> A = [1 2;4 3]+j*[2 3;-1 1] A = 1.0000 + 2.0000i

2.0000 + 3.0000i

4.0000 - 1.0000i

3.0000 + 1.0000i

>> [X,D] = eig(A) X = 0.7518 + 0.0000i

0.4126 + 0.4473i

-0.6068 + 0.2580i

0.7935 + 0.0000i

D = -1.6435 + 0.2649i

0.0000 + 0.0000i

0.0000 + 0.0000i

5.6435 + 2.7351i

2.15 THE EIGENVALUE PROBLEM

135

>> [Y,B] = eig(A.’) Y = 0.7935 + 0.0000i -0.4126 - 0.4473i

0.6068 - 0.2580i 0.7518 + 0.0000i

B = -1.6435 + 0.2649i 0.0000 + 0.0000i

0.0000 + 0.0000i 5.6435 + 2.7351i

Note that the eigenvalues given by B and D are identical but the right and left eigenvalues given by X and Y are different. We can verify, for this problem, the relationships given by (2.66) and (2.67) as follows: >> p = Y.’*X p = 0.9623 + 0.1650i -0.0000 + 0.0000i

-0.0000 - 0.0000i 0.9623 + 0.1650i

>> X(:,1) = X(:,1)/p(1,1); X(:,2) = X(:,2)/p(2,2); >> Y.’*X ans = 1.0000 - 0.0000i -0.0000 + 0.0000i

-0.0000 - 0.0000i 1.0000 + 0.0000i

>> Y.’*A*X ans = -1.6435 + 0.2649i 0.0000 - 0.0000i

-0.0000 - 0.0000i 5.6435 + 2.7351i

2.15.1 EIGENVALUE DECOMPOSITION Consider the eigensolution of a Hermitian (or real symmetric matrix) given by (2.51). Then, postmultiplying by X−1 gives A = XDX−1 This is not a particularly useful decomposition of A since it involves a matrix inversion. However, from (2.52), we have X−1 = XH and hence A = XDXH

(2.68)

Thus A has been decomposed into the product of two orthogonal or unitary matrices and a diagonal matrix. This is useful decomposition. For example, to compute the square of a Hermitian matrix. We

136

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

have A2 = AA = (XDXH )(XDXH ) and thus A2 = XDIDXH = XD2 XH

(2.69)

since XH X = I. Computing D2 is trivial since D is diagonal and so each element on the diagonal is squared. In fact, this result can be generalized to Ap = XDp XH

(2.70)

Another useful result is obtained by take in inverse of (2.68) to give A−1 = (XDXH )−1 = X−H D−1 X−1

(2.71)

Since X−1 = XH , on inverting X = X−H so (2.71) becomes A−1 = XD−1 XH

(2.72)

This is an efficient method of inverting a matrix since inverting the diagonal matrix D is trivial; it is obtained by taking the reciprocal of each element on the diagonal. Of course, this assumes that the eigenvalue problem has been solved. Example 2.9. Eqs. (2.68) to (2.72) are illustrated in M ATLAB by continuing Example 2.7 >> A = [1 2;2 3]+j*[0 2;-2 0]; >> [X D]= eig(A); >> AA = X*D*X’ AA = 1.0000 + 0.0000i

2.0000 + 2.0000i

2.0000 - 2.0000i

3.0000 + 0.0000i

>> Ainv = X*inv(D)*X’ Ainv = -0.6000 + 0.0000i

0.4000 + 0.4000i

0.4000 - 0.4000i

-0.2000 + 0.0000i

>> inv(A) ans = -0.6000 + 0.0000i

0.4000 + 0.4000i

0.4000 - 0.4000i

-0.2000 + 0.0000i

2.15 THE EIGENVALUE PROBLEM

137

Consider the decomposition of a general, complex matrix. By post-multiplying (2.64) by X−1 we have A = XDX−1

(2.73)

Similarly, pre-multiplying the second expression of (2.65) by Y−T we have A = Y−T DYT

(2.74)

Neither of these are particularly useful transforms since they both require matrix inversion. However, from (2.66), X−1 = YT . Substituting in (2.73) gives A = XDYT

(2.75)

This decomposition avoids the need to invert a matrix, but it requires the solution of two eigenvalue problems for A and AT . Example 2.10. Consider the matrices of Example 2.8 >> A = [1 2;4 3]+j*[2 3;-1 1]; >> [X,D] = eig(A); >> [Y,B] = eig(A.’); >> p = Y.’*X; >> X(:,1) = X(:,1)/p(1,1); X(:,2) = X(:,2)/p(2,2); >> A1 = X*D*inv(X) A1 = 1.0000 + 2.0000i

2.0000 + 3.0000i

4.0000 - 1.0000i

3.0000 + 1.0000i

>> A2 = inv(Y.’)*D*Y.’ A2 = 1.0000 + 2.0000i

2.0000 + 3.0000i

4.0000 - 1.0000i

3.0000 + 1.0000i

>> A3 = X*D*Y.’ A3 = 1.0000 + 2.0000i

2.0000 + 3.0000i

4.0000 - 1.0000i

3.0000 + 1.0000i

Eigenvalue decomposition will fail if the matrix being decomposed is defective. A defective matrix is a square matrix that does not have a full set of linearly independent eigenvectors and consequently can not be diagonalized. A is a normal matrix if AH A = AAH . A normal matrix (which includes Hermitian and real symmetric matrices as special cases) is never defective.

138

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

2.15.2 COMPARING EIGENVALUE AND SINGULAR VALUE DECOMPOSITION Eigenvalue decomposition (EVD) can only decompose square matrices whereas singular value decomposition (SVD) can decompose rectangular or square matrices. There is, however, a relationship between singular values and eigenvalues. Consider a general, n × m matrix A. The SVD of A is A = USVT

(2.76)

From this SVD we can write AAT = (USVT )(VST UT ) = USST UT because VT V = I. Let the EVD of AAT be AAT = XDXT

(2.77)

Comparing (2.76) and (2.77), we see D = SST and U = X. That is U is equal to the eigenvector of AAT and SST is equal to the eigenvalues of AAT . From the SVD of A we can also write AT A = (VSUT )(UST VT ) = VSST VT since UT U = I. Let the EVD of AT A be AT A = YDYT Then V = Y and again D = SST . The fact that the eigenvalues of AT A are the square of the non-zero singular values of A has implications for the condition number A. The condition number of AAT is the condition number of A squared. Thus if A is ill-conditioned, working with AT A will make the ill-conditioning worse. Thus, in difficult cases singular values are computed more accurately than eigenvalues. Example 2.11. Determine the eigenvalues and condition number of AAT and the SVD and condition number of A. Note that A is a 3 × 4 matrix. >> A = [4 -2 0 7;-2 7 -5 3;0 -5 12 -2]; >> [X,D] = eig(A*A’) X = 0.0620

0.9979

-0.0163

0.8321

-0.0607

-0.5512

0.5511

-0.0206

0.8342

20.0386

0

0

0

69.1435

0

0

0

239.8179

D =

2.15 THE EIGENVALUE PROBLEM

139

>> [U,S,V] = svd(A) U = -0.0163 -0.5512 0.8342

0.9979 -0.0607 -0.0206

0.0620 0.8321 0.5511

S = 15.4861 0 0

0 8.3153 0

0 0 4.4764

0 0 0

0.4947 -0.2787 0.0068 0.8231

-0.3164 0.6580 0.5478 0.4084

-0.8067 -0.4719 -0.1423 0.3261

0 69.1435 0

0 0 20.0386

V = 0.0670 -0.5164 0.8244 -0.2219 >> S*S’ ans = 239.8179 0 0 >> cond(A) ans = 3.4595 >> cond(A*A’) ans = 11.9678

We note that S2 is equal to D (but the eigenvalues are listed in a different order) and the condition number of AAT is the square of the condition number of A. There is a special case when A is a square, symmetric, positive definite matrix, illustrated in the following example. Example 2.12. The matrix A is square and symmetric. Determining its eigensolution we have >> A = [4 -2 0;-2 7 -5; 0 -5 12]; >> [X,D] = eig(A)

140

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

X = -0.7055 -0.6312 -0.3224

-0.7022 0.5611 0.4382

0.0957 -0.5355 0.8391

2.2107 0 0

0 5.5980 0

0 0 15.1913

D =

Since all the eigenvalues of the square matrix A are positive, it must be a positive definite matrix. Determining its SVD we have >> [U,S,V] = svd(A) U = -0.0957 0.5355 -0.8391

-0.7022 0.5611 0.4382

0.7055 0.6312 0.3224

15.1913 0 0

0 5.5980 0

0 0 2.2107

-0.0957 0.5355 -0.8391

-0.7022 0.5611 0.4382

0.7055 0.6312 0.3224

S =

V =

This example has illustrated the fact that the eigenvalues and singular values of a square, symmetrical, positive definite matrix are identical. Because the eigenvalue and singular values have been computed in different orders it isn’t immediately obvious that X is equivalent to V. However, if we rearrange X so that it becomes equal to [−x3 x2 −x1 ] then X does equal V. Note that the columns of X have been reordered in line with the order of the singular values, and two columns have been multiplied by −1. Recall that each eigenvector can be multiplied by a constant, in this case −1, and still be a solution to the eigenvalue problem.

2.16 ITERATIVE METHODS FOR SOLVING THE EIGENVALUE PROBLEM The first of two simple iterative procedures described here determines the dominant or largest eigenvalue. The method, which is called the power method or matrix iteration, can be used on both symmetric and non-symmetric matrices. However, for a non-symmetric matrix the user must be alert to

2.16 ITERATIVE METHODS FOR SOLVING THE EIGENVALUE PROBLEM

141

the possibility that there is not a single real dominant eigenvalue value but a complex conjugate pair. Under these conditions simple iteration does not converge. Consider the eigenvalue problem defined by (2.40) and let the vector u0 be an initial trial solution. The vector u0 is an unknown linear combination of all the eigenvectors of the system provided they are linearly independent. Thus: n  αi xi (2.78) u0 = i=1

where αi are unknown coefficients and xi are the unknown eigenvectors. Let the iterative scheme be u1 = Au0 , u2 = Au1 , ..., up = Aup−1

(2.79)

Substituting (2.78) into the sequence (2.79) we have u1 =

n 

αi Axi =

i=1

n 

since Axi = λi xi

αi λi xi

i=1

u2 =

n 

αi λi Axi =

i=1

n 

αi λ2i xi

i=1

................................. up =

n 

p−1

αi λi

Axi =

i=1

p

αi λi xi

i=1

The final equation can be rearranged, thus:  p up = λ1

n 

(2.80)

α1 x1 +

n 

 αi

i=2

λi λ1



p xi

(2.81)

It is the accepted convention that the n eigenvalues, λi of a matrix are numbered such that |λ1 | > |λ2 | > ... > |λn | Hence



λi λ1

p

tends to zero as p tends to infinity for i = 2, 3, ..., n. As p becomes large, we have from (2.81) p

up ⇒ λ1 α1 x1

(2.82)

Thus up becomes proportional to x1 and the ratio between corresponding components of up and up−1 tends to λ1 . The algorithm, often referred to as the Power Method, is not usually implemented exactly as described previously because problems could arise due to numeric overflows. Usually, after each iteration,

142

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

the resulting trial vector is normalized by dividing it by its largest element, thereby reducing the largest element in the vector to unity. This can be expressed mathematically as ⎫ vp = Aup ⎬   1 p = 0, 1, 2, . . . up+1 = vp ⎭ max(vp )

(2.83)

where max(vp ) is the element of vp with the maximum modulus. The pair of Eqs. (2.83) are iterated until convergence is achieved. This modification to the algorithm does not affect the rate of convergence of the iteration. In addition to preventing the build up of very large numbers, the modification described before has the added advantage that it is now much easier to decide at what stage the iteration should be terminated. Post-multiplying the coefficient matrix A by one of its eigenvectors gives the eigenvector multiplied by the corresponding eigenvalue. Thus, when we stop the iteration because up+1 is sufficiently close to up to ensure convergence, max(vp ) will be an estimate of the eigenvalue. The rate of convergence of the iteration is primarily dependent on the distribution of the eigenvalues; the smaller the ratios |λi /λ1 | where i = 2, 3, ..., n, the faster the convergence. The following M ATLAB function eigit implements the iterative method to find the dominant eigenvalue and the associated eigenvector. function [lam u iter] = eigit(A,tol) % Solves EVP to determine dominant eigensolution by the power method % Sample call: [lam u iter]=eigit(A,tol) % A is a square matrix, tol is the accuracy % lam is the dominant eigenvalue, u is the associated vector % iter is the number of iterations required [n n] = size(A); err = 100*tol; u0 = ones(n,1);

iter = 0;

while err>tol v = A*u0; u1 = (1/max(v))*v; err = max(abs(u1-u0)); u0 = u1;

iter = iter+1;

end u = u0;

lam = max(v);

We now apply the power method to find the dominant eigenvalue and corresponding vector for following eigenvalue problem. ⎡ ⎡ ⎤ ⎤⎡ ⎤ 1 2 3 x1 x1 ⎣ 2 (2.84) 5 −6 ⎦ ⎣x2 ⎦ = λ ⎣x2 ⎦ x3 x3 3 −6 9 >> A = [1 2 3;2 5 -6;3 -6 9]; >> [lam u iterations] = eigit(A,1e-8)

2.16 ITERATIVE METHODS FOR SOLVING THE EIGENVALUE PROBLEM

143

lam = 13.4627 u = 0.1319 -0.6778 1.0000 iterations = 18

The dominant eigenvalue, to eight decimal places, is 13.46269899. The power method can also be used to determine the smallest eigenvalue of a system. The eigenvalue problem Ax = λx can be rearranged to give A−1 x = (1/λ) x Here the iterations will converge to the largest value of 1/λ, that is, the smallest value of λ. However, as a general rule, matrix inversion should be avoided, particularly in large systems. We have seen that direct iteration of Ax = λx leads to the largest or dominant eigenvalue. A second iterative procedure, called inverse iteration, provides a powerful method of determining subdominant eigensolutions. Consider again the eigenvalue problem of (2.40). Subtracting μx from both sides of this equation we have (A − μI) x = (λ − μ) x

(A − μI)−1 x =



 1 x λ−μ

(2.85)

(2.86)

Consider the iterative scheme that begins with a trial vector u0 . Then using the equivalent of (2.83) to (2.86) we have ⎫ vs = (A − μI)−1 us ⎪ ⎬   s = 0, 1, 2, . . . (2.87) 1 ⎪ us+1 = vs ⎭ max(vs ) Iteration will lead to the largest value of 1/(λ − μ), that is the smallest value of (λ − μ). The smallest value of (λ − μ) implies that the value of λ will be the eigenvalue closest to μ and u will have converged to the eigenvector x corresponding to this particular eigenvalue. Thus, by a suitable choice of μ, we have a procedure for finding subdominant eigensolutions.

144

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

Iteration is terminated when us+1 is sufficiently close to us . When convergence is complete 1 = max (vs ) λ−μ Thus, the value of λ nearest to μ is given by λ=μ+

1 max (vs )

(2.88)

The rate of convergence is fast, provided the chosen value of μ is close to an eigenvalue. If μ is equal to an eigenvalue, then (A − μI) is singular. In practice this seldom presents difficulties because it is unlikely that μ would be chosen, by chance, to exactly equal an eigenvalue. However, if A − μI is singular then we have confirmation that the eigenvalue is known to a very high precision. The corresponding eigenvector can then be obtained by changing μ by adding a small quantity and iterating to determine the eigenvector. Although inverse iteration can be used to find the eigensolutions of a system about which we have no previous knowledge, it is more usual to use inverse iteration to refine the approximate eigensolution obtained by some other technique. In practice (A − μI)−1 is not formed explicitly; instead, (A − μI) is usually decomposed into the product of a lower and an upper triangular matrix; explicit matrix inversion is avoided and is replaced by two efficient substitution procedures. In the simple M ATLAB implementation of this procedure shown next, the operator \ is used to avoid matrix inversion. function [lam u iter] = eiginv(A,mu,tol) % Determines eigenvalue of A closest to mu with a tolerance tol. % Sample call: [lam u] = eiginv(A,mu,tol) % lam is the eigenvalue and u the corresponding eigenvector. [n,n] = size(A); err = 100*tol; B = A-mu*eye(n,n); u0 = ones(n,1); iter = 0; while err>tol v = B\u0; f = 1/max(v); u1 = f*v; err = max(abs(u1-u0)); u0 = u1; iter = iter+1; end u = u0; lam = mu+f;

We now apply this function to find the eigenvalue of (2.84) nearest to 4 and the corresponding eigenvector. >> A = [1 2 3;2 5 -6;3 -6 9]; >> [lam u iterations] = eiginv(A,4,1e-8) lam = 4.1283

2.17 SOLUTION OF THE GENERAL EIGENVALUE PROBLEM

145

u = 1.0000 0.8737 0.4603 iterations = 6

The eigenvalue closest to 4 is 4.12827017 to eight decimal places. The functions eigit and eiginv should be used with care when solving large-scale eigenvalue problems since convergence is not always guaranteed and in adverse circumstances may be slow.

2.17 SOLUTION OF THE GENERAL EIGENVALUE PROBLEM There are many algorithms available to solve the eigenvalue problem. The method chosen is influenced by many factors such as the form and size of the eigenvalue problem, whether or not it is symmetric, whether it is real or complex, whether or not only the eigenvalues are required, whether all or only some of the eigenvalues and vectors are required. We now describe the algorithms that are generally used to solve eigenvalue problems. The alternative approaches are as follows. If A is a general matrix, it is first reduced to Hessenberg form using Householder’s transformation method. A Hessenberg matrix has zeros everywhere below the diagonal except for the first sub-diagonal. If A is a symmetric matrix the transform creates a tridiagonal matrix. Then the eigenvalues and eigenvectors of the real upper Hessenberg matrix are found by the iterative application of the QR procedure. The QR procedure involves decomposing the Hessenberg matrix into an upper triangular matrix and a unitary matrix. The following script, e4s213.m uses the M ATLAB function hess to convert the original matrix to a Hessenberg matrix, followed by the iterative application of the QR decomposition using M ATLAB function qr to determine the eigenvalues of a symmetric matrix. Note that in this script we have iterated 10 times, rather than use a formal test for convergence since the purpose of the script is merely to illustrate the functioning of the iterative application of the QR procedure. % e4s213.m A = [5 4 1 1;4 5 1 1; 1 1 4 2;1 1 2 4]; H1 = hess(A); for i = 1:10 [Q R] = qr(H1); H2 = R*Q; H1 = H2; p = diag(H1)’; fprintf(’%2.0f %8.4f %8.4f’,i,p(1),p(2)) fprintf(’%8.4f %8.4f\n’,p(3),p(4)) end

Running script e4s213.m gives

146

1 2 3 4 5 6 7 8 9 10

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

8.3636 9.4940 9.8646 9.9655 9.9913 9.9978 9.9995 9.9999 10.0000 10.0000

6.2420 5.4433 5.1255 5.0329 5.0084 5.0021 5.0005 5.0001 5.0000 5.0000

2.3944 2.0627 2.0099 2.0016 2.0003 2.0000 2.0000 2.0000 2.0000 2.0000

The iteration converges to the values 1, 10, 5, and 2 which are the correct. This QR iteration could be applied directly to the full matrix A but in general it would be inefficient. We have not given details of how the eigenvectors are computed. When the arguments in the M ATLAB function eig are two real or complex matrices, the QZ algorithm is used instead of the QR algorithm. The QZ algorithm (Golub and Van Loan, 1989) has been modified to deal with the complex case. When eig is called using a single complex matrix A then the algorithm works by applying the QZ algorithm to eig(A,eye(size(A))). The QZ algorithm begins by noting that there exists a unitary Q and Z such that QH AZ = T and QH BZ = S and hence both T and S are both upper triangular. This is called generalized Schur decomposition. Providing skk is not zero then the eigenvalues are computed from the ratio tkk /skk , where k = 1, 2, ..., n. The script e4s214.m demonstrates that the ratios of the diagonal elements of the T and S matrices give the required eigenvalues. % e4s214.m A = [10+2i 1 2;1-3i 2 -1;1 1 2]; B = [1 2-2i -2;4 5 6;7+3i 9 9]; [T S Q Z V] = qz(A,B); r1 = diag(T)./diag(S) r2 = eig(A,B)

Running script e4s214.m gives r1 = 1.6154 + 2.7252i -0.4882 - 1.3680i 0.1518 + 0.0193i r2 = 1.6154 + 2.7252i -0.4882 - 1.3680i 0.1518 + 0.0193i

Schur decomposition is closely related to the eigenvalue problem. The M ATLAB function schur(a) produces an upper triangular matrix T with real eigenvalues on its diagonal and complex eigenvalues in 2 × 2 blocks on the diagonal. The matrix A decomposed to give A = UTUH

2.17 SOLUTION OF THE GENERAL EIGENVALUE PROBLEM

147

where U is a unitary matrix such that UH U = I. The output of the script e4s215 shows the similarity between Schur decomposition given by the M ATLAB function schur and the eigenvalues determined using eig, of a given matrix. % e4s215.m A = [4 -5 0 3;0 4 -3 -5;5 -3 4 0;3 0 5 4]; T = schur(A), lam = eig(A)

Running script e4s215 gives T = 12.0000

0.0000

-0.0000

-0.0000

0

1.0000

-5.0000

-0.0000

0

5.0000

1.0000

-0.0000

0

0

0

2.0000

lam = 12.0000 + 0.0000i 1.0000 + 5.0000i 1.0000 - 5.0000i 2.0000 + 0.0000i

We can readily identify the four eigenvalues in the matrix T. The real parts of the eigenvalues lie on the diagonal of T and any non-zero elements in the sub- and super-diagonal of T are the imaginary parts of the eigenvalues. The script e4s216.m compares the performance of the eig function when solving various classes of problem. % e4s216.m disp(’

real1

realsym1

real2

realsym2

comp1

comp2’)

for n = 100:50:500 A = rand(n); C = rand(n); S = A+C*i; T = rand(n)+i*rand(n); tic, [U,V] = eig(A); f1 = toc; B = A+A.’; D = C+C.’; tic, [U,V] = eig(B); f2 = toc; tic, [U,V] = eig(A,C); f3 = toc; tic, [U,V] = eig(B,D); f4 = toc; tic, [U,V] = eig(S); f5 = toc; tic, [U,V] = eig(S,T); f6 = toc; fprintf(’%12.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n’,f1,f2,f3,f4,f5,f6); end

This script, e4s216.m gives the time taken (in seconds) to carry out the various operations. The output is as follows:

148

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

real1 0.188 0.053 0.196 0.100 0.127 0.162 0.202 0.253 0.300

realsym1 0.002 0.009 0.011 0.011 0.013 0.021 0.027 0.038 0.046

real2 0.027 0.060 0.138 0.153 0.244 0.391 0.653 0.886 1.250

realsym2 0.019 0.107 0.120 0.151 0.249 0.375 0.635 0.913 1.270

comp1 0.033 0.148 0.151 0.182 0.236 0.338 0.463 0.610 0.792

comp2 0.142 0.190 0.209 0.276 0.423 0.709 1.066 1.696 2.220

Note that each row is for a different size of n × n matrix, from 100 to 500 in steps of 50. In some circumstances not all the eigenvalues and eigenvectors are required. For example, in a complex engineering structure, modeled with many hundreds of degrees of freedom, we may only require the first 15 eigenvalues, giving the natural frequencies of the model, and the corresponding eigenvectors. M ATLAB provides the function eigs which finds a selected number of eigenvalues, such as those with the largest amplitude, the largest or smallest real or imaginary part, etc. This function is particularly useful when seeking a small number of eigenvalues of very large sparse matrices. Eigenvalue reduction algorithms are used to reduce the size of eigenvalue problem (for example, Guyan, 1965) but still allow selected eigenvalues to be computed to an acceptable level of accuracy. M ATLAB also includes the facility to find the eigenvalues of a sparse matrix. The script e4s217.m compares the number of floating-point operations required to find the eigenvalues of a matrix treated as sparse with the number of floating-point operations required to find the eigenvalues of the corresponding full matrix. % e4s217.m % generate a sparse triple diagonal matrix n = 4000; rowpos = 2:n; colpos = 1:n-1; values = ones(1,n-1); offdiag = sparse(rowpos,colpos,values,n,n); A = sparse(1:n,1:n,4*ones(1,n),n,n); A = A+offdiag+offdiag.’; % generate full matrix B = full(A); tic, eig(A); sptim = toc; tic, eig(B); futim = toc; fprintf(’Time for sparse eigen solve = %8.6f\n’,sptim) fprintf(’Time for full eigen solve = %8.6f\n’,futim)

The results from running script e4s217.m are Time for sparse eigen solve = 0.511121 Time for full eigen solve = 0.591209

Clearly a significant saving in time.

2.18 THE GOOGLE ‘PAGERANK’ ALGORITHM

149

2.18 THE GOOGLE ‘PAGERANK’ ALGORITHM An interesting application of eigenvalue analysis is the Google algorithm for determining the most important pages found in a web search. Search engines use crawlers to continuously search the web, finding new pages and adding them to the vast index of pages. When you undertake a web search, the search engine looks at the index and finds the pages containing the key words used in your search. For example, a Google search for “eigenvalue” finds about 1,670,000 results in 0.54 seconds. The next step is to decide in what order to present these pages: clearly the most relevant pages should be shown first. Different search engines use different ways of ranking the pages. What makes the Google search engine particularly attractive is the ‘PageRank’ algorithm that decides which are the most important pages. The Google ‘PageRank’ algorithm was first described by Brin and Page (1998). The pages containing the key search words can be viewed as a network where a large number of nodes, representing topics or pages, are connected by edges. Such a network is often called a graph. Certain nodes have more connections than others, depending on their importance and relevance. A node which has many other node connections pointing to it, and thus referencing it, may be considered more important than others with fewer references. When searching on the web for specific information amongst the huge numbers of items, it would clearly be useful to rank the pages in terms of their relevance to the user, since this is what the user requires. The Google algorithm provides an objective and efficient page ranking process. It ranks the pages by relevance based on the number of citations a page receives, assigns a page rank to each page, and lists the pages in order of relevance to the user. A key feature of page rank is that a page has high rank if many other pages of high rank point to it, so it is not only the issue of the number of pages which point to a page, the rank of these pages is also important. So how can an algorithm be built to implement this process? First a model of the network of potentially interesting pages must be devised and then we must produce an objective mathematical process to establish the required page ranking in numerical form. From our previous discussion it is clear the model of the pages could be a connected network or graph and this is the approach adopted. Since we wish to use a numerical process we must then represent the graph numerically. This can be achieved by representing the graph as a matrix. If we think of the matrix having elements at positions (i, j ) where i is the row number of the matrix and j is the column of the matrix then the elements of A may be denoted by Aij . In our model of the internet, we may define the elements of the matrix as the probabilities that information will be sent from page j to page i. We can set the probability values for each node i which has a connection to node j , to 1/nj assuming each node j has nj outward pointing connections to other nodes i. If there is no connection to any other node then the probability is of course zero. Thus the sum of probabilities in each column must equal 1 and probabilities lie in the range zero to one. Such a matrix is a column stochastic matrix, defined as one whose elements are less than or equal to one and the columns of which sum to 1. A stochastic matrix has some specific properties, one of which is that its largest or dominant eigenvalue is equal to one. The question now arises, how can we use this model of the internet to rank the internet pages? It can be shown that the dominant eigenvector of the stochastic matrix A which we have used to model the internet links provides the page ranking we require. The numerical ordering of the individual values of the eigenvector rank the relevance of the specific page on the web to the searcher; the largest being

150

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

FIGURE 2.9 Connections of different strengths between five pages of the internet.

the most relevant and the smallest being the least relevant. However, this dominant eigenvector of the stochastic matrix A is not unique. To ensure uniqueness the Google matrix G was introduced as follows: G = αA + (1 − α)UPT where 0 ≤ α < 1 where U is a matrix of 1s and P is sometimes called a personalization vector such that P ≥ 0 and ||P|| = 1. This may be taken as the initial probabilities each given by 1/n where n is the size of the matrix. This matrix, G, the Google matrix, provides a unique dominant eigenvector and is also stochastic. The Google page rank vector r may now be defined by the equations: Gr = r where r ≥ 0 and ||r|| = 1 This is the standard eigenvalue problem, see (2.40), with the dominant eigenvalue equal to one. The Google page rank vector r may now be used to rank the internet pages and display them in order, largest value first and then in decreasing order of magnitude. The final question is: how can we find the dominant eigenvector of G? Fortunately it can be found by using the relatively simple power method, as described in Section 2.16. In the description of the power method in Section 2.16 it is explained that in order to avoid a numeric overflow caused by effectively multiplying the eigenvector by the eigenvalue raised to a power equal to the number of iterations performed, the iteration algorithm must be modified as described by (2.83). This is implemented in the function eigit. In determining the eigensolution of a stochastic matrix, since the dominant eigenvalue is close to and tending to unity, this difficulty does not arise and the algorithm can be simplified. In fact, since we iterate the equations r = Grtrial and rtrial = r this iteration is effectively r = Gp rtrial , where p is suitably large. Thus irrespective of rtrial the product Gp rtrial is the same, but scaled differently. This is because every column of Gp is identical. Consider the graph shown in Fig. 2.9. The graph representing the connections between five scientific papers on the internet with various links connecting each node. Then this diagram may be represented by the stochastic matrix:

2.19 SUMMARY

⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Then the Google matrix is ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

0.0300 0.3133 0.3133 0.3133 0.0300

0 0 1/3 0 1/3 0 1/3 1 0 0

0.0300 0.0300 0.0300 0.8800 0.0300

0 1/2 0 1/2 0

⎤ 0 1/2 0 1/2 ⎥ ⎥ 0 0 ⎥ ⎥ 0 0 ⎦ 1 0

0.0300 0.4550 0.0300 0.4550 0.0300

0.0300 0.0300 0.0300 0.0300 0.8800

151

(2.89)

0.4550 0.4550 0.0300 0.0300 0.0300

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

(2.90)

Note that the sum of each column is one and for each element Gij in the matrix the condition 0 ≤ Gij ≤ 1 is satisfied. Then G50 is ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

0.1468 0.2188 0.0716 0.2880 0.2748

0.1468 0.2188 0.0716 0.2880 0.2748

0.1468 0.2188 0.0716 0.2880 0.2748

0.1468 0.2188 0.0716 0.2880 0.2748

0.1468 0.2188 0.0716 0.2880 0.2748

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

(2.91)

If rtrial = [1 1 1 1 1]T then G50 rtrial = [0.7340 1.0940 0.3580 1.4400 1.3740] In contrast, suppose rtrial = [1 −1 1 −1 1]T then G50 rtrial = [0.1468 0.2188 0.0716 0.2880 0.2748] In both cases, the vectors can be normalized to give r = [0.5097 0.7598 0.2486 1.0000 0.9542] Thus the most important web location or paper is number 4, closely followed by number 5. To implement these procedures as a formal algorithm a test must be included to determine if the procedure has converged to a sufficiently accurate estimate for r.

2.19 SUMMARY We have described many of the important algorithms related to computational matrix algebra and shown how the power of M ATLAB can be used to illustrate the application of these algorithms in a revealing way. We have shown how to solve over- and under-determined systems and eigenvalue problems. We have drawn the attention of the reader to the importance of sparsity in linear systems and

152

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

demonstrated its significance. The scripts provided should help readers to develop their own applications. In Chapter 10, we show how the symbolic toolbox can be usefully applied to solve some problems in linear algebra.

2.20 PROBLEMS 2.1. An n × n Hilbert matrix, A, is defined by aij = 1/(i + j − 1) for i, j = 1, 2, ..., n Find the inverse of A and the inverse of AT A for n = 5. Then, noting that (AT A)−1 = A−1 (A−1 )T find the inverse of AT A using this result for values of n = 3, 4, ..., 6. Compare the accuracy of the two results by using the inverse Hilbert function invhilb to find the exact inverse using (AT A)−1 = A−1 (A−1 )T . Hint: Compute norm(P − R) and norm(Q − R) where P = (AT A)−1 and Q = A−1 (A−1 )T and R is the exact value of Q. 2.2. Find the condition number of AT A where A is an n × n Hilbert matrix, defined in Problem 2.1, for n = 3, 4, ..., 6. How do these results relate to the results of Problem 2.1? 2.3. It can be proved that the series (I − A)−1 = I + A + A2 + A3 + ..., where A is an n × n matrix, converges if the eigenvalues of A are all less than unity. The following n × n matrix satisfies this condition if a + 2b < 1 and a and b are positive: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

a b .. . 0 0

b a .. . 0 0

0 b .. . 0 0

... ... ... ...

0 0 .. . b 0

0 0 .. . a b

0 0 .. . b a

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Experiment with this matrix for various values of n, a, and b to illustrate that the series converges under the condition stated. 2.4. Use the function eig to find the eigenvalues of the following matrix: ⎡ ⎤ 2 3 6 ⎣ 2 3 −4 ⎦ 6 11 4 Then use the rref function on the matrix (A − λI), taking λ equal to any of the eigenvalues. Solve the resulting equations by hand to obtain the eigenvector of the matrix. Hint: Note that an eigenvector is the solution of (A − λI)x = 0 for λ equal to specific eigenvalues. Assume an arbitrary value for x3 .

2.20 PROBLEMS

153

2.5. For the system given in Problem 2.3, find the eigenvalues, assuming both full and sparse forms with n = 10 : 10 : 30. Compare your results with the exact solution given by λk = a + 2b cos {kπ/(n + 1)} , k = 1, 2, . . . 2.6. Find the solution of the over-determined system that follows using pinv, qr, and the \ operator. ⎡ ⎡ ⎤ ⎤ 2.0 −3.0 2.0 ⎡ 1.01 ⎤ ⎢ 1.9 −3.0 ⎢ 1.01 ⎥ 2.2 ⎥ ⎢ ⎢ ⎥ x1 ⎥ ⎢ 2.1 −2.9 ⎥ ⎥ ⎦ ⎣ 2.0 ⎥ x2 = ⎢ ⎢ ⎢ 0.98 ⎥ ⎣ 6.1 ⎣ ⎦ x3 2.1 −3.0 4.94 ⎦ −3.0 5.0 2.1 4.10 2.7. Write a script to generate an n × n matrix E, where E = {1/(n + 1)}C where cij = i(n − i + 1) = ci,j −1 − i = cj i

2.8.

2.9. 2.10.

2.11.

if i = j if j > i if j < i

Having generated E for n = 5, solve Ex = b where b = [1 : n]T by (a) using the \ operator; (b) using the lu function and solving Ux = y and Ly = b. Determine the inverse of E of Problem 2.7 for n = 20 and 50. Compare with the exact inverse which is a matrix with 2 along the main diagonal and −1 along the upper and lower sub-diagonals and zero elsewhere. Determine the eigenvalues of E defined in Problem 2.7 for n = 20 and 50. The exact eigenvalues for this system are given by λk = 1/[2 − 2 cos {kπ/(n + 1)}] where k = 1, ..., n. Determine the condition number of E of Problem 2.7 using the M ATLAB function cond, for n = 20 and 50. Compare your results with the theoretical expression for the condition number which is 4n2 /π 2 . Find the eigenvalues and the left and right eigenvectors using the M ATLAB function eig for the matrix ⎡ ⎤ 8 −1 −5 A = ⎣ −4 4 −2 ⎦ 18 −5 −7

2.12. For the following matrix A, using eigit, eiginv, determine: a. the largest eigenvalue b. the eigenvalue nearest 100 c. the smallest eigenvalue ⎡ 122 41 40 26 25 ⎢ 40 170 25 14 24 ⎢ A=⎢ 27 26 172 7 3 ⎢ ⎣ 32 22 9 106 6 31 28 −2 −1 165

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

154

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

2.13. Given that

⎡

⎤ ⎡ 2 2 2 6 −2 ⎦ and B = ⎣ 4 −1 0 1

1 A=⎣ 5 1 and defining C by

C=

A B B A

⎤ 0 1 −5 1 ⎦ 0 0



verify using eig that the eigenvalues of C are given by a combination of the eigenvalues of A + B and A − B. 2.14. Write a M ATLAB script to generate the matrix ⎡ ⎢ ⎢ ⎢ ⎢ A=⎢ ⎢ ⎢ ⎣

n n−1 n−2 n−1 n−1 n−2 n−2 n−2 n−2 .. .. .. . . . 2 2 2 1 1 1

... 2 ... 2 ... 2 .. .. . . ... 2 ... 1

1 1 1 .. . 1 1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Taking n = 5 and n = 50 and using the M ATLAB function eig, find the largest and smallest eigenvalues. The eigenvalues of this matrix are given by the formula

 1 (2i − 1)π −1 λi = , i = 1, 2 . . . n 1 − cos 2 2n + 1 Verify your results are correct using this formula. 2.15. Taking n = 10, find the eigenvalues of the matrix ⎡ ⎢ ⎢ ⎢ ⎢ A=⎢ ⎢ ⎢ ⎣

1 0 0 .. . 0 1

0 1 0 .. . 0 2

0 0 1 .. . 0 3

... ... ... ... ... ...

0 1 0 2 0 3 .. .. . . 1 n−1 n−1 n

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

using eig. As an alternative, find the eigenvalues of A by first generating the characteristic polynomial using poly for the matrix A and then using roots to find the roots of the resulting polynomial. What conclusions do you draw from these results? 2.16. For the matrix given in Problem 2.12, use eig to find the eigenvalues. Then find the eigenvalues of A by first generating the characteristic polynomial for A using poly and then using roots to find the roots of the resulting polynomial. Use sort to compare the results of the two approaches. What conclusions do you draw from these results?

2.20 PROBLEMS

155

2.17. For the matrix given in Problem 2.14, taking n = 10, show that the trace is equal to the sum of the eigenvalues and the determinant is equal to the product of the eigenvalues. Use the M ATLAB functions det, trace, and eig. 2.18. The matrix A is defined as follows: ⎡ ⎤ 2 −1 0 0 ... 0 ⎢ −1 2 −1 0 ... 0 ⎥ ⎢ ⎥ ⎢ 0 −1 2 −1 . . . 0 ⎥ ⎢ ⎥ A=⎢ .. .. .. .. .. .. ⎥ ⎢ . . . . . . ⎥ ⎢ ⎥ ⎣ 0 0 . . . −1 2 −1 ⎦ 0 0 ... 0 −1 2 The condition number, c, for this matrix is given by c = pnq where n is the size of the matrix and p and q are constants. By computing the condition number for the matrix A for n = 5 : 5 : 50 using the M ATLAB function cond, fit the function pnq to the set of results you produce. Hint: Take logs of both sides of the equation for c and solve the system of over-determined equations using the \ operator. 2.19. An approximation for the inverse of (I − A) where I is an n × n unit matrix and A is an n × n matrix is given by (I − A)−1 = I + A + A2 + A3 + ... This series only converges and the approximation is valid only if the maximum eigenvalue of A is less than 1. Write a M ATLAB function invapprox(A,k) which obtains an approximation to (I − A)−1 using k terms of the given series. The function must find all eigenvalues of A using the M ATLAB function eig. If the largest eigenvalue is greater than one then a message will be output indicating the method fails. Otherwise the function will compute an approximation to (I − A)−1 using k terms of the series expansion given. Taking k = 4, test the function on the matrices: ⎡ ⎤ ⎡ ⎤ 1.0 0.3 0 0.2 0.3 0 ⎣ 0.3 0.2 0.3 ⎦ and ⎣ 0.3 1.0 0.3 ⎦ 0 0.3 1.0 0 0.3 0.2 Use the norm function to compare the accuracy of the inverse of the matrix (I − A) found using the M ATLAB inv function and the function invapprox(A,k) for k = 4, 8, 16. 2.20. The system of equations Ax = b, where A is a matrix of m rows and n columns, x is an n element column vector, b an m element column vector, is said to be under-determined if n > m. The direct use of the M ATLAB inv function to solve this system fails since the matrix A is not square. However, multiplying both sides of the equation by AT gives: AT Ax = AT b AT A is a square matrix and the M ATLAB inv function can now be used to solve the system. Write a M ATLAB function to use this result to solve under-determined systems. The function should allow the input of the b vector and the A matrix, form the necessary matrix products

156

CHAPTER 2 LINEAR EQUATIONS AND EIGENSYSTEMS

and use the M ATLAB inv function to solve the system. The accuracy of the solution should be checked by using the M ATLAB norm function to measure the difference between Ax and b. The function must also include the direct use of the M ATLAB \ symbol to solve the same underdetermined linear system, again include a check on the accuracy of the solution by using the M ATLAB norm function to measure the difference between Ax and b. The function should take the form udsys(A,b) and return the solutions given by the different methods and the norms produced by the two methods. Test your program by using it to solve the under-determined system of linear equations Ax = b where:



 1 −2 −5 3 −10 A= and b = 3 4 2 −7 20 What conclusions do you draw regarding the two methods by comparing the norms that the two methods produce? 2.21. An orthogonal matrix A is defined as a square matrix such that the product of the matrix and its transpose equals the unit matrix or AAT = I Use M ATLAB to verify that the following matrices are orthogonal: ⎡ ⎢ B=⎢ ⎣

C=

√1 3 √1 3 √1 3

√1 6 −2 √ 6 √1 6

cos(π/3) − sin(π/3)

− √12

⎤

⎥ 0 ⎥ ⎦

√1 2

sin(π/3) cos(π/3)



2.22. Write M ATLAB scripts to implement both the Gauss–Seidel and the Jacobi method and use them to solve, with an accuracy of 0.000005, the equation system Ax = b where the elements of A are aii = −4 aij = 2 if |i − j | = 1 aij = 0 if |i − j | ≥ 2 where i, j = 1, 2, ..., 10 and bT = [2 3 4 ... 11] Use initial values xi = 0, where i = 1, 2, ..., 10. (You might also like to experiment with other initial values.) Check your results by solving the system using the M ATLAB \ operator.