ON THE INVERSE OF A MATRIX WITH SEVERAL RANK ONE UPDATES Romesh Saigal Department of Industrial and Operations Engineering The University of Michigan Ann Arbor, MI 48109-2117 Technical Report 93-41. December 1993

On the inverse of a matrix with several rank one updates Romesh SAIGAL Department of Industrial and Operations Engineering, The University of Michigan, Ann Arbor, Michigan 48109-2117, USA December 28, 1993 Abstract In this short note, we consider the problem of computing the inverse of a symmetric positive-definite matrix after several rank one updates have been subtracted from a similar matrix. We give a generalization of the Sherman-Morrison-Woodbury formula for a single rank one modification. We also present a scheme which generates this inverse in m2r + mr2 + (-1) multiplications where m is the size of the matrix and r is the number of rank one updates. In case these updates are of the identity natrix, the resulting multiplications is mr2 + r(r-1) As an application, we present 2 an alternative implementation of the simplex method that uses only positive definite symmetric matrices, and may be suitable for parallel/vector/distributed computing environments. This also results in a unification of the implementations of interior point methods and simplex method. Key words: Symmetric positive-definite matrix, rank one updates, Sherman-MorrisonWoodbury formula, Interior point methods, Simplex method, Vector/parallel/distributed computing. Abbreviated title: Inverse of a matrix 0

1 Introduction In this paper we derive the following formula for a symmetric positive definite matrix: (A - EET)-1 = A-1 + GDGT (1) where E, and G are m x r matrices, A and A - EET are m x m symmetric positive-definite matrices; and, D is a diagonal matrix with positive diagonal entries. The matrix appearing on the left hand side of formula (1) can also be viewed as r rank one modifications of A, with r A- EET = A- EjEj j=1 where E = (El, E2,*, Er). The formula (1) is a generalization of the following ShermanMorrison-Woodbury formula, Hager [4], when r = 1: (A- uuT)-' = A-1 + vDvT (2) where At = u and D = (1 - uTA-lu)-l which is positive if the resulting matrix after the rank one modification has a positive determinant. Matrices of the form A - LUD-V are called Schur complements, and an excellent review andl applicationls of such matrices can be found in Cottle [2]. We deal here with situations tlhat generate a symllnetric and positive definite Schur complement. Examples of such forms arise (luring Cholesky factorizations of large symmetric positive definite systems, George awtd Liu [3], and partitioning techniques for interior point methods, Saigal [6]. In this paper, we derive formula (1) for symmetric positive definite case. We do this by an algorithmn which, inductively, uses the formula (2) and generates the required matrices D and; using m2r + mr2 + ( multiplications. When A is the identity matrix, this numl)er is i7mr2 + r(-1) In addition, we present an implementation of the simplex method that uses only symmetric positive definite matrices. This is done in an attempt to unify the implemenltations of interior point methods and simplex method. This implementation contrasts with the standard implementation of Cholesky factorization for simplex method, see for example Bartels, Golub and Saunders [1]. This paper consists of four sections. The introduction section 1 is followed by section 2 1

which presents the inductive scheme. Section 3 gives an alternative formula A is the identity matrix and r > m. Section 4 presents the alternative implementation of the Simplex method. 2 The Inductive Scheme We will now show that there are matrices G and D which satisfy the formula (1). For this purpose we introduce the following: Define [ A if k = O Bk+l - Bk- EkET if 1 <k < r and note that Br+1 = A - EET. Also, for each k = 1,, r, define BkE(k) = Ek. We can then prove: Theorem 1 There exist matrices G and D such that D is diagonal with positive diagonal entries, and (A- EE)-1 = A-1 + GDGT where, for each k = 1,., r, Gk = E(k) and Dkk = (1- < Ek) Ek >)-1 Proof: Since A - EET is positive definite, for each k = 1,, r Bk is positive definite. Also, det(Bk+l) = det(Bk)(1- < Bk1EkEk >), and so for each k = 1,.,r, Dkk > 0. From the formula (2), the theorem holds for r = 1. Thus, assume that the theorem holds some for r = 1 - 1 1 > 2. We now show the result for r = 1, and thus the theorem follows by induction. Now Bl + = Bi - El Er thus from the Sherman-Morrison-Woodburv formula B^1 = Brl +B 1EI(BiIEl)rT 1- < El, B- El > From the induction hypothesis, B-1 = A-1 + I — j —-E —E) E = 12- < E, Ej > I

