Karmarkar's Algorithm, in depth: the canonical form, projective scaling, the potential function and a tested implementation

The canonical form and how to reach it

Karmarkar's algorithm solves LPs in a canonical form: minimise c^Tx subject to Ax = 0, e^Tx = 1 and x ≥ 0, where e is the all-ones vector. Two extra assumptions come with it: the centre of the simplex, e/n, is feasible, which means Ae = 0, and the optimal objective value is exactly 0. The feasible region is the intersection of the probability simplex with a subspace, and the target value is known in advance; what is unknown is where it is attained.

Getting a general LP, minimise c^Tx with Ax = b and x ≥ 0, into this shape takes three steps. First, bound the variables: the input size L gives a bound on the sum of the coordinates of a basic solution, so add a slack to make the sum equal that bound, and scale to make it 1. Second, homogenise: since the coordinates now sum to 1, b = b(e^Tx), so Ax = b becomes (A - be^T)x = 0. Third, add an artificial variable whose column makes the centre feasible, with a large cost so it is driven to zero. The remaining assumption, optimal value 0, is handled either by merging the primal and dual into one problem whose objective is the duality gap, which is 0 at optimality, or by Todd and Burrell's sliding-objective variant that updates a lower bound on the optimum from dual estimates as it goes.

The projective transformation

The key difficulty for any interior method is that the current point x_k may sit close to a face of the polytope, where only a tiny step is possible before some coordinate turns negative. Karmarkar's insight was to change coordinates so the current point becomes the centre. Let D = diag(x_k) and define the projective map T(x) = D^-1x / (e^TD^-1x). It maps the simplex to itself, sends x_k to e/n, and sends the subspace Ax = 0 to the subspace ADy = 0. Because the objective is homogeneous and its optimum is 0, minimising c^Tx is equivalent to minimising the linear function (Dc)^Ty in the new coordinates as far as reaching zero is concerned.

At the centre of the simplex the largest inscribed ball has radius r = 1/√(n(n-1)), and every point of that ball intersected with the feasible subspace is feasible. So the step is: project the transformed cost Dc onto the null space of B = [AD; e^T], move from e/n a distance αr against that direction, and map back with x = Dy / e^TDy. Karmarkar used α = 1/4.

The potential function and the complexity bound

The objective alone does not decrease monotonically under this step, because the map back distorts distances. Progress is measured instead by the potential f(x) = n ln(c^Tx) - Σ_j ln x_j. The first term rewards a small objective, the barrier term penalises approaching the boundary, and the function is invariant under the projective map up to a constant, so progress in y-space is progress in x-space. Karmarkar proved that each step reduces f by at least a fixed constant independent of n.

That one inequality gives the complexity bound. Starting at the centre, getting c^Tx below 2^-O(L), from where an exact vertex can be recovered by rounding, needs the potential to fall by O(nL), so O(nL) iterations suffice. Each iteration is dominated by solving a linear system with BB^T, O(n³) naively, giving O(n⁴L). Karmarkar then observed that D changes little per step, so the inverse can be maintained with rank-one updates of only the entries that changed noticeably, bringing the total to O(n^3.5L) arithmetic operations.

A tested implementation

The implementation is short because NumPy does the linear algebra. It forms BB^T explicitly and solves with least squares for clarity; a serious implementation would factor it with a Cholesky decomposition and exploit sparsity.

import numpy as np

def karmarkar(A, c, alpha=0.25, eps=1e-8, max_iter=10_000):
    """Minimise c @ x s.t. A @ x = 0, sum(x) = 1, x >= 0.
    Canonical form: A @ ones = 0 and the optimal value is 0."""
    m, n = A.shape
    x = np.full(n, 1.0 / n)
    r = 1.0 / np.sqrt(n * (n - 1))           # inscribed-ball radius at the centre
    for k in range(max_iter):
        if c @ x < eps:
            return x, k
        D = np.diag(x)
        B = np.vstack([A @ D, np.ones(n)])   # constraints after scaling
        cd = D @ c                           # cost after scaling
        w = np.linalg.lstsq(B @ B.T, B @ cd, rcond=None)[0]
        cp = cd - B.T @ w                    # projection onto null(B)
        norm = np.linalg.norm(cp)
        if norm < 1e-14:                     # cost constant on the face: done
            return x, k
        y = np.full(n, 1.0 / n) - alpha * r * cp / norm
        x = (D @ y) / (D @ y).sum()          # inverse projective map
    return x, max_iter

