Ellipsoid Algorithm, in depth: central cuts, the volume argument, certified bounds and why polynomial did not mean fast

The problem it solves

Start with the feasibility version. You are given a convex set K in n dimensions and asked for a point inside it. You do not see K directly; you can only ask the oracle about a point x. The oracle replies either that x is in K or with a vector a such that every point of K satisfies a·y ≤ a·x, a hyperplane through x with K entirely on one side. For a polytope {y : Ay ≤ b}, the oracle checks each row and returns any violated one.

Two assumptions make the problem well posed: K lies inside a ball of radius R, and if K is non-empty it contains a ball of radius r. The ellipsoid method then either finds a point of K or certifies that K is empty after O(n² log(R/r)) oracle calls. For LPs with integer data of total bit length L, both radii can be bounded in terms of L, after a small perturbation to give the feasible region volume, and that is how polynomial time in L arises.

The central-cut update

Represent an ellipsoid by its centre x and a symmetric positive definite matrix P: E = {y : (y − x)^T P⁻¹ (y − x) ≤ 1}. A ball of radius R is P = R²I. At each step ask the oracle about the centre. If it returns a cut a, K lies in the half ellipsoid {y in E : a·y ≤ a·x}. Replace E with the smallest ellipsoid containing that half. With g = Pa / √(a^TPa), the closed form is:

x_next = x - g / (n + 1)
P_next = n**2 / (n**2 - 1) * (P - 2 / (n + 1) * outer(g, g))

The centre moves along −g by 1/(n + 1) of the distance from the centre to the boundary, away from the cut. The matrix shrinks along g, the direction of the cut, and expands slightly in every other direction by the factor n²/(n² − 1). The figure shows an exact step in two dimensions: an ellipse with semi-axes 4 and 2, cut through its centre along x + y = 0, is replaced by a rotated ellipse with semi-axes 3.35 and 1.84 that just covers the kept half. The formula needs n ≥ 2; in one dimension the method reduces to bisection.

Why the volume shrinks

The convergence proof rests on one ratio. For the central-cut update, vol(E_k+1) / vol(E_k) = (n/(n + 1)) · (n²/(n² − 1))^(n−1)/2, which is less than e^{−1/(2(n+1))}. In the two-dimensional example the actual ratio is 0.770 against a bound of e^−1/6 = 0.846. The ratio does not depend on the data or on which cut was chosen.

After k steps the volume has fallen by at least e^{−k/(2(n+1))}. The starting ball has volume proportional to Rⁿ, and if K is non-empty it contains a ball of volume proportional to rⁿ that every ellipsoid must contain, because every cut keeps all of K. Once the ellipsoid is smaller than that, K must be empty. Solving for k gives about 2n(n + 1) ln(R/r) steps. Notice that the shrinkage per step is tiny in high dimension: at n = 100 each step removes only about half a percent of the volume. That is the root of the method's practical slowness.

From feasibility to optimisation

To minimise c·x over K, use the same loop with a second kind of cut. When the centre is infeasible, cut with the violated constraint. When the centre is feasible, record it as the best point if it improves on the incumbent and then cut with the objective itself, keeping {y : c·y ≤ c·x}. Every optimal point survives every cut, so the optimum stays inside every ellipsoid.

That invariant gives a certificate for free. The minimum of c·y over the ellipsoid is c·x − √(c^TPc), and since the optimum lies inside the ellipsoid, this is a valid lower bound on the optimal value whenever the initial ball contained an optimum. Stop when the best feasible value minus the best lower bound is below ε. You then hold a feasible point and a proof that it is within ε of optimal.

A numpy implementation

import numpy as np

def ellipsoid_lp(c, A, b, R, eps=1e-6, max_iter=100_000):
    """Minimise c@x subject to A@x <= b, assuming an optimum lies in the ball |x| <= R.

    Returns (best_x, best_value, lower_bound, iterations). best_x is None when no
    feasible centre was found before the ellipsoid degenerated.
    """
    n = len(c)
    x = np.zeros(n)
    P = (R * R) * np.eye(n)
    best_x, best_val, lower = None, np.inf, -np.inf
    for k in range(max_iter):
        viol = A @ x - b
        i = int(np.argmax(viol))
        if viol[i] > 0:
            a = A[i]                      # feasibility cut: keep a@y <= a@x
        else:
            val = c @ x
            if val < best_val:
                best_x, best_val = x.copy(), val
            a = c                         # objective cut: keep c@y <= c@x
            lower = max(lower, val - np.sqrt(c @ P @ c))
            if best_val - lower <= eps:
                return best_x, best_val, lower, k
        Pa = P @ a
        width = np.sqrt(a @ Pa)
        if width < 1e-12:
            break
        g = Pa / width
        x = x - g / (n + 1)
        P = (n * n / (n * n - 1.0)) * (P - (2.0 / (n + 1)) * np.outer(g, g))
        P = 0.5 * (P + P.T)               # keep it symmetric against round-off
    return best_x, best_val, lower, k

