〕尺cx 一一ーー一でーー～ 当 &&”『！11-------1--111-----1--Jl--1-1-]--1-J---- 1 SLC013377 DRC- 2 T J VES Efficient Updating of the Basis Inverse in Linear Programming via Partitioning by Johannes Bisschop and Alexander Meeraus VggS AT16NAV, BAI Folk 'NI RZCo14-'IRUC'IIr'- - ANtl I)E"'I I OPME A G 211989 Technical Note No. 2 Research Project 670-24 July 1976 Development Research Center World Bank 1818 H Street, N.W. Washington, D.C. 20433 T 74 SLCO1 3377 J. Bisschop and A. Meeraus, Development Research Center, World Bank, Washington, D.C. EFFICIENT UPDATING OF THE BASIS INVERSE IN LINEAR PROGRAMMING VIA PARTITIONING. Matrix modifications, such as column interchanges, row inter- changes or column and row additions, are viewed as matrix augmentations. Primal simplex iterations are treated as a special case, and an extremely compact updating procedure based on partitioning methods is developed. The method is especially suited for large sparse linear programs that are solved entirely in core. The procedure is not related to any of the methods based on the LU decomposition or the product form of the inverse. The new method has the distinct advantage that the growth of additional nonzero elements is not related to the size of the problem but only to the number of iterations following reinversion. It is shown that the representa- tion of the updated inverse does not grow monotonically in size, and that it may actually contract during certain simplex iterations. Implementation of the new procedure is straightforward, and does not require any involved data manipulation activities. Computational requirements differ from iteration to iteration. There is every indication, however, that the method requires on the average more computations per simplex iteration than other existing methods. Stability properties of the new method are considered. Compari- sons are made on the basis of real-world problems using a commerical code designed for large sparse linear programming problems. Efficient Updating of the Basis Inverse in Linear Programming via Partitioning by Johannes Bisschop and Alexander Meeraus ABSTRACT: A compact and flexible updating procedure using matrix augmentation is developed. It is shown that the representatio of the updated inverse does not grow monotonically in size, and may actually decrease during certain simplex iterations. 1. Introduction. There are basically two methods which are being used for the inversion and updating of the basis matrix in linear programming. Until recently the product form of the inverse (PFI) was used predominantly because of its simplicity to implement. Since then the usage of triangular (LU) decomposition has become widely recognized because of the superiority of this form of the inverse over the original product form in terms of speed, compactness, and accuracy (see [1], [9]). Research in this area is still extensive, and several sparsity-exploiting variants of the basic LU decomposition are under study (see [3], [6], and [8]). In this paper we propose yet another method for updating the inverse basis using augmentation and partitioning. Our main interest lies in solving large sparse systems that are kept entirely in core. The proposed updating method compares very favorably to LU decomposition in terms of compactness and flexibility, and is relatively straightforward to implement. - 2- Although actual comparisons between codes have not been made, the computa- tions required in both methods are comparable. We start out showing how general modifications of a reinverted basis B can be handled using augmentation and partitioning. Primal simplex iterations are then viewed as a special case, and studied in detail. 2. Matrix Modifications The study of matrix modifications or matrix tearing dates back to the work of Kron [5]. The updating of a basis B in linear programming is an example of a matrix modification where several columns previously not in B are interchanged for columns of B . We will consider various modifications of the matrix B The basic ideas in the following theorem stem from the work of Kron t5], and can also be found in Rose and Bunch [7]. Theorem 1: The matrix B is mx m , U is mx p , V is p x m -1 and S is p x p . Assume that S is available. Let B B + USV be a modification of the matrix B . Then i) x is a solution of Ex b if and only if x solves the augmented system B _U_ x b V -S y o ii) x is a solution of xB = c if and only if x solves the augmented system B U [x y -1 = ' c o1 - 3- Proof: By construction, the vector x solves Bx = b if and only if it solves the augmented system -1 = . (2.1) V -S y o To solve (2.1) is equivalent to solving B U_ x b(2) [ -1 [= L] (2.2) y -S y o One verifies this by applying the non-singular transformation [I -S (2.3) p to both sides of (2.1) to obtain the system (2.2) . Similarly, the vector x solves xB = c if and only if it solves the augmented system x y [A-S- = [c o . (2.4) To solve (2.4) is equivalent to solving [x y] V= [c o . (2.5) One verifies this by applying the non-singular transformation m 1(2.6) -SV I I to both sides of (2.4) to obtain the system (2.5) . This concludes the proof of the theorem. -4- Corollary: If B B + USV , and S1 is not readily available, one may use the above theorem for B = B + (US)I V = B + UI (SV) p p and consider either of the matrices B _US an B U L:i and . ~ V -I S_V _- In actual applications the collorary to Theorem 1 is of interest as the inversion of S can be avoided when inverting the augmented matrix. If one adds the assumption that B-1 exists, Theorem 1 provides us with a representation of B Without loss of generality, assume that B B + UV . Then using Kron's partitioning scheme, x = B b can be computed via i) y - Q-1 VB-l b ii) x = B-1 (b - Uy) and x = cB can be computed via i) y = - c B-1 U Q-1 ii) x = (c - yV) B-1 where -1 Q = ( -I - VB U) p The following theorem concerning the expansion and contraction of matrices will be used in the sequel. Theorem 2 Consider the partitioned matrices A = [ ] and B = [, ] A 3 A 4B 3 B4 whAr1 -1 -1 where B =A , and A1 and B4 are assumed to exist. Then knowledge -5- of A 1, A2 , A3 , and Q1= (A 4- A3A1i A2)1 is sufficient to determine the inverse of the larger matrix A , and knowledge of B -1, B 2 , and B is sufficient to determine the inverse of the smaller matrix Al ' Proof: It is straightforward to verify that AB I whenever -1 -1 -1 -1 B = A +A A Q AA 1 1 1 2Q A3A1' B -' - B2 = 1 A A2Q -1 -1 B = - Q A3 A ' B = Q-1 -1. and Q = (A - A3 A1 A2) -1 -1 Then, by substitution, one can verify that A = B - B B B 1 1 2 4 3 This concludes the proof of the theorem. Limiting out study of modification to linear programming, there are four basic modifications of a basis matrix B that result in an augmented system of the form [ ] (2.7) IV -I 1. Column replacement. Let the column vector at of the m x m matrix B be interchanged for the column vector a. not yet in B . Then B = B + UV where the 1 x m vector V is a unit vector with the one in the ith position, and the m x 1 vector U contains the vector difference a. - a . 2. Row replacement. Let the row vector a of B be interchanged -6- for the row vector a . not yet in B . Then B B+ UV where the 1 x m vector V contains the vector difference a. - a. , and the m x 1 vector U is a unit vector with the one in the ith position. 3. Row and column addition. A modification adding a row and a column to the previous basis results immediately in an augmented system. 4. Row and column deletion. After a row and column have been added during one stage of the solution procedure, a subsequent deletion will result in a reduced system of the form (2.7) . It is of interest to point out that whenever a sequence of column interchanges has taken place and either the PFI or LU decom- position has been used, one cannot perform any of the modifications 2 , 3 and 4 , unless one reinverts. When partitioning is used, any of the modifications 1 through 4 can be done without the immediate need for reinversion. In the next section we will study column replacements as they occur in the revised simplex method. The same methodology, however, applies also to row replacements, and column/row additions or deletions. 3. Primal Simplex Iterations after Reinversion. It is assumed that the sparse matrix A representing the LP problem is stored in core, together with the most compact representation of the basis inverse B-1 . In [2] we show that a combination of the Hellerman-Rarick preassigned pivot procedure [4] and a recursive partitioning scheme result in a more compact inverse representation than either the PFI -7- or the LU decomposition. By Theorem 1 , finding an inverse representation for the modified matrix B = B + UV is equivalent to finding an inverse representation for the augmented system (2.7) . This involves the matrices B-1 , U , V , and Q-1 , where Q = (-I - V B-1 U) . The size of Q is bounded from above by the number of simplex iterations following the latest reinversion of the basis matrix. As we will see, the dimensions of Q may in some simplex iterations be unaltered, or even diminished. It is important to note that only Q-1 needs to be stored explicitly, and that both U and V are stored with pointers. Recall that U contains vector differences involving columns of A , and that V contains unit vectors. The size of the matrix Q effectively determines the build-up of additional non-zero elements in the simplex method. It is therefore necessary to distinguish between four situations that may arise during simplex iterations following the latest reinversion of a basis B . Let the column vectors ai and a belong to B , while the column vectors a. and a do not belong to B . The subscripts i, p, j, and q are used to j* q distinguish between the vectors, and do not necessarily correspond to column positions. Case 1 : The non-B vector a is interchanged for the B vector a. The dimensions of the Q matrix associated with the augmented system (2.7) will increase by one. Case 2 : The non-B vector a is interchanged for the non-B vector a. which was entered as a nonbasic vector during a previous iteration in exchange for the B vector ai . The net effect is that the non-B vector a enters in exchange for the B vector ai . The dimensions of the Q matrix will be unaltered during such interchange of vectors not belonging to B . - 8- Case 3 The B vector ai , which left the basis during a previous iteration in exchange for the non-B vector aj , is re-entered into the basis in exchange for the B vector a . The net effect is that the non-B vector a. has entered in exchange for the B vector a . As in J P Case 2 , the dimensions of the Q matrix will be unaltered during such interchange of vectors belonging to B . Case 4 : The B vector a , which left the basis during a previous iteration in exchange for the non-B vector a , re-enters the basis in exchange for the non-B vector a. which has entered into the basis for the B vector ai in some other previous iteration. The net effect is that the non-B vector a is interchanged for the B vector a, , and the dimensions of the matrix Q will be decreased by one. Each of the above cases will be considered separately when discussing the modification of the Q matrix in a simplex iteration. Assume that k iterations have passed since the latest reinversion, and that a total of k < k non-B column vectors have been interchanged for B column vectors. The resulting augmented system is [B j k(3.1) 1 2 1 1 The mxZ matrix U = U - U , where U contains all columns that 2 have entered into the basis, and UT2 contains all columns that have left the basis. As part of an inductive argument (the case 9. = 1 is straight- -1 2 forward) we assume that V B U = I after k iterations. It will be shown for each of the above cases that this assumption will also be -9- satisfied after the k + 1 st iteration. Note that Q = (-I -V B-1 U ) -1 1 -1 2 = (- -V B U + VB U ) -1 = (-V B1 U ) At the end of the k th iteration after a reinversion during the revised simplex method we have the matrices/vectors B1 Q1 , X (the solution k~ k 1 vector), X k kB1 (a partially updated pricing vector), and a k(the entering B vector). Step 1 = -1 -1 k -1 k Compute w = Q (-V B a) ; (save B a ) y = B-1 (a - U w) = B1 a k B U1 w + B1 U2 w k z. In this step one needs to access the matrix B-1 twice when computing B-1 k -11k B a and B U w . In the event that a belongs to B , the matrix B-1 needs to be accessed only once. This occurs whenever a B vector -1 2 leaves and re-enters the basis. The matrix B U is a set of I Unit vectors where the location of the ones is known from the set of pointers determining U2 Step 2 Use x to determine the vector a , the i th column of B , that will 0 leave the basis. Establish which of the four cases discussed previously has occurred. k tbeoe k+1 k Update x to become x using the vector y . o o Step 3 Modify Q into Q (-V B U ) . This can be done in - 10 - all four cases with information created in steps 1 and 4. Each case will be discussed separately. Step 4 Let c denote the vector of cost coefficients associated with B the basic variables of iteration k . Then compute k+1 -1 ^-1 w = (-cB B U)Q k+l k+1 ^ -1 t (c - w V) k+1 -1 ^k+1 -1 ^1 k+1 -1 ^2 Here cB B U = c B U - cB B U , where -1 ~2 i) the matrix B U is a set of unit vectors where the location of the ones is known from the set of pointers deter- mining U2 k+i -1 k -1 k+l -1 ii) cB B = c B + d e. B ; B B 1 k k+l -1 = X d e k+1 with d denoting the difference between the cost coefficients associated with a and a respectively, and ei denoting a unit vector with the one in the i th position, where i is the number of the position,of the leaving vector a in B (the latest modified basis). In this step -1 k+1 -1 one needs to access the matrix B when computing (cB - w )B- and ei B-1 (a partial access). Some of the information used in steps 1 and 4 is also needed when Q Y is modified into Q . Each of the four cases discussed previously will be considered separately. Case 1. The non-B vector a. is interchanged for the B vector 3 ai . The matrix U is equal to the matrix U k plus an added column vector containing the difference a - a. The matrix V is equal to - 11 - the matrix Vz plus an added row vector containing an m-dimensional unit th vector e. with a one in the i position, where i is the number of the position of the leaving vector a. in B (the latest modified basis). As we will see later, any B vector a. in the basis always occupies its original position in B . By construction, the matrix Q = (-I - V B U) is such that V B U 2 +1 Here Q zqC R I' q q where qC = V B-1 a and [qR qI = ei B U1 . Both B-1 a. and ei B-1 are needed in each simplex iteration so that no extra sweep through B is required for the updating of Q Let = Q2 Q3 Q4 Then, using Theorem 2 -1 -1 C R -1 1 Qz + (1q)z q q Q ' -1C 2 =(-1/q) QZ q R -1 3 = (-q) q I R -1C and q = (q - q Q q ) Case 2. The non-B vector a is interchanged for the non-B vector a which occupies the jth column of U 1 and which was entered during a previous iteration in exchange for the B vector a.. The matrix U is equal to the 1 th^ difference a - a . The j Column of U - 12 - is a - a. in this case. The matrix V is exactly equal to V . The resulting matrix = - V B U ) haj the same dimensions as Q , and is such that V B U = V B-1 = = . Here Q = Q+Q1 Q2 where Q1 q-q (Zx 1) Q = e. (1 x i) , ^ -1 -1 with q = -V B a = -V B a q £ q q , the j th column of Q , e. , the (1 x R) unit vector with a one in the 3 jth position. Note that q is available from step 1 in the simplex method. Using Kron's tearing formula (5], one can compute -1 -1 -1 -1 -l -1 Q Z Q Q 1 -Q2 k 1 2 1 Note that (- -1 Q -1 -q (- 2 91 1) = - 2 1 - Q-1 + = (-1 - 2 k +1 (-Q2 Q ' and that - -1 -1 -1 -1-l -1qj -1 whr Q1 + -1Q2 -Z th where Q-1 q Q2 Q1 is the zero matrix with the exception of its j th -1 row which is equal to the j row of Q . - 13 - Case 3. The B vector a. , which left the basis during a previous 1 iteration in exchange for the non-B vector a (presently occupying the j th column of U1 ),is re-entered into the basis in exchange for the B vector a p * Pt The matrix U is equal to the matrix U with the exception of the j h column th^ of U containing the difference a - ai . The j column of U is a - a in this case. The matrix V is equal to the matrix V with the exception of th th the j row of V containing a unit vector with a one in the i position. The j,t row of V is a unit vector with a one in the p thposition in this case. The resulting matrix Q = (-I - V B U ) has the same dimensions as -1 ^2 and by construction, has the property that V B U 1 . Here Q = Q+Q1 Q2 where Q1 = e (x 1), Q2 = q q (1 x Z) -1 ^1 with q = -eB U , P q , the j th row of Q , ej , e , unit vectors oftiength k tflnd m with a one in the j and p position respectively. -1 Note that e B is used in step 4 in the simplex method. Using Kron's tearing formula [5], one can compute ^-1 -1 -1 -1 -1 -1 Q YQ + Q Q1 (--Q2 k 1 2 k ' Note that -1 -1 (-1- Q2 1 = (-1- (q - q ) Q Q1 -1 = (- - Qk Q1 + 1), = (-q Q 1 Q' - 14 - an ha -1 = -1 + - -1 1 T'-1 QjQ-1 Q-1 Q'-1 QqQ-1 andthtQ = QA1+( A1A\ 1qQE )\E 1j where -1 -1 th where Q1 Q1 q Q1 is the zero matrix with the exception of its j th -1 column, which is equal to the j column of Q Case 4. The B vector a , which left the basis during a previous iteration in exchange for the non-B vector a , re-enters the th basis in exchange for the non-B vector a (presently occupying the j column of U ) which was entered into the basis for the B vector a in some other previous iteration. Let P be a k x k permutation matrix which equals the identity -matrix with columns j and q interchanged. Note that P is symmetric, and therefore P1 P = P . Let R be an Z x k reduction matrix which th equals the identity matrix minus the one in the j column. Let U , then ^1 1 U = U R ^22 U2 = 2 PR V = RPVX; I = RI R The resulting matrix Q = (-I - V B U ) has its jth column and j th row equal to the zero vector, and contains a smaller matrix with dimensions (2-1) x (k-l) . One may verify this as follows: Q = -RI R - RP V B U R + RP V B UPR 2Z k.2.2 2. -1 1 -1 2 = R [-I k -PV X B U 9 +PV k B Uj P] R -1 1 = R [-I -PV B- U 1+P I P] R = R [-P V B-1 U ] R - RP Q Z R - 15 - This shows that no new information is needed to update the basis inverse in this case. For the sake of convenience, redefine Q to be the (Z-1) x (Z-1) matrix contained in R P Q R . To determine -1 -1 Q we use the larger inverse (P Q ) , and rely upon Theorem 2 Let q be the intersection of the j th row and j th column of (P Q)-1 which is equivalent to the intersection of the j th row and th -1 R q column of Q . Let the 1 x (Z-1) vector q be equal to the th -1 I j row of Q without the element q . Similarly, let the (k-1) x 1 C beqatoteqth -1l vector q be equal to the q column of Q without the element q -1 th th Also, let Q1 be the matrix Q without the j row and q column. Then, by Theorem 2 , -1I C R Q1- (1/q,) q q This completes the discussion for each of the four cases that can occur in the updating of the basis in the simplex method. Next we will consider some storage requirements associated with the above updating procedure. Let A be M x N , and let K denote the maximum size that -1 Q-1 is allowed to assume (usually K < 50) . Then, using previously developed notation, the basic core requirements involve storage for B-1 A (both in super-sparse form), an N-dimensional array for c (usually 2 -1 packed), K storage locations for Q , one K-dimensional working array and three M-dimensional working arrays to perform steps 1 through 4 of the simplex method. Total core requirements beyond that needed for B , A and c are therefore roughly K2 + K + 3M between every reinversion. The quantity K2 compares extremely favorable with the build- up of new nonzero elements in either the PFI or LU decomposition. The following table is extracted from the paper of Tomlin [9] in which he compares the two methods. - 16 - Table 1-: Growth of Nonzero Elements Iterations Problem A Problem B Problem C Partitioned After Latest (822 rows) (2978 rows) (3496 rows) Form of Reinversion PFI LU PFI LU PFI LU Updated Inverse* 10 4340 450 10807 972 3692 409 100 20 8572 868 21453 1792 22193 1975 400 30 12965 1118 33400 2971 40677 4357 900 40 17232 1874 45792 4443 59297 6113 1600 50 21582 2691 57903 5662 77935 8240 2500 */ It should be noted that core allocation is fixed in advance, so that every number in this column is actually the constant 2500 . In this table the partitioned form of the updated inverse is superior to either the PF1 or LU 'decomposition from a storage viewpoint. This storage advantage only adds to our previous observation that B-1 can be stored in a recursive form much more compactly than either the PFI or LU decomposition of B (see 12]). Comparing the LU decomposition and the partitioned form of the updated inverse from a computational viewpoint is not a straightforward matter. At the time of this writing, no actual comparisons on real-world problems using efficient implementations of the two methods have been made. In the LU decomposition one forward sweep, one backward sweep, and one partial backward sweep through the updated inverse is required (see [9]). The computations involved are comparable to steps 1 , 2 , and 4 as -1l- described previously where U. , V, Q1 are accessed twice, and B1 three times fully and one time partially. In the event that the incoming - 17 - vector is a B vector, this reduces to two complete sweeps and one partial sweep through B-1 . The updating of the Q matrix in the partitioned form is straightforward, and does not require many computations. This superficial analysis is an indicator that the two methods are not far apart when considering the computational cost for each simplex iteration. 4. Summary and Conclusions. In this paper we address the general problem of matrix modifications, and the specific problem of updating a basis inverse in linear programming. Our interest is in large sparse systems to be solved entirely in core. We develop a partitioned form of the updated inverse which we found to be more compact than either the PFI or the LU decomposition. This parti- tioned form involves an interesting mixture of sparse and full matrix technology. The method is flexible in that one is not locked into just column operations as is the case with either the PFI or the LU decom- position. Any of the four modifications discussed in Section 2 can be made in subsequent iterations. Relative to LU decomposition the method is straightforward to implement, and no excessive work needs to be done in the actual updating of the inverse. A thorough comparison on real- world problems using well implemented codes is needed to make any definite conclusions regarding computational efficiency. Realizing that the method has clear advantages over other methods in terms of compactness and implementation, we feel that the partitioned form of the updated inverse deserves definite attention when developing linear programming codes for large sparse problems. References [1] Bartels, R.H., G.H. Golub, "The Simplex Method of Linear Programming Using LU Decomposition", Comm. ACM. 12 (1969), pp. 266-268. [2] Bisschop, J. and A. Meeraus, "A Recursive Form of the Inverse of General Sparse Matrices", DRC Technical Note # 1 1976. 13] Goldfarb, D., "On the Bartels-Golub Decomposition for Linear Programming Bases", AERE Report CSS,Harwell, Didcot: Oxfordshire, 1975. [4] Hellerman, E., and D. Rarick, "Reinversion with the Preassigned Pivot Procedure", Math. Programming 2 (1971), pp. 195- 216. [5] Kron, G., Diakoptics, MacDonald: London, 1956. (61 Reid, J.K., "A Sparsity-Exploiting Variant of the Bartels-Golub Decomposition for Linear Programming Bases", AERE Report CSS, Harwell, Didcot: Oxfordshire, 1975. [7] Rose, D.J., and J.R. Bunch, "The Role of Partitioning in the Numerical Solution of Sparse Systems", Sparse Matrices and Their Applications, edited by D.J. Rose and R.A. Willoughby, Plenum Press: New York, 1972. [8] Saunders, M.A., "The Complexity of LU Updating in the Simplex Method", The Complexity of Computational Problem Solving, edited by R.S. Anderssen and R.P. Brent, University Press: Queensland, 1975. [9] Tomlin, J.A., "Modifying Triangular Factors of the Basis in the Simplex Method", Sparse Matrices and Their Applications, edited by D.J. Rose and R.A. Willoughby, Plenum Press: New York, 1972.