revised 14
The simplex method
maximize cT x = z subject to Ax = b x 0,
Consider once again the standard equality-form linear program
Corresponding to any basis B is a unique tableau. Theorem 12.2 gave a formula for that tableau, which we can rewrite slightly as follows: z -
jN
(cj - y T Aj )xj = y T b A-1 Aj xj = xB , ^ B
jN
xB +
^ where the vector y solves the equation AT y = cB , and the vector xB (whose B components are the current values of the basic variables) is A-1 b. Recall B that N is a list of the nonbasic indices, the vector Aj is the jth column of the matrix A; as usual, AB is the basis matrix, and analogously, the vector cB has entries ci as the index i runs through the basis B. Using this notation, consider how an iteration of the simplex method proceeds. We begin by choosing an entering variable xk (where k N ) with positive reduced cost ck -y T Ak . To do this, we must first calculate the vector y. Step 1: Solve the equation AT y = cB for the vector y. B Step 2: Choose an entering index k N such that ck > y T Ak . The next step is the ratio test, for which we need the column vector whose components are the coefficients of the entering variable xk in the body of the tableau. We can see from the formula for the tableau that this vector, which we call d, is just A-1 Ak . B Step 3: Solve the equation AB d = Ak for the vector d. Step 4: Calculate xi ^ : i B, di > 0 , di and choose a leaving index r from among those indices i attaining this minimum. t min 68
Since the ratio test at the next iteration will involve the new values of the basic variables, it is helpful to calculate these values. Step 5: Update the current values of the variables, xk t, xB xB - td, ^ ^ ^ and then the basis: Replace r by k in the ordered list B. We write B B {k} \ {r}. This description of the simplex method, which we call the revised simplex method, gives a more accurate idea of the calculations behind computer implementations like CPLEX. One obvious advantage over the hand calculations we have used until now is that, in the revised simplex method, we only calculate those elements of the tableau that we need for the current iteration. To illustrate, consider the following simple example. We wish to solve the linear program maximize x1 subject to x1 - x2 1 x1 + x2 3 x1 , x2 0. After introducing slack variables, the vectors and matrix defining our standard equality-form problem are 1 0 1 -1 1 0 1 c = , A = , b= . 0 1 1 0 1 3 0 The revised simplex method then performs the following sequence of iterations. Iteration 1 We have an obvious initial basis corresponding to the two slack variables, x3 , x4 , so B = [3, 4],
69
and the corresponding values of these two basic variables are given by xB = ^ 1 3 .
The basis matrix and the objective vector indexed by the basis are AB = 1 0 0 1 , cB = 0 0 .
Step 1: We first compute the vector y. 1 0 0 1 y1 y2 = 0 0 y= 0 0 .
Step 2: To choose an entering index k, we note j cj 1 1 2 0 AT y j > 0 0
so we can choose k = 1.
Step 3: We next compute the vector d. 1 0 0 1 d3 d4 = 1 1 d3 d4 = 1 1 .
Step 4: To choose a leaving index by the ratio test, we note t min (i = 1 1 3 , , 3 1 4) = 1,
so we can choose r = 3.
Step 5: We update the basis and the current values of the basic variables: x1 t = 1 ^ x4 x4 - td4 = 2 ^ ^ B B {k} \ {r} = [1, 4].
70
Iteration 2 From the previous iteration, the new basis and current values of the basic variables are 1 B = [1, 4], xB = ^ . 2 The basis matrix and the objective vector indexed by the basis are AB = 1 0 1 1 , cB = 1 0 .
Step 1: We compute the vector y. 1 1 0 1 y1 y2 = 1 0 y= 1 0 .
Step 2: For the entering index k, note j cj 2 0 3 0 AT y j > -1 0
so we can choose k = 2.
Step 3: We compute the vector d. 1 0 1 1 d1 d4 = -1 1 d1 d4 = -1 2 .
Step 4: For the ratio test, note t min (i = - , 1 , 2 2 4) = 1,
so r = 4.
Step 5: Again we update the basis and the current values of the basic variables: x2 t = 1 ^ x1 x1 - td1 = 2 ^ ^ B B {k} \ {r} = [1, 2]. 71
Iteration 3 Now the new basis and current values the of basic variables are B = [1, 2], xB = ^ 2 1 .
The basis matrix and the objective vector indexed by the basis are AB = 1 -1 1 1 , cB = 1 0 .
Step 1: We compute the vector y. 1 1 -1 1 y1 y2 = 1 0 y= 1/2 1/2 .
Step 2: When we try to find an entering index k, we note j cj 3 0 4 0 AT y j 1/2 1/2
so the current solution is optimal.
To conclude, the optimal solution is x = [2, 1, 0, 0]T . One book-keeping matter to keep in mind in calculations like the one above is that the order of the indices in the basis B is important. We should fix the order of B at the beginning of each iteration, because the numerical linear algebra we do involves the matrix AB and the vectors cB , xB , d: each ^ of these vectors depend on the order of the list B. At first sight, on a small example like the one above, the revised simplex method seems to be more involved than the standard simplex method that we have considered until now. Indeed, for small examples using hand calculation, it is usually easier simply to use the standard method. However, practical implementations of the simplex method are invariably in the form of the revised simplex method. To appreciate why, we must first understand which steps in each iteration involve the most computational work. The two most substantial steps consist of solving the two linear systems involving the basis matrix AB : AT y = cB and AB d = Ak . B 72
If the constraint matrix A has m rows, these are m-by-m linear systems. As a general guideline, solving such a system by Gaussian elimination typically involves of the order of m3 /3 multiplications/divisions. Hence the computational overhead of solving these two systems typically dominates by far the remaining calculations in the algorithm. The revised simplex method is numerically efficient for several related reasons. Basis factorization. Practical implementations factor the basis matrix AB as a product of simple matrices, giving two advantages: The factorization is easy to update after each pivot; The two linear systems involving AB are easy to solve using the factorization. Indeed, the work per iteration even with dense matrices becomes a multiple of m2 , not m3 . Problem shape. Linear programs in practice (when written in standard equality form) typically involve many more variables than constraints. Hence the basis matrix AB is typically much smaller than the full constraint matrix A, again making linear systems involving AB relatively easy to solve. Sparsity. Linear programs in practice are highly structured. The resulting constraint matrices A have many zeros, again making linear systems involving AB easy to solve. In addition, sparse constraint matrices make the computation of y T aj for j N very cheap. To illustrate the first and last points, suppose 5 4 1 -1 -1 2 -1 4 4 , with inverse A-1 = 3 3 -1 2 -1 AB = B 2 2 -1 2 -1 1 1 -1 2
3 3 3 2 1
2 2 2 2 1
1 1 1 1 1
.
(Here blank entries in AB are zeroes.) Note that AB is quite sparse while its inverse is dense. So performing operations explicitly with A-1 is quite costly. B
73
However, AB can be expressed as the product of two sparse triangular matrices L and U (for lower and upper triangular respectively): 1 1 -1 -1 1 1 -1 , U = . -1 1 1 -1 L= -1 1 1 -1 -1 1 1 To solve the linear system AT y = U T LT y = cB , we can solve very cheaply B U T z = cB and then LT y = z, and to solve AB d = LU d = Ak , we can solve Lf = Ak and then U d = f . While L and U are very sparse, their inverses are not: 1 1 1 1 1 1 1 1 1 1 1 1 -1 1 1 1 , U -1 = 1 1 1 . L = 1 1 1 1 1 1 1 1 1 1 1 1 While this is a very special case, the ideas are not: the basis matrix will usually be very sparse, while its inverse is usually much more dense; there is usually a very sparse factorization of the basis matrix (possibly after permuting its rows and columns) into triangular (or other very simple) matrices for which linear systems can be cheaply solved; and such factorizations can be efficiently updated after a simplex iteration.
74