38 ASSURING NUMERICAL STABILITY IN THE INFORMATICA 4/89 PROCESS OF MATRIX REFACTORIZATION WITHIN LINEAR PROGRAMMING PACKAGE ON PC Keywords: linear programming, HR matrix, matrix Jane2 Barle and Janez Grad factorization, supersparsity, microcomputers Ekonomska fakulteta Borisa Kidriča, Ljubljana ABSTRACT! One of the most challenging tasks of those who develop linear programming software is development of quick, efficient and reliable matrix refactor ization subroutine. The paper describes the implementation of this subroutine within the PC-LIP programming package, which we developed for the IBM-PC personal computers. ft major design criterium for PC-LIP was to combine storage economy with numerical stability. The former was achieved using data structures which exploit super-sparsity and the latter implementing state of the art algorithms for basis matrix refactor ization. These algorithms were combined with different tools for improving numerical stability. The resulting subroutine performs satisfactory even on a badly scaled data which are quite common in practice. ZAGOTAVLJANJE NUMERICNE STABILNOSTI MED REFAKTORIZACIJO BAZNE MATRIKE V PROGRAMSKEM PAKETU ZA LINEARNO PROGRAMIRANJE NA OSEBNEM RAČUNALNIKU: Razvoj hitrega, učinkovitega in zanesljivega podprograma za faktor izaci jo bazne matrike spada med najbolj zahtevne naloge pri izgradnji programske opreme za linearno programiranje. ¿Slanek podaja opis implementacije tega podprograma v okviru programskega paketa PC-LIP, ki smo ga razvili za IBM kompatibilne osebne računalnike. Pri načrtovanju programskega paketa PC-LIP je bil glavni cilj vskladitev ekonomične izrabe pomnilnika z numericno stabilnostjo. To je bilo doseženo predvsem z uporabo podatkovnih struktur, ki izrabljajo hiperrazpršenost, in najbolj učinkovitih sodobnih algoritmov za izvajanje refaktor izaci je bazne matrike. Ti algoritmi so bili kombinirani z različnimi postopki za zagotavljanje numericne stabilnosti. Razviti podprogram za refaktorizacijo je bil uspešen tudi na slabo pogojenih problemih, ki so v praksi dokaj pogosti. Introduction A major concern of those who develop linear programming software is how to produce efficient, reliable and numericaly stable computational procedures for solving large-scale problems. When microcomputer software is considered, the problem of fitting algorithms and data structures within the limited storage is also very important. Contemporary literature on computational linear programming offers a plethora of different methods for achieving these goals. Roughly speaking these algorithms and techniques can be divided into following groups: i) Data structures which are designed for exploiting sparsity in LP data. They can be also tailored for utilization of structure and distribution of nonzero elements contained in LP matrix. ii) Revised simplex algorithm with product form factorization of basis matrix. State of the art implementations of this method are based on LU factorization. iii) Subroutines for refactorization of basis matrix. They are designed for controlling size and accuracy of product form factorization of basis matrix during the LP solving process. Each of these groups offers a great choice of 39 more or less elaborate techniques or algorithms. Sometimes it is not practical to use only the most advanced methods. For example, it is often desirable to sacrifice some computational speed in favour of reliability or storage economy. But in any case refactorization subroutine have to do its job correctly. Refactorization subroutine is also very important part of PC-LIP package which we developed for the IBM-PC personal computers. Practical experiences on real life problems show as that this package is capable to solve even very badly scaled problems. We attribute such a performance mainly to the careful implementation of matrix refactorization subroutine. Our implementation can be described as a successful combination of several state of the art numerical methods with tools for controlling numerical stability. That is why description of our refactorization subroutine may be of interest. Analysing Structure of Basis Matrix Re-factorization subroutine starts with analytical phase, where the structure of matrix nonzero elements is analysed. In our implementation Hellerman and Rarick algorithm (Hellerman, Rarick, 1971) is used for this purpose. This is quite a famous algorithm which is de facto standard for analytical phase implementations. Results of this algorithms are row and column permutations which transform matrix into form of so called HR matrix. Every HR matrix can be, after suitable rearangement of rows and columns, represented in the form of block lower triangular matrix: F1 o F2 H1 H2 j J£ IL where F 1 (i = l,...,k) are square matrices and H1 are rectangular matrices. HR matrices are distinguished from general block lower triangular matrices by structure of matrices F1. These matrices are always nonsingular. When their dimension exceed 1*1 they are called bumps or external bumps and must have structure similar to one on a next picture, where it is symbol for element which must be different from zero and * symbol for element which can be both zero or nonzero. F1 = tt # *» # * **# * » **»#* * *»**» * Columns with, at least one nonzero above the main diagonal are called spikes. There are two rules concerning spikes: i) Nonspike columns within bump must have nonzero diagonal elements, ii) Last column within the bump is a spike having a nonzero uppermost element. Typical overall number of spikes is much smaller than the number of columns within the matrix. This is an important fact which can be exploited -for the economical storing of the factor i 2ed basis matrix. It is easy to prove that when product form factorization is formed, only those elementary matrices which correspond to spikes'must be actually computed and stored. Other matrices from the product can be replaced by pointers to the non-spike columns (Chvatal, 1983). It is easy to check that matrix B with described structure can be represented as a product of matrices having a following form: Is 0 F * H1 where F* and H' are situated in the same rows and Ip are unit It and columns as in matrix B. I5 matrices of dimension s and p respectively, is assumed that s and p are numbers of columns of matrix F1, Therefore to the left and to the right which is of dimension r (s+r+p »a .Bk m) , B = nl B Bc < 1 ) Adequate definition for this identity can Bi be "generalized product form of matrix B". B1 can be defined as generalized elementary matrices (ordinary elementary matrices, which are contained in product form factorization, can differ from identity matrix only in one column). For matrices structured in such a way the following factorization formula is valid: Is 0 Fi 0 H* !P I s C 0 F1 0 XP Is c Ir 0 Hi (2) This identity explains why elementary matrices within particular bump can be computed 40 completely independent -from other parts of matr i x . If submatrix F1 is a bump, then factorization (2) is called splitting the bump (Helgason, Kennington, 19B2). Main purpose for its usage is to reduce number of nonzero elements in the product form factorization of B. Fill in (creating of new nonzero elements) during the factorization of Bl is restricted to r rows which belong to external bump F1 . It is an improvement if compared with the usual product form factorization, where creation of new nonzero elements is possible in rows belonging to submatrix H1 as well. Experiences show that this approach saves computer storage in spite of some overhead which is necessary to store additional r elementary matrices (Hellerman, Rarick, 1971, Helgason, Kennington, 19B2). Numerical Phase of the Algorith» Our algorithm for the numerical phase of refactorization which includes splitting the bump is presented in the continuation. This algorithm is modification of recent algorithm (Helgason, Kennington, 19BH) in which we included additional techniques for assuring numerical stability. It was necessary due to the fact that in the mentioned algorithm well scaled matrix was assumed. This assumption is in general quite a realistic one because many mainframe programming packages use some procedures for automatic scaling of data prior to applying the revised simplex algorithm. For example, this is true when MPSX/370 is considered (Benichou et al, 1977). However, such kind of procedure is not included in the PC-LIP. Ue avoided this for the sake of storage economy. The use of scales for rows and columns requires additional storage and, what is more important, practicaly prevents employing of such data structures which take advantage of supersparsity. This is a characteristics of large scale problems which means that the number of distinct numerical values in the problem matrix is usually much smaller than the number of nonzero coefficients (Greenberg, 1976). In the PC-LIP supersparsity is exploited in a standard way: nonzero values within problem matrix are represented by pointers to the table of all distinct nonzero' values (Bar 1e, Brad, 19B7). When basis matrix is not well scaled, automatically or by means of proper problem formulation, preassigned pivot can appear to be too small and for this reason inadequate. Two cases, which must be treated differently ares I. Inadequate pivot is situated within the external bump. In such a case its corresponding column can be treated in a same way as a spike. This means that such columns are included in the process of "spike swapping" (Helgason, Kennington, 1980). In fact this procedure is a variant of partial pivoting which is restricted to spikes within the same bump. 2. Inadequate pivot belongs to column which is outside the bumps (column from the triangle part of the matrix). In this case the only solution is to permute this pivot to the right bottom of the matrix. Such pivots will be referred to as "unstable p i vats". In the continuation of the paper we describe our implementation, which includes handling of above cases. Algorithm S CProduct form factorization for a HR matrix including splitting the bump] Preassigned sequence of pivots is represented with vector C, consisting of column indices, and vector R, consisting of row indices. It is assumed that these sequences are the results of Hellerman and Rarick's algorithm. Othor information obtained with this algorithm can be included into R and C using the fallowing method: indices of spike columns are stored in C with opposite (negative) sign, as well as components of R where external bumps are beginning. Algorithm's input is also basis matrix B, which is of dimension m and parameter TPIVR ("pivot relative tolerance"). All pivots yr, for which the inequality jyr| < TpIVR«yma>( holds, where ymax is the largest absolute value of available pivots, are counted as inadequate. Typical values for TPIVR are 0.001 or 0.01. SO: [Divide the pivots into stable and unstable] a) Set m, go to SI. c) If Rj<0, go to SO g). d) Set (i ) r = Rj ( i i) , 1 = Cj (iii) ymax = mg* |Bkl| e) If |Brj| < TPIVR*ymax, set (i) for i = i..... n-1: Rj = Rj + 1 snd cj " cj + l ( i i ) Rn = r and 0^ = 1 (iii) n = n-1 f) Set i = i+1 and go to SO b> g> Set (i ) t =■ number of columns in this external bump ( i i ) i = i + t hi Go to SO b) 41 SI J [Initialize] a) Set (i) i = 1 (pivot counter) (i i > 1 = Cj (pivot column index) (iii) r = |Rj| (pivot row index) (iv > j = 1 (counter for- elementary matrices > b> Alocate storage for (i) y — real vector of dimension m (ii) ETA—fi1e (data structure for storing matrices Ej>. c) Set y = B^j (l»" column of the matr i x B). d> If RjiO, go to S5. S2: [Obtain new lower triangular elementary matr i x] a) Set z = y. Vector z is the rtH column of the elementary matrix Ej . b) Set j = j+1. S3: CTest for the last stable pivot] a) If i=n, go to S16. < i i ) (iii) 1 = lRi| If i 0, go to S2. S5: [Initialize for Bump] a) Set (i) d = j (current length of ETA file) (ii) p = i (first column in this external bump) (iii) b = number of columns in this external bump (iv) t = 1+b-l (last column in this external bump) b) Set k - p. c) If k>t, go to S6. d) If Ck<0, set k = k+1 and go to S5 c ) e) Set (i) u = Ck ( i i > v = Rk (iii> Ymax = |Bqu| f) If |Svu| < TPIVR*ymax, set Ck = - Ck g) Set k = k+1 and go to S5 c). S6: [Obtain new lower triangular elementary matrix] a) Set fyt< ' f°r k = lRsl Set j = j + 1. S7: CTest for end of bump] a) If i = t, go to SI 1. b ) Set ( i ) i = i + 1 SB: CTest for spike] If ci>0, set y = B#J and go to S6. S9: CSpike update] Solve system of linear equations Ed---Ej-lV = B*1 S10: [Swap spikes if |yr| < TPIVR»ymax] a) Compute ymax = max |yk| for k £ CR j,...,Rt> b) If |yr| 1 TPIVR*ymax, go to S6. c) Obtain new l; for which 1 = |CS| Of set i = i — 1 and go to S12. b) If C^<0, solve system of equations Ed-"Ej-lY " B*1 and go to Sll. S14: [Obtain new lower triangular elementary matrix] a) Set il , for k = r yk , for k = 0 , otherwise Vector z is the r»n column of the elementary matrix Ej. b) Set j = j+1 S15: CTest for end of bump] a) If i=t, go to S3. b) If i i =* i+1 ( ii ) 1 = |Cj| (iii) r = |R j j (iv) y = B»! c) Go to SI*». S16: CPartial pivoting] a) If m=n, go to S17. 42 b) Set i = n+1 c) Sort columns having indices from Cj to Cm in such a way that the number of their nonzero elements is increasing. d) In the submatrix containing the rows from Rj to Rm and the columns from C^ to Cm perform Gaussian elimination with partial pivoting. S17 5 CEnd3 Product form of B including splitting the bump is obtained. At the termination of algorithm S, matrix B is represented as a product of elementary matrices E After obtaining this new factorisation of B, its accuracy must be tested. State of the art method for doing this is the so called Aird-Lynch estimate (Rice, 1985). If this estimate shows that C 3) is not accurate enough, algorithm S must be repeated with larger value of TPIVR. In our implementation of algorithm S within the PC-LIP, old value of TPIVR is multiplied by 10. Partial pivoting within step S10 c) can be performed by using subroutine BTRAN, which can be restricted to only those elementary matrices which belong to the current bump (Saunders, 1976). BTRAN (Backward TRANsformation) is a historical name for the subroutine which solves systems of the form BTy = d, where B is in a product form. The systems of the form Bz = u, which appear within several steps of the S algorithm, can be solved by using subroutine FTRAN (Forward TRANsformation). BTRAN and FTRAN are also important subroutines within the revised simplex algorithm. If row and column permutations determined by R and C are taken into account, all matrices Ej (i=1,2,...,j-1) are either upper triangular or lower triangular, but they are intermixed. That is why (3) is not LU factorization of B. For this reason the revised simplex algorithm with ordinary product form of basis matrix must be applied after performing the refactorization (as it is in PC-LIP). It is not possible to use those algorithms which use and maintain LU format of the basis matrix, for example Forrest and Tomlin method (Forrest, Tomlin, 1972). If one wishes to use LU format of the basis matrix, splitting the bump can be used only partially on an overall bump or kernel. Kernel is that part of HR matrix which is obtained after the lower and upper triangles have been removed from the matrix. Recently an algorithm was proposed (Helgason, Kennington, 19B2) which performs splitting the bump while maintaining the LU format. Ule briefly sketch how our methods for handling the unstable columns can be incorporated in this algorithm. By rearangement of rows and columns, the HR matrix may be placed in the following formj U V W o L T 0 M N L and U are lower and upper triangular matrix respectively, 0 are zero matrices of suitable dimensions. We use T instead of 0 which is used in the mentioned algorithm (Helgason, Kennington, 1992). This enables transfer of nonstable columns to the rightmost part of the matrix. It is easy to check that for matrix B* the following factorization is valid: I 0 o U V u 0 L T * 0 I 0 O M N 0 0 I LU factorization can be performed in usual way for the first matrix at the right hand side. The second matrix is already upper triangular. Due to the fact that a product of two upper triangular matrices is also an upper triangular matrix, LU factorization of B' is obtained. Conclusion« The matrix refactorization subroutine as described in the paper has been included in the PC-LIP linear programming software package. Our main contribution was that we have combined the already known methods for "splitting the bump" with some methods for assuring numerical stability. The algorithm was tested on many real life problems and proved to be stable even on a very badly scaled data. Algorithm satisfies also with respect to the computational speed. Unfortunately we have not yet had an opportunity for comparising it» performance with some other algorithm performance. It is possible however to measure the amount of reinversion computational time in overall run time. Another interesting test is to examine the effect of inversion frequency on the solution time. We performed these two test on a real life problem with 342 constraints, 454 structural variables and 204B nonzero elements. With the inversion frequency 20, the optimal solution was obtained after 221 iterations and 565.1 seconds of elapsed time. During the process 12 refactorizations were 43 performed in 149.7 seconds. This means that refactorizations amounts 26.49'/. to the overall computational time. We used the same problem for the test with the inversion -frequency 50. The effect of inversion frequency on solution t i me: Inver s ion f requency (i terat ions) Solution t i me (sees) I terat ions Time per i teration (sees) 20 50 ■ 565.1 544 221 224 2.55 2 .43 Results show that higher inversion frequency does not effect much the overall solution time. This can be explained with relatively slow execution of the product form variant of revised simplex method. With the use of the Forrest-Tomlin method the performance could be slightly improved (Ashford, Daniel, 1988). References 1. Ashford R.W., R.C. Daniels " A note on evaluating LP software for personal computers", European Journal of Operations Research,' 35(1989), pp. 160-164. 2. Benichou M., J.M. Gauthier, G. Hentges, G. Ribiere: "The efficient solution of large-scale linear programming problems - some algorithmic techniques and computational results", Mathematical Programming, 13(1977), pp. 280-322. 3. Chvatal V.: Linear Programming, New York — San Francisco, W.H. Freeman and Company 1983. 4. Forrest J.J.H., Tomlin J.A.: "Updated triangular factors of the basis to maintain sparsity in the product form simplex method", Mathematical Programming, 2(1972), pp. 263-278. 5. Greenberg H.J.s "A Tutorial on Matricial Packing", Design and Implementation of Optimization software, Urbino (Italy), (Ed. Greenberg H.J.), Alphen aan den Rijn (Netherlands), Sijthoff and Nordhoff 1978, pp. 109-142. 6. Helgason R.V., Kennington J.L.: "Spike swapping in basis reinversion", Naval Research Logistics Quarterly, 27(1980), pp. 697-701. 7. Helgason R.V., Kennington J.L.s "A note on splitting the bump in an elimination factorization", Naval Research Logistics Quarterly, 29(1982), pp. 169-178. 8. Hellerman E., Rarick D. s "Reinversion with the preassigned pivot procedure", Mathematical Programming, 1(1971), pp. 195-216. 9. Rice J.R.s Numerical Methods, Software, and Analysis, New York, McGraw-Hill 1985. 10. Saunders M.A.s "A fast, stable implementation of the simplex method using Bartels-Golub updating", Sparse Matrix Computations, (Eds. Bunch J.R., Rose D.J.), New York, Academic Press 1976, pp. 213-226. 11. Tomlin J.A.s "On scaling linear programming problems", Mathematical Programming Study, 4(1975), pp. 146-166. 12. Bar 1e J., Grad J.: "PC-LIP: A Microcomputer Linear Programming Package", (program description), Ljubljana, 1987.