Interior Point Methods, in depth: the central path, primal-dual Newton steps, Mehrotra's predictor-corrector and real LP solvers

Optimality conditions and the central path

Write the LP in standard form: minimize c.x subject to A x = b and x >= 0, with A an m by n matrix. Its dual is: maximize b.y subject to A^T y + s = c and s >= 0, where s holds the dual slacks. By LP duality (see the Lagrangian duality article), a pair (x, y, s) is optimal exactly when:

A x = b,  x >= 0                 primal feasibility
A^T y + s = c,  s >= 0           dual feasibility
x_i * s_i = 0  for every i       complementary slackness

The first two blocks are linear. All the difficulty sits in the third: for every variable, either x_i or its reduced cost s_i must be zero. Simplex handles this combinatorially, by guessing which x_i are zero (the basis) and repairing the guess. An IPM relaxes it smoothly. Replace x_i * s_i = 0 with x_i * s_i = mu for a positive parameter mu. For each mu there is a unique solution with x and s strictly positive; the curve of these solutions as mu falls to zero is the central path, and it ends at an optimal solution.

On the central path the duality gap is c.x - b.y = x.s = n * mu, so the primal objective is within n mu of optimal. That gives an exact stopping rule: drive mu below the accuracy you need divided by n.

The barrier view

The same path appears from the primal side. Replace the constraints x >= 0 with a logarithmic barrier and minimize t * c.x - sum(log x_i) subject to A x = b. The log term goes to infinity at the boundary, so the minimizer stays strictly inside, and its optimality conditions are exactly the perturbed system above with mu = 1/t. As t grows, the barrier's pull weakens and the minimizer slides toward the optimal vertex.

Figure 1 shows this for the two-variable example used below, with its five inequalities as barrier terms. The guarantee is a bound, not an equality: at t = 10 the point (1.99, 5.97) has objective 35.80, within 0.2 of the optimum 36, inside the bound 5/t = 0.5. The early barrier methods followed this path by solving one barrier problem per t; primal-dual methods instead take a single Newton step per mu, updating x, y and s together, which is far more efficient.

Newton steps and Mehrotra's predictor-corrector

At a current point with x, s greater than zero, write the residuals rp = b - A x and rd = c - A^T y - s. Linearizing the perturbed KKT system gives the Newton equations for the step (dx, dy, ds):

A dx                = rp
A^T dy + ds         = rd
S dx  + X ds        = sigma*mu*1 - X S 1      (X = diag(x), S = diag(s))

Eliminating ds and dx leaves one symmetric positive-definite m by m system, the normal equations: (A D A^T) dy = rhs with D = diag(x / s). Every iteration forms this matrix, factorizes it with Cholesky and back-substitutes. The step is then shortened so x and s stay strictly positive: compute the largest step that keeps them non-negative and take 99 percent of it.

The centering parameter sigma between 0 and 1 trades progress against safety. With sigma = 0 the step aims straight at optimality (the affine-scaling direction) but quickly runs into the boundary. With sigma = 1 it aims back at the central path and makes no progress on mu. Mehrotra's predictor-corrector (1992) picks sigma adaptively: first solve with sigma = 0 to see how far a pure step could go, measure the mu it would reach, and set sigma = (mu_aff / mu) ** 3. Then solve again with that target plus a second-order correction term, reusing the same Cholesky factor, so the second solve is cheap.

A tested implementation

This is a complete infeasible-start Mehrotra solver for dense standard-form LPs. It starts from x = s = 1 and y = 0, which need not satisfy A x = b; the residuals are driven to zero along the way.

import numpy as np

def ipm_lp(A, b, c, tol=1e-8, max_iter=100):
    m, n = A.shape
    x, s, y = np.ones(n), np.ones(n), np.zeros(m)
    for k in range(max_iter):
        rp, rd = b - A @ x, c - A.T @ y - s
        mu = x @ s / n
        if (mu < tol and np.linalg.norm(rp) < tol * (1 + np.linalg.norm(b))
                and np.linalg.norm(rd) < tol * (1 + np.linalg.norm(c))):
            return x, y, s, k
        L = np.linalg.cholesky((A * (x / s)) @ A.T)        # A D A^T

        def solve(rxs):
            rhs = rp + A @ ((x * rd - rxs) / s)
            dy = np.linalg.solve(L.T, np.linalg.solve(L, rhs))
            ds = rd - A.T @ dy
            return (rxs - x * ds) / s, dy, ds

        def step(v, dv):
            neg = dv < 0
            return min(1.0, (-v[neg] / dv[neg]).min()) if neg.any() else 1.0

        dx, dy, ds = solve(-x * s)                         # predictor
        ap, ad = step(x, dx), step(s, ds)
        mu_aff = (x + ap * dx) @ (s + ad * ds) / n
        sigma = (mu_aff / mu) ** 3
        dx, dy, ds = solve(sigma * mu - x * s - dx * ds)   # corrector
        ap, ad = 0.99 * step(x, dx), 0.99 * step(s, ds)
        x, y, s = x + ap * dx, y + ad * dy, s + ad * ds
    raise RuntimeError("no convergence")

Tested on 200 random feasible, bounded LPs (m from 5 to 39, up to about 100 variables), it matched the optimum reported by SciPy's HiGHS to a relative 1e-6 on every one, in 6 to 11 iterations with a median of 8. On larger random problems the count grew slowly: 9 iterations at m = 50, 12 at m = 200, 13 at m = 800. That slow growth is the practical signature of IPMs; the theoretical worst case is O(sqrt(n) log(1/epsilon)) iterations for short-step variants.