and our result follows. U Given E'j) for j = 1,..,, Ek, +) can be readily computed by using Theorem 1, which gives Bk+1 = A-1 + G(k)D(k)G(k) with G(k) = (E(1), E(2)., Ek)) and D(k) = (1- < Ej, Ej) >)1, and E(k+l) = A-Ek+l + G(k)D(k)G(k)TEk+. (3) We can then prove: Theorem 2 Given A-', G and D can be generated by m2r + mr2 + r(-1) multiplications. In case A is the identity matrix, the number of multiplications is mr2 + 2r() Proof: Since B1 = A, generating E(1) requires m2 multiplications; and none when A is the identity matrix. It takes m multiplications to generate Dll. Given E(1) through E(k), using the formula (3), it is readily seen that it takes m2+ 2mk + k multiplications. In case A is the identity matrix, this number is 2mk + k. Also, it takes m multiplications to generate Dk+l,k+. Thus, the number of multiplications required is r-l r(r - 1) nr2r +,L(2mk + k) + mr = m2r + mr2 + 2 k=1 Tile result when A is the identity matrix follows by the same argument.! The inductive procedure used in the proof of the Theorem 1 may not exist when the matrix A - EET is indefinite. As an example consider I - uuT - vvT where uTu = 1 and 7' r = 1 with UTr $ 0. In this case both the matrices I - uuT and I - vvT are singular, but tle lIatrix after two rank 1 updates is non-singular and thus has an inverse. It may have the form of the inverse stipulated by formula (1), but it cannot be generated by the procedure developed in this section. 3 Rank One Updates of the Identity Iri nlany applications, like Saigal [5] and [6], several rank one updates of the identity matrix are involved, i.e., the matrix to be inverted is I - EET 3

where E is an m x r matrix. In case r < m, the result of Theorem 2 shows that the inverse of this matrix should be of the form given in formula (1). Otherwise, one can use the following standard identity, Hager [4], (I - EET)-1 = I + E(I - ETE)-ET. Since I - EET is positive definite if and only if I - ETE is, the result of Theorem 1 gives, in rm2 + m(,M-1) multiplications, the following identity (I - ETE)-1 = I+ GTDG where G is an m x r matrix. Thus we obtain an alternative formula for the inverse that follows (I - EET)-1 = I + EET + EGTDGET. 4 An Application An application, of the formula developed here, to interior point methods for solving the linear programming problem, can be found in Saigal [5] and [6]. In this section we show how the formula can be used to implement the simplex method. This can be a framework for unifying the implementation of these two methods. We consider the following linear programming problem: minimize cTz Ax = b (4) x > 0 where A is an 17 x n matrix of full rank im, and, b and c are appropriate vectors. We will assume that a feasible basis B is available, and familiarity with the simplex method will be assumed. Let N = {ji,j2,'',n-m} be the set of indices of the non-basic variables. We now generate an implementation that uses the Cholesky factorization of the matrix AAT. For large and sparse matrices, a sparse factor can be generated by using the techniques of George and Liu [3]. 4

Let L be a lower triangular matrix which is the Cholesky factor of AAT, i.e., AAT = LLT. Then, for the given basis matrix B, we have BBT = AAT - S AjAT. (5) jEN where Aj is the jth column of the matrix A. Thus, BBT = L(I- E AjAj)LT jEN where LA = Aj for each j E N. Defining D= I- E AjAj (6) jEN we can write BBT = LDLT. (7) Our aim is to use the decomposition (7) during the steps of the simplex method, and to update this decomposition during the iterations of the method. This involves updating D when one column in B is replaced by one indexed in N. During the application of the simplex method, two linear systems, the primal system Ba = a (8) for some give vector a, and the dual system B'! = CR (9) where CB is the cost vector of the basic variables, are solved. In addition, the reduced costs for each nonbasic index j C N Cj = c, - yTAj (10) are also computed. It is readily confirmed that, using the decomposition (7), system (8) is solved by the following sequence of steps: Step 1 Compute a such that Da = L-a 5