The oracle here is the loop over rows, A @ x - b, and choosing the most violated row is one reasonable policy among many. Replace those two lines with any routine that returns a separating vector and the rest of the function is unchanged, which is the whole point of the method. Each iteration costs O(mn) for the oracle and O(n²) for the update.

Worked example: a two-variable LP

Take the classic two-variable LP: maximise 3x + 5y subject to x ≤ 4, 2y ≤ 12, 3x + 2y ≤ 18, x ≥ 0, y ≥ 0. Its optimum is 36 at (2, 6). Pass c = (−3, −5) to minimise, encode the five constraints as rows of A, and start from a ball of radius 20 around the origin, which contains the whole feasible region.

The early iterations alternate between feasibility cuts, when the centre has wandered outside the polygon, and objective cuts, which push the centre toward higher values of 3x + 5y. The run stopped after 141 iterations at (1.99999996, 5.99999999) with value −35.9999998 and a certified lower bound of −36.0000005: the gap is under 10⁻⁶, so the returned point is provably within one millionth of optimal. The simplex method solves the same problem in two or three pivots, which already hints at the comparison to come.

Measured: iterations grow as n squared

To see how iterations grow, the harness generated random LPs with n variables, 3n random constraints with right-hand sides between 1 and 2, and a box |x_i| ≤ 5. Each was solved to ε = 10⁻⁶ and checked against scipy's linprog; all 25 matched to 10⁻⁵, and every reported lower bound was below the true optimum.

Iterations grow almost exactly as n², with a constant of about 32. The volume argument explains that number: shrinking the ellipsoid's width in every direction by a factor F needs about 2n² ln F steps, and here F, the ratio of the starting width to ε, is around 10⁷, so 2 ln F is about 32. Halving ε adds only 2n² ln 2 steps; doubling n quadruples the count. With O(n²) work per step, total cost grows as n⁴ for fixed accuracy, before counting the oracle.

Separation oracles: the lasting result

The result that kept the ellipsoid method important is due to Grötschel, Lovász and Schrijver, in a 1981 paper and a 1988 book: for well-described convex sets, optimisation and separation are polynomially equivalent. If you can separate over a polytope in polynomial time, you can optimise a linear function over it in polynomial time, however many facets it has. Some consequences:

• Polytopes with exponentially many facets, such as the matching polytope, can be optimised over directly, because their separation problems reduce to polynomial-time cut computations.
• The first polynomial-time algorithm for minimising a general submodular function came from this framework, before combinatorial algorithms were found.
• The Lovász theta function, and with it maximum-weight stable sets in perfect graphs, can be computed in polynomial time.
• In approximation algorithms, an LP relaxation with exponentially many constraints is often solved via its separation oracle before rounding.

In practice, most of these problems are now solved with cutting-plane loops around the simplex method or specialised algorithms; the ellipsoid method supplies the proof that a polynomial algorithm exists.

Failure modes

• Numerical drift. P loses positive definiteness through round-off in long runs, especially once some axes become tiny. Symmetrise every step, as the code does; for exact guarantees Khachiyan's analysis inflates the ellipsoid slightly each step and keeps enough bits of precision. In floating point, detect a negative or zero a @ P @ a and stop.
• A starting ball that misses the optimum. The lower bound is valid only if an optimum lies inside the initial ball. Derive R from bounds on the variables, never guess.
• Flat feasible sets. An equality constraint gives K zero volume, so the method can only certify approximate feasibility. Eliminate equalities or relax them by a tolerance first.
• Slow convergence mistaken for a bug. Tens of thousands of iterations at n = 40 is the expected behaviour, as the table shows.

Trade-offs against simplex and interior point

Use the ellipsoid method to prove that a problem is polynomially solvable, to teach how cutting-plane methods work, or for small problems where only an oracle exists. Use interior point methods or simplex for real LPs, and read LP duality to see why the lower bound here is a special case of a broader certificate.

What to do next

1. Implement the update and check the two-dimensional step by hand: semi-axes 4 and 2 should become 3.35 and 1.84.
2. Run the solver on the worked LP and confirm the 141-iteration answer and its lower bound.
3. Replace the row loop with a different oracle, such as one constraint per subset, to see that the method never needs the full list.
4. Plot iterations against n on your own random LPs and confirm the n² trend.
5. Read about separation and optimisation equivalence, then compare with cutting-plane loops built on simplex for the same problems.