Worked example: the Wyndor LP

Take the classic Wyndor problem: maximize 3 x1 + 5 x2 subject to x1 <= 4, 2 x2 <= 12, 3 x1 + 2 x2 <= 18, x non-negative. In standard form, add three slack variables and minimize -3 x1 - 5 x2, so A is 3 by 5. The solver's log:

iter  mu        |rp|     |rd|     objective
 0    1.00e+00  1.5e+01  7.4e+00   -8.000000
 1    4.32e+00  1.3e+01  3.7e+00  -13.695093
 2    3.09e+01  1.3e-01  3.7e-02  -35.742104
 3    3.51e-01  1.3e-03  4.1e-04  -35.987755
 4    3.52e-03  1.3e-05  4.1e-06  -35.999843
 5    3.52e-05  1.3e-07  4.1e-08  -35.999998
 6    3.52e-07  1.3e-09  4.1e-10  -36.000000
 7    3.52e-09  1.3e-11  4.1e-12  -36.000000

Two phases are visible. In iterations 0 to 2 the method is mostly fixing infeasibility: the start x = s = 1 violates A x = b by 15, and the primal residual falls by a factor of 100 in iteration 2. Mu actually rises to 30.9 there, because the large step makes x and s bigger; that is not divergence, just the infeasible start being repaired. From iteration 3 on, the point is essentially feasible and mu falls by a factor of 100 per iteration, the fast local convergence Mehrotra's corrector is known for.

The result is x = (2, 6, 2, 0, 0): x1 = 2, x2 = 6, slack 2 on the first constraint and zero on the other two, objective 36 for the maximization. The dual for the minimization form is y = (0, -1.5, -1); negated for the original maximization, the shadow prices are (0, 1.5, 1). One more unit of the second constraint's capacity is worth 1.5, and the first constraint, which has slack, is worth nothing, exactly as complementary slackness requires.

What dominates cost at scale

On a real problem, almost all the time goes into the Cholesky factorization of A D A^T. The matrix has the same sparsity pattern every iteration (only D changes), so solvers compute a fill-reducing ordering once, with approximate minimum degree or nested dissection, and reuse the symbolic factorization. Three things then decide whether a problem is fast or painful:

• Dense columns. One column of A with nonzeros in every row makes A D A^T completely dense. Solvers detect such columns and handle them separately.
• Ill-conditioning. As mu goes to zero, entries of D head to zero or infinity, and the matrix becomes badly conditioned. Practical codes survive with careful pivot handling, and this is why IPMs usually stop around 1e-8 rather than machine precision.
• Crossover. An IPM ends near the centre of the optimal face, not at a vertex. When you need a basic solution (for warm starts, sensitivity reports or integer programming, see the integer programming article), the solver runs a crossover phase that pushes to a vertex using simplex-style pivots.

Using interior point in practice

Do not use the code above in production; it is dense and has none of the safeguards below. Call a solver. In SciPy, linprog(c, A_ub, b_ub, A_eq, b_eq, method='highs-ipm') selects HiGHS's interior-point code and method='highs' lets HiGHS choose. The legacy method='interior-point' still ran on SciPy 1.17.1 in our test but printed a DeprecationWarning, so do not write new code against it. Commercial solvers expose the same choice, usually as a method parameter selecting barrier, primal simplex or dual simplex, plus a switch to enable or disable crossover.

Beyond LP, the same machinery solves convex quadratic programs, second-order cone programs and semidefinite programs, using self-concordant barriers for each cone. That is why modelling tools such as CVXPY, covered in the convex optimization article, hand most problems to interior-point conic solvers. Robust implementations also use a homogeneous self-dual embedding, which turns infeasibility and unboundedness into a detectable property of the solution rather than a solver that never converges.

Failure modes

• Infeasible or unbounded input. A basic infeasible-start method like the one above simply fails to converge. Check the solver's status code, never just the returned x.
• Rank-deficient A. Redundant equality rows make A D A^T singular, and Cholesky fails. Presolve removes them; do not switch presolve off.
• Bad scaling. Coefficients spanning many orders of magnitude worsen conditioning and inflate iteration counts. Rescale units so values sit roughly between 1e-3 and 1e3.
• Expecting a vertex. With crossover disabled the answer may be a mix of several optimal vertices, with fractional values where simplex would give a clean basis.
• Tolerance confusion. A reported optimum is optimal to about 1e-8 relative. Rounding a nearly integral value and treating it as exact can silently produce an infeasible plan.

Trade-offs against simplex

The rule most practitioners use: try barrier first for a large one-off LP, use dual simplex when you re-solve many slightly changed problems, and benchmark both on your own instances. See the simplex article for the other half of the comparison.

What to do next

1. Add a print of k, mu and the residual norms inside the loop, run the Wyndor problem, and reproduce the log and the duals (0, -1.5, -1).
2. Verify complementary slackness on the result: x_i s_i should be below about 1e-7 for every i (the largest product in our run was 1.7e-8).
3. Generate random LPs as described and compare objectives and iteration counts with linprog(method='highs-ipm').
4. Replace sigma with fixed values such as 0.1 and 0.5 and count how many more iterations the solver needs; that is the value of Mehrotra's heuristic.
5. On a real model, run your solver with barrier and with dual simplex, and with crossover on and off; record time and solution type for each.
6. Turn on the solver log and watch the primal and dual residuals and the gap; a stall there is your cue to check scaling and redundant rows.