Homework 6. First draft due Mar 30.
Final submission due Apr 6. There is no class Mar 16 or Mar 23, but I
will read email I am away.
This homework is long, but it should be manageable,
as long as you are following what is happening in class and
becoming reasonably familiar with Matlab.
If you find you are stuck and wasting too much time on it,
PLEASE SEND ME EMAIL. We can arrange a time to talk when I return.
The original Matlab program simplx.m uses the "\" operation to solve the
necessary systems of linear equations from scratch at every iteration,
using Gauss elimination (which is equivalent to LU factorization).
This does not take advantage of the fact that the basis matrix B is changing
only by one column at each iteration.
Modify simplx.m to implement three EFFICIENT versions
of the primal simplex method, as follows. In all cases assume that
b > 0 so a primal feasible point is available (original variables x=0,
slack variables w=b) and that the initial basis matrix B is therefore
the identity matrix.
The FIRST method to be implemented is explicit updating of the inverse
of the basis matrix, using the Sherman-Morrison-Woodbury (SMW)
formula. All you have to do is to change the code to remove the "\"
operations, used to get either xB and yN or Delta-xB and Delta-yN (depending
on which version you are using), and
replace these with multiplications by the inverse, say "Binv"
and its transpose. Do NOT compute all of "Binv*N", which is not needed.
The "Binv" matrix is then updated using the SMW formula, as explained in
class (see also p.127 of the text).
The SECOND method to be implemented is
the "product form of the inverse", avoiding both explicit computation of
the inverse and LU factorizations. In this method, we simply keep track of
all the E matrices (the "eta file") since the beginning of the simplex method,
which means simply keeping track of the "delta-xB" vectors and an associated
integer for each (see p.130 of the text). Each "solve" operation to get
xB and yN or Delta-xB and Delta-yN requires a sequence of "solves" with
the E matrices, each of which is trivially done using the formula on p.129,
as also explained in class.
The THIRD method to be implemented uses the LU factorization of the
basis matrix. The initial basis matrix is B=I which has L=I, U=I, and the
method starts by doing exactly the same as the second method, namely just
keeping track of the "eta file". The difference is that after a certain
number, say K, iterations, the "eta file" is abandoned and the LU factorization
of the current basis matrix is computed from scratch by simply calling
the Matlab "lu" function: use "[L,U] = lu(B)": this produces factors L and U
such that B = LU, with U upper triangular
and L lower triangular with permutations.
We don't want K to be too small because it's too expensive to compute
L and U factors frequently, but a larger K may pay off.
When the LU factorization is computed, the equation for
xB (or Delta-xB) requires only solving equations with L and U, like this:
xB = L\(U\b). Since L and U are triangular (or triangular with permutations),
these are solved very efficiently by Matlab, as long as it recognizes
the structure: to test this, experiment with "flops" (type "help flops"); in
my experience both on workstations and PC's, Matlab does recognize that
L and U can be solved efficiently and so does NOT do Gauss elimination all
over again, as it would do with x = B\b where B is not triangular.
Likewise the system for yN involves solves with the transpose matrices
L' and U'. Subsequent iterations require repeated solves with these
triangular matrices together with the eta matrix operations using the
"eta-file" saved since the most recent LU factorization was done. A new
LU factorization is then computed after K more iterations, etc.
This can be controlled in the main loop by a statement beginning
"if mod(iter,K) == 0".
There is another important aspect to the THIRD method: it must work efficiently
when the input matrix A is stored using Matlab's sparse data structures.
Type "help sparse" for more information on this.
Most Matlab operations, including both the LU factorization and the
solve operations using "\", work with both dense and sparse data.
As far as I can see, the only addition you have to make to this program
to accommodate sparse data is that, if A is sparse,
you must use speye(m) (sparse identity) instead of eye(m) when initializing
the basis matrix B. You can check whether A is sparse by "if issparse(A)..."
(type "help issparse").
When you are writing your programs remember to use PLENTY OF COMMENTS,
and use other good programming habits, such as indenting loops sensibly.
using meaningful variable names, and generally writing readable code.
The grader will be reading the source code.
Start by debugging your programs on the small examples you have used
before. Continue until you are sure the programs are working correctly.
If you are using a version of Matlab with the optimization toolbox installed
(which is available at Courant but probably not outside), you can check
your answers by calling an interface for the Matlab LP solver:
lpmatlab.m .
Otherwise you will have to check them against simplx.m, remembering that
my original version produces a permuted solution without keeping track
of the permutations (you were supposed to have fixed this).
Once you are sure your programs are working correctly, you should suppress all
output by adding semicolons and turn your script files into function
subprograms (type "help function" for information on how to do this).
These functions should have the calling sequence
[x,y] = LPname(A,b,c,maxit), and they should all return the optimal
original vectors x and y (often called x_tilde and y_tilde in class)
for the primal-dual pair
max c'x subject to Ax <= b, x >= 0
min b'y subject to A'y >= c, y >= 0
with the zero variables set accordingly and all vectors in correctly
permuted order. You don't need to return the slack variables w and z
since these are easily computed by w = b - Ax, z = A'y - c
if they are needed. Your functions should generate user-friendly error messages
if any of the following errors occur: the data dimensions don't match, the
iteration count exceeds maxit, the initial b has a negative entry or the
problem is found to be unbounded. See "help error".
Once everything is working correctly, test your programs on the following
examples:
rand('state',0); m = 40; n = 50;
A = full(sprand(m,n,0.2)) + 10*eye(m,n); b = rand(m,1); c = rand(n,1);
This example generates a sparse matrix A with about 20% nonzero entries,
adding 10 to the diagonal entries, and passing it to your codes as a full
matrix, i.e. NOT taking account of the sparsity. Do "help sprand" and
"help full" for more
information. The first statement initializes the random number generator
so we are all using the same example (assuming you are using Matlab 5).
rand('state',0); m = 40; n = 50;
A = sprand(m,n,0.2) + 10*speye(m,n); b = rand(m,1); c = rand(n,1);
This generates the SAME example except that now A is passed to your codes
as a sparse matrix. This should make the third method far MORE efficient,
since the LU factorization will take advantage of the sparsity in the matrix,
as will the "solves" using the \ operator. Don't forget that your code must
use eye, not speye if the data is sparse. You can also try passing the
sparse data to the codes for the first and second methods. What happens?
rand('state',0); m = 100; n = 200;
A = sprand(m,n,0.2) + 10*speye(m,n); b = rand(m,1); c = rand(n,1);
A = hilb(8); b = A*ones(8,1); c = b;
This example has been chosen to be very difficult numerically. The Hilbert
matrix A (type "help hilb") is notoriously ill-conditioned, meaning that
the solution of Ax=b and Ap=q will have very different solutions x and p
even if b is nearly equal to q. Another way of saying this is that the
norm of the inverse of A is huge, or that A is nearly singular. The LP
has been set up so the solution is x = y = A\b = A'\c = ones(8,1).
Be sure to use "format long" so you can see all digits of your computed
solution. How do your codes do on this difficult problem?
Compare the results of your codes with each other and with either lpmatlab.m
or the fixed-up simplx.m on these problems in a nice table,
showing number of ITERATIONS (should be same for all methods),
FLOP count (do "help flops"; remember
"flops(0)" resets the flop count), and ACCURACY (comparing with lpmatlab.m
or the fixed-up simplx.m for all problems but the last, and comparing
with the exact solution of all ones in the last case).
EXPERIMENT with various choices for K, which should affect BOTH the efficiency
AND the accuracy (in the case of the Hilbert LP) of the third method.
Try values of K in the range 5 to m for the random problems and 1 to 8
for the Hilbert LP.
You should also experiment with some other dimensions for the examples,
including at least one other in your table.
Finally, PROVE that the solution to the Hilbert LP is indeed x and y all ones.
Hint: show that the optimality conditions hold for these x and y, by
showing feasibility and complementarity.