ELo ELSEVIER Applied Mathematics and Computation 92(1998)49 58 The generalized Cholesky factorization method for saddle point problems Jinxi Zhao Department of Mathematics. Nanjing Unicersity. Nanjing 210 008. China Abstract Over the past 10 years. a variety of iterative methods for saddle point problems have been proposed. In this paper, we present a class of direct methods, the so-called gener- alized Cholesky factorization method. for solving linear systems arising from saddle point problems or discretization of the Stokes equations. Numerical results illustrate the efficiency of new methods given in this paper. 1998 Elsevier Science Inc. All rights reserved. Keywords: Saddle point problem: Generalized Cholesky factorization: Condition I: Generalized positive definite 1.Introduction Consider the linear system of equations )-( where A is symmetric positive definite, B of full row rank, and C symmetric positive semi-definite. Problems in this class arise frequently in the context of minimization of quadratic forms subject to linear constraints [8]. An important example arises from the numerical discretization of the Stokes equations. In particular, we are concerned with the discretization of Project supported by the Science Foundation of the State Education Commission of China and the "863" plan of China. 0096-3003/98/S19.00 1998 Elsevier Science Inc. All rights reserved. PII: S0096-3003(97)10040-6
50 .wf,.1th.(au.2190g)495N -A4-下pf. Vu -0).on 2. 4=0.oni2. (2) /P=0 where is a simply connected bounded domain in R.s==2 or 3.This system of the Stokes equations is a fundamental problem arising in computational fluid dynamics [4].Discretization of Eq.(2)by finite difference or finite element techniques leads to a linear system of equations of (1). In recent years.a variety of iterative algorithms have been devised for solving saddle point problems [1.3.4.6.7].In this paper.we have developed the generalized Cholesky factorization for four typical matrices arising in numerical optimization and computational fluid dynamics.Using the matrix factorization.we establish a class of direct methods for solving the corre- sponding linear system.New methods proposed in this paper remain main advantages of the classical Cholesky factorization for positive definite sys- tems.Hence the new method is referred to as the generalized Cholesky fac- torization method. In the following we always assume that matrices 4E R"".BR",and C"satisly the following condition. Condition I.A is symmetric positive definite.B is of full row rank and C is symmetrie semi-positive definite. 2.Symmetric indefinite case Let us assume that G is an (m+n)*(m+n)matrix and express it as 3 where A.B and C satisfy Condition I.It is easy to see that G is symmetric in- definite.The purpose of this section is based on the matrix factorization of G to give a new algorithm for solving the linear system (1). Firstly.we can prove the following theorem. Theorem 1.Ler G he an (m+*im+n)matrix expressed hy Ey.(3)and A.B.C sutisfr Condition I.Then there ahuvs exists the factorizution form G=B B (4)
J.Zhuo Appl.Muth.Comput.92 (1998)49.58 51 where andL1∈Rmwi8 low triangular,.Ln∈is low triangular..Lg∈R×m. Proof.Since 4 is symmetric positive definite,there always exists the Cholesky factorization A=LiLi where LE Ris nonsingular low triangular.Take La B(LI) (51 so that B=LgLi It follows that Le is full row rank because B is full row rank.Because C is sym- metric semi-positive definite,the matrix C+LaL must be symmetric positive definite.Hence we have the Cholesky factorization LuL=C+LaLg (6) also let 4=(2)=(日) Thus we have L1L=（ Remark 1.From Theorem 1.we set that matrices Li and L can be obtained conveniently as long as submatrices L.Lg and L have been computed.This is why our method is as fast as the classical Cholesky factorization for symmetric positive definite. We now discuss the realization of the generalized Cholesky factorization (4) of G.Let A=a∈Rmxm.L=l, Using the Cholesky factorization of A,the elements of L can be computed from
52 J.Zhao Appl.Math.Comput.92 (1998:49 58 ii月 i.j=1: Set LB=lg,=B(L)∈R"m Then all elements g;can be obtained from the following: (8) Since LyLT=C+LaLg let Lw=比∈E”“ be low triangular. Hence i1: i=1:n.j=1:n (9) From the above discussion.the triangular factor of G can be obtained pro- vided that submatrices L.L and L have been formed.i.e.. 4=(位) and =(日)》 Using the factorization expression G=LiL
J.Zhuo Appl.Muth.Comput.92 (1998;49 58 53 it is easy to give the solution procedure of the linear system (1).Here we des- cribe the following. Algorithm 1. 1.Given 4=[d E R*",B=b,ER"and C=Ci E R"satisfying Con- dition I and given f∈Rm,g∈R”. 2.L==chol(A)or using expression (7). 3.Computing L8 =[g from LaL=B or using expression (8). 4.L =[v]chol(C+LgL)or using expression (9). 5.Computing from ,2)0)() 6.Obtained the final solution of the linear system (1)from (日))-) As a special case,if we have C=0 in the linear system (1).i.e.. -(8) then we can obtain the analogous factorization form. Using matrices La and La obtained in Theorem I,and Lg full row rank.we have the Cholesky factorization of LaL.i.e.. L,LX=LBL牙 where LER"*"is low triangular.Hence the following theorem is proved Theorem 2.Let A G- B B O and matrices A and B sutisfy Condition I.then there always exists a factorizution G:=L2Lg, (10) where
54 J.Zhuo Appl Muth Comput.92 (1998,49 58 LE is low triangular and the Cholesky factor of A.L8=B(L) L.ER""is low triangular and the Cholesky factor of LaLg. 3.Generalized semi-positive definite case Consider the linear system (8))() (11 where matrices A,B and C satisfy Condition I. Let 4-B G=(8c (12) and obviously.the matrix G3 is nonsymmetric.But the symmetric part of G3. i.e..(G+G3)/2.is symmetric semi-positive definite [2.5].So we call G3 the generalized semi-positive definite. As in the above discussion,we can always form the low triangular matrix L from symmetric positive definite matrix A and La B(LI) Since C'is symmetric semi-positive definite and La full row rank.there always exists the Cholesky factorization,i.e.. L LI =C+LaLg where L,ER""is low triangular.Thus we have proved the following result. Theorem 3.Assume thut Ga is of where A,B and C satisfv Condition I.Then there always exists the generalized Cholesky factorization. G3 =L3Lg. (13) where =(2) (日) La ER"is low triangular and the Cholesky factor of A.L.ER""is low trian- gulur and the Cholesky factor of C+LaLg.Ln B(LI)E R
J.Zhuo Appl.Math.('omput.93 1998.49 58 55 Since the realization procedures of submatrices L.La and L,are the same as in Section 2.Hence we can propose an algorithm for solving the linear system Algorithm 2. 1.Given 4=lanER".B=huR"xm and C'=cul E R""satisfying Condition I and given f∈R".g∈R" 2.L=]=chol(4)or using expression (7). 3.ComputingL=gl from Ll.B or using expression (8). 4.L=[r=chol(C+LaL)or using expression (9). 5.Computing =( from ())(） 6.Obtained the final solution of'the linear system (1)from (:）)() Consider the special case where C=0 in Eq.(12).i.e.. Obviously.we can prove the following Theorem 4./f G= und A.B sutisfr Condition 1.then we hate Ga =LaL. 14) 17r 4-(伦2)=（日） Here matrices L and L&are the same as in Theorem 3.and L-is the Cholesky fac- tor ofL Since L is of full row rank.there always exists the Cholesky factor L. of LaLk
56 .Zhuw/1ppl.Mh.C.921199然！495 4.Discussion The method described in this paper retains the main advantages of the clas- sical Cholesky factorization.During the course of the computation.N =m+n square roots must be taken.Condition I assures us that the arguments of these square roots will be positive.About N/6 flops are needed beyond n square roots.Finally.because 4 is symmetric positive definite,the elements of L will be controllable.In fact,we have the following relation, 1a≤a.i=1:m.k=1:i (15) Since L8=B(LI). it follows that ca+ If we take M=,max（g+·+gim) 5成 then ≤cw+M.i=1:n.j=1:i. That is,the elements of L (or Ly.L.L:cannot become too large. 5.Numerical results We now present the results of numerical experiments for solving Eqs.(1) and (11).All experiments were performed in MATLAB on a PC-386 computer. Take Am=a=Hn+Inm∈Rm" where Hm= [1 i+i-1 is an m x m Hilbert matrix and/is an m x m unit matrix. Also take B=h,l=maxi.i川∈R"m and Cw=c=UG,U∈R
J.Z0/pl.Mh.Cr.9211998)49.58 57 where 2 Un=1- ww w=(:n wTw and om=diag(1,2...n-l.(0. Thus.matrices A.B and C satisfy Condition I. Example 1.Solve the linear system (。))() where =+2xm+=m 8=∑xj-xm+,1=m The vectors =(=( denote the computed solution and the exact solution.respectively Example 2.Consider the linear system (g))() Table 1 The result of Example 1 Order Algorithm 1 Gauss elimination 名 n Flops Norm (- Flops Norm (x-x') I 10 3998 9.4259×101 8207 2.8856×101: 30 10 11533 3,4882×1031 24125 1.0046×10u 3(0 20 50)321 4.7859×10-0 100187 2.2251×10m 5(0 30 188938 6.1818×109 384447 1.7993×10" 50 40 267548 1.7401×10* 540643 4.4792×10" 50 50 344524 2.0480×10 733965 9.4886×104 The Gauss elimination is provided by MATLAB
58 J.Zhoo Appl.Muth.Comput.92 (1998:49 5& Table 2 The result of Example 2 Order Algorithm 2 Gauss elimination U Flops Norm (- Flops Norm (-x 10 10 3998 6.7242×108 8207 1.8781×101日 20 10 11533 2.5209×101 24125 1.0121×100 30 20 50321 5.2676×100 10187 1.9155×10" 50 30 18938 6.3810×109 384447 1.3570×10" 50 40 267548 8.712s×109 540643 8.1077×10“ 50 50 344524 1.0074×10s 733965 1.5418×0× The Gauss elimination is provided by MATLAB. where: R=2,-m+从1=n From the above result we can see that the generalized (holesky method pre- sented in this paper will be efficient enough also for practical application.In fact,it would be still efficient when C=0 in linear systems (1)or (11). References [1]R.E.Bank.B.D.Welfert.H.Yserentant.A class of iterative methods for solving saddle point problems.Numer.Math.56 (1990)645 666. [2]A.Bjorck.Least squares methods.in:P.G.Ciarlet.J.L.Lions (Eds.)Handbook of Numerical Analysis vol.I:Solution of Equations in R".Elsevier.Amsterdam.1990. [3]J.H.Bramble.J.E.Pasciak.A preconditioning technique for indefinite systems resulting from mixed approximations of elliptic problems.Math.Comp.50 (1988)I 17. [4]H.C.Elmain.G.H.Golub.Inexaet and preconditioned Uzawa algorithms for saddle point problems.Technical Report NA-93-012.Computer Science Depariment.Stanford University. 1993. [5]G.H.Golub.C.F.Van Loan.Matrix Computations.2nd ed..The Johns Hopkins University Press.Baltimore.1989. [6]G.H.Golub.M.L.Overton.The convergence of inexact Chebyshe and Richardson iterative methods for solving linear systems.Numer.Math.53 (1988)571 593. [7]1.Gustafsson.A class of first order factorizations.BIT 18 (1978)142 156. [8]Y.Yuan.Numerical methods for nonlinear programming.Modern Mathematics Series. Shanghai Scientific and Technical Publishers.1993