To test it, random canonical instances were built by choosing a sparse optimal point x*, generating A orthogonal to both e and x*, and giving positive costs only to coordinates where x* is zero, so the optimum is 0 by construction. On 200 such instances with 5 to 60 variables, the result agreed with SciPy's linprog using HiGHS to within 1.0e-8, the stopping tolerance, and every iterate stayed strictly positive and feasible.

Worked example and measured iterations

Take n = 3, the single constraint x₁ + x₂ - 2x₃ = 0 and cost c = (1, 0, 0). With the simplex constraint this forces x₃ = 1/3 and x₁ + x₂ = 2/3, so the optimum is (0, 2/3, 1/3) with value 0. The centre (1/3, 1/3, 1/3) is feasible and the potential starts at 0.

The first step moves to (0.2612, 0.4055, 0.3333), the second to (0.1950, 0.4717, 0.3333) and the third to (0.1401, 0.5266, 0.3333): the method slides along the feasible segment, keeping x₃ fixed as the constraint demands. The objective goes 0.333, 0.261, 0.195, 0.140 and the potential 0, -0.68, -1.42, -2.19, falling by roughly 0.7 to 0.8 per step on this instance, well above the worst-case constant the proof guarantees. With α = 1/4 it takes 44 iterations to reach an objective below 1e-8; with α = 0.9, still inside the ball, it takes 7.

On the random instances the median iteration counts at α = 1/4 were 159, 278 and 393 for n = 20, 60 and 120, and at α = 0.9 they were 45, 78 and 111. Both grow with n, as the O(nL) bound allows, and that slow, short-step behaviour is exactly what practical codes abandoned.

From Karmarkar to modern interior-point solvers

Within two years the field understood the method better. Gill, Murray, Saunders, Tomlin and Wright showed in 1986 that Karmarkar's step is a projected Newton step for a logarithmic barrier function, connecting it to Fiacco and McCormick's barrier methods from the 1960s. Several groups, including Barnes and Vanderbei, Meketon and Freedman, proposed affine scaling, which drops the projective map and the known optimum: scale by D, step against the projected cost, and move a large fraction, often 90 to 99 percent, of the way to the boundary along that direction, using a ratio test on the coordinates that are shrinking. That long step is the real practical change; the α = 0.9 run above is still a short step inside the inscribed ball. Affine scaling has no known polynomial bound, yet it usually needs far fewer iterations. It later emerged that Dikin had published the same idea in 1967.

AT&T patented the method (US 4,744,028, granted in 1988) and sold it in the KORBX system, a commercial episode that fed the long debate about patenting algorithms. Meanwhile the theory converged on primal-dual path-following methods and Mehrotra's predictor-corrector, which take long steps, need no canonical form, and typically finish in tens of iterations almost regardless of size. That is what modern solvers run, alongside the simplex method, which remains strong for warm starts and re-optimisation. Karmarkar's lasting contributions are the idea of re-centring by scaling, potential-function analysis, and the proof that interior methods can be both polynomial and fast.

Failure modes

• Ill-conditioned normal equations. As coordinates approach zero, D squashes columns of AD and BB^T becomes nearly singular; use a stable factorisation and watch its condition number rather than inverting.
• Wrong optimal value. If the optimum is not exactly 0, the potential argument fails and the iterates stall; the canonical form must be built honestly.
• Infeasible centre. If Ae is not 0, the first iterate is infeasible and everything that follows is meaningless; assert it before iterating.
• Stopping without a vertex. The method converges to the optimal face, not a vertex; if you need a basic solution, add a crossover or rounding step.
• Dense linear algebra. Forming BB^T densely costs O(m²n) per step and erases the sparsity real LPs depend on.
• Pushing α to 1. At α = 1 the step reaches the edge of the inscribed ball, and rounding can make a coordinate zero or negative, after which the logarithms in the potential are undefined; keep α strictly below 1 and assert x > 0 after every step.
• Absolute stopping tests. A fixed threshold such as 1e-8 on c^Tx means different things for different cost scales; normalise c, or stop on a relative duality gap as production solvers do.

Trade-offs

What to do next

1. Run the implementation on the n = 3 example and check the first three iterates against the values above.
2. Add the potential function to the loop and assert it falls every iteration.
3. Convert a small LP of your own to canonical form by hand, including the primal-dual merge, and solve it.
4. Replace lstsq with a Cholesky solve and measure how the condition number grows near the optimum.
5. Compare the iteration counts with linprog's interior-point and simplex options on the same instances.
6. Continue with the simplex method, LP duality and convex optimization.