Homework 8, due May 4.

Email me or phone me if you have problems.

Implement several EFFICIENT versions of the primal dual interior point algorithm given in Figure 17.1 on p. 274 of the text. The difference from Homework 7 is that we will reduce the size of the linear system to be solved at each step, and that we will take advantage of sparsity.

The first step is to get Homework 7 fully debugged and tested, and enhanced as follows. There are 3 possible ways the algorithm might terminate: (1) the optimality conditions are satisfied within the given tolerance; (2) the norm of x, y, z or w exceeds a given user-defined parameter M, indicating that the LP is unbounded (this parameter could have the value 1.0e+6); (3) a maximum number of iterations is exceeded (this could have the value 50). If more than 50 iterations are required, it probably means there is a bug in the program, or that the tolerance is too small (it should be substantially bigger than the machine precision 1e-16), or that the line search parameter r is too small (it should be at least 0.9), or that the centering parameter eta (called delta in the text) is too big (it should be less than or equal to 0.5). The algorithm should be coded as a function whose parameters are: the termination tolerance; the bound M; the maximum number of iterations; the steplength parameter r; and the centering parameter eta. If your code is not working, test it on very small problems and display all the variables to see what is going on. You must get it debugged before proceeding.

Once your code is working and fully tested you are ready to go on to the more efficient implementations. There are four of these:

The first version continues to factor the big matrix J, but exploits its sparsity. Change the construction of J to

J = sparse(2*(n+m), 2*(n+m)); % initialize J as a zero SPARSE matrix
J(1:m,1:n) = A; % top left block (A must be sparse too)
J(1:m,m+n+1:m+n+m) = speye(m); % sparse identity block
J(m+n+1:m+n+n,1:n) = sparse(diag(z)); % sparse Z block

etc. Observe how many flops this saves.

The second version gets dx and dy from solving the "reduced KKT system" at the end of Section 18.1 of the text (1/3 of the way down p.286). The matrix on the left-hand side is symmetric but not positive definite. Store it as a sparse matrix, so as to exploit the sparsity in both A and the diagonal matrices. Solve the system using \ (which will invoke Gauss elimination). Then dw and dz are trivially obtained from dx and dy. Do NOT store the diagonal matrices X, Y etc. Store them as vectors x, y. Use .* and .\ to operate with these: thus XZe is the vector x.*z, X^{-1}e is obtained from x .\ ones(n,1), and X^{-1}Z is sparse(diag(x.\z)).

The third version gets dy from the primal normal equations (18.9) (also p.286). The matrix on the left-hand side is symmetric and positive definite (move the minus sign to the right-hand side), and should be solved using the Cholesky factorization (type "help chol"), which is more efficient than the LU factorization (Gauss elimination). Let's call the left-hand side matrix B. The Cholesky algorithm computes an upper triangular matrix R such that B = R'R. Then solving the system Bp = q is equivalent to solving R'R p = q, so you need to type R = chol(B), followed by p = R\(R'\q). After getting dy in this fashion, you get dx from (18.8). When computing the matrix B, make sure it is stored as a sparse matrix.

The fourth version gets dx from the dual normal equations (18.10) (p.287). Again the matrix on the left-hand side is symmetric and positive definite and should be factored using Cholesky factorization. Again make sure the matrix is stored as a sparse matrix. You then get dy from the equation above (18.10).

Test these on a small test problem before you go on to larger ones. If they don't generate almost exactly the same iterates as the original version using J, they are buggy, and you need to debug them before going further.

Run all 4 versions on the same test problems that you used in Homework 6, plus an additional test problem: add a dense column (say ones(m,1)) to the largest sparse example. Prepare a table as you did before, showing the number of ITERATIONS (should be same for all versions), the FLOP count, the ACCURACY of the final solution, and the approximate CONDITIONING of the final matrix which is factored (J for the first version, the KKT matrix for the second, and the matrices passed to "chol" for the third and fourth versions). The condition number is estimated by CONDEST, but do NOT do this inside the main loop and do NOT include this in the FLOP count, as computing the condition number is more work than solving the linear system. Just do it once at the end. Also compare the results with those for the simplex method in Homework 6. Prepare several tables experimenting with different values of the steplength parameter r (e.g. 0.9, 0.999 and one other) and eta (e.g. 0.25, 0.1 and at least one other). Summarize the results in a written discussion.

Important Note: I will expect you to be able to discuss the results at the oral final exam, so it is important to get both Homework 6 and this Homework working properly. Come to see me if you have trouble. Now for some good news: this is the last programming assignment.