Working with the Basis Inverse over a Sequence of Iterations

Working with the Basis Inverse over a Sequence of Iterations 1 1 Equations Involving the Basis Matrix At each iteration of the simplex method, we have a basis consisting of an index of variables: B(1), . . . , B(m) , from which we form the basis matrix B by collecting the columns AB(1) , . . . , AB(m) of A into a matrix: B :=
        AB(1)  AB(2)  . . .  AB(m−1)  AB(m)  . In order to execute the simplex method at each iteration, we need to be able to compute: x = B −1 r1 and/or pT = r2T B −1 , (1) for iteration-specific vectors r1 and r2 , which is to say that we need to solve equation systems of the type: Bx = r1 and/or pT B = r2T (2) for x and p. 2 LU Factorization One way to solve (2) is to factor B into the product of a lower and upper triangular matrix L, U : B = LU , and then compute x and/or p as follows. To compute x, we solve the following two systems by back substitution: • First solve Lv = r1 for v • Then solve U x = v for x. 2 To compute p, we solve the following two systems by back substitution: • First solve uT U = r2T for u • Then solve pT L = uT for p. It is straightforward to verify that these procedures yield x and p that satisfy (2). If we compute according to these procedures, then: Bx = LU x = Lv = r1 3 pT B = pT LU = uT U = r2T . and Updating the Basis and its Inverse As the simplex method moves from one iteration to the next, the basis matrix B changes by one column. Without loss of generality, assume that the columns of A have been re-ordered so that B := [ A1 | . . . | Aj−1 | Aj | Aj+1 | . . . | Am ] at one iteration. At the next iteration we have a new basis matrix B̃ of the form: B̃ := [ A1 | . . . | Aj−1 | Ak | Aj+1 | . . . | Am ] . Here we see that column Aj has been replaced by column Ak in the new basis. Assume that at the previous iteration we have B and we have computed an LU factorization of B that allows us to solve equations involving B −1 . At the current iteration, we now have B̃ and we would like to solve equations involving B̃ −1 . Although one might think that we might have to compute ˜ that is not the case. Herein we describe how an LU factorization of B, the linear algebra of working with B̃ −1 is computed in practice. Before we describe the method, we first need to digress a bit to discuss rank-1 matrices and rank-1 updates of the inverse of a matrix. 3 3.1 Rank-1 Matrices Consider the following matrix:   −2 2 0 −3  −4  4 0 −6  . W =  −14 14 0 −21  10 −10 0 15 W is an example of rank-1 matrix. All rows are linearly dependent and all columns are linearly dependent. Now define:   1  2   u=  7  −5 and v T = ( −2 2 0 −3 ) . If we think of u and v as n × 1 matrices, then notice that it makes sense to write:  1  2 W = uv T =   7 −5   −2 2   −4 4  × ( −2 2 0 −3 ) =    −14 14 10 −10 0 0 0 0  −3 −6   . −21  15 In fact, we can write any rank-1 matrix as uv T for suitable vectors u and v. 3.2 Rank-1 Update of a Matrix Inverse Suppose we have a matrix M and we have computed its inverse M −1 . Now consider the matrix M̃ := M + uv T ˜ −1 for some rank-1 matrix W = uv T . Then there is an exact formula for M −1 based on the data M , u, and v, which is called the Sherman-Morrison formula: ˜ is invertible if and only if v T M −1 u = −1, in which case Property. M ˜ −1 = I − M M −1 uv T M −1 . 1 + v T M −1 u 4 (3) Proof: Let Q= I− M −1 uv T M −1 . 1 + v T M −1 u ˜ Q = I, which we now compute: Then it suffices to show that M

˜ Q = M + uv T × I − M

