White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI

Подождите немного. Документ загружается.

362 CHAPTER 9. KRYLOV METHODS FOR AX = D

nonzero components. If the matrix has a small number of nonzero components,

then one can use a small sequence of Givens transformations to ﬁnd the QR fac-

torization. Other methods for ﬁnding the QR factorization are the row version

of Gram-Schmidt, which generates more numerical errors, and the Householder

In order to formulate the modiﬁed (also called the column version) Gram-

Schmidt method, write the

D = TU in columns

· · · d

· · · t

· · · u

= t

+ t

= t

+ t

+ · · · + t

First, choose t

= d

where u

= (d

)

= Second, s ince t

= 0 for all

n A 1> compute t

= 1u

+ 0= Third, for n A 1 move the columns t

the left side, that is, update column vectors n = 2> ===> p

 t

= t

 t

= t

+ · · · + t

This is a reduction in dimension so that the above three steps can be repeated

on the

q × (p  1) reduced problem.

Example 2. Consider the 4 × 3 matrix

D =

1 1 1

1 1 0

1 0 2

1 0 0

= (d

)

= (

1 1 1 1

)

= 2=

@2 1@2 1@2 1@2

= u

= 1 and d

 t

@2 1@2 1@2 1@2

= u

= 3@2 and d

 t

@4 3@4 5@4 3@4

This reduces to a 4 × 2 matrix QR factorization. Eventually, the QR factoriza-

tion is obtained

D =

1@2 1@2 1@

@2 1@2 1@

1@2 1@2 2@

@2 1@2 2@

2 1 3@2

0 1 1@2

0 0

10@2

transformation, see [16, Section 5.5].

9.4. LEAST SQUARES 363

The modiﬁed Gram-Schmidt method allows one to ﬁnd QR factorizations

where the column dimension of the coe

!cient matrix is increasing, which is

the case for the application to the GMRES methods. Suppose

D> initially

q × (p  1)> is a matrix,whose QR factorization has already been computed.

Augment this matrix by another column vector. We must ﬁnd t

so that

= t

1>p

+ · · · + t

p1

p1>p

+ t

p>p

If the previous modiﬁed Gram-Schmidt method is to be used for the q×(p1)

matrix, then none of the updates for the new column vector have been done.

The ﬁrst update for column p is d

t

1>p

where u

1>p

= t

= By overwriting

the new column vector one can obtain all of the needed vector updates. The

following loop completes the modiﬁed Gram-Schmidt QR factorization when

an additional column is augmented to the matrix, augmented modiﬁed Gram-

Schmidt,

= d

for l = 1> p  1

l>p

= t

 t

l>p

endloop

p>p

= (t

)

if u

p>p

= 0 then

stop

else

= t

p>p

endif.

When the above loop is used with

= Dt

p1

and within a loop with respect

p> this gives the Arnoldi algorithm, which will be used in the next section.

In order to formulate the Givens transformation for a matrix with a small

number of nonzero components, consider the 2 × 1 matrix

D =



The QR factorization has a simple form

D = T

(TU) = (T

T)U = U



By inspection one can determine the components of a 2 × 2 matrix that does

this

= J



f v

v f

364 CHAPTER 9. KRYLOV METHODS FOR AX = D

where

v = e@u

> f = d@u

and u

+ e

= J is often called the Givens

rotation because one can view

v and f as the sine and cosine of an angle.

Example 3. Consider the 3 × 2 matrix

1 1

1 2

1 3

Apply three Givens transformations so as to zero out the lower triangular part

of the matrix:

D =

2 1@

2 0

1@

2 1@

2 0

0 0 1

1 1

1 2

1 3

2 3@

0 1

1 3

D =

3 0 1@

0 1 0

1@

3 0

2 3@

0 1@

1 3

3 2

0 1

and

D =

1 0 0

0 1

3@2



3@2 1@2

3 2@

0 1

3 2

0 0

This gives the "big" or "fat" version of the QR factorization where

T is square

with a third column and

U has a third row of zero components

D = J

U =

5774 =7071 =4082

=5774 0 =8165

=5774 =7071 =4082

1=7321 3=4641

0 1

=4142

0 0

The solution to the least squares problem in the ﬁrst example can be found by

solving

