From: "Alan Miller" Subject: Re: computation of the inverse of a matrix Date: Sat, 20 Jan 2001 00:27:47 GMT Newsgroups: sci.math.num-analysis Summary: Solving linear systems robustly after rank-1 updates Laurent Martin wrote in message <873deffadd.fsf@195.132.222.185>... >Victor Eijkhout writes: > >> Laurent Martin writes: >> >> > I'm inverting a sparse symmetric definite positive matrix. >> >> Are you sure you want to do that? > >Yes I'm sure (well I think I am :-). I want to solve Ax=b but I also >need the results of (A+dA)x=b where dA are slight modifications of A >(all dA are matrices of rank 1). To do so, I use the >Morrison-Sherman-Woodbury formula which needs the inverse of A. > >Laurent. But rank 1 updates of the Cholesky factorization of A are fast and don't have ill-conditioning problems. The Morrison-Sherman-Woodbury update is slow and can lose lots of digits of accuracy. There is no comparison. Also, with luck, the update of the Cholesky factorization may still be moderately sparse. See Morven Gentleman's papers on using planar (Jacobi, Givens) rotations for updating Cholesky factorizations directly - i.e. you don't first update A. See also his Applied Statistics algorithm, I think it is AS75. My algorithm AS274 extends it. You can get AS75 from statlib: http://lib.stat.cmu.edu Get the latest version of my AS274 from my ozemail web site. -- Alan Miller, Retired Scientist (Statistician) CSIRO Mathematical & Information Sciences Alan.Miller -at- vic.cmis.csiro.au http://www.ozemail.com.au/~milleraj http://users.bigpond.net.au/amiller/ ============================================================================== From: sci.math.num-analysis@elsner.org (Ulrich Elsner) Subject: Re: computation of the inverse of a matrix Date: 20 Jan 2001 14:13:24 GMT Newsgroups: sci.math.num-analysis According to Laurent Martin : >Victor Eijkhout writes: > >> Laurent Martin writes: >> >> > I'm inverting a sparse symmetric definite positive matrix. >> >> Are you sure you want to do that? > >Yes I'm sure (well I think I am :-). I want to solve Ax=b but I also >need the results of (A+dA)x=b where dA are slight modifications of A >(all dA are matrices of rank 1). To do so, I use the >Morrison-Sherman-Woodbury formula which needs the inverse of A. Luckily you don't need the inverse. Let dA=F*G^T, with F,G \in \R^(n,1) (i.e., dA has rank 1). Then the Sherman-Morrison-Woodbury formula says: (A+dA)^{-1} = A^{-1} - A^{-1} F (I - G^T A^{-1} F)^{-1} G^T A^{-1} You only want h x = (A+dA)^{-1} * b. But this is x = A^{-1}b - (A^{-1} F (I - G^T A^{-1} F)^{-1} G^T A^{-1})b So you calculate y= A^{-1}b z= A^{-1}F w= (I - G^T z )^{-1} G^T y // Note that (I - G^T w) is of size 1 x= y - z w Hence you only need two solutions of linear systems with A. So you only calculate the Cholesky factorization once and then solving each modified system is cheap. Another possibility is to update the Cholesky factorization directly. Ulrich Elsner -- Ulrich Elsner, Fak. fuer Mathematik, TU Chemnitz, 09107 Chemnitz, Germany