ITERATIVE ALGORITHMS FOR THE INVERSION OF MATRICES ON DIGITAL COMPUTERS by ARTHUR DORIAN SHAW HARRIS B. A., U n i v e r s i t y o f B r i t i s h Columbia, 1957 A THESIS SUBMITTED I N PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF ARTS i n t h e Department o f MATHEMATICS We a c c e p t t h i s t h e s i s as c o n f o r m i n g to the required THE standard UNIVERSITY OF BRITISH COLUMBIA April, I960 In the presenting requirements of British it freely agree for Columbia, Department copying gain shall f o rr e f e r e n c e and study. I f o rextensive or publication copying of this n o t be a l l o w e d w i t h o u t A p r i l 13. I960 Columbia, make further of this by t h e Head thesis o f my I t i s understood thesis f o r financial my w r i t t e n Mathematics The University of British V a n c o u v e r 8, Canada. Date shall p u r p o s e s may b e g r a n t e d of fulfilment of the Library o r by h i s r e p r e s e n t a t i v e s . that Department i n partial I agree that permission scholarly thesis f o r an advanced degree a t t h e U n i v e r s i t y available that this permission. ii ABSTRACT A f t e r a g e n e r a l d i s c u s s i o n o f m a t r i x norms and d i g i t a l o p e r a t i o n s , m a t r i x i n v e r s i o n p r o c e d u r e s based on power s e r i e s e x p a n s i o n s a r e examined* The g e n e r a l c l a s s o f methods o f w h i c h t h e D i a g o n a l and G a u s s - S e i d e l i t e r a t i o n s a r e i l l u s t r a t i v e i s s t u d i e d i n some d e t a i l w i t h bounds f o r t h e e r r o r m a t r i x b e i n g o b t a i n e d assuming, b o t h e x a c t and d i g i t a l operations. The c o n c e p t o f t h e c o n d i t i o n o f a m a t r i x and i t s b e a r i n g on i t e r a t i v e i n v e r s i o n procedures i s looked i n t o . A s i m i l a r d e r i v a t i o n and examination i s t h e n made f o r H o t e l l i n g ' s a l g o r i t h m , H o t e l l i n g ' s i t e r a t i o n i s f u r t h e r examined w i t h a v i e w t o m o d i f y i n g Higher-order f o r m u l a e a r e o b t a i n e d and c r i t i c i z e d and a new v a r i a t i o n o f the a l g o r i t h m c a l l e d t h e Optimized on. it. H o t e l l i n g method i s d e r i v e d and commented Some schemes f o r c o n s t r u c t i n g i n i t i a l a p p r o x i m a t i o n s i n connection with H o t e l l i n g ' s i t e r a t i o n ( a n d s i m i l a r methods) a r e d i s c u s s e d and a new m o d i f i c a t i o n o f a p r o c e d u r e p r o p o s e d by B e r g e r and S a i b e l i s c o n s t r u c t e d . The f i n a l p a r t o f t h e t h e s i s d i s c u s s e s a c l a s s o f f i n i t e - s t e p i t e r a t i v e i n v e r s i o n s based on an i d e n t i t y o f H o u s e h o l d e r ' s . Three members o f t h e c l a s s , namely J o r d a n - t y p e C o m p l e t i o n , t h e Symmetric method and t h e Q u a s i optimum method a r e d e f i n e d and b r i e f l y d i s c u s s e d . The Quasi-optimum method i s t h e n examined i n f u r t h e r d e t a i l and some o f i t s p r o p e r t i e s d e r i v e d f o r t h e s p e c i a l case w i t h t h e u n i t m a t r i x f o r a n i n i t i a l approximation. I hereby c e r t i f y t h a t t h i s a b s t r a c t i s s a t i s f a c t o r y . iii TABLE OF CONTENTS Page CHAPTER I INTRODUCTION CHAPTER I I MATRIX AND VECTOR NORMS CHAPTER I I I 1 1. D e f i n i t i o n o f a v e c t o r norm 3 2. Some u s e f u l norms and t h e i r p r o p e r t i e s ....... 3 3e N a t u r a l m a t r i x norms; some p r o p e r t i e s . 4 4. I n e q u a l i t i e s i n v o l v i n g n a t u r a l norms ......... 6 5. D i g i t a l o p e r a t i o n s ; d e f i n i t i o n s and r e l a t i o n s 7 6. An i l l u s t r a t i o n o f t h e u s e o f d i g i t a l c o n c e p t s 9 ITERATIVE METHODS BASED ON MATRIX POWER SERIES EXPANSIONS ..o. 1. P r e l i m i n a r y remarks 2. Total-step (Diagonal) i t e r a t i o n 14 3. Single-step (Gauss-Seidel) i t e r a t i o n 17 4. The g e n e r a l i t e r a t i o n ; some p r o p e r t i e s 19 5. S t r i c t bounds f o r t h e e r r o r m a t r i x 22 6. Remarks on t h e c o n d i t i o n o f m a t r i c e s ........... 14 ( a ) Bounds f o r t h e norm o f t h e change i n t h e i n v e r s e r e s u l t i n g f r o m a change i n the matrix 26 (b) C o n d i t i o n numbers o f a m a t r i x ........ 28 ( c ) Remarks on an a s s u m p t i o n made i n o b t a i n i n g a bound f o r t h e e r r o r m a t r i x i n (5) 29 7. Speed o f convergence 8. The s i z e o f t h e r e s i d u a l f o r t h e i t e r a t i o n d e f i n e d i n (4) u s i n g d i g i t a l c o n c e p t s .... 29 30 iv Page CHAPTER I I I (Continued) 9. The error i n solving a system of l i n e a r equations by i n v e r t i n g the c o e f f i c i e n t matrix using d i g i t a l concepts ............. 10 » A remark on double precision procedures as applied to the i t e r a t i o n i n (4) 11. 35 • ••• 37 Hotelling's i t e r a t i o n 39 (a) Derivation from Euler*s i d e n t i t y * 42 (c) S t r i c t bounds f o r the error matrix ... 42 (b) The r e s i d u a l (d) Analysis of the r e s i d u a l i n terms of d i g i t a l concepts............... (e) Remarks on double p r e c i s i o n e 47 procedures CHAPTER IV 44 HOTELLING'S ALGORITHM: MODIFICATIONS AND INITIAL APPROXIMATIONS 1. Three possible approaches to the algorithm ..... 4# 2. Related formulae suggested by R^ 50 3. The Optimized Hotelling algorithm 4» Suggestions f o r i n i t i a l approximations +1 = R^ • (a) Various simple devices 56 (b) Berger and Saibel's method 57 (c) An improved method CHAPTER V 52 RANK-REDUCTION METHODS ~ •••• 57 A CLASS OF FINITE-STEP ITERATIVE INVERSION PROCEDURES 1. The basic Rank-reduction formula 62 ' V Page CHAPTER V (Continued) 2 0 Jordan-type Completion ( a ) Development of. f o r m u l a ••* •••*• 63 (b) R e l a t i o n t o J o r d a n i n v e r s i o n 64 (diagonalization) 3. 4. Symmetric Method (a) Banachiewicz f a c t o r i z a t i o n ......... 66 (b) Development o f Symmetric f o r m u l a ... 67 ( c ) H. S. W i l f ' s p r o c e d u r e 69 Quasi-optimum Method • (a) H e u r i s t i c c o n s i d e r a t i o n s 70 73 (b) Development o f f o r m u l a e ( c ) S p e c i a l c a s e : columns o f A - S Q as • the b a s i c s e t (d) S p e c i a l c a s e : t h e u n i t v e c t o r s the b a s i c s e t 74 as 75 (e) Symmetrical p r o p e r t i e s f o r t h e case C Q = I 78 ( f ) O p e r a t i o n count f o r Quasi-optimum 79 algorithms 5. A p r o p e r t y o f t h e r e s i d u a l f o r t h e Symmetric method BIBLIOGRAPHY • 80 82 vi ACKNOWIEDGMENTS The a u t h o r w i s h e s t o e x p r e s s h i s a p p r e c i a t i o n o f t h e a s s i s t a n c e r e c e i v e d from Dr. T. E. H u l l a n d Dr. B. N. M o y l s i n c o n n e c t i o n w i t h this thesis. E s p e c i a l t h a n k s a r e due t o D r . H u l l f o r h i s c o n s t a n t i n t e r e s t and encouragement. Thanks a r e a l s o due t h e Computing C e n t r e o f t h e U n i v e r s i t y o f B r i t i s h Columbia f o r t h e u s e o f t h e i r computing f a c i l i t i e s . Much o f t h e m a t e r i a l o f t h e t h e s i s was s u g g e s t e d b y r e s u l t s o b t a i n e d by t e s t i n g v a r i o u s a l g o r i t h m s and c o d i n g schemes. The adequate evaluation and c h e c k i n g o f v a r i o u s schemes s u g g e s t e d i n t h e t h e s i s would have been v i r t u a l l y i m p o s s i b l e w i t h o u t t h e u s e o f t h e ALWAC I I I - E computer. F i n a l l y , t h e a u t h o r w i s h e s t o acknowledge h i s i n d e b t e d n e s s t o Mr. H. Dempster o f t h e Computing C e n t r e f o r numerous s u g g e s t i o n s , comments and d i s c u s s i o n s i n c o n n e c t i o n w i t h t o p i c s o f t h e t h e s i s and t o Mrs. D. H a m i l t o n f o r h e r u n f a i l i n g c a r e and p a t i e n c e i n c o n n e c t i o n w i t h t h e t y p i n g and p r e p a r a t i o n o f t h e t h e s i s . 1. CHAPTER I Introduction I n r e c e n t y e a r s s e v e r a l s u r v e y s o f known methods f o r s o l v i n g l i n e a r equations and i n v e r t i n g m a t r i c e s have been made by v a r i o u s Perhaps the most complete of t h e s e i s t h e one authors. c o m p i l e d by G. E. Forsythe e n t i t l e d "A T e n t a t i v e C l a s s i f i c a t i o n of Methods and B i b l i o g r a p h y S o l v i n g Systems o f L i n e a r ; E q u a t i o n s " Mathematics S e r i e s No. o f about 450 and may titles 29. — p u b l i s h e d i n N. B. S. This l i s t i n g — Applied which includes a b i b l i o g r a p h y i s based on F o r s y t h e ' s be t a k e n as e x h a u s t i v e on own system o f c l a s s i f i c a t i o n as o f t h e t i m e of c o m p i l a t i o n (1951). The p r e s e n t p a p e r , on t h e o t h e r hand, makes no attempt t o c o v e r a v e r y wide range o f i t e r a t i v e i n v e r s i o n methods. dealt with: F i v e main p o i n t s a r e a u n i f i e d treatment of a c l a s s of i t e r a t i v e i n v e r s i o n s w h i c h i n c l u d e some w i d e l y u s e d methods, an a n a l y s i s of t h e limiting e r r o r i n h e r e n t i n such methods i n terms o f d i g i t a l o p e r a t i o n s ( i n v o l v i n g t h e concept o f r o u n d - o f f ) , t h e development of t h e O p t i m i z e d H o t e l l i n g i t e r a t i o n ( a new v a r i a n t of t h e f a m i l i a r a l g o r i t h m t h a t o b v i a t e s need f o r a good i n i t i a l a p p r o x i m a t i o n ) and a l s o o f a new the scheme f o r c o n s t r u c t i n g i n i t i a l a p p r o x i m a t i o n s t h a t p r o v i d e s a good s t a r t i n g p o i n t for a c l a s s of matrices and, f i n a l l y , t h e d e f i n i t i o n and e x a m i n a t i o n o f a c l a s s o f f i n i t e - s t e p i t e r a t i v e i n v e r s i o n p r o c e d u r e s some -of w h i c h are a l t e r n a t i v e f o r m u l a t i o n s o f known methods and some o f w h i c h have n o v e l f e a t u r e s t o recommend them. While developing t h e main p o i n t s l i s t e d , v a r i o u s r e l a t i o n s h i p s p e r t a i n i n g t o i n v e r s i o n p r o c e d u r e s i n g e n e r a l are i n t r o d u c e d . t h e s e a r e not new Many o f b u t , b e i n g , f o r the most p a r t , s c a t t e r e d t h r o u g h o u t t h e l i t e r a t u r e , are p r o f i t a b l y g a t h e r e d t o g e t h e r i n a paper such as the 2. p r e s e n t one w h i c h attempts t o p r e s e n t a c o m p u t e r - o r i e n t e d approach throughout * A l t h o u g h t h e r e i s no emphasis h e r e on t h e comparative e v a l u a t i o n of s p e c i f i c a l g o r i t h m s o r on t h e r e l a t i v e m e r i t s o f d i r e c t i n v e r s i o n schemes as opposed t o i t e r a t i v e p r o c e d u r e s , be drawn f r o m t h e t e s t . c e r t a i n c o n c l u s i o n s may reasonably E x c e p t i n t h e case o f " s p e c i a l " m a t r i c e s ( e . g . , s p a r s e o r v e r y s t r o n g l y d i a g o n a l m a t r i c e s ) t h e r e seems t o be r e a s o n f o r n o t c h o o s i n g H o t e l l i n g ' s i t e r a t i o n over both a l g o r i t h m s and t h i r d and h i g h e r - o r d e r f o r m u l a e . no first-order Provided t h a t storage c o n s i d e r a t i o n s and t h e number o f o p e r a t i o n s a r e not o f o v e r r i d i n g importance, t h e O p t i m i z e d H o t e l l i n g a l g o r i t h m would seem t o be an coded and w i d e l y a p p l i c a b l e i n v e r s i o n p r o c e s s . easily- T h i s method h a s , moreover, s t o o d up w e l l under n u m e r i c a l t e s t i n g o f i l l - c o n d i t i o n e d m a t r i c e s on t h e computer* I n a l l c a s e s , t h e use o f d o u b l e - l e n g t h a c c u m u l a t i o n o f i n n e r products before rounding to sing]e-length i s h i g h l y d e s i r a b l e . I f , however, t h e r e i s some r e a s o n t o b e l i e v e t h a t t h e m a t r i c e s t o be d e a l t w i t h a r e r e a s o n a b l y w e l l - c o n d i t i o n e d , t h e r e w o u l d seem t o be no r e a s o n f o r not s e l e c t i n g a f i n i t e - s t e p procedure i t e r a t i o n s of the H o t e l l i n g t y p e . i n preference t o The economy i n t h e number o f o p e r a t i o n s a c h i e v e d i n t h e f i n i t e - s t e p methods may w e l l be a d e c i d i n g f a c t o r . Add t o t h i s the f a c t t h a t warnings of p o s s i b l e i l l - c o n d i t i o n i n the m a t r i x b e i n g i n v e r t e d and p i v o t - s e a r c h t e c h n i q u e s a r e a v a i l a b l e f o r s e v e r a l o f t h e s e methods and t h e case f o r u s i n g such methods seems f a i r l y s t r o n g . Arguments f o r u s i n g , i n p a r t i c u l a r , t h e Quasi-optimum f i n i t e - s t e p i n v e r s i o n i n p r e f e r e n c e t o one o f t h e f a m i l i a r d i r e c t i n v e r s i o n s i n c l u d e c o n s i d e r a t i o n s o f symmetry, s u i t a b i l i t y f o r d o u b l e - p r e c i s i o n a r i t h m e t i c , and amenability to a d i g i t a l error analysis. 3 CHAPTER I I Matrix and Vector Norms In any discussion of i t e r a t i v e processes involving matrices, conciseness of treatment i s often aided by the use of vector and matrix norms. For t h i s reason, a brief statement of pertinent norm properties and associated notation follows. I n t u i t i v e l y described as a ''measure of size", a norm on a linear vector space i s defined as a real-valued function on the space S of elements x written as [| x || or N(x) and satisfying the following requirements: (2.1.1) for a l l x I S > 0 (2.1.2) x I - 0<^x = 0 for a l l x £ S for a l l scalars (2.1.3) (2.1.4) x + 7 M ll ll x for a l l x, y £ S (triangle inequality) l|7 + X Since, i f /Jx// i s a norm, ||x//' = CL6X # ||Ax||also defines a norm for any fixed non-singular matrix A, i t i s easy to construct vector norms which do not give equal weight t o a l l the components of the vector. A class of vector norms having wide application i s , however, symmetrical i n the foregoing respect. This i s the class of norms defined by: n llxll fc (2.2) def 1, 2, 3, n i=l Minkowsky's inequality ensures that these norms satisfy the triangle inequality} the remaining properties of (2.1) are t r i v i a l l y s a t i s f i e d . Of t h i s class of norms, three i n particular are easily computed, namely f x •I2 -t/2r i = x 8X1(1 l i m i t i n g member of the class, which i s obtainable by inspection. Max • over i |xJL x 4. It i s not always necessary to specify a particular norm i n the context since i t i s well known that a l l vector norms are equivalent i n the following sense: for any two norms (£ and ^3 [[xfj* and |(xjj" there exist positive numbers (which depend only on the norms) such that Thus, for example, $ (2.3.1) l\x\\ (2.3.2) llx/l^^ flxll! < ||x|| ^ y^" 2 llxll n HXW^ Since "equivalence" as used above i s a true equivalence relation, the property mentioned i s most conveniently proved by showing that an arbitrary norm i s equivalent, i n particular, to || x f j ^ . Of the various possible ways of defining a matrix norm, the following, which defines the matrix norm i n terms of a vector norm, has useful theoretica l properties: (2.4) Ax 11A || -over Sup.a l l ||x// x, x 4 0 The norm symbol [ || |/ or N( )] applied i n the same context to both matrices and vectors w i l l here be understood to refer to associated norms. That i s to say, the matrix norms w i l l be "natural norms", (to use the terminology of F. JOHN [Ref. 1]) associated with (or "subordinate to") their respective vector norms by definition (2.4). Using this notation, one can establish the following basic relations for natural matrix norms: 5. (2.5.1) | | l | =1 (2.5.2) || Ax|| < II A||' || x l| ( f o r matrices A and vectors x) (2.5.3) II ABIj ^ llAll • II Bll ( f o r matrices A, B) (2.5.4) l|A + Bll i (where I i s the unit matrix) + || Bll llAll It i s c l e a r from definition (2.4) that i f yx. i s the (numerically) <) largest eigenvalue of A, l^ul constitutes a lower bound f o r a l l natural norms of A. A and any One can, however, also make the statement that f o r any matrix € > 0, there e x i s t s a n a t u r a l norm, IIAII , s a t i s f y i n g the inequality tuU IUII 4 \jx\ (2.6) 6 + Proofs of t h i s u s e f u l r e s u l t use f i r s t the f a c t that an a r b i t r a r y A can be t r i a n g u l a r i z e d by pre- and p o s t - m u l t i p l i c a t i o n by a non-singular matrix t o obtain T = B AB" -, 1 pre- and p o s t - m u l t i p l i c a t i o n by A further s i m i l a r i t y operation (namely D ~ l a Diag ( ) and D respectively, where € i s s u f f i c i e n t l y small) then reduces a l l off-diagonal elements of the t r i a n g u l a r matrix T t o a r b i t r a r i l y small quantities, i . e . , ; P » D"! TD = y \ + E, where J\. i s the diagonal matrix of eigenvalues and a l l elements of E are small. that llxll \\(T~ BxH 2 of A F i n a l l y , using the f a c t , mentioned before, defines a v a l i d norm as long as G'^B i s non- singular (which i n t h i s case i t i s ) , one can show that the natural matrix norm based on the vector, norm just defined has the required properties. For a general matrix A, the natural norm square norm llxl^ s a t i s f i e s the r e l a t i o n z (2.7) lUH^ . = max \ \ \ ± II AD2 based on the vector 6 where \^ i s an e i g e n v a l u e o f AA* ( t h e s t a r d e n o t i n g t h e t r a n s p o s e d complex conjugate). T h i s f o l l o w s e a s i l y from t h e d e f i n i t i o n (2.4) u s i n g t h e f a c t t h a t t h e non-negative H • U* A u H e r a d t i a n m a t r i x H = AA* c a n be e x p r e s s e d i n t h e form. where U i s u n i t a r y and A . i s t h e d i a g o n a l m a t r i x o f e i g e n v a l u e s ( n e c e s s a r i l y r e a l and non-negative) since the eigenvalues w i l l Biave from Hermitian, o f AA* w i l l t h e n be the squares o f t h o s e o f A, we II A L = max IjijJ' i = 1, 2, t I f , further, A i s i t s e l f (2.7) (2.S) where ja. ± o f H. n, are t h e e i g e n v a l u e s o f A. m a t r i c e s one c a n a c t u a l l y t a k e €» Thus f o r H e r m i t i a n 0 i n (2,6). S i n c e t h e e i g e n v a l u e s o f an a r b i t r a r y A a r e n o t r e a d i l y a v a i l a b l e , t h e n a t u r a l norm i s o f t e n more c o n v e n i e n t t o u s e t h a n IIAII H A I ^ . I t c a n be seen from (2.4) t h a t one has " A II (2.9) i.e., H A I I ^ » max. W l £ 1 I i s t h e " l a r g e s t a b s o l u t e row sum". One can a l s o make u s e o f G e r s h g o r i n ' s i n e q u a l i t i e s which r e s t r i c t t h e e i g e n v a l u e s o f A t o l i e w i t h i n t h e u n i o n o f a l l c i r c l e s i n t h e complex p l a n e w i t h c e n t r e s a ^ n and r a d i i r^• £1 j»l l iil a • F o r e x ampl * e l a r g e s t e i g e n v a l u e o f AA*, i t f o l l o w s from T r AA* - Zl Uil i,j (2.10) 2 i f ^max n ^ m i a x s t l a e } numerically T r AA* ^ and that: 3 ( L |aJ I « lUlL 2 V* i , j ^ ) 2 £ ( TL i , j laj ) * . 2 ^ J 7. A f u r t h e r i n e q u a l i t y can be obtained by applying an extremum property of IIAII2 • I t can be shown that |1AI| * (2.11) 2 ^ „ L* 4yll 2 \ ±j\ =|e| Ae^j where Since a » (0, . . . , 0, 1 , 0, . . . , 0) (2.12) 2 e^ i s the k^* unit vector 1 one has immediately from t h i s and from 1 (a^U Ulo ^ (jL kJ )^ ^ n max 2 over x , j a J * J X max overi,j (2.10) | a, -1 1 3 The foregoing properties and relations involving matrix and vector norms have t o be modified when operations with truncated or rounded numbers are introduced. In c o n t r a - d i s t i n c t i o n t o "exact operations on exact numbers" ( i . e . , the ordinary processes of arithmetic) we w i l l define "pseudo-operations on d i g i t a l numbers" as arithmetic operations (and processes based on these) on numbers rounded t o a f i n i t e number of places (assuming a/2-based number system) and r e s u l t i n g i n a s i m i l a r such number by a process of truncation with or without rounding. More restrictive d e f i n i t i o n s of d i g i t a l numbers are often used but w i l l be introduced here only as required. Following (with minor v a r i a t i o n s ) the notation used by J . von NEUMANN and H. H. GOLDSTINE [Ref. 2 ] , "a" w i l l represent a d i g i t a l number. More s p e c i f i c a l l y , i t w i l l represent the d i g i t a l number a r i s i n g from a number "a" by any process of truncation implied by the context. S i m i l a r l y , v and A w i l l denote the " d i g i t a l vector" and " d i g i t a l analogously obtained. written as Q , (+), @, The four arithmetic pseudo-operations (x). F i n a l l y , matrix" w i l l be the pseudo-product of two matrices 8. w i l l be represented as i n "AoB" with an inner product always being expressed i n the form x^ o y . If one assumed that addition were exact (as would be the case i f a l l numbers were rounded to s places but were unrestricted i n magnitude), then for the natural norm II II2 one would have the following inequalities as (2.5.2) and (2.5.3): modifications of (2.13.1) (2.13.2) where £ = ~ — II IoxH 4 2 I1AH 2 Hxll + n^£ 2 II AOBII ^ IIAI 'IIBH 2 o r 2 ~ S + n € 2 2 according as the numbers are rounded to s places or are truncated without rounding. These inequalities follow by use of since i f a i s the largest element of either a - b I v n £ where b i s the corresponding (2.12), Aox or AoB, one has element of Ax or AB, respective- iy. Relationships (2.13) are relevant to fixed-point operation. I f , however, floating-point scaling i s used, the operations + and (+) are no longer equivalent and further round-off error i s introduced. Suppose that the d i g i t a l quantities considered have the form ± / 3 m.lawhere p i s an integral power of the b a s e , ^ , and m i s a fraction rounded to s - ^ j - < m < 1. places and lying i n the interval exponent i n the sum A + B . One obtains: (2.14.1) II A Let q be the largest / 0Bll ^ lllll + HBI1 2 2 where €. i s defined as before. plication: i s derived next. 2 The inequality for (square) matrix multi- The largest round-off error i n an element of 9. the product i s bounded as follows: where i f the max element of A«B - A B i s g over i , j then Igf q < (2n - l ) £ s / q i s the floating point exponent of the largest element i n A B . xi <<p This statement w i l l hold provided here. which w i l l always be assumed B The factor (2n - l ) i n the above i s composed of from multiplication of elements and (n - l ) & y 3 additions necessary to- form an inner product. I|AOBI| 2 - llloB - A B + ^ llloB - A B I 2 + A B| III q n & fi^ resulting from the (n - l ) One then has 2 ll| 2 and by (2.12) this gives rise to (2.14.2) I A . B I 4 Illll2 2 IIB1 + n(2n - l ) £ 2 £ q and similarly f o r matrix-vector multiplication: (2.14.3) l|A-x|| ^ 2 HAII 2 l\xll + Vn (2n - l ) € ^ 2 , A theorem - basic for the types of iterative process examined i n what follows - states that a necessary and sufficient condition for the limit of Ak to be the zero matrix i s that the eigenvalues of A be a l l less than unity i n absolute value. This result, which i s awkward to establish for a general matrix A unless norms are used, w i l l be proved i n detail to show a typical application of the theorem involving relationu(2i6).(The result w i l l then be modified to incorporate the use 10. of pseudo-operations and d i g i t a l quantities. X of the stated condition, l e t of of A and l e t | X j j ^ l ^.1. a (natural) norm of A \nJ^ To e s t a b l i s h the s u f f i c i e n c y be the numerically largest eigenvalue Inequality (2.6) then assures the existence such that M 1. < The norm properties then immediately give IIA^ Hence A^ approaches, the zero matrix i n t h i s norm and - by the equivalence of norms - i n any other norm a l s o . vector corresponding to A k i d ||x'|| . m a x ^ For the necessity, l e t x* be the eigen- Siince we have H XM. ||x'l| Kki 9 I \k x I * I - IAl i t follows that f o r any natural norm: over x^O but, by hypothesis, I I A H I —> ||x|| 1 » y 0, whence I Xl < 1. The simple c r i t e r i o n f o r convergence to zero established above i s , unfortunately, not d i r e c t l y applicable when one i s confronted with a d i g i t a l matrix A obtained as a r e s u l t of previous computations matrix with accompanying round-off. or by input of a Even i f IXjug^l < 1 i s known to hold f o r A, one does not, i n general, have such assurance for the condition that the l a r g e s t eigenvalue of Th e Furthermore, A be numerically less than unity would not ensure that the i t e r a t e d pseudo-product would approach the zero matrix. A. pseudo-product of A with i t s e l f 5? Io ••• «A (with lie k factors) i s , i n f a c t , not even defined, since matrix pseudo-multiplication i s not a s s o c i a t i v e . A ® w i l l be used t o denote In spite of t h i s , the symbol such a pseudo-product with the understanding that the symbol stands f o r any k—2 one member of the set o f 2 be formed from the " k possible "pseudo-powers" that can i n general factors, I t i s of i n t e r e s t to determine a good upper bound for We make the assumption that (2.15.1) JlAll ^ A i s such that , Y - n(2n - 1) G where y\ - Y Lim IIA k-»~> » n(2n - 1) Inequalities (2.15.1) and (2.12) imply I max A < 1 and hence lover i , j the largest exponent of any element of A i s l e s s than or equal t o zero (assuming From (2.15.1) and (2.5.3) we also have ^> l ) . I I A 2 I 2 * I I A I | ^ «I. Inequality (2.14.2) can therefore be applied with q • 0 since i t has just been shown that a c t u a l l y does not i n v a l i d a t e the i n e q u a l i t y . q ^ 0 and increasing q Thus we have by (2.14.2): l|A®t,*iA-|.n U ( 2 „-l) £ (by 2.15.1) . 12. Proceeding by induction one can show that for a k =» 1, 2, 3, • 'over i , j n. Having established this, one can repeatedly q - 0 . This gives the sequence of PPly (2.14.2), always taking inequalities and i n general: Summing the series gives (2.15.2) «A©II S I j j k 2 from which, since || AII Lim k-^oo + n ( 2 n . l ) £ 1-IUlL 2 2 (2.15.3) 1 2 1, one obtains | l ® L <C n ( 2 n - l ) f> 2( 1 - || A || 2 ) 13. As immediate corollaries to the foregoing one can make the statements: 1 (2.15.4) If n max. over x S , 1 A J then II A ^ l l ^ ( l - n(2n-l) A-3 * 2 ) ^ ll A II f o r a l l k - 1, 2, 3, where A = (a, ,) ^-J n. The above condition being met, the largest element of A ® , for sufficientl y large (2.15.5) k, w i l l be bounded above by I max A ® | <: »<f " Hf* I over i , J ' ^ 2(1 - | I | L ) where "sufficiently large" i s conservatively indicated by 1 (2.15.6) k > 1 l o g * (~ max . 7 T aU A) P n over i v-^ ~ J' n a which insures that I] A | | ^ Z 3 " 3 • 14 CHAPTER I I I I t e r a t i v e Methods Based on Matrix Power-Series Expansions Although the problems of matrix inversion and of solution of systems of l i n e a r equations are c l o s e l y related, the l a t t e r has received more e x p l i c i t attention over the years. The reason f o r t h i s i s probably that the majority of investigators (many of whom were not, primarily, mathematicians) were confronted with some physical problem necessitating the s o l u t i o n of a system of l i n e a r equations of n o n - t r i v i a l order and, consequently, sought a s o l u t i o n i n the most d i r e c t manner. As an additional consequence of t h i s predominantly pragmatic history, much of the material that has accumulated on the subject concerns i t s e l f with a great v a r i e t y of s p e c i f i c algorithms based on a r e l a t i v e l y small number of d i s t i n c t approaches t o the problem. The more sophisticated methods (of which the Conjugate Gradient f i n i t e - s t e p i t e r a t i o n [Refs. 3, 4 , 53 i s representative) are, f o r the most part, quite recent developments stemming, more or l e s s d i r e c t l y , from the increased i n t e r e s t i n t h i s f i e l d aroused by current s c i e n t i f i c and i n d u s t r i a l trends and from the coining of age of high-speed, general-purpose computing devices. With the possible exception of some schemes f a l l i n g i n the last-named category, most equation-solving procedures are t r i v i a l l y adaptable to matrix inversion. However, some methods are e s s e n t i a l l y better suited f o r such adaptation than others. Two of the oldest, simplest, and most f a m i l i a r i t e r a t i v e methods f o r solving l i n e a r equations have immediate applications t o the matrix i n v e r s i o n problem. Since they both give rise t o inversion schemes which are members of the general class of inversion methods which w i l l be examined next, t h e i r nature w i l l be b r i e f l y r e c a l l e d here. The f i r s t 15. method i s variously referred to as 'simple', 'ordinary', 'block', 'total-step 1 or 'diagonal' iteration. an approximation, x ^ , .to the vector obtained by dividing each element element of (3.1) x i n the equation Ax = b is b by the corresponding diagonal The i ' t h element of A. Tv solution, x ^ , In the f i r s t step of the process, of the next approximate x ^ i s then obtained by (1) xj x ; ^ (0) n = - ± - [b, - a., x, ] i - 1, 2, Successive steps of the iteration are analogous to the one for n x^. The use of matrix notation, besides lending conciseness, brings out general features which are often obscured by verbal descriptions and relations given i n terms of matrix elements. Thus, for the process Just given, one has the following development. (3.2) Consider the equation (3.3) Let (3.4) whence x .- (D - W)" b (3.5) or x = (I - D " W) A= D- W Ax = b where D = Diag(a ), ii A= 1 1 _1 D" b 1 The iterative process given by (3.1) can be expressed as (3.6) x( k+1 > = D- (b + Wx X (k) ) (a^) 16. which, with (3.7) X Q= D x ( k + l ) b, gives r i s e by induction t o the series expansion = [ I + ( D ^ M D - V + ... + (D" W) ]D" b + ( D W ) 1 k 1 _1 k+:L XQ I t e r a t i o n (3.6) i s thus, e s s e n t i a l l y , an algorithm f o r evaluating successive p a r t i a l sums of the series expansion f o r ( I - D**%)"^" i n (3.5). The algorithm f o r f i n d i n g these sums i s not unique, but before describing other approaches, a few comments should be made on the foregoing method. The adaptation t o matrix inversion follows t r i v i a l l y from the matrix equation AC = I (where C = A"*"*") which replaces (3.2). the same l i n e s as before, but w r i t i n g C for x and I Following f o r b, one obtains: (3.3) C k + 1 - irHl + WC ) k It i s apparent from (3.7) that a necessary condition f o r convergence of the i t e r a t i o n defined by (3«6) i s that fX „ a 1 the absolutely largest eigenvalue of D~^¥. < 1 X where IucUC ' is luaX I t follows that the method i s suitable f o r matrices A with a numerically dominant diagonal, D. Provided A i s p o s i t i v e d e f i n i t e , the domination w i l l c e r t a i n l y be strong enough f o r convergence i f | cL^ I > £1 j=l with D = Diag(d^) f o r i = i 2, ... n ' ' *' and W = («ij)j since i n t h i s event one has II D^wllg $ II D^wllj^ < 1. D jwjj | Nothing i s e s s e n t i a l l y a l t e r e d by replacing by the diagonal matrix obtained a f t e r permuting A so as to place as many numerically dominant elements as possible on the main diagonal,,, I f the process described by (3»l) i s used, the permutation need not 17. a c t u a l l y be c a r r i e d out. Two concluding remarks concern the computational aspects of the algorithm. F i r s t , an i t e r a t i v e procedure such as (3.6) i s preferable t o a straightforward formation and summation o f successively higher powers o f the matrix since the i t e r a t i v e procedure i s s e l f correcting f o r minor errors ( i . e . , errors that do not increase the relevant matrix norm beyond permissible bounds). Second, the accuracy of the inverse obtained i n the l i m i t by d i g i t a l operations o f the kind defined e a r l i e r , i s of major concern to the numerical analyst. terminal accuracy w i l l be dependent, i n general, upon the This particular algorithm considered, the choice of the sequence of operations w i t h i n any one algorithm, the mode o f computation used throughout or at d i f f e r e n t 3tages ( i . e . , fixed-point, double-precision, e t c . ) , the choice of the i n i t i a l approximation, CQ, and, f i n a l l y , upon the nature of the matrix A being inverted. The f a c t o r s l i s t e d w i l l , however, enter i n with varying 'weight* and the r e l a t i v e importances w i l l vary appreciably with the algorithm under consideration. The second example of i t e r a t i v e schemes of the type under considera- t i o n i s u s u a l l y named the Seidel o r Gauss-Seidel method. Since, as OstrowskL and Forsythe have mentioned, Gauss never r e f e r r e d t o the method and S e i d e l advised against i t i . : use, i t i s perhaps better r e f e r r e d to as the single-step (as opposed t o the t o t a l - s t e p ) i t e r a t i o n . As the name suggests, the single-step procedure consists of applying a t y p i c a l step o f the t o t a l - s t e p method (as given by (3.1)) t o each equation i n t u r n ( i n a c y c l i c manner) using the most recently a v a i l a b l e vector elements i n the approximating vector at each stage. The matrix expression f o r -this i t e r a t i v e procedure i s therefore obtained as f o l l o w s . 18 (3.9) where Let A = L - U L i s the "lower t r i a n g l e of A", including the main diagonal. As i n (3.5) we have (3.10) x = (I - L ^ U T 1 L _ 1 b The equation giving the r e l a t i o n s h i p between the relevant matrices and vectors at the end of the (k + l ) s t cycle o f n (3.11) Lx(k+D - ux(k) „ steps i s b which i s o f the same form as (3.6) namely (3.12) x< k + 1 ) = L" (b 1 + Ux^) Once again, the algorithm obtains the series expansion of ( I - IT^U) i n the form: (3.13) x ( k + l ) « [ I - ( L U ) • ... + _1 (iThjfl IT * 1 + (L^U)^ The convergence condition therefore, i s that the absolutely eigenvalue o f L~~TJ be l e s s than one i n absolute value. converges i f A XQ largest The series i s p o s i t i v e d e f i n i t e (PIZZETTI [Ref. 6]) and i n f a c t , f o r a l l r e a l , symmetrical ness o f A 1 A with a l l a ^ > 0, the positive d e f i n i t e - i s both necessary and s u f f i c i e n t f o r convergence (REICH [Ref. 7 ] ) , (OSTROWSKE [Ref. 8 ] ) . I f , instead of taking the equations i n a c y c l i c order, one selects at each step the one that would give a reduction of the largest element of the residual vector, b - A x ^ , one no longer has an expansion of the form (3.13) and the method then merges into the 19. more f l e x i b l e i t e r a t i o n schemes known as Relaxation methods i n which elements of the r e s i d u a l vector are reduced either s i n g l y or i n groups by judicious ad hoc changes i n the approximation vector based on experience and aided by a "relaxation t a b l e " of changes and e f f e c t s . The i n v e r s i o n algorithm derived from (3.12) i s of course = I T +• I T U C (3.14) 1 1 k which i s formally the same as ( 3 . 8 ) . Having given two examples, the general algorithm of form (3.15) C k + 1 -P+Q C k w i l l be examined, the corresponding equation-solving algorithm - which w i l l be returned t o l a t e r - being (3.16) x - Pb + Qx<k) ( k + l ) The i t e r a t i o n defined by (3.15) i s of the kind c a l l e d "stationary" by several authors such as FORSYTHE [Ref. 12] and HOUSEHOLDER [Ref. The term r e f e r s here to the fact that and not functions of k. we can define (3.17) C We s h a l l assume S = P""^ and SQ = T = S' + S 1 k + 1 P _ 1 TC and P Q C HE A""*" one has C = S" and (S - T) C = I 1 + S" 1 T C are f i x e d matrices i s non-singular so that i n order t o write (3.15) i n the form K Since the algorithm should y i e l d an i d e n t i t y f o r where 13]. = C k = C 20. whence (3.18) A = S-T Thus the algorithm (3.15) with non-singular P gives r i s e t o the series expression (3.19) C K + 1 - [ I + ( S ^ T ) + ... + ( S - ^ k j S " 1 + (S^T)^ 1 C Q and i s therefore, as previously i l l u s t r a t e d , a device f o r obtaining p a r t i a l sums of the expansion o f ( I - S'^T)"^ (3.20) C = A - 1 = (I - S " ? ) 1 S - 1 i n the equation _ 1 The k'th l e f t r e s i d u a l of an i t e r a t i v e inversion process w i l l be defined as R ^ ) ^ ^1 - G^A and the right residual by The d i s t i n g u i s h i n g subscritps r or & choice i s c l e a r from the context. Ejj = A" - = RjjA . 1 (3.21.1) -1 R k w We also define the k'th error matrix as i t h t h i s notation we have f o r (3.17): = I - C A k 1 1 k 1 - S ^ T C l - C _ (S - T ) ] K : • *** \-i whence by induction: R k ( k ) r dif. * ~ ^^k* may be dropped i f the required = i - ( s " + sr ? c _ ) ( s - T) (3.22) R = ( S ^ T ^ d - CQA) 21. Taking norms, we have therefore ( f o r any norms) (3.23) II ^ II ^ H s ^ f l l - C A|| 0 jX which, once again, shows that of m a x l < 1 where X i s an eigenvalue S~^T, i s a s u f f i c i e n t condition f o r convergence. shown by (3.19)]. chose (3.25) C Q « S" 1 S-h Since S i s by assumption non-singular, we can A = C Q ^ - T, so that [The necessity i s We can then say - C - f C J - A ) «. I - C Q A 1 0 so that (3.17) assumes the form (3.26) C k + 1 = Co + ( I - C A ) Q and (3.23) becomes (3.27) K I U I I I - C O A I I ^ 1 From t h i s point of view, the matrix P i n (3.15) i s most naturally interpreted as being a " s u f f i c i e n t l y good" i n i t i a l approximation t o A ~ ^ » I f i n p a r t i c u l a r one takes (3.28) f o r summing A" 1 CQ = I one has an algorithm » I + B + B + ... 2 where B = I -A. This series - the Neumann expansion f o r A""^ - i s frequently used f o r the Leontieff matrices of economics. These matrices have the property that the off-diagonal elements of each row t o t a l l e s s than one and the diagonal elements are a l l equal t o one. 22 Apart from the question of obtaining a s u f f i c i e n t l y good points are of i n t e r e s t i n connection with (3.26): and l i m i t i n g accuracy of the approximations. C Q , two speed o f convergence The former w i l l be influenced by the d i g i t a l nature of the actual operations used and the l a t t e r , of course, arises s o l e l y on that account. Postponing, temporarily, an analysis i n terms of d i g i t a l operations, bounds f o r the e r r o r matrix at the k'th step w i l l be derived assuming exact operations. The notation Max I B I w i l l designate the absolutely largest element of a matrix B , max over i , j hjj[. To determine bounds f o r - \[\ A- ! - | | ( I - C one may say: 1 2 (3.29.1) and since 0 A ) K ( A - - C0)l2 1 A~ = E^. + Ci. t h i s gives IIBgl ^ 2 1(1- C A) Q k \ + (I - C A) (C k 0 k - C )ll 0 2 whence (3.29.2) II|| 2 ^ II I - CQA|| 2^11 2 + )| ( I - C Q A ) ( C but the l a s t term i s , from (3.19) and (3.25): K - Cq)\\ 2 {?.V).l) ||(I - C A ) [ { l + ( I - C A) + ... • ( I k 0 » Q I (I - C A ) G A) }c k 0 0 - C ] Q [ ( I - C A) + ... + ( I - C A ) ] C | | k k Q Q Q 0 2 whence from ( 3 . 2 9 . 1 ) and 3 . 2 9 . 2 ) (3.29.4) Max [ E l < k 1 1 1 - 1k - CQA | I For lower bounds one has from (2.12) (3.30.1) Ejj = R But i K A \\\l\ ^ 2 - 1 so that Max I R K I = E^. A RQ* = E ^ and also 1 by ( 3 . 2 2 ) . This immediately gives the two bounds: R - °- (3 3 2) n K jiTlfT * ^ l E k ' and (3.30.3) n IIAll 2 Unfortunately, ,-1 RQ -11| k+1 bo I 2 < Max IU I K w i l l probably be unavailable. However, use of 24. (2.10) gives another version: k+1 2 k+1 [Tr (R* R ^ " 1 ] 2 n ^ k+1 2 ^ ( K Q RQ) min whence ^ ( R / R Q ) (3.30.4) ^ Max I E J > n* IUII 2 Inequality (3.29.4) i s useful for estimating the error after obtaining C k but i s not suitable for obtaining lower bounds for the number of iterations required to produce an approximation with a maximum specified error. For this purpose a poorer bound i s available. Since (from (3.29.1) and (3.29.2) once more) we have k Cn = [ f i (I - C.A) ] C 1 f i=0 u 1 i t follows that Max I E J ^ HEJI $ j|£ 4 r s 1=0 Cl-CoA)*l U 2 ; l i - v l S l M 2 2 25. whence , , , I I I - C Q A I I ^ | Max I E J ^ (3.31) IIc il 1 0 2 -—2 1- IU-C A.. 0 2 2 Using this relation we can find a lower bound for Max I E J 4 0( . The result comes out to guarantee that k >y (3.32) k which w i l l log + log (1 - || I C - lOg A II Q ) - III - C Q A If on the other hand one assumes that log (IICQL CQAJ 2 Hc^^lls ^ 1^0 ^2 more loosely speaking, one assumes that II C Q l| - °* ^> can be substituted 2 for llCj_-|. 1^2 ^- (3.29.4) without introducing too serious an error, one n c can obtain the bound ( 3 b 3 3 ) k > l6gQC-log ( l l c f l 0 log lll~C A|l 0 2 II i C A|| - 0 2 + 0() 2 If the best bounds obtainable from (3.32) and (3.33) are called and that k , respectively, one can make the following comparison. Assume 2 CQ i s a very poor approximation, i.e. ||RQ l l > > 0, but that 2 the iteration does converge, i . e . , able assume k 2 k^ jjRQ || 1. < 2 Since we can reason- c( < ^ I | C Q f f , and since the denominators of k^ and 2 are equal we need only write 1 - ||RQ 1 2 • «£ > 0, where € i s by hypothesis, small, and say that for sufficiently small numerator of k, £=J log 1 „ „ '/c » 0 < 0 2 ) 26. whereas 1, numerator of k_ ~ l o g —rrr,—s— 2 II 0 L <o c The presence of the £ f a c t o r i n the numerator of shows that, f o r the case considered, (3.33) gives a much better bound. The assumption necessary f o r the v a l i d i t y of (3.33), namely || C _^ I g K ^ || CQ I i s quite reasonable i f A matrix since i n t h i s case A~^ CQ and i s an i l l - c o n d i t i o n e d cannot d i f f e r too r a d i c a l l y from (or, hence, from each other) without the l i k e l i h o o d that the convergence condition flI - CQA 1 2 < 1 will.be violated. To support t h i s assertion, a few r e l a t i o n s involving matrix norms and the condition of a matrix w i l l be developed next. Let. A C be the change i n the inverse from a change A A (3.34.1) C of a matrix i n A. By d e f i n i t i o n one has (A - A A ) ( C - A C) = I whence (3.34.2) ( A A - A)( A C) = ( A A) C One can re-arrange (3.34.2) t o obtain various bounds, thus: L\ A = ( A A - A ) ( A C) A (3.34.3) llACj2 > A ' A IIA |L 2 )| A gives - ~ rr - A A l| A resulting ( A A). C - ( A A - A ) ( A C ) whereas A Cl > (3.34.4) 2 gives rr-S -n U - A A L X For an upper bound, one can say A C = (A A - A ) " (A 1 A) C - ( AC - C ) ( AA) C I U C L (3.34.5) Suppose now that S whence lie - A c I L A A and lie I L U A A C are small compared to A and C more precisely, suppose that I A- A A | | 2 ~ lUt, and |c- A C ^ Then from (3.34.5) and (3.34.3) we obtain |o|* (3.34.6) » - l ^ i l |AA|| 2 , 1 Ull* Thus i f p i s the von Neumann condition number of A defined by z X 3 def. 'max X where X > max ' X • are the absolutely largest and 'min * & min smallest eigenvalues of A , reference to ( 2 . 8 ) shows that for real symmetric A one has: 28. (3-34.7) ^ ^ 2 ? >y A :A|I 2 Inequalities (3.34.7) show, by the following considerations, why the von Neumann condition number i s i n f a c t a reasonable measure of the condition of a matrix. By the condition of a matrix one means, loosely, speaking, the extent to which a small perturbation of i t s elements does or does not r e s u l t i n a large perturbation of the elements of i t s inverse, A reasonable t r a n s l a t i o n of t h i s concept i n t o terms of norms i s obtained by defining a condition number p* norm (J || (3.34.8) and A by s e l e c t i n g some defining A * = r f o r a matrix Lim def. k - * 0 .. Max, over // A A || I A .. One can then say that than 1. i s ill-conditioned i f For the u n i t matrix, always be taken to be possible A A. 1 p * p* i s much l a r g e r can, for s u i t a b l y defined norms, which w i l l then be a lower bound of Since the upper bound, p =» IIA \\^ ll C ll^ a c t u a l l y be attained (namely when the column vectors of A p * for a l l , i n (3.34.7) can A are a l l 29. eigenvectors belonging to that ft \jax)> o n e s e e s (3.34.7) f r o m and (3.34.8) i s the particular condition number, ft *, resulting from selecting |( | 2 as n | . ( * With this preliminary discussion, we return to the plausibility of the assumption & Jl C _ 1| k 1 || C H for ill-conditioned Q 2 A, under the restrict- 2 ion that II I - C«A II < 1. Consider any approximate inverse C^ 11 BJ 2 relation (3.34.9) =• ft || I - C A l| 2 « (AJ I - 2 * then which, by the previously used , gives . 2 2 I II ( A C) A || ||C|| A C || and set C^ = C - A C , IIRJI . C II ^ 2 II A C ^ ft || 2 2 c 2 If one makes the reasonable assumption that for randomly chosen matrices I I R ^ || A C the resulting values of * 2 w i l l not be clustered near the lower end of the possible range given by ill-conditioned matrices (large ft ), (3.34.9), |j A C Jl 2 ——-r——— fl C || i t follows that for must be kept very 2 small for i f i t i s not. HR^ || 2 w i l l , under the above assumption, probably become larger than 1. The speed of convergence of appears only implicitly i n (3.26) also depends on ft . This fact (3.29.4), (3.31) i n the following guise: if A 30. i s i l l - c o n d i t i o n e d , i t i s correspondingly li I - C A 1) i s very small. 2 N u * A i s positive definite. > W - n ' > f maxl> >Si-l'—' X A eigenvalues of converges. S p e c i f i c a l l y , i f o( = k (3.35) where be Then A . I | I - C ft 0( . 0 A | I I I I - CQA||^ (3.36) 1 < | ^ m a x \ SO - Let C Q «= i I being the that the i t e r a t i o n (3.26) , k > 1, we have | - 1 - ^ - 2 i s as defined previously. J| A 1. C Q such that To bring out the dependence more e x p l i c i t l y , suppose, f o r convenience, that w h e r e hard t o obtain a In p r a c t i c e , a suitable oC would or, i f A were not p o s i t i v e d e f i n i t e , 0C- i|A A T li with C = ^ - A T . For a consideration o f the l i m i t i n g accuracy of C K obtained from (3.26) using d i g i t a l operations, one has to decide on the mode of computation and, i n some cases, on the exact sequence i n which the operations are t o be performed. (3.26) i s to be computed i n the Suppose that the i t e r a t i o n of form (3.37.1) C K + 1 = C © (I© C 0 Q o J) o c k where a l l matrices are now d i g i t a l ( i n the sense of Chapter I I ) . The method of round-off i s not s p e c i f i e d but i s i m p l i c i t i n the d i g i t a l operations (+), etc. I t i s assumed that A" (rather than some matrix A 31. from which i t i s possibly obtained by round-off or other approximations) i s the matrix being inverted. As a measure of the accuracy of compute the norm of the l e f t r e s i d u a l , (3.37.2) R ( k + 1 ) £ - = C -^ , we k+ 2 I - C ^ A I - [ CQ © ( I © C Q A) O C^] A O whence also (3.37.3) R( D? " " C © ( I © C 1 • A) • C ] A C k+ 0 0 k +• Cc + ( i © c 0 - C^A - [ ( 1 0 + C(I0C - C + ( T - (CQA) k . A) . C Q k o 0 K (0)£ R (k)0 (3.37.4) ' t R ( h K e f o r e + 1 ) E O A) C j 0 A) C A - C A + (CQ k A+ CQA C K ( I A) O k Since the "surplus" terms may be written as R ] A k • A) C ] A - ( I © C O C C O I) o c ] A 0 CJA A - CQA)(I - C A) K or S ° l g has the form: n - [ ( C 0+ M) L ( C 0 © i^)] A+ + [ ( I - M ) - (l©M )j ( I - R 3 + [CQ o A - 3 CQA](I - R \ . \] [M^- ( k ) g ) * ( K R ) ) G ( 0 A ) E R ( k ) f t 32. where Mg normalized and M^ are constant throughout the i t e r a t i o n . Hence f o r f l o a t i n g point, taking norms and using (2.14.1), (2.14.2) one has: IIH/^JI (k+l)t (3.37.5) $ n e / Hill 1 + 2 n(2n-l)£ Hill (3^ 2 2 + VST € /S* (1 + llR( )£ll ) 3 k • n(2n -1) 6 p + 1 1 where in W 2 (l + II R k 1 1 W ( k ) g H) 2 2 q^, q , q^, q^ are the l a r g e s t f l o a t i n g point exponents occurring 2 CQ + Mj,, MgCj^, I - M^ and CQA respectively. Replacing q^ i = 1,, by the largest of these, q, gives: (3.37.6 U R( k + 1 )^" 2 ^ C2n 2 III +C H S || ( 0 ) + 2 n 11^,^1 a - n + vn] e p q + ( 2 n - n + Vn) € £ ] ||R g/l 2 Q i which we write, f o r convenience as 0.37.7) 2 « ( ( e * "Vie + |l5 0) 11 2 S) 11 33 For one can therefore set CC + ( (3.37.8) S - 1) H % H / ( R + H 0 2 This condition w i l l hold i f 2 < 0 2 l i e s i n the following i n t e r v a l (neglecting terms of l e s s e r importance): l5)ll (3.37.9) where 6 (|[l-2n e/3 -Y], 2 2 q |[l-2n £ /3 2 q + V]) V i s defined by V -.Sl (3.37.10) + 4n 4 £ 2 /3 2 q - (12 + 8 III II ) n Furthermore, since s a t i s f a c t i o n of (3.37.8) ensures that 2 6 £ ll R- II X s a t i s f i e s the condition - provided that hold f o r k = 1 q q also 2 i s chosen large enough to also - (3.37.8) i s a s u f f i c i e n t condition f o r convergence of i t e r a t i o n (3.37.1) u n t i l the " l i m i t i n g accuracy" of the i t e r a t i o n i s attained. Referring back to (3.37.4) and (3.37.5) i t i s seen that q^, i - 1, 4 are determined e s s e n t i a l l y by 0^^, ^Qp-^t RQ and CQA respectively (or, more accurately by quantities d i f f e r i n g from these only through e f f e c t of round-off). Thus a suitable "q" which would hold throughout the i t e r a t i o n would be of the order of magnitude of the l a r g e r of 0 (since ,-H CnA <«l) and l o g j Max A J. P lover i , j . ! 1 value of q one has from (3.37*7) w i t h such a 34. llR )gil (3.37.H) < H*bH (k ( I I R " + 8 ) k+ c - §( I M C 0 + S ) i 2 or (3.37.12) H H C L 0 T 4 H IIR II (IIEJI 0 + a ) k — + ^ — 1-11^)11 - <5 2 since for convergence + 3 < HRLII (3.37.13) II , H H f o o M S \°° W I we have i n the limit as k 1, 2 0 II- i ( 1 * , y , , 2 - n ^ ii II - 1 - II R ) ( 0 ) g ^ 2 ~ 2 ll - n q-s 2 where terms i n n and Vn" have been neglected and € , the basic round1 —s off error, i s taken to be - p , Since ( 3 . 3 7 . 1 3 ) also holds, i n particular, when R ^ j g i s substitute for one obtains with (3.37.8) R^oo)^ R( ){? 0 °n the right side once more with the inequality < replaced by ^ and replacing R(Q)6» tThus an upper bound for II R ^ ^ ^ g (| in f u l l : . (3.37.14) where, as i n (3.37.15) II R(00)ell ^ - (2n 2 - n +• n ) € ^ q - V] (3.37.10) V - / ( 8 » [1 + - l) 2 4 * - 4(n -n3 n 4 + 5 / 2 )^V* - 4 ( 3 + 2 IIAll^n 1 + 2(nJHT)€/S l /2 i 1 2 6/5 q is - 35. t h e exponent or of q b e i n g now t h a t o f t h e l a r g e s t element o f C^^^ ^ A A «• I , w h i c h e v e r i s t h e l a r g e r . (3.26) r a t h e r t h a n Iteration C K + « 1 C + Q 0^(1 was ACQ) - (3.37.13), s i n c e i t l e a d s n a t u r a l l y t o t h e bound f o r t h e l e f t r e s i d u a l (3.37.14). A bound f o r f( fc rather t h a "^(r)^ n i s d e s selected i r a D l e since i t l e a d s most d i r e c t l y t o a bound f o r t h e norm o f t h e e r r o r v e c t o r o b t a i n e d when s o l v i n g t h e l i n e a r system inverse by C^g v = C^ ^g k (3.38.1) I n t h i s case, i f ^ ( ^ j ) ^ ° b (where i t i s assumed t h a t x - v = Cb - C C b R c (k)g i s of the form X ^ P ^ A = A x and i s approximated b = b) one h a s : 00* - (k)£ • (C-C u by t h e use o f t h e approximate . - where Ax = b ( k ) £ + b + ) b c (k)* + b " (k)£ c b u U (Mow - Mw). U s i n g (2.9) as t h e most c o n v e n i e n t norm, I, one has t h e n (3.38.2) o Max i ( x - v)| = ||x- vL ^ llR (k)e L UxlU ||u|| + whence (3.38.3) Max (x - v) £ p Max x over i lover i + (2n - 1) G £ ( 36. p - Max where, by ( 2 . 9 ) , 0 72 ± Ir^ \ , ± R( )£ = ( r ^ j ) , and q i s the k f l o a t i n g p o i n t exponent, o f t h e l a r g e s t element o f v . I t might seem a t r i v i a l p o i n t t o w o r r y o v e r - as i n t h e f o r e g o i n g whether a good l e f t i n v e r s e o r a good r i g h t i n v e r s e I s being computed s i n c e one h a s (3.39.1) whence R (R(o)r) ( 0 ) p = A - ( ^(0)0A ||i (3.39.2) I-AC || P (0)r 0 ) r = A d - C ^ A " 1 311(1 n e n c e i | | A || 2 2 |U (0)6 l| P HA" !! 1 2 However, as MENDELSOHN [ R e f . 9] has shown, f o r any p o s i t i v e exist matrices X and A k and € t h e r e such t h a t Max (I-XA) I >. k over i , j 1 j Max (I-AX) j S £ ' over i , j but That such b e h a v i o u r i s , a s one would s u s p e c t , c h a r a c t e r i s t i c o f i l l o o n d i t i o n e d m a t r i c e s i s shown, by t h e f o l l o w i n g r e l a t i o n due t o H o u s e h o l d e r . Letting by 1 Max . ( I - AX) ^ k over j ' (2.12): (3.40 1) 8 k ^ ||l - AX || 2 and Max ( I - XA) < £ I over i , j • one c a n s a y 37 whence b y (3.39.1) k S || A I I J|l-XA|| 2 - I U I I S n € gives 1 2 one c a n t h e r e f o r e c o n c l u d e A A * X X (3.40.2) max where '2 /(A" !! 2 From (2.7) 1 2 but b y use o f (2.12) once more, t h i s k H A " A" 1 A" * 1 n€ max * here i n d i c a t e s the transposed conjugate. Hence f o r H e r m i t i a n one c a n s a y X (3.40.3) max x . n£ A min w i t h t h e l e f t s i d e b e i n g t h e vo^Neumann c o n d i t i o n number o f A A . f i n a l comment on i t e r a t i o n (3.26<)' a s computed by (3.37.1) w i l l be made b e f o r e p a s s i n g on t o a n o t h e r t y p e o f i t e r a t i o n . I f i n (3.37.1) one assumes t h e u s e o f a double " p r e c i s i o n procedure r a t h e r t h a n t h e ordinary normalized f l o a t i n g p o i n t method p o s t u l a t e d i n t h e f o r e g o i n g , one l e a v e s unchanged t h e r o u n d - o f f error associated with the operation (+) b u t r e d u c e s t h e e r r o r f o r d i g i t a l m a t r i x m u l t i p l i c a t i o n ("o") same s i z e a s t h a t f o r (+) s i n c e r o u n d - o f f accumulation o f a double-length assumed e x a c t . t o the i s now o n l y p e r f o r m e d a f t e r i n n e r p r o d u c t w h i c h c a n be reasonably One must t h e r e f o r e r e p l a c e i n e q u a l i t y (3.37.5) b y A 38 (3.41.1) A ^ n £ /S^ + v^H e I I I II + n € / S T ftA II 2 ^ (1 +l|R( k )ell ) + n 6 / 4 (l+lfe(k)£ w h i c h r e s u l t s , i n p l a c e o f (3.37.6), i n (3.41.2) || R ( k + 1 )(,l| ^ ( 2 n IIA II + i/TT+ n ) € H(o)e- 'I c + + R + f ^ ^ "(k)J' 6 R 2 w h i c h g i v e s a s t h e analogue t o (3.37.14), (3.37.15)s (3.U.3) llH ( o 0 ^ll § [1 - (n + vTT )• c 4 -Y3 where (3.41.4) ^ - Cl - 2 {(4 UA II + 3)n + V i r } £ p 1 / 9 « / + ( n + 2n 3 \ 2 2q , 2 + n) € £ ] 2 2 q h a v i n g t h e same s i g n i f i c a n c e as i n (3.37.14) e t c . o f l e s s e r s i g n i f i c a n c e , assuming £ ft^ By d r o p p i n g terms i s s m a l l , and a p p r o x i m a t i n g 39. fpr ' V , one c a n make t h e f o l l o w i n g c o m p a r i s o n : f o r ordinary normalized f l o a t i n g point (3.42.1) " H 4 0- 2 * ll A" n £/S> 2 q whereas f o r t h e d o u b l e p r e c i s i o n procedure d e s c r i b e d , (3.42.2) '^(oo)^ ^ 2(| + llll )ne/3 q 2 A l l f o r e g o i n g e s t i m a t e s c o u l d be improved b y u s i n g p r o b a b a l i s t i c r a t h e r than s t r i c t i n e q u a l i t i e s . However, t h i s a s p e c t w i l l n o t be d e v e l o p e d here. The m a i n f a u l t w i t h a l g o r i t h m (3.26) i s , as s u g g e s t e d by r e l a t i o n (3.35), t h e s l o w n e s s o f convergence f o r v a l u e s o f p> l i a b l e t o a r i s e i n practice. V a r i o u s m o d i f i c a t i o n s o f (3.26) a r e n o t h a r d t o c o n s t r u c t b u t one i n p a r t i c u l a r g i v e s a marked g a i n i n t h e speed o f convergence. This i s t h e a l g o r i t h m recommended b y HOTELLING [ R e f . 10] and sometimes r e f e r r e d t o a s Newton's method f o r f i n d i n g t h e i n v e r s e . The method i s b a s e d on t h e i d e n t i t y (named v a r i o u s l y a f t e r C a n t o r a n d E u l e r ) w h i c h t h e i n f i n i t e sum f o r ( l - x ) (3.43.1) (1 - x ) " = ZI x k k=0 1 as an i n f i n i t e J | (1 + x k=0 2 expresses product: ) for \x| < 1 The p r o o f i s s t r a i g h t f o r w a r d and a p p l i e s e q u a l l y t o m a t r i c e s . 40. Letting A be any m a t r i x , one h a s : I - A " * = (I+ A^Xl (3.43.2) 2 3 - A **) 1 2 (I • A °) — ( I + A ) ( I + A ) ( I + A) 2 j TT 2 k (I+ A 2 ) k=0 I f now b o t h ( I- A (3.43.3) (I- >1 A ) " ) and ( I - A) = (I- A 1 " )- 2 0 1 (I - A) a r e n o n - s i n g u l a r , one c a n w r i t e TT 1 (I+ A 2 K ) k=0 I f we s p e c i f y ||A II < 1 then 2 2 and b o t h ( I - A) and |U ? 0+l ) ( I- A ll ^ 1 f o r j - 1, 2, 3, ... w i l l be n o n - s i n g u l a r ( f o r a l l j ) s i n c e t h e i r v o n Neumann s e r i e s then n e c e s s a r i l y converge t o l i m i t s which are t h e i r i n v e r s e s . a l s o i n t h e l i m i t as (since [U 2J+1 (3.43.4) || Thus (3.43.3) must h o l d f o r a l l j and hence )| A || 4. 1. j — o o , provided We t h e r e f o r e o b t a i n -?0) -1 ( I - A) TT ( i + A * ) k=0 Next we a p p l y (3.43.4) t o (3.44.1) A " = iH 1 i=0 for Ul < 1. t h e summation o f t h e s e r i e s i n t h e e x p r e s s i o n ( I - C A)] C Q 0 w h i c h r e s u l t s f r o m (3.20), (3.24), (3.25). A - = 1 [I - (I - whence b y (3.43.3), s i n c e (3.44.2) A -1 -1 CQA)]" C Q | | l - C A \\ 2 < 1 Q » [ I- ( I- CQ 1 ? C 0 We say k + 1 A ) ^ ] by a s s u m p t i o n , -1 2 k CI + ( I- T T X 1 CQA)*" ] i=0 for k - l , 2, 3,, ... I f we d e f i n e t h e k ' t h a p p r o x i m a t i o n t o t h e i n v e r s e a s (3.44.3) Cj, = C T T [ I+ (I- C A) ]] C 2 0 0 i=0 Then f r o m (3.44.2) (3.44.4) = \ .-1 „ rr~ ^2 A - C j = [ ( l - OQA) „ „ „ . 2 0.(1 - C Q A ) k + 1 k + 2 N + c ] 2 k + 1 Thus t h e e r r o r f o r the^k'th a p p r o x i m a t i o n i s o f o r d e r ( I - C Q A ) k+1 (I- CQA) rather than of as f o r i t e r a t i o n (3.26). To o b t a i n an a l g o r i t h m f o r f o r m i n g t h e s u c c e s s i v e a p p r o x i m a t i o n s (3.44.3) we s a y , d i r e c t l y from (3.44.3): k+l 9 (3.44.5) \ = [ I+ ( I- +1 CQA)'' ] and f r o m (3.44.2) and (3.44.3) t h i s i s (3.44.6) C K + 1 = [ I- [ (A" 1 C " 1 ) " 1 - I} ] C K 42. or (3.44.7) = [21 - C A] C k Algorithm (3.44.7) expression (3.45.1) k w i l l be r e f e r r e d t o a s H o t e l l i n g ' s i t e r a t i o n . f o r t h e r e s i d u a l i s e a s i l y obtained. R( )£ k An We have: 1 - 0 ^ = = I- C ^ A J C ^ A [21 or (3.45.2) R (k)£ = R (k-l)e whence b y i n d u c t i o n (3.45.3) R , . „ ^"(0)g xLrt-w/ \k)t = T h i s g i v e s , t a k i n g norms, t h e analogue o f (3.45.4) l R (k)JI 2 $ ' I 1 ' C 0 A 1 I (2.27) namely: 2 S i m i l a r l y one o b t a i n s (3,45.5) N (k)r R A bound f o r t h e e r r o r of (3.46.1) $ 11 III - A C r F ^ and hence f o r t h e a b s o l u t e l y l a r g e s t element i s o b t a i n e d from (3.44.4). We have: (I - C A) F^ - Q 1=1 ,k+l ] i j 2 43V or (3.46.2) 2 = (I - C A ) r K + 1 k + 1 2 Z_> { ( I " G A ) i=0 ^ 2 0 ? 2 0 J whence, t a k i n g norms (3.46.3) IIE || _ 4 I I I - C k A ( | A ? 2 k + £{lli-c Al\ 1 2 k + 1 0 i=0 ] llcjl % and summing t h e G.P. g i v e s _k+l (3.46.4) Max JE l k $ H F ^ II $ I - C Al, 0 a ,k+l 1 - , || I - C A II T h i s i s t h e analogue o f (3.29.A). Once a g a i n we need a bound i n t e r m s o f C Q T h i s i s o b t a i n e d , a f t e r LONSETH [ R e f . 1 1 ] by w r i t i n g r a t h e r t h a n C^. (3.47.1) \ = R -1 k A' - a f A-H whence Max Q [I \\\ i $ 1 - ( I C A)] W\ H III 0 - C Q A II f _ 1 C Q I|C II 0 HE, i=0 (I - C Q A ) 1 |i 44. or Max I E J ^ (3.47.2) I - C A )| llc ll 1 - ||l- CQA lj Q \\\ 11 Q 2 , z (3.47.2) i s t h e analogue o f (3.31). U s i n g t h i s r e l a t i o n one can, as b e f o r e , f i n d a l o w e r bound f o r a I'EjJ Max . T h i s comes out t o b e : QC + l o g ( l - log log 2 lc H He i 0 k * (3.49) log , t h e n t h e analogue o f (log log log 2 II III - CnA II) - l o g l l C n I - CQA II oC CQ o r RQ i ss t i l l , (3.33) namely (3.33) i s f o u n d t o be - log (oc + 2 log s i t u a t i o n as r e g a r d s l o w e r bounds f o r E ^ as f o r t h e s i m p l e i t e r a t i o n . of t h a t w i l l ensure t h a t on t h e o t h e r hand one makes t h e same assumption a s f o r that The Of k :> l o g (3.48) If ^ k HCQH)} II I - CQA l | i s s u b s t a n t i a l l y t h e same (3.30.2) s t i l l h o l d s and a bound i n terms effectively, unavailable. F i n a l l y , b e f o r e making some c o m p a r i s o n s and summarizing comments, we need t h e l i m i t i n g a c c u r a c y o b t a i n a b l e by (3.44.7) u s i n g d i g i t a l operations on.digital quantities. As i n t h e case o f a l g o r i t h m (3.26) we examine l | R ( k + l ) ^ 2 as a measure o f t e r m i n a l a c c u r a c y . We c a n w r i t e 45. (3.50.1) R T - C ( k + l ) r I k + 1 = I - [(2I©G k o A) o C ] A = I - C(2I©C k o A) o C ] A + 0 ( 2 i e C - (21 © C k k k . A) + (21 - ^ k k . A) C A k - 2 C A + ( C o A) C A - ( C ^ ) k o A) C ] A k + k S i n c e t h e " l e f t - o v e r " terms c a n be w r i t t e n as I - 2C A + C A C j l = (I - C A ) k k 2 k , t h e above has t h e f o r m (3.50.2) 5 (k+l)f" W o C ] A + [(21 - Mg) - (21 e ^ ) ] " % l C A k + [ \ o I - CjX] C^A + ( I - C A ) k 2 k Thus u s i n g (2.14.1) and (2.14.2) and w r i t i n g (3.50.3) " f(kk ii))e ? J R l + + 2 ^ N ( 2 N - ^ o n e h a s ! k "A € q + = I - R( )g 2 "3 n(2n - 1)6/3 (1+ Ill^H ) + l|R(k)g " 2 2 46. where q^, q , q^ a r e exponents ( o f t h e f l o a t i n g p o i n t base 2 p ) d e t e r m i n e d - e s s e n t i a l l y - b y t h e l a r g e s t e l e m e n t s o f C ^ ^ , 2 1 and C^A respectively. C a l l i n g the l a r g e s t o f these q and s e t t i n g , f o r brevity, (3.50.4) Of = { ( 2 n - n ) III (3.50.5) 5= ( 2 n - n + V~n) e )l 2 z + 2 n - n + V~n ] € 2 f 2 q we c a n s e t f o r t h e improvement upon t h e k ' t h s t e p : (3.50.6) II R ( k + 1 ) e a II ^ 2 • 5 llR ( k ) t II || R 2 + J| ( k ) 2 ^ 15 ( k ) t II 2 or (3.50.7) d + ( 5 - 1 ) llR(k)e,U + H (k)&l' < 0 G £ Q R 2 This r e l a t i o n s h i p holds i n t h e i n t e r v a l (3.50.8) II\\\ 2 I (| [ 1 - ( 2 n - n + VU) 2 | [ 1 - ( 2 n ^ - n + V n ) € (3 q - Y , + T) where V = i /(5- l ) 2 - 4oC q (3.50.9) ^ CI * 4 ( n * - n +> [ ( 6 3 + n ^ )C ^ * - III II )n - 1 /llj^) - 4(3+2 6/n} 6 ^ ] 1 / 2 n - £ /3 2 47 T h i s i s i d e n t i c a l w i t h (3.37,9)> (3.37.10) i f terms o f l e s s e r s i g n i f i c a n c e a r e n e g l e c t e d w i t h t h e e x c e p t i o n t h a t t h e t e r m (6 replaces the term 2n £ i n (3.37.15). f| A \\^) n 6. £ The s i z e o f Q "q" i s governed, r o u g h l y s p e a k i n g , b y t h e l a r g e r o f 2 and t h e l a r g e s t element o f C Q O ^ r a t h e r t h a n o f 1 and t h e l a r g e s t element o f & A™ 1 A as f o r (3.37.14), (3.37.15). Use o f d o u b l e - p r e c i s i o n accumulation i n t h e t r e a t m e n t o f a l g o r i t h m (3.26), t e r m i n a l accuracy of inner products parallels that I n view of the roughly equivaleht o b t a i n a b l e b y u s e o f i t e r a t i o n f o r m u l a s (3.26) and (3.44.7) and t h e s u p e r i o r speed o f convergence p r o p e r t i e s o f t h e l a t t e r , t h e r e seems t o be no p o i n t i n u s i n g (3.26) except when t h e m a t r i x so s t r o n g l y d i a g o n a l ( o r , more g e n e r a l l y , when so good a f i r s t tion CQ e x i s t s as t o make II RQ l| approxima- very small) that 4 i t e r a t i o n s o r l e s s o f (3.26) w i l l g i v e a s u f f i c i e n t l y good a p p r o x i m a t i o n A \ Ai s t o the inverse Use o f t h e d o u b l e - p r e c i s i o n p r o c e d u r e f o r both a l g o r i t h m s i s h i g h l y d e s i r a b l e from the p o i n t o f view of t e r m i n a l accuracy. 48. CHAPTER IV Hotelling's Algorithm; Modifications and I n i t i a l Approximations The foregoing sections of Chapter III developed and analysed the Hotelling algorithm (3*4472 )• Since the successful application of the iteration depends upon the availability of an appropriate Cg, methods for constructing approximate inverses w i l l be examined here. Various modifications of (3»44"»7) designed to make i t converge faster or with more easily available starting values, w i l l also be of interest. The derivation i n terms of Euler's identity, while suitable for obtaining error bounds and for linking the algorithm with the basic series expansion (3.19)* i s , from a heuristic point of view, inappropriate as a departure point for some modifications which w i l l be introduced. Accordingly, three other approaches are listed each of which can lead to different lines of development. The f i r s t and simplest derivation of (3.44.7) follows from the stipulation that the residual at any step be the square of the one 2 immediately preceding. llR (f < k+1 II Rol| R^^ » R k , then gives directly and (3.44.7/) then follows from 2k I - Cfc The resulting equation +1 A - (I - C k A) 2 . The second approach i s more of a recapitulation of earlier material developed from a slightly different point of view - one which leads more naturally to a modification which w i l l be examined later i n some detail. Since i n (3«l6) P i s non-singular but otherwise arbitrary, we can select 49. Q o I and write the iterative formula as (4.1.1) x( k + 1 ) = x< > + k where z^ ) i s as yet arbitrary. The "perfect." choice of z^ ) i s k k therefore, z(k) , ( k ) v where v(k) i s the error vector defined by v ( ) » x - x ( ) . If rOO, the residual vector, i s defined as k k r ( ) = b - Ax^ ) k we have v^ ) - A ^rC ) so that the error vector k k - k satisfies the same equation as does x i t s e l f , provided that b is replaced by r ( ) . This fact provides the basis for a common method of k obtaining improved solutions to Ax = b with the use of any suitable direct method for solving linear equations. The direct method i s applied repeatedly to solve (approximately) Ax = b, A v ^ = r ^ ) , Av^ ) = r ^ ) , etc., with a l l quantities being appropriately scaled at each 2 2 step so as to gain the maximum advantage from the diminishing of the r^'s. The sum x + v ^ ) + v^2) + ,,, ± then an improved approximation to the 9 s true solution, x. Here, however, we are concerned only with a particular way of writing an iterative formula and not i n any iterative procedure as such. Taking the identity x = x(k) + y(k) - x( ) • A"lr( ) k and replacing A" 1 k with an approximation which, looking ahead, we can c a l l CQ, we are led to the iterative formula (4.1.2) x(k+D « x(k) + Cor(k) 50. This formula i s seen to be merely a re-arrangement of the simple iteration x (k+l) x (k)^ a c x ^ + ( i _ CQA) X ^ examined earlier. b (k+l) . Replacing b by I and Cjj+i, respectively, gives us the matrix inversion by form of the iteration, namely C K + 1 =» + C Q R ^ ^ • If we now approximate A"^ by the best available inverse, C , at each step rather than by our k fixed approximation, CQ, we have (4.1.3) C k + 1 =C C^ k + ik)r or, equivalently, (4.1.4) ^1 k* Hktfk mG This last formula i s (3.44.7) once more. Third, and last, (3.44.7) can be taken as a simple special case of the general class of iterations defined by (4.2.1) where C ^ => f ( C A ) C K k f ( ) i s some function which approximates the inverse, i . e . , f(M) ^M" for a class of M's. In the case of (3.44.7), the function 1 f(C A) i s given by )? (4.2.2) f(C A) = I + (I - C^A) k that i s to say, the series expansion for (Cj^A)""''" i s approximated by the f i r s t two terms of the series. The modification of (3.44.7) suggested by the f i r s t derivation listed 51 above i s the obvious one of insisting on the relation (4.3.1) where m \ + 1 = if i s some fixed integer greater than 2. For example, m = 3 results i n the recurrence formula (4.3.2) - C3I - 3 0 ^ + (C^A) ] 2 which converges according to the scheme »R k + 1 JU llEoll 3k In general, such formulae w i l l have the form (4.3.3) where - tfijl -yS (AC ) 2 k + — WP? ^?* ! 1 + 1 Cfe ^ i s the binomial coefficient of x^ i n the expansion of ( l - x ) . m If i n particular m » 2 q where q i s a positive integer, one can easily obtain formulae of the type exemplified by (for the case q = 3)s (4.3.4) C k + 1 = [I + (I - C A)][I + (I - C A)2][I + (I - 0^)4] C k k k where the powers of (I - C^) are of form 2*- and are hence economically obtained by repeated squaring. At f i r s t sight i t would appear that these formulae give an increased speed of convergence over (3.44.7). In terms of the gain per iteration, this i s , of course, true but when rate of convergence i s measured i n terms of operations - specifically i n terms of matrix multiplications - this 52. advantage i s seen to be illusory. k'th iteration the next m For formulae of type (4«3.l). at the matrix multiplications decrease R^ to Rjj+1 m and Hfc+1 k» however, i f for convenience, we suppose a R Hotfelling iteration (with m = 2) takes tions and R^ p = R^ P + R k • Since 2P > m = 2p to Rj^p with for a l l m =« 2p then the m matrix multiplica- p > 1, any selection of m, m > 2 merely decreases the efficiency of the iteration formula i n terms of "convergence per matrix multiplication". are equal - for a l l m = 2 Formulae of type (4.3.4) - to the Hotelling iteration but their added q complexity (and increased round-off error per iteration - a factor that adversely affects the limiting accuracy obtainable) enables one to dismiss a l l formulae of both types (4.3.3) and (4.3.4) as worthless for m>2. The second approach mentioned leads to a modification of Hotelling's algorithm which might be called the Optimized Hotelling method. Returning to (4.1.1), and writing i t i n the form (4.4.1) x(k+D =• x(k) + A&) (k) z (where X^k) i s a scalar as yet undetermined which will vary with the well-known fact that for arbitrary (non-zero ) k) we use z ( ) there exists a k such that v(k l) + i s a better approximation to the solution x than i s To formulate this more exactly, one can specify that (using r(k) as defined before) (4.4.2) || r ( k + D |j ^ | | ( k ) || \ Thus we want r ||r(k)||- - ||p(k+l)|.- > 0 53. | ( r ( M l 2 - | r M - X< > A ( k ) | | k z but since for any vector 2 > o p we have ||p|| » p»p, this i s 2 2 2 2 or 2[r< W >] k (4.4.3) k X< > k -llAz^ll A< > 2 k 2 > 0 2 This function of \(k) A 2r > » A z " —77T~5 (k 4 4 ( t) 0 and ( k ) ; l | A « . ( . . has roots at k ( 1 ani has a (positive) maximum at % - ) x Az « This value of " 2 (the inverse Rayleigh quotient) i s therefore such as make iteration (4.4.1) converge (for a l l z(^) that do not make r ^ ) • Az^ ) zero). k is k o A " M k \ The "best" direction for z ^ ) , as previously mentioned, k For this choice X ^ however, even i f only an approximation one can s t i l l compute the length for the direction becomes 1 and xC i s used i n place k + 1 ^« xj of X" , 1 corresponding which gives the optimal C^r^). If i n (4.1.3) the column-vectors 54. of C j ^ , Cjj and pk^(k)r regarded as components of a set of n vector a r e equations of form (4.1.1), one can think of the corresponding set of n optimized equations of form (4.4.1) with column-vectors of vectors of form column of C^R^^ given by (4.4.4), can, moreover, be regarded as a set z^ ) = C^r^^. n To multiply each such vector (i.e., k CjjR^jp) by i t s particular matrix by a diagonal matrix of The /\ k one has to post-multiply the /v where • Diag( k i«*l,2,...,n. Furthermore, A.^ must, by (4.4.4), be the product of two diagonal matrices (4.4.5) (4.4.6) and D D^ * hPk i = ^ D 2 where 2 k( = (ACkRk) where the symbol * 1 * (ACkR ) k here indicates that only the diagonal elements of the "star-product" are computed to form D^, D . 2 i n (4.1.3) then gives, with re-arrangement: (4.4.7) C j ^ = C (l + A - AC k k A k ) k where (4.4.8) /V = [(I - A C ) k T * A C ( l - AC )3 k k CACjjCl-A^)]*' * [AC (l-AC )] k where the "quotient" represents DjD"^* k Making the above change 55. From the nature o f necessarily converges. the Optimized Hotelling i t e r a t i o n This r e s u l t s from the fact that at each step a l l the column-vector norms and hence, i n p a r t i c u l a r , the largest columnvector norm, decreases. matrix C k Since the largest column-vector norm of the i s a norm f o r C^ i t s e l f , convergence follows. Moreover, ( 4 . 4 . 8 ) shows that J\^ —>I as k—><=*=>. i n the l i m i t as C^—>A The Optimized and the standard H o t e l l i n g formulae are thus equivalent i n the l i m i t as k grows l a r g e . The Optimized version of the formula converges regardless of the size of III - Over t h e i r common ACQII. range of convergence the Optimized version must necessarily be the faster-converging of the two since /\- , being optimum, pSo/iOfver-be k worse than the unit matrix as a choice. From (4.4.7) one can e a s i l y find (4.4.9) R( l ) r =" ( k ) r C * R k + so that, i f we define (4.4.10) R ( k + 1 ) r - 1 H k A c k^k) = (I -AC /Yk)* R O H Q H ^ . . . , k w e ^ a v e ^ Results obtained i n practice with the Optimized Hotelling i t e r a t i o n were quite good, p a r t i c u l a r l y i n the case of very i l l - c o n d i t i o n e d matrices. Convergence was slow (of geometric character) i n the range HR^II & 1 but was extremely f a s t (of exponential character) both above and below t h i s middle p o i n t . The good convergence c h a r a c t e r i s t i c s , however, are obtained at the cost of three matrix m u l t i p l i c a t i o n s and the approximate equivalent o f three matrix-vector m u l t i p l i c a t i o n s . Furthermore, three 56. entire matrices must be stored throughout the process. an i n i t i a l approximation In the absence of CQ which i s known to be good, the choice 0. = A? i s a good choice since this gives the familiar "Steepest Gradient" step on the i n i t i a l cycle [Refs, 12 and 13]. Some methods of obtaining a starting approximation CQ are discussed next. For positive definite A, the Diagonal and Seidel types of approxima- tion w i l l necessarily give convergence. mentioned earlier, take CQ =» D""^A^ where the elements of D are the absolute row-sums of A^A. Convergence i a l l these methods unless to some matrix algorithm, of B For a general A, one can, as S j however, apt to be slow with A i s quite well-conditioned. If A i s "close" of special form which can be inverted by some special B"-'- w i l l be a suitable i n i t i a l approximation for the inversion A. The term "close" as used here refers to closeness element-by- element rather than i n any norm. Special inversion algorithms exist for such matrices as the triangular and tri-diagonal types. Although any direct method of inversion can be used to provide a good starting approximation for further improvement (provided that the matrix i s not so i l l conditioned as to make the result of the direct inversion meaningless) i t would be preferable to use some simpler method that would provide a poorer i n i t i a l approximation for much less work. Such an "approximate inversion" i s available for matrices with inverses which have a l l offdiagonal elements i n any given column (or any given row) of the same sign and - very roughly - of the same size. This type of matrix often occurs among the Leontieff matrices mentioned earlier. The method to be developed leads i n a natural way to a class of inversion methods to be discussed i n Chapter V. These methods have the interesting characteri s t i c of being finite-step iterative inversions. 57. If p and are a r b i t r a r y non-zero vectors and matrix with a known inverse, the inverse of the matrix (4.5.1) Z =» S + p q S i s any Z where T i s e a s i l y found t o be given by (4.5.2) where Z" » S" - - i 1 S" p q S" 1 OC = q^ S"" p. 1 1 T 1 This formula i s a s p e c i a l case of the general r e l a t i o n s h i p given by HOUSEHOLDER [Ref. 14] (4.5.3) (A + USV ")" - A " - A " US(S + S V A " U S ) " SV A " 1 1 1 1 T 1 1 T 1 where A, U, S, V are l i m i t e d only by the r e s t r i c t i o n that a l l terms i n the above be defined. Z" 1 One therefore thinks of approximating of form (4.5.2) where S A"^- by a i s selected from a class of e a s i l y i n v e r t i b l e matrices such as, f o r example, the diagonal matrices. BERGER and SAIBEL [Ref. 15] obtain t h e i r approximation scalar matrix A - Z where D* = d l and a vector Z = D 1 + pv , v T norm i s the one defined by by (4.5.2) p so as t o minimize the norm of i s the vector ( l , 1, ||M|L«» (YL l,j 11,2 by s e l e c t i n g a I)* . 2 l ) , and the They then i n v e r t Z 3 A more d i r e c t approach i s to minimize the norm of I - AZ" 1 since t h i s quantity i s the one which a c t u a l l y determines the convergence of the H o t e l l i n g i t e r a t i o n . (4.5.4) where Z" M - I - AZ" 1 We define 1 i s of the form (4.5.2) namely I - D'""•'* - j - " ^ D 1 1 1 and f p q D"*^" 1 7 (D* non-singular, diagonal) 1 1 «= 1 - Y+aL m. . ) • 2 D —1 p 01 x » —1 T q D For convenience we l e t T y s o t h a t w e D'"*" a D and a c t u a l l y work with the equation (4.5.5) Z" - x y + D 1 T 1 so that the problem becomes one of minimizing f 2 a ||M| + or 1 (4.5.6) f 2 7 j D ) „ UJ _ Az-lll^ || j _ A(xy = Since t h i s i s equivalent t o minimizing work with f + D)l|* . 7 f (x,y,D) « Jf I - AX we 1 2 rather than f-- . For the p a r t i c u l a r norm being used, the above i s equivalent t o f(x,y,D) a T [(I - ? AZ^Kl - AZ ) ] _ 1 (4*5.7) T = T r [ I - AZ- - Z" 1 T A 1 T T + AZ' Z" 1 1 T A ] Considering the lastterm, we have T Tr AZ" Z 1 _ 1 (4.5.8) T A = Tr [ A ( x y T + D)(yx + D) A ] T = Tr (AxyT yx? A?) + Tr(ADyxT A?) + Tr(Axy T DA?) + TrCAD^A^) S i m i l a r l y , the other terms i n (4.5.6) give (4.5.9) T Tr ( Z ~ l T A ) - Tr ( y x A ) + Tr DA T T T T 59. and (4.5.10) Tr (AZ" ) = Tr(Axy ) + Tr(AD) 1 T I f , f o r conciseness, we define \ \ \ Ax = <t , ADy = , we can re-write o ,Xr insert (4.5.8), (4.5.9), (4.5.10) and^the new expressions i n terms of i n t o (4.5.7) t o obtain: f(x,y,D) = \\^\j2 j|y||2 (4.5.11) + + + Tr(AD&D ) - *n, y - Tr(AD) T T - y V - Tr(AD) + n 2 - \\M 2 l(y| + 2<A & - 2-UVy + Tr[(AD)(AD) T T T - (AD) T - (AD) + I ] or (4.5.12) 2 f(x,y,D) « ]|y|l 2 2 |M| + 2 $ T 2 -y )^ + T Now, d i f f e r e n t i a t i n g (4.5.12) with respect to ^ p a r t i a l derivatives equal t o zero f o r a maximum: (4.5.13) (4.5.14) (4.5.15) T 2 - 2 / M l y + 2 DA tfc -2-^ = 0 dy 2 T 2 - ~ 4 » 2lly11 6 A}]/ ^ f .2 D J R +2$ - 2y = 0 2 D A T A - 2 D. • 2 A V A D =0 2 II I - AD 11 * , y, D and setting the 60. The notation used i n equations following interpretation. "^u *^u~^ • to ( 4 . 5 . 1 5 ) has the -\ * u, -ff-=- designates the vector (4.5.12) For any vector a S i m u i l a r l y , for any diagonal matrix, G => Diag (g^), ~b f indicates the matrix M • O^j), % vector Diag (^ ) . Finally, for any (square) matrix, designates the diagonal matrix u = (u^) D t u Diag (m^) and for any represents the diagonal matrix, It i s impracticable to solve explicitly for x (4.5.13) and to but i f x (4.5.15) D as follows. From Diag(uj_). from equations i s assumed given, one can solve for y one can obtain the two equivalent (4.5.13) expressions for y: (4.5.16) y = — ± x (4.5.17) I L = - y - i and likewise from D„ B y Solving (4.5.19) D Ax T (D - D D A Ax (4.5.15) (D - = (x A Ax T one has A A (4.5.18) T D D~l A Ax T for D gives - A ) ^Ax. D D m, ) A ( 4 . 5 . 1 7 ) , - AD) A Ax x (4.5.18) (I D D m Ax A Ax T )(x A Ax T T D T A A - D 2 m A Ax T ) " 1 61. whence substitution i n (4.5.20) y » (D A T A X (4.5.16) gives AX - D a A T Ax) ( x A T I f one sets x = v • ( l , 1, in Z Ax - D * . ^ ) " 1 1) i n (4.5.19) and (4.5.20), one -1 - I T obtains an approximation T a xy + D , to A . The "denominators" (4.5.19), (4.5.20) being i d e n t i c a l , the number of operations necessary to produce Z " * by the method just proposed can, by suitable grouping, be reduced t o 1 3n 2 + 5n additions, 2 n 2 + 7n m u l t i p l i c a t i o n s and 2n divisions. Since the Berger and S a i b e l approximation i s included i n the class of matrices obtainable by (4.5.19), (4.5.20), the proposed method w i l l n e c e s s a r i l y produce a Z " " which i s at l e a s t as good as t h e i r s i n the sense of reducing the norm 1 III could be used i n which the vector - x A Z - l l l . . An i t e r a t i v e technique obtained by the Berger and Saibel method could be taken as the given vector for the present procedure i n place of the choice x =» v = ( l , 1, l ) suggested i n the foregoing. 62. CHAPTER V Rank-reduction Methods A Class of F i n i t e - s t e p Iterative Inversion Procedures The r e l a t i o n (5.1) (B uvV + - B" *f f - 1 1 + v Bi ± 1 u used previously i n Chapter Iv" to obtain i n i t i a l approximations f o r use i n power series expansions i n ( I - A C Q ) , can also be used to obtain the inverse of a matrix i n a f i n i t e number of steps. matrix A Expressing a non-singular i n the form m (5.2) A = S + 0 u k vj k-1 where S 0 i s a matrix with a known inverse, one can, by repeated use of (5.1) i n v e r t the succession;of matrices (5.3) S = S k + 0 £T u v£ k k-1, 2, m This procedure i s expressed by the i t e r a t i v e r e l a t i o n T (5.4) C h + 1 = C* - ^ 1 U + k v ^ V k C, - S" k k 1 and C = A" m 1 ° 1 k l °k + Where + u k k = 0, 1, (m-1) k l + . Any method based on the above procedure w i l l be c a l l e d a rankreduction method. The name derives from a theorem stating that 63.. m TZ u v k k has rank m (m $ n) i f { u j , | v J are two linearlyk k independent sets. Such i t e r a t i v e procedures are conveniently c l a s s i f i e d according to the corresponding decomposition of type (5.2). Three types of rank-reduction inversion w i l l be considered here and w i l l be termed Jordan-type Completion, Symmetric, and Quasi-optimum methods respectively. Since a l l such methods could be termed direct inversion procedures (some i n particular being alternative formulations of familiar direct algorithms), some support should be given at this point for classifying them as f i n i t e step i t e r a t i o n s . The most obvious j u s t i f i c a t i o n i s the i t e r a t i v e character of C5e4) and the consequent appropriateness of applying an error analysis t y p i c a l of i t e r a t i v e methods. Another j u s t i f i c a t i o n i s the fact that some rank-reduction methods have the defining property of an i t e r a t i v e inversion, namely that they produce successive approximations, A" 1 C , to k which give residuals (i; - AC^) which decrease monotonically i n some norm. F i n a l l y , a t y p i c a l feature of finite-step iterations i s that they may often be profitably continued beyond the stage where they would y i e l d the exact inverse i f exact computations were performed throughoutj rank-reduction methods also have this property. The first-named rank-reduction method, Jordan-type Completion, results from perhaps the most obvious decomposition of A: m n (5.5.1) A - S Q + B ^ Z ^ where ^ k»l t h i s gives (5.5.2) A = SQ + 21 C o l nk=l k B ) k ej k B -A - S Q some 64e with the related iteration (5.5.3) C - C k +1 Cu C o l - -3S ^ h l + [Roi° ] k) ( B ) ^ k [Ro T k - 0, 1, ( b ) W k + 1 ] (n-1) Col^ Discussion of this method w i l l be limited here to justifying the designation "Jordan-type" by showing that, at each step, the C^'s obtained by (5.5.3) with SQ = I are identical with certain matrices obtained at the corresponding stages of ordinary Jordan inversion ("Diagonalization"). Jordan inversion consists of forming the successive products JgJ-jA, ... J]A J «J _]_ n n a and, concurrently, the related products I J,, J„J,, .... J J . . . . J , =« A" 1 2 1 n n-1 1 product J ^ J j ^ i ... J]A. and premultiplying J _-^ ••• where 1 J, , i s obtained from the k+1 J-^ i s obtained from by k J^A, A, The effect of J^ i s to "zero the k'th column" leaving a l l previously zeroed columns unaltered. The form of J ^ i s that of a unit matrix with the k'th column replaced by the column vector . (s «; n ( a i k a kk . t i ... a 2 k \k ^ a kk a kk where the primes indicate that the elements i n question are obtained from the product J k J _]_ ... JjA k rather than from itself. The inverse of has a similar from but with the k'th column replaced by the k'th column of the product matrix J _-j_ ••• k A only i n columns J -j. ••• ^ i ^ . k<-> k+1 to n Thus, inclusive. differs from It w i l l be shown 65. that = J _ ^ ••• J k f o k = 1, 2, r k obtained by (5.5.3) and S n where = Sj^" is i s defined by k k (5.5.5) S. = I + (A - I) H e, ej i=l 1 1 It i s clear from (5.5.2) and the foregoing description of J ^ that 1 = J ^ ." On.the assumption that 1 S k = J l 1 *2 J k 1 we now establish that S. , = J T ... J , J. ^_ • S, J„?", k+1 1 k k+1 k K+1 1 1 From the definition of J""\ we have k+1 k+1 J = ( J k k-1 — J J l " A Z ) e k+l k+l e - [(A" S,)- - I] e 1 whence, since S 1 (A - S ) e k k + 1 k + 1 T e k + 1 e e 1 e ^ + I = (A - i ) e k k+1 - <A " I) k+l k+l J + + f c + 1 T e k + 1 , gives \ but, directly from (5.5.5)* the right side i s S -j_ j this gives the k+ desired result. 66. It i s well known that i n the Jordan inversion i t i s possible to meet zero divisors even though A i s non-singular, the process being "safe" only i f A i s positive definite. The Jordan-type completion defined by (5.5.3) suffers from the same defect. In conclusion, one may note that Jordan-type completion i s particularl y suitable for finding the inverses of matrices differing only i n a few columns (or rows) from a matrix with a known inverse. In such cases, only a few applications of (5.5*3) are necessary with each cycle requiring 2 2n 3 + n multiplications as opposed to a minimum of vr ' multiplications for a direct inversion of the entire new matrix. The decomposition A = UDL of a non-singular matrix into an upper triangular, a lower triangular and a diagonal matrix, i s unique. Various direct inversion algorithms make implicit or explicit use of such related factorizations as A =» (D + L ) ( l + U) which characterizes the Doolittle and Banachiewicz procedures. More generally, l e t A = (D + L)(D' + U) where one of D, D* may be chosen arbitrarily. (5.6.1) or (5.6.2) k=l Setting D + L.= X, D + U= Y one has 1 67. The elements of k = 1, 2, ..., n R° w (Y) ' k can be successively computed f o r by the r e l a t i o n s /.•j X T k [Row ] = e A = S A; k k [e (yr) (y) T Col^CRow^] T k which gives, since the f i r s t f a c t o r i n each term i s zero f o r Ito^r (5.6.3) = ( x ) k e where (Y) Row. 3 T k fx) C o l ^ ' • d^ 6 lei - Col^ ][Row^ ] X ) ) T k i s known (D being assumed known), and W*> - g coxW t ) [Row where T fx) Col^ j < k. , Similar- i s obtained from X . . ] A ) k are assumed to have been previously obtained f o r l y , Col£ ^ ( 5 |[Row 1 T j > k [Row^^]^ e k ] e k ' CHof V e cof) k •" != L i s obtained from (5.6.2). The process described above, with Banachiewicz f a c t o r i z a t i o n of A. D 1 » I, i s e s s e n t i a l l y the However, the form of (5.6.1) suggests the existence of a c l o s e l y r e l a t e d c l a s s of rank-reduction inversion methods. The c h a r a c t e r i s t i c of decomposition (5.6.2) i s that the k'th term of the sum i s a matrix of order and above by (n - k + l ) bordered on the l e f t (k - l ) zero columns and zero rows r e s p e c t i v e l y . The symmetric method i s a generalization of a process based on such a 68. decomposition. It should be noted that the name "symmetric method" i s not intended to imply that a l l symmetric rank-reduction inversions are of t h i s type. A succession of matrices, Z , are defined by k as usual, Sg i s any matrix with a known inverse, and (5.7.1) k - 0, 1, - Z - % *+l *** " Z n Z K -^u^ j The set ZQ = A - Sg , where, of n (n-1) vectors, which w i l l be referred t o as the basic set, i s required to be a l i n e a r l y independent s e t . Substitution i n T (5.7.1) shows that more, i f more, ZJJ Uj_ Z Z k u^ =» u^ » 0 f o r i = 1, 2, 1^=0 k + 1 0 provided that f o r i = 1, 2, a 0 f o r i a 1, 2, n. k u-^Zg which, since (5.7.2) Z n (k + l ) . I t follows that However, since From the d e f i n i t i o n of - 0 A = Sg gives n-1 + H k=0 k + 1 ) which from (5.7.1) i s n (5.7.3) A = S 1 T + £ k % + l k + l Zk k=0 T Z Zg n k u Vl \ V l £ kJ u ^ s a ll n e a r l y inde> 2^ v » 0 f o r a l l v, whence A = SQ + (Zg - Z ) (Z - Z Further- then, from (5.7.1) once pendent set, t h i s l a s t r e l a t i o n implies that 2^ = 0 i d e n t i c a l l y . 0. u^ ^ we have A • Sg + Zg whence we can write 69. provided that u -^ ^ u k+ k + ^ 4 0 for k = 0, 1, U^.l writing ^ u ^ for and u (n-l). Zi. * k l 2k k l u + *<* v Now, f k + 1 + i n the basic formula (5.4) one obtains the iteration relation (5.7.4) C | K - 0 1 - R °lc k 1 » 1 4* Z t u k+l ^ u k+l + u Z f c °k k+l haW u k+l k-1 However, by definition S • S k whence, substituting (5.7.5) (Zj - Z ^ ) +1 for k = 1, 2 , . . . , ( n - l ) SQ » A « ZQ one obtains = A - S \ + Q k « A - Cj^ for 1 Inserting this last expression for k = 0, 1, n Zj^ i n (5.7.4) one has, after simplification: (5 1 5 .77 .66) ; C C a = c C k + 1 + + k ( I " -° k V l A ) — - k—+ l ( I u "V which defines the symmetric method. for k a 0, 1, i s clear from (5.7.6). u k —+l & - A C k> " V l ' The retention of symmetry i n n i n the case of a symmetric A and CQ This feature has the computational advantage of all.owing a considerable saving i n the amount of computation (although not as much as 5 0 $ ) . For the case j^fc^ = | kj ^ e v , n the unit vectors as the basic set, (5.7.6) becomes a procedure for performing a method of H. S. WILF's [Ref.16] 70. His proposed method involves computing the successive with the basic set £ (5.7.7) Z k + 1 = Z u k ^ taken as the unit vectors and with Col ^ [Eowft k ) * ^ k Z^'s z ] ) k±l k+l, of (5.7.1) SQ = I , T — k+1 He then uses the basic formula ( 5 . 4 ) i n conjunction with ( 5 . 7 . 7 ) . procedure thus involves computing and storing both k =• 1, 2, ..., n. Z^ and His for The labour involved i s reduced by the presence (as shown below) of zeros i n the successive From ( 5 . 7 . 7 ) i t i s seen that Furthermore, i f Z. has m Z-^ Z^'s. has a zero first-row-and-column. zero rows and columns then Z. retains a l l these and adds a zero j * t h row and column provided that the index j has not been used before. columns of the successive Thus i f , i n p a r t i c u l a r , the rows and Z^'s are taken i n the natural sequence, the Z *s are of the type previously shown to characterize the Banachiewicz k decomposition of A - I. The l a s t type of rank-reduction method examined here i s the Quasi-optimum one. considerations. we define The name derives from the following h e u r i s t i c Considering the basic decomposition ( 5 . 2 ) once more, Z^ = A - S and ( 5 . 7 . 5 ) . Z^ K i n agreement with the notation used i n ( 5 . 3 ) i s no longer, however, t o be determined by (5.7.1). 71. It i s natural to examine choices of x ^, ^k+1 wki * w i l l minimize 0 k+ T some convenient norm of A - [Sg + Z - k iyk+l * + (5.8.1) T * " * * k+l^k+l^ x 1 , e * °^ Accordingly, we define x k 1 f - HZj^ll » » Z - xy H T k and use subscripts on x, y to denote components of these vectors. The norm here w i l l be defined as (5.8.2) f = II i=l II j=l (z i X i y.) 2 r J where Z k = 3 (z^) . differentiating, (5.8.3) ^j- = -2II (5.8.4) 1 ^ « -2:2 (z ± j -x i y^vy. - x. y.) ^ i - 1, 2, n j « 1, 2, n equating to zero and using vector notation gives (5.8.5) [Row^f y - x y y » 0 i = 1, 2 , n (5.8.6) (Zj T T [Colj ] x - y^ x x = 0 j = 1, 2, n T 72. which i n matrix notation gives (5.8.7) \ 1 = y yx (5.8.8) z£ x . x x y T T Since these l a s t two equations imply (5.8.9) zj Z (5.8.10) Z x are eigenvectors of and y y = (x x)(y y) y T K T z£ x = ( x x ) ( y y ) x T K the common eigenvalue T T Z £ and T Z K respectively, with T T (x x ) ( y y ) . T Suppose, therefore, that i s an eigenvector of Z From (5.8.8) and the d e f i n i t i o n of Z F C K Z K . follows: T (5.8.11) Z K + 1 = ^ - x^ + 1 i t i i *k+l *k+l Taking S Q so that Z Q = A - S Q i s ( r e a l and) non-singular, i s p o s i t i v e d e f i n i t e and has n orthogonal s e t . Further, i f ^ Z Q T Z Q d i s t i n c t , r e a l eigenvectors forming an i s an eigenvector of Z^ Z£ with T eigenvalue C(^, then i t i s also an eigenvector of Z K + ^ with the 73. same eigenvalue. for This follows immediately from (5.8.11) since one has, i jf k + 1, \ l * i ' \ £ *1 - ^ * L l CZfc ^ *k+l ^+1 + ^) T T - (z. . k ^ V 7 i ^ i Z k x + T A *k+l ^+1 T W V J H L ) T ^+1 ; ^+1 T *k+l i X T since x^ =» 0 . Suppose were a basic set f o r the rank-reduction i t e r a t i o n (5.4) i . e . expansion (5.2) (ufcj" C^k} ^ { k] V = ( k] y terminated with i *) T t h e n n t h terms i n the sum f o r e i t e r a t i o n based on t h i s basic set would be optimal i n the sense of reducing the norm step. It i s easy t o show that ^ ' | x k orthogonal set of n a basic set; f o r l e t from (5.8.11), (b) similarly, (c) x? » 0 implies T x. Z i n = 0 f o r a l l i = 1, 2, .... n, hence ' i d e n t i c a l l y zero. at each be any vectors then we have: (a) It follows that f x£ x£ + 1 = 0 - 0 xT Z f c + 1 f o r k » 0, 1, 2, » 0 for (n-l) i / k + 1 Z n is The remainder of the development i s i d e n t i c a l with 74. that followed f o r (5.7.2) etc., i n the case of the Symmetric method, namely, n-1 A -o s + Z 0 " 0 S + ( Z 0 * n> ' 0 Z S + ZT (Zfc " Z k«0 k + 1 ) which by (5.8.11) gives (5.8.12) A = S S * . Q + H f X ^ +1 V l Making the s u b s t i t u t i o n V l » A - C^ 1 , as f o r the Symmetric case, gives the corresponding i t e r a t i o n of type (5.4) to be (5.8.13) W a Cfc • ° k ^ V l U A C k " ° A k ) V l Relation (5.8.13) defines the Quasi-optimum method. gives an inverse, A - 1 , in n The i t e r a t i o n steps with any orthogonal set of n vectors as the basic set, as i s c l e a r from the d e r i v a t i o n . In p a r t i c u l a r , i f the basic set i s the set of eigenvectors of T (A - S Q ) ( A - S Q ) , the i t e r a t i o n i s a c t u a l l y optimal i n the sense already described. Two s p e c i a l cases of the Quasi-optimum method are of i n t e r e s t . The f i r s t i s unusual i n that the basic se6 i s not constructed as an (ZJ orthogonal set at a l l but i s defined instead by » G o \ i + * 75. Since the derivation given for (5.8.12) does not, therefore, apply, the proof that Since Z H as thus defined i s really a basic set follows. ^Xfcj n-1 can be expressed as Z » £ T i=0 h ( \ Col. e Z T , we can write (5.8.11) K 1 + 1 1 i n the form: n-1 n-1 (5.9) Z h t l , > ^ e = £5) £ Ckd'T 1 + 1 1 C61.k+1 V [Col.K+1 V ] f-ji ^ CdL ej . i+1V i+1 T j + " 1 (Zk) T C C o l Now the (k + l ) ' s t column of 2^.^ summations with k+l ] ( (Zfc) °Cl C i s given by the term i n both i = k but this i s a zero column since the denominator of this term cancels i n the second summation leaving i t equal to the corresponding term of the f i r s t summation. Also, i f has a zero column, then (5.9) shows that Z^^ has one i n the same place. Hence Z n i s identically zero and (5.8.12) follows as before. Unfortunately, no analogue of (5.8.13) exists for this choice of the basic set and the iteration can only be performed as i n Wilf's method (mentioned earlier) i n which both the and the C^'s are successively computed using (5.4). The other Quasi-optimum iteration of interest i s the simple one with £ k | ° { k| * x e t h e urr ^ ' Sectors. (5.8.13) then assumes the form: 76. /, (5.10) c C C W (Cft.) £ T , k+l k+l = C + C o l _ Ro e c m h [Row r W k+1 With Col (A),T °k V k+1 CQ » I , the above iteration has two useful features. number of operations involved i n computing A" 1 First, the i s materially reduced by the presence of zeros i n determined locations as w i l l be shown later. Second, the iteration acquires certain symmetrical properties. A further property i s that the denominator w i l l remain close to unity i f the i n i t i a l CQ i s sufficiently good. These properties w i l l be gone into after the residuals, (I - AC^), associated with the Quasi-optimum iteration have been discussed. The Quasi-optimum method lacks an annoying feature of the Symmetric method. Even when the expansion (5.7.3) of A i s not defined, the associated iteration (5.7.6) does not, necessarily, break down. In this event, the Symmetric method produces a completely spurious "inverse". It i s , moreover, hard to guard against this possibility since the Z^'s are not actually computed. The expansion (5.8.12), on the other hand, i s necessarily defined and iteration (5.8,13) therefore gives a v a l i d inverse (in the absence of round-off) provided the iteration does not break down at any stage. Working with the "right-residual", R ( ) k • I - AC^, we have from r (5.8.13), for any Quasi-optimum method: (5.11.1) R ( K t l ) r a I - A C K - A ^ X k+l °k A X k+1 77. or T i - ^ L AC, (5.n. ) 2 . w whence, taking the natural norm H J (5.U.3) [j 4 Vl) x ) V k+1 AC r of (2.7) II II g I - \ X we have: lH l x k V l ( k ) r l l 2 2 T AC, x, . x. . However, the matrix ( I ) K K + J V l K A C idempotent and i t s eigenvalues + i 1 s e a s i l y v e r i f i e d to be V l k From (2.7) therefore a l l u n i t y . one can therefore write HR( (5.11.4) 4 k+l)l * llR( k ) r ll 2 which shows that f o r any Quasi-optimum i t e r a t i o n , the norms of the successive residuals are monotone non-increasing. the l e f t r e s i d u a l , (5.11.5) R R ( k + 1 ) £ ( i)£ = I - C ^ A k + « R(i) r or (5.11.2) R(i)£ f° and l e f t r e s i d u a l s . ^ k V which, with: one obtains k + fl - ° x k+l A C I f one now works with \ k R ( k ) t kV l ' shows that i f a zero column occurs i n e i t h e r r s o m e 1> i t i s retained by a l l subsequent right |tor the p a r t i c u l a r i t e r a t i o n (5.10) using the unit 78. vectors as the basic set, a stronger statement then (with can be made. One has P « ( p ^ ) « AC^)' (P ) (H k (5.U.6) [ 8 IL... -IL., (k+l)r (k)r "i»l ( k ) r ) T 1 ( ) k P k+1, This r e l a t i o n shows two new properties. i s zero f o r k = 0, 1, R.(k+l)r Rk are retained by subsequent k+1 F i r s t , the (k+l)st row of (n-l); R^ s, i >k. These l a s t r not only show how the r e s i d u a l diminishes as "symmetric p r o p e r t i e s " of i t e r a t i o n i n i t i a l approximation From CQ = I k—>n statements but lead to the (5.10) mentioned e a r l i e r when the i s used. (5.10) i t i s c l e a r that i n a l l rows a f t e r the k'th. second, zero rows of CQ = I must be i d e n t i c a l with One can,therefore, p a r t i t i o n C K i n the form where P k i s a. square matrix of order k and j) , i s a zero matrix (rectangular) and I the unit matrix of order (n - k)» previously, R( ) p a r t i t i o n e d as k r Since, as shown has zero rows f o r i t s f i r s t k rows, i t can be 79 Also, we can p a r t i t i o n where A compatibly as (j) stands f o r a zerb matrix of appropriate type. We then have the following k sub-matrix r e l a t i o n s as a r e s u l t of the equation A C k - 1 * (k)r (5.11.7) R 1 E P = I, E Q K K = -F, GP = -U , K K G^- + H - I - W k The f i r s t two equations give (5.11.8) P (5.11.9) 0^ » - E " F . -P F K = FT 1 1 FC from which follows the "symmetrical property" i n question, namely that for sypaetrical A, the square sub-matrices of order k, P , are symmetK r i c a l , being i n fact the inverses of successive square along the main diagonal of A. Since only sub-matrices 0^. and h a l f o f P K need t o be computed at each step, considerable savings i n the number of operations to obtain • A - 1 can be effected. The actual number of operations (excluding additions and m u l t i p l i c a t i o n s i n v o l v i n g quantities known t o be zero and m u l t i p l i c a t i o n s by q u a n t i t i e s known to be i n v e r s i o n by ( 5 . 1 0 ) with CQ = I i s : + l ) f o r the entire 80 multiplications r? - n divisions 2?. + n total n 3 additions: + T n^ - | n No allowance for possible symmetry of A 2 + 0(n) i s taken i n the above count. 3 If a general CQ i s used the operation count i s increased to multiplications, n 2 divisions and 2n3 + n 2n additions. A similar examination of the residuals can be made for the Symmetric iteration defined by (5.7.6). From (5.7.6) and the fact that R ( k ) r A » A R( ).g k one can write and (5.12.2) R . (K » E l)r » W /1 - ^ 1 "M* Since similar expressions can be written for R ( k ) g that both R( ]_)r occurring i n R j u ^ l = | k^> * e R k+ ( ) k h e y and (k+l)6 R ( )£ k reta: *- n a n y z e r o respectively. o n e c a n conclude rows and zero columns With the case for which vectors, (5.12.1) and the analogous expression 81. for (}£+l)£ show that the (k + l j s t row of R( ]_) R k+ column of R( +l)2 ^ * ^ence ^ an< k y t n e r ( + l) k previous remarks, a l l preceding columns and rows, respectively, as well) are necessarily zero. s - t 82 BIBLIOGRAPHY [I] JOHN, F, Advanced Numerical Analysis, New York University, 1956, p. 9. [2] von NEUMANN, J . and GOLDSTINE, H.H. Numerical Inverting of Matrices of High Order, B u l l e t i n of the American Mathematical Society, v o l . 53, no. 11, Nov. 1947, pp. 1033-1047. [3] STIEFEL, E. Uber einige Methoden der Relaxationsrechnung, Z e i t s c h r i f t f u r Angewandte Mathematik und Physik, v o l . 3, 1952, pp. 1-33. [4] LANCZOS, C. Solution of Systems of Linear Equations by Minimized Iterations, Journal of Research of the National Bureau of Standards, v o l . 49, 1952, PP. 33-53. [5] HESTENES, M.R. Iterative Methods f o r Solving Linear Equations, NAML Report 52-9, National Bureau of Standards, Los Angeles, 1951. [6] PIZZETTI, P. S u l l a Compensazione d e l l e Osservazioni Secondo i l Metodo dei Minimi Quadrati, Rendiconti d e l l a Reale Accademia Nazionale dei L i n c e i , Nota I, I I , Rome, (4) v o l . 3, 1887, pp. 230-235 and 288-293. [7] REICH, E. [8] OSTROWSKI, A. On the Convergence of C y c l i c Linear Iterations f o r Symmetric and Nearly Symmetric Matrices, I I , NBS Report 1759, National Bureau of Standards, Los Angeles, June 26, 1952. [9] MENDELSOHN, N.S. Some Elementary Properties of Ill-Conditioned Matrices and Linear Equations, American Mathematical Monthly, v o l . LXIII, no. 5, May 1956. [10] HOTELLING, H. P r a c t i c a l Problems of Matrix Calculation, Proceedings of the Berkeley Symposium on Mathematical S t a t i s t i c s and P r o b a b i l i t y , 1946, pp. 275-293. [II] LONSETH, A.T. An Extension of an Algorithm of Hotelling, Proceedings of the Berkeley Symposium on Mathematical S t a t i s t i c s and P r o b a b i l i t y , 1946. On the Convergence of the C l a s s i c a l Iterative Method of Solving Linear Simultaneous Equations, Ann. Math. S t a t i s t . , v o l . 20, 1949, pp. 448-451. 83. [12] FORSYTHE, G.S. Solving Linear Algebraic Equations Can Be Interesting, Bulletin of the American Mathematical Society, v o l . 59, no. 4, July 1953, pp. 310-313. [13] HOUSEHOLDER, A.S. Principles of Numerical Analysis, McGrawH i l l Inc., 1953, pp. 47, 48. [14] HOUSEHOLDER, A.S. Principles of Numerical Analysis, McGrawH i l l Inc., 1953, p. 79. [15] BERGER, W.J. and SAIBEL, E. Power Series Inversion of the Leontieff Matrix, Econometrica, v o l . 2 5 , no. 1, Jan. 1957, pp. 154-165. [16] WILF, H.S. Matrix Inversion by the Annihilation of Rank, SIAM Journal, vol. 7, no. 2 , pp. 149-151.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Iterative algorithms for the inversion of matrices...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Iterative algorithms for the inversion of matrices on digital computers Harris, Arthur Dorian Shaw 1960
pdf
Page Metadata
Item Metadata
Title | Iterative algorithms for the inversion of matrices on digital computers |
Creator |
Harris, Arthur Dorian Shaw |
Publisher | University of British Columbia |
Date Issued | 1960 |
Description | After a general discussion of matrix norms and digital operations, matrix inversion procedures based on power series expansions are examined. The general class of methods of which the Diagonal and Gauss-Seidel iterations are illustrative is studied in some detail with bounds for the error matrix being obtained assuming, both exact and digital operations. The concept of the condition of a matrix and its bearing on iterative inversion procedures is looked into. A similar derivation and examination is then made for Hotelling's algorithm. Hotelling's iteration is further examined with a view to modifying it. Higher-order formulae are obtained and criticized and a new variation of the algorithm called the Optimized Hotelling method is derived and commented on. Some schemes for constructing initial approximations in connection with Hotelling's iteration (and similar methods) are discussed and a new modification of a procedure proposed by Berger and Saibel is constructed. The final part of the thesis discusses a class of finite-step iterative inversions based on an identity of Householder's. Three members of the class, namely Jordan-type Completion, the Symmetric method and the Quasi-optimum method are defined and briefly discussed. The Quasi-optimum method is then examined in further detail and some of its properties derived for the special case with the unit matrix for an initial approximation. |
Subject |
Calculators Matrices |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-12-14 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080608 |
URI | http://hdl.handle.net/2429/39666 |
Degree |
Master of Arts - MA |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1960_A8 H25 I8.pdf [ 4.16MB ]
- Metadata
- JSON: 831-1.0080608.json
- JSON-LD: 831-1.0080608-ld.json
- RDF/XML (Pretty): 831-1.0080608-rdf.xml
- RDF/JSON: 831-1.0080608-rdf.json
- Turtle: 831-1.0080608-turtle.txt
- N-Triples: 831-1.0080608-rdf-ntriples.txt
- Original Record: 831-1.0080608-source.json
- Full Text
- 831-1.0080608-fulltext.txt
- Citation
- 831-1.0080608.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080608/manifest