U{ = T



1=7321 3=4641

0 1

=4142



{



5774 =5774 =5774

=7071 0 =7071



10=9697

=8284

9.5. GMRES 365

The solution is

{

= 2=0000 and {

= 2=3333> which is the same as in the ﬁrst

example. A very easy computation is in M

ATLAB where the single command

A\d will produce the least squares solution of D{ = g! Also, the MATLAB

command [q r] = qr(A) will generate the QR factorization of A.

9.4.1 Exercises

1. Verify by hand and by MATLAB the calculations in Example 1.

2. Verify by hand and by M

ATLAB the calculations in Example 2 for the

mo diﬁed Gram-Schmidt method.

3. Consider Example 2 where the ﬁrst two columns in the Q matrix have

been computed. Verify by hand that the loop for the augmented modiﬁed Gram-

Schmidt will give the third column in Q.

4. Show that if the matrix D has full column rank, then the matrix U in

the QR factorization must have an inverse.

5. Verify by hand and by M

ATLAB the calculations in Example 3 for the

sequence of Givens transformations.

6. Show

T = L where T is a product of Givens transformations.

9.5 GMRES

If D is not a SPD matrix, then the conjugate gradient method cannot be directly

used. One alternative is to replace

D{ = g by the normal equations D

D{ =

g> which may be ill-conditioned and subject to signiﬁcant roundo errors.

Another approach is to try to minimize the residual U({) = u({)

u({) in place

M ({) =

{

D{{

g for the SPD case. As in the conjugate gradient method,

this will be done on the Krylov space.

Deﬁnition. The generalized minimum residual method (GMRES) is

{

p+1

= {

l=0



where u

= g  D{

and 

5 R are chosen so that

U({

p+1

) = min

U(|)

| 5 {

+ N

p+1

and

p+1

= {} | } =

l=0

> f

5 R}=

366 CHAPTER 9. KRYLOV METHODS FOR AX = D

Like the conjugate gradient method the Krylov vectors are very useful for

the analysis of convergence. Consider the residual after

p + 1 iterations

g  D{

p+1

= g  D({

+ 

+ · · · + 

)

 D(

+ 

+ · · · + 

)

= (

L  D(

L + 

D + · · · + 

))u

Thus

p+1



p+1

(D)u

(9.5.1)

where

p+1

(}) = 1  (

} + 

}

+ · · · + 

}

p+1

)= Next one can make appro-

priate choices of the polynomial

p+1

(}) and use some properties of eigenvalues

Let

D be an q × q invert-

ible matrix and consider D{ = g=

1. GMRES will obtain the solution within q iterations.

2. If

g is a linear combination of n of the eigenvectors of D and D =

Y Y

where Y Y

= L and  is a diagonal matrix, then the GMRES

will obtain the solution within

n iterations.

3. If the set of all eigenvalues of

D has at most n distinct eigenvalues and

D = Y Y

1

where  is a diagonal matrix, then GMRES will obtain

the solution within

n iterations.

The Krylov space of vectors has the nice property that

 N

p+1

= This

allows one to reformulate the problem of ﬁnding the 

D({

p1

l=0



) = g

p1

l=0



l+1

= g

p1

l=0



l+1

= u

= (9.5.2)

Let bold K

be the q × p matrix of Krylov vectors

· · · D

p1

The equation in (9.5.2) has the form

 = u

where (9.5.3)

= D

· · · D

p1

and

 =



· · · 

p1

and matrix algebra to prove the following theorem, see C. T. Kelley [11, Chapter

3].

Theorem 9.5.1 (GMRES ConvergenceProperties)

9.5. GMRES 367

The equation in (9.5.3) is a least squares problem for

 5 R

where DK

q × p matrix.

In order to e

!ciently solve this sequence of least squares problems, we

construct an orthonormal basis of

one column vector per iteration. Let

= {y

> y

2>···

> y

} be this basis, and let bold V

be the q × p matrix whose

columns are the basis vectors

· · · y

Since DN

 N

p+1

, each column in DV

should be a linear combination of

columns in V

p+1

= This allows one to construct V

one column per iteration

by using the modiﬁed Gram-Schmidt process.

Let the ﬁrst column of V

be the normalized initial residual

= ey

where e = ((u

)

is chosen so that y

= 1= Since DN

 N

> D times

the ﬁrst column should be a l inear combination of

and y

= y

+ y

Find k

and k

by requiring y

= y

= 1 and y

= 0 and assuming

 y

is not the zero vector

= y

}

= Dy

 y

= (}

})

and

= }@k

For the next column

= y

+ y

Again require the three vectors to be orthonormal and Dy

 y

not zero to get

= y

and k

= y

} = Dy

 y

= (}

})

and

= }@k

368 CHAPTER 9. KRYLOV METHODS FOR AX = D

Continue this and represent the results in matrix form

= V

p+1

K where (9.5.4)

· · · Dy

p+1

· · · y

p+1

· · · k

0 k

· · · k

0 0

0 0 0 k

p+1>p

l>p

= y

for l  p> (9.5.5)

} = Dy

 y

1>p

· · ·  y

p>p

6= 0>

p+1>p

= (}

})

and (9.5.6)

p+1

= }@k

p+1>p

= (9.5.7)

Here

D is q × q> V

is q × p and K is an (p + 1) × p upper Hessenberg matrix

(

= 0 when l A m + 1). This allows for the easy solution of the least squares

problem (9.5.3).

Theorem 9.5.2 (GMRES Reduction) The solution of the least squares prob-

lem (9.5.3) is given by the solution of the least squares problem

K = h

e (9.5.8)

where

is the ﬁrst unit vector, e = ((u

)

and DV

= V

p+1

Proof. Since u

= ey

> u

= V

p+1

e= The least squares problem in (9.5.3)

can be written in terms of the orthonormal basis

 = V

p+1

Use the orthonormal property in the expression for

bu() = V

p+1

e  DV

 = V

p+1

e  V

p+1

K

(bu())

bu() = (V

p+1

e  V

p+1

K)

p+1

e  V

p+1

K)

= (

e  K)

p+1

e  K)

= (h

e  K)

e  K)=

Thus the least squares solution of (9.5.8) will give the least squares solution of

(9.5.3) where K

 = V

=

9.5. GMRES 369

} = Dy

 y

1>p

· · ·  y

p>p

= 0> then the next column vector y

p+1

cannot be found and

= V

K(1 : p=1 : p)=

Now K = K(1 : p=1 : p) must have an inverse and K = h

e has a solution.

This means

0 =

 DV



= g  D{

 DV



= g  D({

+ V

)=

If } = Dy

 y

1>p

· · ·  y

p>p

6= 0, then k

p+1>p

= (}

})

6= 0 and

= V

p+1

K= Now K is an upper Hessenberg matrix with nonzero compo-

nents on the subd iagonal. This means K has full column rank so that the least

squares problem in (9.5.8) can be solved by the QR factorization of

K = TU=

The normal equation for (9.5.8) gives

K = K

e and

U = T

e= (9.5.9)

The QR factorization of the Hessenberg matrix can easily be done by Givens

rotations. An implementation of the GMRES method can be summarized by

the following algorithm.

GMRES Metho d.

let

{

b e an initial guess for the solution

= g  D{

and Y (:> 1) = u

@((u

)

for k = 1, m

Y (:> n + 1) = DY (:> n)

compute columns

n + 1 of Y

n+1

and K in (9.5.4)-(9.5.7)

(use modiﬁed Gram-Schmidt)

compute the QR factorization of

(use Givens rotations)

test for convergence

solve (9.5.8) for



{

n+1

= {

+ Y

n+1



endloop.

The following M

ATLA B code is for a two variable partial dierential equation

with both ﬁrst and second order derivatives. The discrete problem is obtained

by using centered di

erences and u pwind dierences for the ﬁrst order deriv-

atives. The sparse matrix implementation of GMRES is used along with the

The code is initialized in lines 1-42, the GMRES loop is done in lines 43-

87, and the output is generated in lines 88-98. The GMRES l oop has the

sparse matrix product in lines 47-49, SSOR preconditioning in lines 51-52, the

SSOR preconditioner, and this is avariation of the code in [11, chapter 3].

370 CHAPTER 9. KRYLOV METHODS FOR AX = D

mo diﬁed Gram-Schmidt orthogonalization in lines 54-61, and Givens rotations

are done in lines 63-83. Upon exiting the GMRES loop the upper triangular

solve in (9.5.8) is done in line 89, and the approximate solution

{

+ V

n+1

 is

generated in the loop 91-93.

MATLAB Code pcgmres.m

1. % This code solves the partial dierential equation

2. % -u_xx - u_yy + a1 u_x + a2 u_y + a3 u = f.

3. % It uses gmres with the SSOR preconditioner.

4. clear;

5. % Input data.

6. nx = 65;

7. ny = nx;

8. hh = 1./nx;

9. errtol=.0001;

10. kmax = 200;

11. a1 = 1.;

12. a2 = 10.;

13. a3 = 1.;

14. ac = 4.+a1*hh+a2*hh+a3*hh*hh;

15. rac = 1./ac;

16. aw = 1.+a1*hh;

17. ae = 1.;

18. as = 1.+a2*hh;

19. an = 1.;

20. % Initial guess.

21. x0(1:nx+1,1:ny+1) = 0.0;

22. x = x0;

23. h = zeros(kmax);

24. v = zeros(nx+1,ny+1,kmax);

25. c = zeros(kmax+1,1);

26. s = zeros(kmax+1,1);

27. for j= 1:ny+1

28. for i = 1:nx+1

29. b(i,j) = hh*hh*200.*(1.+sin(pi*(i-1)*hh)*sin(pi*(j-1)*hh));

30. end

31. end

32. rhat(1:nx+1,1:ny+1) = 0.;

33. w = 1.60;

34. r = b;

35. errtol = errtol*sum(sum(b(2:nx,2:ny).*b(2:nx,2:ny)))^.5;

36. % This preconditioner is SSOR.

37. rhat = ssorpc(nx,ny,ae,aw,as,an,ac,rac,w,r,rhat);

38. r(2:nx,2:ny) = rhat(2:nx,2:ny);

39. rho = sum(sum(r(2:nx,2:ny).*r(2:nx,2:ny)))^.5;

9.5. GMRES 371

40. g = rho*eye(kmax+1,1);

41. v(2:nx,2:ny,1) = r(2:nx,2:ny)/rho;

42. k = 0;

43. % Begin gmres loop.

44. while((rho

A errtol) & (k ? kmax))

45. k = k+1;

46. % Matrix vector product.

47. v(2:nx,2:ny,k+1) = -aw*v(1:nx-1,2:ny,k)-ae*v(3:nx+1,2:ny,k)-...

48. as*v(2:nx,1:ny-1,k)-an*v(2:nx,3:ny+1,k)+...

49. ac*v(2:nx,2:ny,k);

50. % This preconditioner is SSOR.

51. rhat = ssorpc(nx,ny,ae,aw,as,an,ac,rac,w,v(:,:,k+1),rhat);

52. v(2:nx,2:ny,k+1) = rhat(2:nx,2:ny);

53. % Begin modiﬁed GS. May need to reorthogonalize.

54. for j=1:k

55. h(j,k) = sum(sum(v(2:nx,2:ny,j).*v(2:nx,2:ny,k+1)));

56. v(2:nx,2:ny,k+1) = v(2:nx,2:ny,k+1)-h(j,k)*v(2:nx,2:ny,j);

57. end

58. h(k+1,k) = sum(sum(v(2:nx,2:ny,k+1).*v(2:nx,2:ny,k+1)))^.5;

59. if(h(k+1,k) ~= 0)

60. v(2:nx,2:ny,k+1) = v(2:nx,2:ny,k+1)/h(k+1,k);

61. end

62. % Apply old Givens rotations to h(1:k,k).

63. if k

64. for i=1:k-1

65. hik = c(i)*h(i,k)-s(i)*h(i+1,k);

66. hipk = s(i)*h(i,k)+c(i)*h(i+1,k);

67. h(i,k) = hik;

68. h(i+1,k) = hipk;

69. end

70. end

71. nu = norm(h(k:k+1,k));

72. % May need better Givens implementation.

73. % Deﬁne and Apply new Givens rotations to h(k:k+1,k).

74. if nu~=0

75. c(k) = h(k,k)/nu;

76. s(k) = -h(k+1,k)/nu;

77. h(k,k) = c(k)*h(k,k)-s(k)*h(k+1,k);

78. h(k+1,k) = 0;

79. gk = c(k)*g(k) -s(k)*g(k+1);

80. gkp = s(k)*g(k) +c(k)*g(k+1);

81. g(k) = gk;

82. g(k+1) = gkp;

83. end

84. rho=abs(g(k+1));