1 Introduction

Our standard linear programming problem is

Here, c and x are n vectors, b is an m vector, and A is an m×n matrix. The dual problem can be written

where y is an M-vector and s is the n-vector of dual slacks.

The primal affine scaling algorithm sometimes makes very slow progress when it gets close to a non-optimal corner of the feasible region. One possible remedy is to try explicitly to center the iterates. We can do this by moving in the direction

where X^k is the diagonal matrix with diagonal entries equal to the components of x^k and A^k = AX^k. Notice that this direction is trying to increase all components simultaneously. Recall the example where we found the primal affine scaling direction,

At the point x^k = (0.8, 0.1, 0.2, 1.9), the centering direction is

Current commercial implementations combine affine steps with centering and duality in predictor-corrector methods. In addition to excellent computational performance, variants of these methods converge in a number of iterations that is polynomial in the size of the problem.

Notation: We will let e = [1,…, 1]^T .

2 Log barrier function

In order to make the iterates more centered, we add a penalty term to the objective of (P), giving

(1)

where μ is a positive constant. Note that if x_i gets close to zero then the function f(x) gets large. This gives a nonlinear program

Given μ > 0, it can be shown that there is a unique optimal solution to the problem (P(μ)).

For problem (LP1), the barrier function is

With μ = 3, the value of the barrier function is 11.33 at the point x^k = (0.8, 0.1, 0.2, 1.9) and it is 13.32 at the point = (0.95, 1.95, 0.05, 0.05), close to the optimal solution of (LP1). Conversely, when μ = 0.1, the barrier function takes the values 1.54 and -2.26 at the points x^k = (0.8, 0.1, 0.2, 1.9) and = (0.95, 1.95, 0.05, 0.05), respectively. Thus, for larger choices of μ, the value f(x^k) of the more central point is smaller than that of f(), but for smaller choices of μ the results are reversed.

The motivation for this problem is that it should be easier to get an approximate optimal solution to (P(μ)) for a large value of μ, and then this solution can help in finding approximate optimal solutions for smaller values of μ. The limiting solution as μ → 0 is an optimal solution to (P).

The gradient of f(x) is c-μX^-1e. This eventually leads to a direction that combines the affine direction and a centering direction, as we will see later.

3 Primal-dual scaling and μ-complementary slackness

The conditions for a point to be optimal for (P(μ)) are closely related to the standard LP optimality conditions. In fact, we have the following theorem:

Theorem 1 Let μ > 0. A point x > 0 is optimal for (P(μ)) if and only if there exist vectors y and s satisfying

This theorem makes it clear that there is no particular reason to favor (P) over (D). Typically, algorithms use iterates (x^k,y^k,s^k) of both primal and dual variables. We will assume the iterates satisfy primal and dual feasibility, and work towards achieving μ-complementary slackness (and complementary slackness in the limit as μ → 0). Further, we assume x^k > 0 and s^k > 0. Because both primal and dual iterates are available, a scaling matrix can be used that combines the primal variables x and the dual slacks s. This matrix is

(2)

We get A^k = AD^k and c^k = D^kc. The directions used to update (x,y,s) are

The dual direction Δy can be calculated on the way to calculating Δx. It can be shown that Δx = -(D^k)²(s - μX^-1e) - (D^k)²Δs.

Note that if x is primal feasible and (y,s) is dual feasible then

Therefore, for a given (y,s), it is equivalent to minimize c^T x or s^T x.

4 The algorithm

1. Initialize with x⁰ > 0 feasible in (P) and (y⁰,s⁰) feasible in (D) with s⁰ > 0. Let k = 0.
2. Let μ^k = ((x^k)^T s^k)∕n². (Other scalings are possible.)
3. Calculate D^k andA^k = AD^k.
4. Calculate Δy, Δs, Δx.
5. Calculate primal and dual steplengths α_P and α_D, ensuring x^k + α _P Δx > 0 and s^k + α _DΔs > 0.
6. Update x^k+1 = x^k + α _P Δx, y^k+1 = y^k + α _DΔy, s^k+1 = s^k + α _DΔs.
7. Let k → k + 1. Calculate (x^k)^T s^k. If small enough, STOP. Else return to Step 2.

5 An example

The problem

has dual

An initial feasible point for (LP1) is x⁰ = (0.6, 1.5, 0.4, 0.5) and an initial dual feasible solution is y⁰ = (-5,-2) giving s⁰ = (4, 1, 6, 3). Note that we have x _i⁰s _i⁰ = 2.4 for i = 1, 3 and x _i⁰s _i⁰ = 1.5 for i = 2, 4. Thus, this point is not centered. The current duality gap is 7.9, so we pick μ = 0.5 (approximately 7.9∕(4²)). We get

Further,

so

and

Now,

. Taking α_P = 2 and α_D = 1 gives

The complementarity values are now

and the new duality gap is 1.45. The individual complementarities are reasonably close to the target value of μ = 0.5, and the duality gap has been noticeably reduced.

6 Computational effort

Each iteration requires a lot of work. The major effort at each iteration is in calculating the projections. It is not necessary calculate (A(D^k)²A^T )^-1 explicitly; Gaussian elimination can be used instead.

It seems that an interior point method will solve almost any problem in at most 40 iterations. The number of iterations required grows very slowly as the size of the problem increases. By contrast, the simplex algorithm seems to need approximately 1.5m iterations to solve a problem with m constraints. The work required at each iteration of an interior point method is far larger than the amount of work needed to compute a simplex pivot, so there is a trade-off: an interior point algorithm needs far fewer iterations, but it takes considerably more time per iteration, when compared to the simplex algorithm.

Commercial implementations of interior point methods use algorithms that are similar to the one presented, but they have some more modifications to make them faster.

7 Further reading

Textbooks on interior point methods include [1, 2, 3, 4, 5].

Current computational results indicate that interior point methods outperform the simplex algorithm on large problems.

References