Math. PRELIMINARIES 609, Fall 2008 (updated: August 25, 2008) Background: Linear algebra. (1) Review matrix multiplication. ABC where 1 2 3 1 A = 1 -1 0 B = 1 1 4 1 4 For example, compute AB, BC and 1 -1 4 4 1 2 3 . 3 1 -1 0
C=
linearly independent? (4) Review span, basis and dimension. Do the above vectors span R4 , if not, then what is the dimension of the linear space which they span? (5) Suppose A is an n n matrix. What is the range of A? What is the kernel (or null space) of A? Suppose A is the 4 4 matrix made by using the 4 vectors above as columns. What is the dimension of the range of A? What is the dimension of the kernel of A? (6) Suppose A is an n n matrix. Given b Rn , which of the following conditions imply that the linear system Ax = b has a unique solution x Rn ? (a) The determinant of A is nonzero. (b) The determinant of A is zero. (c) The kernel of A is trivial, i.e., Ker(A) = {0} where 0 denotes the zero vector in Rn . (d) The dimension of the kernel of A is greater than zero. (e) The range of A has dimension n. (f) The dimension of the range of A is less than n. (g) The rows of A are linearly independent. (h) The rows of A are linearly dependent.
1
(3) Review linear independence, e.g., are the vectors 1 1 1 1 2 2 1 5 1 1 1 3 1 1 -1 1
(2) Review matrix vector multiplication, e.g., compute Av and Cw where 11 2 4 v = 2 w = . 3 3 2
2
(i) The columns of A are linearly independent. (j) The columns of A are linearly dependent. Discussion: One of the main goals of this class is to study iterative solution of linear systems, i.e., to compute x Rn solving (1.1) Ax = b where A is an n n matrix and b Rn is given. This system has a unique solution if and only if any of the conditions (a), (c), (e), (g), (i) hold (they are all equivalent). In this case, we say that the matrix is non-singular. When any of the conditions (b), (d), (f), (h), or (j) hold then the system either does not have any solution or it has infinitely many of them. In this case, we say that the system is singular. We shall assume that A is such that we have a unique solution. An iterative method for computing the solution of (1.1) is a process which generates a sequence of iterates x1 , x2 , . . . in Rn given an initial iterate x0 Rn . Here is a simple example, (1.2) (1, 2)t , xi+1 = xi + (b - Axi ). Let b = x0 = (0, 0)t (for a row vector v t denotes the corresponding column vector) and (1.3) A= 3 0 . 0 1
Using MATLAB run the above iteration for, say, 20 steps. What is happening? Is it converging? Do the same experiment with (1.4) and A= 1 1/2 . 1/2 1 A= 2 1 1 2
In all but the last example, the iterations were DIVERGING. To investigate the convergence or divergence of an iterative method it is useful to derive a recurrence for the error. We first note that the solution x is a fixed point of the iteration, i.e., x = x + (b - Ax). If the solution is a fixed point for the iteration then we say that the iteration is "consistent." Subtracting the above equations, we find that the error ei = x - xi satisfies (1.5) ei+1 = (I - A)ei .
3
Here I denotes the identity matrix (2 2 in this example). We see that the new error is related to the old by a simple matrix multiplication. Of course, multiplication by a matrix is a linear mapping. An iterative method is called "linear" if the new error results from applying a linear map (matrix multiplication) to the old error. Note that the sequence {xi } converges to x if and only if the sequence of errors converges to the zero vector. Let us consider the first example above. We have ei+1 = from which it follows that ei = (-2)i 0 e 0 0 0 -2 0 e 0 0 i
and is clearly divergent (unless the first component of the initial error is zero). Let us generalize this problem and consider A to be an n n diagonal matrix with entries Ajj = dj , j = 1, . . . , n. Applying the above analysis, we again have (1.5) and so the j'th component of the error ei is given by (ei )j = (1 - dj )i (e0 )j . This method will converge (for any b and any choice of starting iterate x0 ) if and only if - |1 dj | < 1, for j = 1, . . . , n. We would like to have methods which converge on a more general class of problems. One way to generalize the iteration is to introduce an iteration parameter and consider the iteration (1.6) xi+1 = xi + (b - Axi ). The above method is called Richardson's method. This iteration is consistent and a simple manipulation (do it!) gives (1.7) ei+1 = (I - A)ei and so (1.6) is a linear iterative method. Since I and A are diagonal, we have that the j'th component of the error is given by (ei )j = (1 - dj )i (e0 )j , j = 1, . . . , n. Thus, this method is convergent for any right hand side and starting iterate if and only if (1.8)
j=1,...,n
max |1 - dj | < 1.
4
Now suppose that 0 dj 1 , for j = 1, . . . , n with 0 > 0. Then if we take = -1 , 1 0 . 1 For our first example (1.3), we have 1 = 3. Run this example (with = 1/3) for 20 iterations and observe the convergence behavior.
j=1,...,n
max |1 - dj | = max (1 - dj /1 ) 1 -
j=1,...,n
Remark 1. This is not the best choice of iteration parameter. Actually, = 1/2 is optimal. Can you see why? We note that the derivation of (1.7) did not assume anything special about A and so could be used for non-diagonal A. Suppose that M is an invertible n n matrix. Set ei = M ei then from (1.7), ~ M ei+1 = M (I - A)M -1 M ei , i.e., (1.9) ei+1 = (I - A)~i where A = M AM -1 . ~ e Now, {ei } converges to the zero vector if and only if {~i } converges to the e zero vector (why?). The transformation M AM -1 is called a similarity transformation. We will be in great shape if we can find such a transformation which results in a diagonal matrix A since, in that case, we can simply apply our earlier analysis for the diagonal case. A matrix is called diagonalizable if there exists an invertible matrix M with M AM -1 diagonal. The matrix (1.4) is diagonalizable, in fact there is a 2 2 matrix M satisfying A = M AM -1 = 3 0 . 0 1
We see from (1.9) that {~i } will converge to the zero vector for any initial e error ei if we choose = 1/3. Note that the matrix M is used only for the ~ analysis. Run the iteration (1.6) using A given by (1.4) and = 1/3. Remark 2. Not all matrices are diagonalizable. There are, however, fairly general classes of matrices which are known to be diagonalizable. This will be the subject of discussion in a later class. Our method with parameter is, unfortunately, not general enough to deal with all problems. Consider the matrix given by A= 3 0 . 0 -1
5
If you try the above method, you will find that it diverges for any choice of . The analysis leading to (1.8) is still valid so for convergence we would need |1 - 3 | < 1 and |1 + | < 1. This is impossible (even if we used complex ). To deal with this situation, we need to change strategy. We note that since A is non-singular so is At (the transpose of A). Thus, multiplying the equation Ax = b by At does not change the solution set. We consider iterative solution of (1.10) At Ax = At b. The above equation is called the "normal" equation corresponding to the system Ax = b. In our simple example, At = A so we could have just multiplied by A. As we shall see in later classes, multiplication by At will work for more general matrices. We can now apply our method (1.6) to (1.10). In this matrix At A and At b as the right hand side. Now At A diagonal entries {1, 9} and we can choose = 1/9. Run times. Note that you get convergence but it is not as fast iteration applied to the case when A was given by (1.3). Problem 1. Consider A given by A= 3 -4 . -4 3 case, we use the has the positive this iteration 20 as the analogous
Use (1.6) applied to the normal equations to iteratively compute the solution to Ax = b (b and x0 as above). Experiment with different values of and try to determine when you get a convergent iteration. Remark 3. We shall also be interested in solving systems involving matrices with complex coefficients. In this case, the transpose in the normal equations is replaced by the conjugate transpose, A where (A )ij = Aji and the bar denotes the complex conjugate.