Description

Standard form LP barrier method
In the following three exercises, you will implement a barrier method for solving the standard form LP
minimize cTx
subject to ,
with variable x ∈ Rn, where A ∈ Rm×n, with m < n. Throughout this exercise we will assume that A is full rank, and the sublevel sets are all
bounded. (If this is not the case, the centering problem is unbounded below.)
1. Centering step. Implement Newton’s method for solving the centering problem
minimize subject to Ax = b,
with variable x, given a strictly feasible starting point x0.
Your code should accept A, b, c, and x0, and return x⋆, the primal optimal point, ν⋆, a dual optimal point, and the number of Newton steps executed.
Plot λ2/2 versus iteration k, for various problem data and initial points, to verify that your implementation gives asymptotic quadratic convergence. As stopping criterion, you can use λ2/2 ≤ 10−6. Experiment with varying the algorithm parameters α and β, observing the effect on the total number of Newton steps required, for a fixed problem instance. Check that your computed x⋆ and ν⋆ (nearly) satisfy the KKT conditions.
To generate some random problem data (i.e., A, b, c, x0), we recommend the following approach. First, generate A randomly. (You might want to check that it has full rank.) Then generate a random positive vector x0, and take b = Ax0. (This ensures that x0 is strictly feasible.) The parameter c can be chosen randomly. To be sure the sublevel sets are bounded, you can add a row to A with all positive elements. If you want to be able to repeat a run with the same problem data, be sure to set the state for the uniform and normal random number generators.
• We recommend computing λ2 using the formula ). You don’t really need λ for anything; you can work with λ2 instead. (This is important for reasons described below.)
• There can be small numerical errors in the Newton step ∆xnt that you compute. When x is nearly optimal, the computed value of λ2, i.e., ), can actually be (slightly) negative. If you take the squareroot to get λ, you’ll get a complex number, and you’ll never recover. Moreover, your line search will never exit. However, this only happens when x is nearly optimal. So if you exit on the condition λ2/2 ≤ 10−6, everything will be fine, even when the computed value of λ2 is negative.
• For the line search, you must first multiply the step size t by β until x + t∆xnt is feasible (i.e., strictly positive). If you don’t, when you evaluate f you’ll be taking the logarithm of negative numbers, and you’ll never recover.
2. LP solver with strictly feasible starting point. Using the centering code from part (1), implement a barrier method to solve the standard form LP
minimize cTx
subject to ,
with variable x ∈ Rn, given a strictly feasible starting point x0. Your LP solver should take as argument A, b, c, and x0, and return x⋆.
You can terminate your barrier method when the duality gap, as measured by n/t, is smaller than 10−3. (If you make the tolerance much smaller, you might run into some numerical trouble.) Check your LP solver against the solution found by cvx, for several problem instances.
The comments in part (1) on how to generate random data hold here too.
Experiment with the parameter µ to see the effect on the number of Newton steps per centering step, and the total number of Newton steps required to solve the problem.
Plot the progress of the algorithm, for a problem instance with n = 500 and m = 100, showing duality gap (on a log scale) on the vertical axis, versus the cumulative total number of Newton steps (on a linear scale) on the horizontal axis.
Your algorithm should return a 2 × k matrix history, (where k is the total number of centering steps), whose first row contains the number of Newton steps required for each centering step, and whose second row shows the duality gap at the end of each centering step. In order to get a plot that looks like the ones in the book (e.g., figure 11.4, page 572), you should use the following code:
[xx, yy] = stairs(cumsum(history(1,:)),history(2,:)); semilogy(xx,yy);
3. LP solver. Using the code from part (2), implement a general standard form LP solver, that takes arguments A, b, c, determines (strict) feasibility, and returns an optimal point if the problem is (strictly) feasible.
You will need to implement a phase I method, that determines whether the problem is strictly feasible, and if so, finds a strictly feasible point, which can then be fed to the code from part (2). In fact, you can use the code from part (2) to implement the phase I method.
To find a strictly feasible initial point x0, we solve the phase I problem
minimize t
subject to
,
with variables x and t. If we can find a feasible (x,t), with t < 1, then x is strictly feasible for the original problem. The converse is also true, so the original LP is strictly feasible if and only if t⋆ < 1, where t⋆ is the optimal value of the phase I problem.
We can initialize x and t for the phase I problem with any x0 satisfying Ax0 = b, and
. (Here we can assume that min 0; otherwise x0 is already a strictly
feasible point, and we are done.) You can use a change of variable z = x + (t − 1)1 to transform the phase I problem into the form in part (2).
Check your LP solver against cvx on several numerical examples, including both feasible and infeasible instances.