-ZJ9... -jwf' L I B RA I^Y OF THE U N IVER.5ITY Of I LLI NOIS 510.84 Iter no. 349-354 cop. Z Digitized by the Internet Archive in 2013 http://archive.org/details/erroranalysisofs351lerm yr\ Report No. 351 AN ERROR ANALYSIS OF SOLUTIONS TO SPARSE LINEAR PROGRAMMING PROBLEMS by Raymond J. Lermit September 17, I969 oct i \ m DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS Report No. 351 AN ERROR ANALYSIS OF SOLUTIONS TO SPARSE LINEAR PROGRAMMING PROBLEMS* By Raymond J. Lermit September 17, 1969 Department of Computer Science University of Illinois at Urbana-Champaign Urbana, Illinois 6l8oi Ibis work was supported in part by the Advanced Research Projects Agency Trc A^w^f ?, ^ the R ° me Mr Devel °P m ent Center under Contract No. US AP 30(602)41^4 and submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science, September, 1969. ii ABSTRACT This thesis discusses the round-off errors incurred in solving Linear Programming problems by the Revised Simplex Method, using the Product Form of the Inverse - the most usual method of solution. The matrices involved are assumed to be sparse as is usually the case. Also discussed are the reinversion of the basis matrix and iterative refinement of the solution obtained. iii ACKNOWLEDGEMENT The author wishes to express his sincere gratitude to his adviser and friend, Professor Robert S. Northcote, for his help in preparing this thesis. He would also like to thank his co-workers Serena Yardley and Teruji Yamamoto. A special word of thanks to Mrs. Kay Flessner for her patient work in typing a very difficult manuscript. iv TABLE OF CONTENTS Pag< 1. INTRODUCTION ± 2. NOTATION AND PRELIMINARY RESULTS 3 3- THE SIMPLEX METHOD 10 k. ERROR ANALYSIS OF THE PRODUCT FORM OF THE INVERSE METHOD 15 5- REINVERSION OF THE BASIS MATRIX 30 6. ERROR ANALYSIS OF THE REINVERSION ROUTINE 38 7- ITERATIVE REFINEMENT l+o 8. CONCLUDING REMARKS 55 LIST OF REFERENCES 57 APPENDIX 58 1 • INTRODUCTION The history of linear programming goes back to the discovery of the simplex algorithm in 1 9 4 7 by George E Dantzig, and its implementation on digital computers in the mid l 95 o»s. Little, however, has been written about the round-off errors produced during the computation, the most notable exception being a thesis by Bartels [l] . A linear programming problem may be defined as : minimize a linear objective function T z = c X subject to the constraints Ax = b, x > 0, where A is an m x n real matrix and x, b and c are real vectors of length n, m and n, respectively. Two approaches to error analysis are possible. The forward error a ^ al 7 s±s meth °d assumes that a solution x + Sx has been obtained and the size of 8x is then estimated. Alternatively, the backward error analysis approach assumes that the solution obtained is the exact solution of the slightly different problem: minimize: z + Sz = (c + 5c) T x subject to: (A + 5A) x = b + 5b, x > and estimates of the size of 5A, 5b and 5c are made. Even the most cursory glance at input data for a linear programming problem will convince the reader that the information is far from exact. For this reason backward error analysis is a much more realistic approach and is the one which will be pursued in this paper. If a bound can be given on the sizes of 5A, 5b and 5c, and these values are within tolerances determined by inaccuracies in the data, then the problem may be considered solved. This places on the user the responsibility of specifying the accuracy of his data. A feature of nearly all linear programming problems is that the matrix A is sparse, with very few, often less than one percent, of its elements differing from zero. This fact has been used almost from the outset to reduce the amount of computation and storage required simply by not storing the zeros and omitting calculations involving them. It is also possible to exploit this sparseness in the error analysis. Control of errors in linear programming has traditionally been done by "tolerances". Quantities that are "very small" are set to zero and numbers greater than some small negative constant are deemed to be positive. Such tolerances are usually set by "the system" but may be changed by the user at his peril. Efficient solution of problems seems to require these tolerances, and techniques for their automatic control have been developed (e.g., Clasen [2]). However, the final solution should be free of any such arti- ficialities and this is easily achieved by iteratively refining an approximate solution to purify it of its tainted past. Chapters 3 and k cover the simplex method and its backward error analysis respectively. Chapters 5 and 6 cover a reinversion routine that concentrates on retaining sparseness and providing error control. If the solution is not good enough after reinversion the refinement techniques of Chapter 7 may be used. Finally, Chapter 8 shows how slight imperfections in the solution may be removed by perturbing the original problem. 2. NOTATION AND PRELIMINARY RESULTS It is assumed that all operations are performed in normalized floating-point arithmetic. Define fl (E) to "be the result of evaluating the expression E. Following Wilkinson [9], the following laws in floating point arithmetic are assumed to hold on a nor- mal machine, providing overflow and underflow do not occur: fl (a + b) = a(l + & 1 ) + b(l + 5 ) fl (a - b) = a(l + 5 ) - b(l + 5^) fl (a X b) = (a x b) (l + 8 ) fl (a / b) = (a / b) (l + 5.) for any two machine numbers a and b, the quantities |&.|, which are relative errors, being bounded by some small constant. Two degrees of precision are utilized: normal precision, in which s il < £ i and extended precision where i' - 2 No assumptions are made about the relative sizes of 6 1 and £g. Normal pre- cision is used throughout, except for iterative refinement. Note, however, that many routines use "double precision" for all calculations. It is assumed that operations involving a zero operand are exact; i.e • ) fl (o + a) = fl (a - o) = a fl (o - a) = - a fl (o X a) -- fl (a Xo) = fl (o / a) = o (afo) In practice, however, careful programming should be used to avoid doing these operations at all, as they are redundant. Upper case letters are used for matrices, lower case letters are used for vectors and for the elements of both matrices and vectors. Where a sequence of matrices or vectors is used, superscripts are employed to dis- tinguish them and their elements are similarly labeled. Matrices and vectors of small terms introduced to account for rounding errors are prefixed by 5 i.e., 5B is a perturbation applied to the matrix B. All vectors are column vectors unless the context implies otherwise. Throughout the analysis, the l^ vector norm v I = max v. is used and A I = max I [ Ax n = max E I a. J=l ij for any n x n matrix A. Since sparse matrices are being used, it is conven- ient to define the function V by V (t) = 1 if t f o o if t = o for any scalar t. Thus the number of nonzero elements in a vector v = (v 1 , v 2 , ..., v ) is Lemma n ■ Z n V (V.) L=l V 2. If |0 | < € k = 1, 2, ..._, i then i n k=l u + y - 1 < i e' where e = iL-feJ- 1 ] and i < n (l) Proof : Let p = J (1+9, ) k=l > J (1 - €) = (1 - e) 1 k=l P < (1 + e) 1 (2) -l p-1 < (1 - e)~ -1 (3) Since [ (1 _ e) i/2 . (1 _ e) "i/2J 1 - (1 - e) 1 < (l - €) _i - l Hence, from (l), < 1 - p < (1 - e) 1 - 1 w Now (2) and (3) together yield -l p-i I < (1 - e) -l < H 1 - For small e, e' is only slightly larger than 6. It is often convenient to use e instead of e in error hounds. In the evaluation of the inner product of two dense vectors, it is necessary to estimate the error in calculating using floating-point arithmetic. Let Then • -0 s. . = f 1 (s. + x.y.) i = 1 2 n i+1 = S. (1 + 0.) + Xi y.(l +Ai ) : where |e | < e, |a. | < € J* 1 ~* and e = either e or e p = s * n+1 n = E x y (1+A f ft (i+e . ) i=l x k=i+l * s Using the lemma, = Z x.y. [ 1 + (n - i + 2) 7± ] i=l where \y A < e ' Z x.y. (1 + b ± ) , |B i l < (n + 1) €' i=l Now consider the case where x and y are sparse vectors. Let n V (> i' Sc = Z V Then = min (q^, q^) n and q' = Z V (x.y.) i=l q' < q If X., Y. are x., y., respectively, where x.y i is the j-th nonzero term in J J the sum, then n I Z x. y .) n/ Z x.y} q' = Z X.Y. (1 + A.) where | A . | <( q' +l)€ f J n Z x.y. (1 + 6.) . , 11 i' i=l where |5. < (q+l)e 1 l — In fact, 5. = A.; however, use of the less stringent bound is justified, since q' is hard to estimate but q is easily obtained. It is sometimes convenient to consider one of the vectors (say y) to be dense, in this case: = n so where 1 = "x ^i?i Xiy v = A x±yi u + 5i) |5.| < (^ + l) e< 10 3- THE SIMPLEX METHOD A brief outline of a version of the Revised Simplex algorithm is given; a much more detailed analysis may be found in Dantzig t3l or Hadley [6] The problem to be solved may be put into the form: minimize T z = c x subject to Ax = b x > o. Let B be the basis matrix [P , ..., P 1 where P y , . . ., P y is a 1 m 1 m set of m linearly independent columns of A. Then a basic solution is given by: * -1 x = B x b, r 1 T where x = L x , • • ■ , x j , 1 m since the inverse exists; this solution is also feasible if x > 0. The suc- cess of the Simplex method depends on the fact that, if an optimal solution exists, there is a basic feasible solution which yields the optimum value of the objective function. It is therefore necessary to consider only the basic feasible solutions in the search for optimality. An algorithm for solving the problem by the Revised Simplex Method is as follows : (a) Find an initial basic feasible solution, or indicate that none exists. 11 (t>) Change the basis until optimality is reached, which is achieved by a series of simplex iterations . The steps re- quired for each iteration, starting with some basis B = [P , P , . . . , p ] are : 1 2 m Step 1: Set y T = [^ , c y , ..., c y ] 12 m Step 2: Compute 7T = 7 T B _1 . T Step 3: Compute d. = c. -tt P. J J J for j = 1 (1) n, J f {v i , v 2 , ..., vj Choose an index s such that d = min {d. : d. < } s 3 J If no d. < 0, the solution is optimal. J Step k: Compute y = B~ P Step 5 : Find an index r such that x mm r y r in I'M If no y is > 0, then the solution is unbounded Choosing the minimum of x./y. assures that the new basic solution x' will be feasible. Step 6: Change the basis to B ' = t p v > "•> p v ' P s > P v ' '"' P v ] 1 r-1 r+1 m 12 Step 7: Compute the new solution, x' , using x" = (B* ) _1 b. Since three sets of equations axe solved using the matrix B, it is advantageous to have B~ J " available. If B " is known, then the inverse of the new basis (B' )" is easily obtained, since (B')' 1 = (B'^BB' 1 = (B" 1 B')" 1 B' 1 . B' is identical to B except that the r-th column, P y , is replaced by P g , so r B~ B* takes on the simple form: B^B' . o • o o o 'm -1 where y = [y 1 , y 2 > ' " ' y m ~ B P s' 13 Such a matrix is called an elementary column matrix , or simply an elementary matrix. It is easily inverted to give: (b"M " = -y 2 /y r V y v ~ y J y i say, and (B 1 )" 1 = T B _1 Glearly, if the starting basis is B , and IT, T , . .., T 5 are the elementary matrices produced during p iterations, then the inverse of the current basis B , say, is given by: " 1 = T 9 n^" 1 ... T 2 T 1 B _1 It is usual to take B = B = I, so that o o ' _1 = T* T 5 " 1 . . . T 2 T 1 . If a unit matrix cannot be found from the columns of A, then one must be manufactured using "artificial" columns. When the corresponding artificial variables are all removed (or reduced to zero) by subsequent iterations, a basic feasible solution has been found. 11+ When the matrix B ~ is not calculated explicitly, but only Or , P rpP" 1 , # . t t 1 are stored, the method is called the product form of the inverse. We will restrict our analysis to this method, since it is used more extensively than any other method in large scale LP algorithms. In order to reduce round off errors, storage requirements and running time, it is desirable from time to time to start afresh and recal- culate the inverse of the current basis. This technique, called re inversion , is discussed in Chapter 5. 15 k. ERROR ANALYSIS OF THE PRODUCT FORM OF THE INVERSE METHOD In this chapter the errors introduced during one iteration of the Simplex method are considered. Suppose that, after p iterations, the approximate inverse of the basis B is given in product form by: T P T 5 " 1 ... T 2 T 1 where each T J is an elementary matrix. Define an error matrix E by: P P T P_1 ... T 2 T 1 (B + E) = I. (1) Suppose also that, in performing one iteration, the r-th column of B is replaced by a non-basis column of A, A q = w say, giving a new basis matrix B'. A new elementary matrix T P+1 is formed which satisfies: T P+1 T P T P-1 m t ^ T 2 T l (B , + E , } = i (2) The round off errors in (l) and (2) are measured by ||e|| and ||e'||, respectively. The problem is to find a bound for | |e' | | in terms of | |e| | . To simplify the notation assume, without loss of generality, that the j-th elementary matrix, T°, is a unit matrix except for its j-th column and that in the j-th iteration, it is the j-th column of B that is replaced. Thus, pivoting occurs down the diagonal. This can be achieved by suitable interchange of rows and columns. In practice the actual pivots may be re- corded by saving the pivot column number corresponding to each component of the basic solution. Thus set the j-th elementary matrix to a unit matrix except for the j-th column which is 16 JlT mj, ti^, ..., ti^] Since T is obtained by inverting an elementary matrix whose j-th column is y say, define the vector T\ to be given exactly by i = l(l)m, i |= j i - 3 Whenever the T is used in calculation, the error introduced by the divisions will be included in the error analysis. Either of two different methods may be used in calculating T] , in each of which t = -y. and y. = -1 are set first. One method is then to divide J J each element of y J by t. The other is to multiply each element by t . The latter method will be faster if division takes longer than multiplication; it does, however, produce two round off errors for each element. It is assumed that the latter method is used. One step in the Simplex algorithm is to compute y = B " w for whi»-h the approximate solution is given by y = f 1 (T P T 9 ' 1 . . . T 2 T 1 w) (3) Setting 1 a = w 17 and calculating successively k+1 k k a* x = f l (tV) (k = 1, 2, ..., P ) (h) yields +1 y = a? Equation (k) is equivalent to fl (a +Tl i O^) ±=lf 2, ..., m ™ k+ l i fl ^k a k) i = k ^) Rearrangement of (5) to include the effect of computations on previous gives : i-1 fl (Z o£ T]*) , ±<£ (6) CC fl K + k£ a k *n±> > ± > ' In particular, setting 1 = £, yields °£ ■ fl < w i + Sl °{ Tlj>- (7) For convenience in notation, let p, = a . Then p i " fl < w i \£ p k ^i) (8) a a:. : 18 r. = «? 1 1 +1 fl ( Li \ 1? '2P (9) n (w i + L \ < ]i> * Now let L be the matrix L = -Tl 41 2 1 3 "'3 < 1 1 2 'm 'm Then (8) may he written in the form: L p = w from which 3 is uniquely given by 3 = L w. Using the relations given in Lemma 1 of the Appendix, it can be seen that 3 exactly satisfies (L+6L)3=w+Sw, where 19 where (I* '" l» y l < (^ + 3)^ «i = i z 1 v cn±) k=i 1 le w ± | < q . | WjL | e£. If A w = L 3 - w then Aw=§w-6L. p 11**11 < ||e w| | + || 5 l|| ||p|| "a l|y H +e l IMfl p. q TT = max / , i > p , i < p (15) and = the maximum number of nonzero elements in any- row of U. "u PS ( l u ijl )# 1>J If 6 P = (8 3 1 , 6 P 2 , ..., 6 ) and A w f = (A w , ..., A w ) then 8 3|| < | |6 7\\ and || A w« | | < || 6 7 \ The perturbation 6 P of 3 gives rise to a perturbation L 6 P in w , ..., Therefore, combining (13) and (15) w y = f 1 (T P T P_1 ... T 2 T 1 w) = T P T 5 " 1 ... T 2 T 1 (w + v) 22 where ||v|| < max {eJ_(q L l|w|| + q L (q L +3) mJ |3| l)+| N | | | 63| U I 67| |) (16) and L < Let stt To simplify (l6), let M = max (|n*|, 1) (17) i,k Then ^ < M, M^ < M. (18) v k q = max ^ V (T^) Then | |v| | < e| ( q | |w| | + q (q+3) M | |p| | ) + £{ qM || U _1 || 1 • 2 £[ q (q+3) M | |U — [q (q+3) M ||y|| + MpII + q||w||] (19) -li The only quantity in the above inequality which is unknown is j |u | \; fortunately approximations to the elements of U " are found during the cal- culation of the T\ vectors. Consider now the calculation of the T] vectors forming the approximate inverse of 3, and consider 3 to be a unit matrix except for the first p columns. If 23 and then ce- ll a 21 a- a 12 1 22 or a lp l 2p B = a a. pl "p2 a. 'PP a ml a m2 a mp k+1 a ±3 for i = 1(1) m , j = i(i) p i f k i = k k = l(l)j-l (20) fl (J d-i :=i < ii> a ij k5i fl (a 1 . i< J i> j (21) and the T\ vectors are given "by *nv fl {l/o?..) 33 if 1 t J 1 = 3 (22) Let r be the matrix 2k r :■■ 11 a a 12 2 22 1 a 13 • Q Q 23 3 33 p+1,1 a J+l,2 a: 'ml a m2 a a ip 2 2p cc x PP p+1, p 1 a mp Then k=i kj 'i (UT).. = < Z n a k . T] k + a 1 £1 k,n i i: ij i < P, J < P i > P, J < P i = J > P otherwise (23) But, for i < j < p L p, A P, J < P otherwise Hence fti'Jsld&^^i' K\) 27 :e n ' .£, la"?. I + e,' f la 1 1 .1^1 i.i '1 j^L ' ij Thus and Now < ^ [q (q+2) M | |g| | + 2 | |h| | + ||b||] Ml < e [ U (q+2) m ||g|| + 2 ||h|| + ||b||] (27) ||r|| < max { ||G||, ||B|| + 1} u" 1 = r (.1 + f)" 1 u'^k llrl (I + F)' 1 || < j |r | I provided ||f||< 1 (28) 28 Substituting the bound for | |u~" into equation (19) gives a bound for It has been shown that 2 _1 y = fl (T 9 T 9 ' 1 . . . T T 1 w) = t 5 T 9 ' 1 . • • T T (w + v) It is now possible to estimate |e'||, where E' satisfies r P+1 T 9 • • • T 2 T 1 (B f + E') = I. The new elementary matrix T 9 has as its p + 1 -st column the vector T| where f " y i / Vi I 1 /Vi i f" P+l i = p+l Thus p zeros )+l T^ 1 y = [ 0, . . . , 0, 1, 0, . . • , 0] = T 9 * 1 T P • • • T T 1 (w + v) since w is the p + 1 -st column of B f , v is the p+l - st column of E' 29 For all other columns B* and B are identical, and, since ^ . . . T 2 T 1 (B + E) = tP +1 I = ^ +1 , (B 1 + E') and (B + E) are the same except for column p+1; i.e., E' is E with the p+1 -st column replaced by v. Hence Ne'II <||e|| + ||v|| Finding a bound for ||v|| is not as difficult as equations (19) and (27) would seem to indicate since the only quantities required that are not recorded are the elements of G. (3 is just the last column of G) . These values are obtained during the calculation of y and may be used to update ||g|| from iteration to iteration. This bound depends heavily on: 1. The sparseness of the T) vectors; 2. The size of their elements. It is, therefore, very important to reduce the magnitude of the T) elements by controlling the size of 'pivots'. Consideration also should be given to maintaining spareseness. 30 5- REINVERSION OF THE BASIS MATRIX During the execution of the simplex algorithm, the choice of pivot is restricted at each iteration by the need to obtain a significant improvement in the value of the objective function and to maintain feasibility. Little can be done to control round-off errors and the sparseness of the inverse. In recalculating the inverse a flexible pivot selection agenda can be adopted, and the errors and density be thereby greatly reduced. Before discussing the general reinversion algorithm, the decomposi- tion of a sparse matrix A by two different methods will be considered. In both cases pivoting occurs down the diagonal. (l) Product form decomposition. 1 2 A = E X E . E 11 " 1 E n where ri E J = e l3 e 2j 1 e e . . 33 e • -i • nj 31 (2) Decomposing A into a product of lower and upper triangular matrices P 11 L = i 21 ^22 '31 ^32 ^33 (jnl ^2 ^n 3 1 nn and U = \2 \ 3 1 U 23 • ' "l» u 2n u h-l,n 1 so that A = LU. Having decomposed the matrix, inversion is easy, s nice is 6 )' 1 -e ./e . . ' V< jo i -e ./e . . . 1 where L" 1 = L n L n " X . . . L 2 L 1 32 L J = 1 • 1 -i i/i,. and similarly for U . Note that the pattern of zeros and nonzeros in E , L J , U J correspond exactly to that of the matrices from which they were derived. The density of the product form decomposition is at least as great as in the LU form, as can be shown thus : ^ = JJ >1,J 1 nj 'id 1 e. . , T J (say). 112? Then A = S T S T s n-l T n-1 s n ^ For i < j , T 1 S J = S J T 1 So the product of the S's and T's may be rearranged to give: A = S 1 S 2 o n-l n 1 2 S S T T T^ 1 T n . .1 „2 S ' S ' * ' • > s are all lower triangular matrices, so L = S 1 S 2 is lower triangular. 1 ? Similarly U = T T n 33 . . T is upper triangular, and LU is the (unique) decomposition into lower and upper triangular parts. L . S 1 S 2 '11 '21 1 'nl 6 23 'n2 .1 'n-l,n-l 'n n-1 1 'I nn '11 e e 21 e 22 e 31 e 32 e 33 e -, e _ e nl n2 n3 . e nn Hence L contains the same number of nonzero elements as in S 1 , S 2 , (excluding unit diagonals which need not be stored). 1 2 U = T T T n .'. U" 1 .^ 11 )- 1 (T 11 " 1 )- 1 . . . ( T 2 )-l ( T l)-1 ., S 3U -e -e In 2n -e -e ln-1 2n-l '12 -e n-ln "U ■ 1 " e l2 " e l3 • • • - e m 1 " e 23 • • • " e 2n 1 • • • ~> • '1 * n-l,n 1 - - -i O yi Therefore the number of nonzero elements in T , T , . . . , T is equal to that in u" (again excluding the unit diagonals). Since the density of U~ will he greater than that of U, due to "fill in" during the inversion procedure, it may be concluded that the decom- position into elementary matrices will, in general, give rise to more non- zero elements than will LU decomposition. (TltiS ignores the possibility of cancellation, which might produce zero elements from sums of nonzero ones; however, this should not often happen in practice). It is therefore desirable to use an LU decomposition method for reinversion, not only to reduce the rounding errors, but also to reduce the number of T| elements representing the inverse. It has been shown that a triangular matrix decomposes readily into elementary matrices. It is therefore desirable to sort the rows and columns of the basis to form a matrix which is as close to a triangular matrix as possible. The sorted basis is then partitioned as follows: 35 M N M N \ o 1 } X X XXX o o o XXX X xxx • xxx XXX xxx o o xxx xxx X xxx xxx xxx xxx xxx xxx x o X X w xxx xxx x o xxx xxx xxx \ o xxx xxx xxx * xxx xxx • xxx o \ In practice the matrix need not be physically sorted, all that is needed are the indices of the pivot positions. The NN "block contains logical vectors in the basis; these will produce no T] records. The RR and CC blocks are triangular- The MM block contains rows and columns which cannot be triangularized by sorting. LU decomposition is applied on this block. Rewrite the matrix, in partitioned form, as : " R S M C N D I 36 where (I and I are zero and identity matrices of appropriate sizes. This decomposes to: R I i — r— I T X M u X I s I N I c I I D I The first and last matrices are lower triangular, pivoting down the diagonal gives T\ vectors with no increase in the number of nonzero elements. If M = LU, L and U being lower and upper triangular respectively, then I M N I I r— I L I I X I u N I I The first matrix is lower triangular and the second can be transformed by row and column interchanges into: I I N I U 37 which is upper triangular. After inversion, the interchanges can be reversed. The only place where elements are added to the inverse is in the formation of L and U. Since the pivoting strategy within the M block may be employed, this may be used to minimize the number of nonzeros created. The exact form of the pivoting strategy depends on such factors as the number of columns of the matrix that can be stored in memory at one time. An important point is to chose the pivots dynamically, at each step choosing that pivot which minimizes the number of nonzeros created. It may be necessary to change the pivot selection method where choosing a particular pivot would result in very large elements in the repre- sentation of the inverse. Discussion of this point is deferred until the next chapter. 38 6. ERROR ANALYSIS OF THE REINVERSION ROUTINE In this chapter the computation of the inverse of the matrix B = R S M N c D I r s t columns columns columns using the product form i-l .2 „1 T 9 T 5 ...T T , p = r + s + t, is considered. It is shown that 2 „1 T 9 T 9 ' . . . T T (B + 6B) = I i.e., T^ . . • T is the exact inverse of some perturbed matrix B + 5B, and a bound for j | &b| | is established. As before B is split into the form 39 B = 1 R r— r- S I I 1 X I r ■ M N I I *R X X I | I c D I M Consider first the inversion of B 11 21 31 rl 22 32 r2 33 b -■ -, b -, n r+1 1 r+1 2 ml m2 rr . b r+1 r mr J by: The T\ vector for the j elementary matrix T ( j=l(l) r) is given T\ 3 . = fl (1/b. .) = JJ ' b..(l+ S ..) 33 33 n J = n (-b. . x r\j)= " b i.i ^ 1+6 i,j^ 10 J >a (^ , f or i > j Uo where I S - - 1 < e, > The terms b (1+6 ) may be regarded as the elements of a matrix; hence ij ij where . . T is the exact inverse of the perturbed matrix (B + SB^) R- ( *Vij u j = r + s + l(l)p i = j(l)m otherwise and hi^h This leaves only the matrix B M' I r~~~ M N I I 1+1 to be inverted. Since the elements of N are simply copied into the 7] vectors with a change in sign, the only error is that in inverting M. Therefore set, M=M 1 = in m 11 1 21 m s2 1 "12 m 22 m s2 1 m i s m 2 s m ss 12 s s+1 and compute a sequence of matrices M , M , . . . , M , M , by forming a new T\ vector from the k-th column of NT using only elements on or below the diagonal. Thus *£- (1+6^) kk J kk ; i fK-4^) = -4(i +8 ^ ) "4 (i+ 4> where Is, ik 1 <*{■ (i) The matrix M* i s then updated by this new elementary matrix giving k+1 v + -i M . If there were no rounding errors, the elements of M would be zero k+1 in the k-th column below the diagonal and el. Would be unity. If we set these to zero and one respectively then: 1+2 i.e. Thus if m. k+1 i. . fi cn. k \ ^Ki^i^) k m. . i = j = k i > k, j = k i = k, j > k 1 > k, j > k otherwise (2) k+1 m =< \ \, + \ \j B k! + T)*mJL +aJd Bid + 2 ^i Vi K* i = j = k i > k, j = k i = k, j > k ij m. iJ i > k, j > k otherwise (3) k k k k+1 Tl it and ET is the matrix with elements 6.., then M k+1 = lV + f? w where U3 k n k * k "kk^k S kk k k L k ~k k m ik 5 ik + ffi kk T1 i *kk k _k k d k k J k k _ k „k ,k i = j = k i > k, j = k i = k, j > k i > k, j > k otherwise (5) and 4's* 6 i*ls4> for all i, j, k. s+1 After s iterations M = U, (say), an upper triangular matrix which has ones down the diagonal. Let F ^ +1 - UV 1 E k - L k M* + L k F k = L k (t^ + P k ) k k But the elements of F are zero except for i > k and j > k (because those of E are) . L 1 F k = F k f or i < k .-. l*" 1 L k - 2 . . . L 2 L 1 ^ - F k kk then U = M s+1 . s s-1 = L L ^L 1 M 1 \ L (L $ I s " 1 . k=l v • L F ) L S L 3 " 1 . . . I 2 L 1 M 1 ^ (L S L 8 " 1 . . . L 2 L 1 F k ) L S L 3 " 1 . . • L 2 L 1 (M + F) (6) where If ■A **• U ^2 ^3 u 23 1 •l U 2s ii u 3s s-ls 1 u" 1 = u 2 U 3 . . . u 8 " 1 u s where "*5 V 1 = ■\: ■ U 2j -U. - . J-lJ from which U 2 U 3 . . . U 8 - 1 U S L S L 5 " 1 . . . L 2 L 1 L L (M+F) = I. A bound on | |f| | is now required. Since F* - (l k )" X ^ , k k m., 8., lk lk ij I k k id ^j k k _ „k k , ,k m. . 8. . + 2 T] . m, . (5 . i > k, j = k i = k, j > k 8 M> i > k, j > k otherwise (7) U6 An important point is that these error terms do not involve the 'pivot' element T|, Summing over k gives F. . j-l m J v / k k „ -k k ,„ ,k A k XN ij ij 1J KJ i > J < i-1 z k k _ ^k k / ,k k N N m ij B ij + k§i (m ij h3 + 2 \\* Kj- "^ i < J (8) Let P M a. . max { m. . max ( |T]_. I ) i,k 4i H i k=l l Z, V (m k .) i > k } Then |F. .1 < (cr. . + kq.) p M £' | |F| | < s (a + l+ q ) p M e| where q = max (q ) i = max ( cr . . ) 1,J ij This gives an a posteriori bound for | |f| | in terms of the sparse- ness and the size of the elements of the inverse and the intermediate matrices M 1 , V? , • • •, $ ' As might be expected, sparse matrices give much smaller errors. 1*7 Combining these results indicates that the computed inverse is the exact inverse, in product form, of the perturbed matrix B + 5B 1 R+6R S+8S M+F N C+8C ' D+6D I It is clear from equation (9) that both sparseness and size of the elements are crucial factors in bounding the error. The previous chapter dealt with preserving sparseness. Consideration now must be given to bounding the size of the elements in the T) vectors and the intermediate matrices that are formed. In forming an T) vector, start with some column vector T v = [v^, v , ..., v ] (say) and 'pivot' on some element v , giving m T[ = fl L-v /v , . . . , -v /v , 1/v , -v _/v , • . ., -v /v ] 1' r r-1' r ' r r+1' r v r where the value of r may be arbitrarily chosen. The usual method is to select the pivot element so that j v | is greater than some 'tolerance' value so that the elements of T] will not be 'too large'. Since, as equation (7) shows, the element l/v does not appear in any error term, the remaining terms may U8 be bounded by c (say) by insisting that r be chosen such that: |v r l > \ |v.| for 1-1, ..., r-1, r+1, ..., p It may be noted that for dense systems, one normally selects the largest element of v for this pivot, which is equivalent to taking c = 1. Thus the above scheme is a compromise between minimizing the number of nonzero elements and reducing their size. h 9 7. ITERATIVE REFINEMENT Suppose that an approximate inverse to the basis B and an approxi- mate solution vector x have been calculated. It is desired to calculate a more accurate solution using iterative refinement. A detailed analysis of iterative refinement using floating point arithmetic is given by Moler [7] • A simple adaption of that work is developed here to cover the sparse matrix case. Each iteration of the refinement requires the following three steps 1) Compute the residuals: r = b - Bx with x = x. 2) Compute the correction: d = B r . ^\ « ^ j_! j_- j> 21+I m .m 3) Add the correction to form: x = x + d . The step where most accuracy is required is in calculating the residuals. It is assumed that extended precision is used in this step. Computing d = B r is done using the approximate inverse to B. Let r m = b - B x m = c m W d m = (B + E 111 ) - 1 r m (2) x m+l= /i +d m +g m (3) mm m where c , E and g are correction terms due to round-off error and are intro- duced to make the equations hold exactly. A bound for ||E || is known, using the results of either Chapter h or Chapter 6. It is required that the inverse be reasonably good so that -1 I I N "I A\ < (Mb" 1 !!) Let F m = B' 1 E" 1 . Then |F m || = Mb' 1 ^\\ < Mb" 1 !! ||E m || + Mfli Estimates of | | c m | | and ||g m || are now required, m = fl (x m + d m ) - (x m + d m ) h = e ±\ + e [^ > ^r |e.|< ei ,|e.|< ei ui Is N <\(UAl + l|d ffi ||) When the process converges, ||d m || < | |x m | I Kg" 1 !! <2e x ||x m || <2e x (|| x || + ||x m -x||). m In calculating r , it is assumed that each nnrrmnnpnt J* assumed that each component r. = fl (b.-Z B. . x m ) is calculated using extended precision and the result then is rounded'to'nor- mal precision. This is denoted by fl (b - Z B x m ) 2 A i i.i .1'' ij J m m m c = r - b + Bx and 52 m c. = n (b. -Zb.. x m ) - (b. -Lb.. x m ) 2,1 v i ij j' i yd . (fi 2 (b t - Z B y *»)) (i + ej -(b. -Eb.. x p = (b.(l^ .) - Z By *" (1 + 7y)l (1 + 9^ - (b. - Z By X^) . (b, - Z By X") , ± where \e \ < e , |£ ± | < e, |7.J < (q + 2) e^ m q = max.L. V (B ) since ||c m || <€^ ||b - Bx m || + ( 1 +€]L ) ^ ||b|| + ( q + 2) (l + e x ) e< ||b|| ||x m || <€' ||B|| ||x m - x|| + (l + 6 X ) ^ ll B H l' X l + ( q + 2) (l + e x ) e 2 ||b|| (||x|| + ||x m - x||) i i nil i M / m \ i i | | x | | = | |x + (x - x)| I I I I I I I m II < |x| I + I |x - x| I . I < [€' + (q + 2) (i + e x ) €'] ||B|| + (q+ 3) (1 + e x ) € 2 H B " " X "' m i x - x m+1 - x| I < I |5B| -1 m i x - x 1 - | l&B + 1 1^11 INI M^M l - I IfiP'l I I |b" 1 1 (e£ + (q + 2) (1 + € 1 ) G 2 ) + (q+ 3) (i + e{) €' I m I x - x + 2 e x (IM m I x - x 53 SB | | | |B | | + k (B) (6' + (q + 2) ( 1 - SB -1] 1 + g 1 ) e£) + 2 O, < 1 1 - 1_ 0H 5h Therefore, provided a < 1, convergence of x to within the tolerance given by * 11-11. l -°i can be guaranteed. cr depends essentially on | |$b||, l B and the measure °f sparse- ness, q. It is essential that B be well conditioned and that ||6B|| be small, i e., an accurate inverse of B is known. The effect of a small value of q is to reduce the need for using double precision in calculating the residuals, since it is always multiplied by ^ 8. CONCLUDING REMARKS The question remains as to what constitutes an optimal feasible solution. Strictly speaking, it is one for which x > and the reduced costs d are also > o for j = l(l)n. In practice, small negative tolerances, 5 , 5 (say) are chosen so that the solution is considered optimal and feasible pro- vided x. > g for all i d j > S Q for all columns j not in the basis An alternative approach is as follows: Determine some sufficiently refined solution x and a basis B, such that Bx = b may be considered exact. If some x. are < 0, define 1 ' x = x + 6x > where Bx. CO , x.> \ X i ' X i < This produces a change in b, 5b (say), where 56 6b = B 6x The vector 5b can either be calculated directly, or its norm may be estimated from ||sb|| < ||B|| ||fix|| If the elements of 8b are sufficiently small that b + 5b is an acceptable right hand side, then x may be considered to be the required solution to the problem. Similarly, in calculating d. = c. -.A. n. a. ., n may be refined in the same way as x, so that the d. are accurately known. If d. < for some j, J J a small change in the corresponding c. will set them to zero, if this new J cost vector is acceptable, then the problem is solved. The techniques described in Chapters 3 through 6 and above enable the exact solution to the slightly perturbed problem; minimize T z + 5z = ( c + 6c) x subject to (A + gA) x = b + 6b x > 0. to be obtained. If | oA| | , | &t> | and ||5c|| are within acceptable limits, then the problem may be considered solved. If more accuracy is required, iterative refinement as described in Chapter 7 may be used. 57 LIST OF REFERENCES [1] Bartels, Richard H. , "A Numerical Investigation of the Simplex Method", Technical report no. CS 104, Stanford University, 1968. [2] Clasen, R. J. , "Techniques for Automatic Tolerance Control in Linear Programming", Comm. ACM, £ (November I966), pp. 802-803. [3] Dantzig, George B. , Linear Programming and Extensions , Princeton University Press, I965. [U] Dantzig, George B., Harvey, Roy P., McKnight, Robert D. , Smith, Stanley S. , "Sparse Matrix Techniques in two Mathematical Programming Codes", Technical report no. 69-I, Stanford University, 1969. [5] Forsythe, George E. and Moler, Cleve B. , Computer Solution of Linear Algebraic Systems , Prentice Hall, I967. [6] Hadley, G. , Linear Programming , Addison Wesley, I962. [?] Moler, Cleve B. , "iterative Refinement in Floating Point", JACM, 1^ (April I967), pp. 316-321. [8] Orchard-Hays, William, Advanced Linear Programming Techniques , McGraw Hill, 1968. [9] Wilkinson, J. H. , Rounding Errors in Algebraic Processes , Prentice Hall, 1963. 58 APPENDIX 59 Lemma 1 ; The solution of a triangular system of equations. Suppose that the components of a vector x are computed in the order \> x 2 f '" ' X n by x, = b 1 1 x. = i-1 fl (b. + JL r M x, ) i k=l ik k' i = 2, 3, ••• , n or, equivalently, by the set of equations ■ r 21 1 ■r -r 31 32 -r -r ml m2 -r ram-1 1 x~ L Xm - m i.e., Rx = b, is solved. Then x is the exact solution of the perturbed tri- angular system: (R + SR) x = b + sb, where, 6R = -r 21 8 21 ' r 31 5 31 " r 32 6 32 ml ml 0. mm-1 mm-1 60 and where 6 b = [b x ai , • • • , \ s;] T , Proof s: I < «# - i-l with ^ ■& V (r ik ) In Chapter 2, it was shown that n fl & vi } -A v± < 1+8 i> Bj. I < (^c + 1) € ' n and c^ -J^ V (x ± ) Similarly, if p = f 1 ( z +2^ x ± y ± ) is calculated from: S l = Z B. . = fl (s. + x i y i ) i = 1, 2, ... , n p = S n+1 then n n fl (z + i§i x i y i } = 2 (1 + ^) + j j l x.y. (1 4 6 .) where A | < q^ e« , I $! I < (%, + 1) €' and q x =.Z 1 V (x.) . Using this result, X i ■ fl (\ + ^l r iK *J " ^ (1 + ^ + li ** ** (1 + ^ where B ik I <(<1 1 + 1)€^ and 8. | < q. el • 6l 62 Hence i-1 X i = b i + \ 6 i \ll (r ik + r ik & ik> \ which completes the proof. Note also that i-1 ||5R|| ^max.^ | r. . 8.. < max (q max |r 8 | ) i x Kj