= M + uv T × M −1 − = I + uv T M −1 − M −1 uv T 1+v T M −1 u M −1 uv T M −1 1+v T M −1 u uv T M −1 1+v T M −1 u = I + uv T M −1 1 −  − 1 1+v T M −1 u M −1  uv T M −1 uv T M −1 1+v T M −1 u − v T M −1 u 1+v T M −1 u =I q.e.d. 3.3 ˜ using M −1 Solving Equations with M Suppose that we have a convenient way to solve equations of the form M x = b (for example, if we have computed an LU factorization of M ), but that we want to solve the equation system: M̃ x = b . Using (3), we can write: ˜ −1 b = I − x=M M −1 uv T M −1 b . 1 + v T M −1 u Now notice in this expression that we only need to work with M −1 , which we presume that we can do conveniently. In fact, if we let x1 = M −1 b and x2 = M −1 u , we can write the above as: ˜ −1 b = I − x=M M −1 uv T x2 v T v T x1 −1 1 1 M b = I − x = x − x2 . 1 + v T M −1 u 1 + v T x2 1 + v T x2 5 ˜ x = b: Therefore we have the following procedure for solving M • Solve the system M x1 = b for x1 • Solve the system M x2 = u for x2 • Compute x = x1 − 3.4 v T x1 x2 . 1+v T x2 Computational Efficiency The number of operations needed to form an LU factorization of an n × n matrix M is on the order of n3 . Once the factorization has been computed, the number of operations it takes to then solve M x = b using back substitution by solving Lv = b and U x = v is on the order of n2 . If we solve ˜ x = b by factorizing M ˜ and then doing back substitution, the number M of operations needed would therefore be n3 + n2 . However, if we use the above rank-1 update method, the number of operations is n2 operations for each solve step and then 3n operations for the final step, yielding a total operation count of 2n2 + 3n. This is vastly superior to n3 + n2 for large n. 3.5 Application to the Simplex Method Returning to the simplex method, recall that we presume that the current basis is: B := [ A1 | . . . | Aj−1 | Aj | Aj+1 | . . . | Am ] at one iteration, and at the next iteration we have a new basis matrix B̃ of the form: B̃ := [ A1 | . . . | Aj−1 | Ak | Aj+1 | . . . | Am ] . Now notice that we can write: ˜ = B + (Ak − Aj ) × ej B T , where ej is the j th unit vector (ej has a 1 in the j th component and a 0 in every other component). This means that B̃ is a rank-1 update of B with u = (Ak − Aj ) and v = ej 6 . (4) If we wish to solve the equation system B̃x = r1 , we can apply the method of the previous section, substituting M = B, b = r1 , u = (Ak − Aj ) and  j v = e . This works out to: • Solve the system Bx1 = r1 for x1 • Solve the system Bx2 = Ak − Aj for x2 • Compute x = x1 − (ej )T x1 x2 . 1+(ej )T x2 This is fine if we want to update the basis only once. In practice, however, we would like to systematically apply this method over a sequence of iterations of the simplex method. Before we indicate how this can be done, we need to do a bit more algebraic manipulation. Notice that using (3) and (4) we can write:
˜ −1 = I − B  = I− B −1 uv T 1+v T B −1 u  B −1 (Ak −Aj )(ej ) B −1 T  1+(ej )T B −1 (Ak −Aj ) B −1 . Now notice that because Aj = Bej , it follows that B −1 Aj = ej , and substituting this in the above yields:  −1 j j ˜ −1 = I − (B ATk −e )(e ) B j −1 (e ) B T Ak  B −1 ˜ −1 = EB where ˜= I− E  −1   T B Ak − ej ej . (ej )T B −1 Ak Furthermore, if we let w ˜ be the solution of the system B w ˜ = Ak , that is, ˜ as ˜ = B −1 Ak , then we can write E w   T  ˜ − ej ej ˜= I− w E (ej )T w̃ We state this formally as: 7 . ˜ is obtained by replacing the j th Property A. Suppose that the basis B column of B with the new column Ak . Let w̃ be the solution of the system B w̃ = Ak and define:   T  ˜ − ej ej ˜= I− w E (ej )T w̃ Then . ˜ −1 . ˜ −1 = EB B (5) ˜ And then we have Once we have computed w ˜ we can easily form E. from above: ˜ −1 r1 = EB ˜ −1 r1 . x=B Using this we can construct a slightly different (but equivalent) method for solving B̃x = r1 : • Solve the system B w ˜ = Ak for w ˜  ˜ j )(ej ) ˜ = I − (w−e • Form and save the matrix E j T T  (e ) w̃ • Solve the system Bx1 = r1 for x1 ˜ 1. • Compute x = Ex Notice that      ˜ E=     1 1 .. .  c̃1 c̃2 . .. c̃j .. . .. c̃m where  c̃ = w ˜ − ej          . 1  (ej )T w̃ . Ẽ is an elementary matrix, which is matrix that differs from the identity matrix in only one column or row. To construct Ẽ we only need to solve 8 ˜ is the n-vector w Bw ˜ = Ak , and that the information needed to create E ˜ and the index j. Therefore the amount of memory needed to store Ẽ is just n + 1 numbers. Also the computation of Ẽx1 involves only 2n operations if the code is written to take advantage of the very simple special structure of Ẽ. 4 Implementation over a Sequence of Iterations ˜ be the basis at this iteration. Now let us look at the third iteration. Let B̃ We have: B̃ := [ A1 | . . . | Ai−1 | Ai | Ai+1 | . . . | Am ] at the second iteration, and let us suppose that at the third iteration we ˜ is of the form: replace the column Ai with the column Al , and so B̃ ˜ := [ A | . . . | A B̃ 1 i−1 | Al | Ai+1 | . . . | Am ] . ˜ be the solution of the system Then using Property A above, let w̃ ˜ ˜ ˜ = Ai . Then Bw ˜˜ B ˜ −1 = E ˜ −1 (6) B̃ where  ˜˜ = I − E  i T  e  . ˜˜ (ei )T w ˜ − ei w̃ ˜˜ B ˜˜ −1 = E ˜˜ EB ˜ −1 = E ˜ −1 , and so: It then follows that B ˜˜ EB ˜˜ −1 = E ˜ −1 . B (7) ˜˜ by forming E ˜˜ ˜ and E Therefore we can easily solve equations involving B and working with the original LU factorization of B. This idea can be extended over a large sequence of pivots. We start with a basis B and we compute and store an LU factorization of B. Let our sequence of bases be B0 = B, B1 , . . . , Bk and suppose that we have computed matrices E1 , . . . , Ek with the property that (Bl )−1 = El El−1 · · · E1 B −1 9 , l = 1, . . . , k . Then to work with the next basis inverse Bk+1 we compute a new matrix Ek+1 and we write: (Bk+1 )−1 = Ek+1 Ek · · · E1 B −1 . This method of working with the basis inverse over a sequence of iterations eventually degrades due to accumulated roundoff error. In most simplex codes this method is used for K = 50 iterations in a row, and then the next basis is completely re-factorized from scratch. Then the process continues for another K iterations, etc. 10