Step 2 Compute a such that a= L-a Step 3 Compute the solution a by a = BTa. Also, the system (9) can be solved by the following: Step 1 Compute y such that D = L-BcB Step 2 Compute the solution y by y = L-TY. And, we note that j = - yTAj = cj - TL-Aj = cj - Aj, and step 2 of the dual system can be eliminated. We note that at Step 1 of both the primal and dual systems, the inverse of matrix D is required. In that sense, we will consider this implementation as an application of the main formula presented in this paper. 4.1 Updating D-1 We now show' how the formula for D-' can be updated when the column Aj, leaves the basic set, and the column Aj, enters. That is, the new non-basic index set is N = {ji,'. js-l js+i',Jrjv}. Thus, if the updated matrix D = I - EET then Ek 1 < k < s-1 Ek = Ek+l s < k < r-1 E, k = r 6

and _kBk 1 < k < s-1 Bk = Bk+l+EsE s < k < r with Br+i = Br+i + EsE[ - EEv. From the above definitions, and section 2, we readily see that for each 1 < k < s -1 -l = Bk-l = I + G(k)F(k)G(k)T Now, consider s < k < r - 1. Defining BkEk = Ek we have (Bk+i + EE )Ek) = Ek+1 or,(k+l) > E(k)- (k+l) < EsI k+l > (k+l ) Ek ) - Ek+l )S- KE i+l, E+) k- i sl 1+ < E,, Ek+l) > We now show how this can be obtained. Define a = GTEs, which has been obtained during the'column updating' phase of the simplex method, and let a(k) = (al,..., ak)T. Then E(k+1) = E, + G(k)F(k)a(k) a 1(1d < Es EEk+1)> = < Es, Es > +ETG(k)F(k)G(k)TE, k = lEsIl2 + ZFiia2 i=1 aIl(l tillls E(k) = E(+) -ak+~E, + G(k)F(k)a(k) k k+1 + IIE| ||2 + e Fiia2 allnd FkkJ = l-< E,E(k)> -k kF - +l 1-C+ 1 + IlEl2 + (< Ek+l,,E > +Ek+lG(k)F(k)a(k)). For the case k = r, we obtain BE(r) = Er 7

or (Br+i + EsET)Er) = E and thus E(r)= Et+1 < EEr+ l > Er+1 1+ <E E, sE where Er+l = Ev + GFGTEV and the formula for Frr = 1- < E,, Er) > can be readily obtained. The above formulae have been developed for updating in a parallel computing environment, where each processor, for k > s works on generating Gk and Fkk. For computing on a single processor, using the above results, a simpler and inductive formula can be developed, i.e., (k+) is generated after k) has been generated. In case the number of non-basic i.e., E~lk+1 is generated afte E indices is larger than in, the above representation is not suitable, and we should use the formula of Section 3. The updating formula for this case is also generated in a similar manner. 8

References [1] R. H. Bartels, G. H. Golub and M. A. Saunders, "Numerical Techniques in Mathematical Programming," in Nonlinear Programming,R. B. Rosen, O. L. Mangasarian, and K. Riter, eds., Academic Press, New York, 1970, 123 - 176. [2] R. W.. Cottle, "Manifestations of Schur complement," Linear Algebra and its Applications, 8(1974), 189-211. [3] A. George and J. W. Liu, Computer solution of large sparse positive definite systems, Prentice Hall, New Jersey, (1981). [4] W. W. Hager, "Updating the inverse of a Matrix," SIAM Review 31 (1989) 221-239. [5] R. Saigal, "An infinitely summable series implementation of interior point methods", Technical Report 92 - 37, Department of Industrial and Operations Engineering, University of Michigan, May 1992, ( To appear, Math of OR ). [6] R. Saigal, "Matrix partitioning methods for interior point algorithms", Technical Report 92 - 39, Department of Industrial and Operations Engineering, University of Michigain, Julne 1992. 9