Part I · trunk A · linear programming · branch

Polynomial-time LP: ellipsoid and interior point

Two ways through the middle of the polyhedron: one that settled a theoretical question and is never used, and one that runs in every modern solver.

needs 02feeds 09~2 h slides + 2½ h lab

This is a branch unit, so the express route skips it. It is worth doing for two reasons. First, unit 09's separation–optimization equivalence runs on the ellipsoid method, and that theorem makes much of Part II's machinery (exponentially many cuts, lazy constraints) theoretically respectable. Second, choosing between simplex and interior point is a real decision you make whenever you configure a solver, and this unit gives you the reasons.

Convention: $\max\{c^\top x : Ax \le b\}$ with $A \in \mathbb{R}^{m\times n}$ and no separate sign constraints. The barrier needs every inequality, sign constraints included, written into $Ax \le b$.

How to read this deck

Both methods are shown on small examples with numbers: one ellipsoid run on a corner of unit 01's polygon, the barrier on a single variable where everything has a closed form, and the central path on the polygon itself. The inequality chains and the bit-size bookkeeping behind the theorems are in the notes under Full proof; nothing in the lab needs them.

04 · overview

The route

1 · is LP in P?

Why the question stayed open after simplex, and what "polynomial" means for an LP.

2 · ellipsoid

Cut an ellipsoid in half, wrap the half in a smaller one, repeat. Polynomial, and impractical.

3 · interior point

A barrier keeps you inside, and following the central path reaches the optimum in a few dozen Newton steps.

getting around

Words with a dotted underline are glossary terms: click one for its definition. G lists every term in the unit, N shows the notes, where the proofs live. Unit 01's polygon is the running example again: a solid underline like that one shows it on hover. Optional extra: interactive figures, linked from the slides they support by a small ▸ chip; skip them freely.

Then: the pure barrier method breaks down in floating point past $t \approx 10^6$. HiGHS's IPM beats its simplex 7× cold at $n = 3000$, and loses 4× on the warm re-solve.

section 1

Is LP in P?

Open from 1947 until 1979, even though LPs were being solved every day.

04 · history

What was known, and what wasn't

1947 Dantzig's simplex method. Fast in practice.
1972 Klee–Minty: exponential on a built family (unit 02).
1970s LP is in NP ∩ coNP: short certificates both ways (units 01, 03). No NP-hardness was expected, but no polynomial algorithm was known.
1979 Khachiyan: the ellipsoid method solves LP in polynomial time.
1984 Karmarkar: a polynomial interior-point method that also looked practical.
1990s– primal–dual IPMs (Mehrotra 1992) in every commercial solver.

"polynomial" for LP

Polynomial in the bit size $L$ of the input: the number of bits needed to write $A, b, c$ as integers.

Example: the constraint $x + y \le 4$ needs about 2 bits per coefficient; $1000x + y \le 4096$ needs up to 13. Same shape of problem, more bits, and a polynomial method may take more steps.

Whether there is a strongly polynomial LP algorithm, polynomial in $m$ and $n$ alone, independent of the numbers, is still open. It is one of Smale's problems for the 21st century.

section 2

The ellipsoid method

Binary search in $n$ dimensions, with ellipsoids as the intervals.

04 · ellipsoid

The idea

In 2-D the new ellipse has semi-axes $\tfrac23$ and $\tfrac{2}{\sqrt3}$ of the old radius, and $\tfrac23 \cdot \tfrac{2}{\sqrt3} \approx 0.77$ of the area.

Start with a ball $E_0$ of radius $R$ that contains the feasible set, if nonempty.
If the centre is feasible, stop.
Otherwise a violated constraint $a^\top z \le \beta$ separates the centre from every feasible point. Keep the half $\{a^\top z \le a^\top \text{centre}\}$.
Replace $E_k$ by the smallest ellipsoid containing that half, and repeat.

Each step shrinks the volume by a fixed factor, so the ellipsoid soon becomes too small to contain the feasible set. The set was therefore empty, or the centre found it.

Step 3 is a central cut: it passes through the centre, so exactly half is kept. The next slide does two steps with numbers.

explore interactively

04 · ellipsoid · example

Two ellipsoid steps, with numbers

Find a point of unit 01's polygon with $x + y \ge 3$, written $-x - y \le -3$. Start: ball of radius $R = 5$ at the origin, $P = 25I$.

$k$	centre	check the constraints	$P$
0	$(0, 0)$	$x + y = 0 < 3$ ✗	$\begin{pmatrix}25 & 0\\0 & 25\end{pmatrix}$
1	$(1.18, 1.18)$	$x + y = 2.36 < 3$ ✗	$\begin{pmatrix}22.2 & -11.1\\-11.1 & 22.2\end{pmatrix}$
2	$(1.96, 1.96)$	every constraint ✓	$\begin{pmatrix}24.7 & -19.8\\-19.8 & 24.7\end{pmatrix}$

step 1, by the formulas

Cut normal $a = (-1, -1)$. $\ Pa = (-25, -25)$, $\ a^\top P a = 50$, so $g = Pa / \sqrt{50} = (-3.54, -3.54)$.

New centre $x' = 0 - g/3 = (1.18, 1.18)$: a third of the radius toward the feasible side. Area kept: $\sqrt{\det P' / \det P} = 0.770$.

Dark strip: the target, the polygon above $x + y = 3$. Dotted: the first cut, $x + y = 0$ through the start centre. Each ellipse holds the half the cut kept.

explore interactively

Step 1 in full. $P' = \frac{4}{3}\left(25I - \frac{2}{3} g g^\top\right)$. Here $g g^\top = 12.5 \begin{pmatrix}1 & 1\\1 & 1\end{pmatrix}$, and $\frac23 \cdot 12.5 = 8.33$, so $P' = \frac43 \begin{pmatrix}16.67 & -8.33\\-8.33 & 16.67\end{pmatrix} = \begin{pmatrix}22.2 & -11.1\\-11.1 & 22.2\end{pmatrix}$. Its eigenvalues are 33.3 along $(1, -1)$ and 11.1 along $(1, 1)$, so the semi-axes are 5.77 (along the cut) and 3.33 (across it): the unit-ball picture of "The idea", scaled by $R = 5$ ($5 \cdot \tfrac{2}{\sqrt3} = 5.77$, $5 \cdot \tfrac23 = 3.33$).

Step 2. The centre still violates the same constraint, so $a$ is again $(-1, -1)$. Now $Pa = (-11.1, -11.1)$, $a^\top P a = 22.2$, $g = (-2.36, -2.36)$, and the centre moves by $-g/3$ to $(1.96, 1.96)$. There $x + y = 3.93$: above 3 and below 4, and $x - y = 0 \le 1$, so every constraint holds and the method stops. The area ratio is again 0.770: in 2-D it is the same at every step, whatever the cut.

The run is the reference ellipsoid_feasible with $R = 5$; it returns after 2 iterations. The figure maps $(x, y)$ to $(190 + 30x,\ 190 - 30y)$, and the ellipses' axes come from the eigenvectors above. Note how thin $E_2$ is across the cut direction: the ellipsoid method learns the direction of the target, not just its location.

04 · ellipsoid

The update, and why the volume shrinks

The ellipsoid update, for the ellipsoid $E = \{z : (z-x)^\top P^{-1}(z-x) \le 1\}$ and cut normal $a$:

$$g = \frac{P a}{\sqrt{a^\top P a}}, \qquad x' = x - \frac{g}{n+1}$$ $$P' = \frac{n^2}{n^2-1}\left(P - \frac{2}{n+1}\, g g^\top\right)$$

An affine map takes $E$ to the unit ball and $a$ to $e_1$, so it suffices to check the unit ball, where the new ellipsoid has the axes on the slide before last.

how slowly, in numbers

$n$	kept per step	steps to halve
2	0.770	2.6
10	0.951	13.8
100	0.995	139

volume ratio $$\begin{aligned} \frac{\mathrm{vol}\,E'}{\mathrm{vol}\,E} &= \frac{n}{n+1}\left(\frac{n}{\sqrt{n^2-1}}\right)^{n-1} \\ &< e^{-1/(2(n+1))} \end{aligned}$$

Why: one axis shrinks by $\tfrac{n}{n+1}$, the other $n - 1$ grow by a little less than that shrinks. The net is a loss, but a small one when $n$ is large.

full proof in the notes (N)

So, with lab step 5: if the feasible set holds a ball of radius $r$ inside the start ball of radius $R$, it can stop after $k \ge 2(n+1)\, n \ln(R/r)$ steps: the iteration bound.

The example table. "Kept per step" is the exact ratio $\frac{n}{n+1}(n/\sqrt{n^2-1})^{n-1}$; "steps to halve" is $\ln 0.5$ divided by its log. At $n = 100$ it takes 139 steps just to halve the volume, and the start ball must shrink to a ball of radius $r$, a factor $(r/R)^{100}$ in volume. That, and not the per-step cost, is what makes the method slow.

Full proof of the volume bound

The inequality: using $1 + u \le e^u$, $\frac{n}{n+1} = 1 - \frac{1}{n+1} \le e^{-1/(n+1)}$ and $\left(\frac{n^2}{n^2-1}\right)^{(n-1)/2} = \left(1 + \frac{1}{n^2-1}\right)^{(n-1)/2} \le e^{(n-1)/(2(n^2-1))} = e^{1/(2(n+1))}$. Multiplying gives $e^{-1/(n+1) + 1/(2(n+1))} = e^{-1/(2(n+1))}$. The first inequality is strict, so the whole thing is.

Lab step 4's test checks both properties numerically. It samples 400 points of the kept half and checks they lie in the new ellipsoid, and it compares the determinant ratio with the bound.

Step 5's bound: $\mathrm{vol}(E_k) < e^{-k/(2(n+1))}\mathrm{vol}(B_R)$, and $\mathrm{vol}(B_r) = (r/R)^n \mathrm{vol}(B_R)$. The ellipsoid is smaller than the $r$-ball once $e^{-k/(2(n+1))} < (r/R)^n$, that is, once $k > 2n(n+1)\ln(R/r)$.

04 · ellipsoid

From feasibility to polynomial-time LP

Three gaps between the lab's version and Khachiyan's theorem, each closed by a fact about bit size ($R$ and $r$ from the data):

How big to start. If any solution exists, one has coordinates at most $2^{O(L)}$, so that is $R$.
How small to stop. A feasible set can be flat, holding no ball at all. Loosening every constraint by a tiny $2^{-O(L)}$ gives it one, without changing whether it is empty.
Optimizing, not just finding. Ask for $c^\top x \ge \gamma$ and search over $\gamma$: the sliding objective.

details in the notes (N)

Khachiyan, 1979

LP feasibility is decided in $O(n^2 L)$ ellipsoid iterations, each with polynomially many bits of precision. So LP ∈ P.

Lab then, part 5 (shifted random polytopes, $R = 100$):

$n$	mean iterations	bound
2	7.8	61
5	29	363
8	74	882

Full details

$R$ from the data. If a solution exists, a vertex solution exists (after adding bounds if needed), and Cramer's rule on its basis gives $\|x\|_\infty \le 2^{O(L)}$. $r$ from the data. Relax each $b_i$ by $2^{-O(L)}$: the system stays feasible exactly when the original is, and now contains a ball of radius $2^{-O(L)}$ (see the relaxation argument below). Optimization. Add $c^\top x \ge \gamma$ and binary-search over $\gamma$, whose range and required precision are again $2^{\pm O(L)}$, or keep a single run and add objective cuts through each feasible centre.

The exact constant in the iteration count depends on how $L$ is defined and on the variant; $O(n^2 L)$ is the order usually quoted. Each iteration must also be carried out with rounding, keeping $O(nL)$ bits per number, and the analysis has to show that the rounding errors don't break the volume argument. That is the technical heart of the proof, in Grötschel, Lovász & Schrijver ch. 3.

The relaxation trick for flat polytopes: if $Ax \le b$ has an integral description with bit size $L$, then $Ax \le b + \varepsilon\mathbf 1$ with $\varepsilon = 2^{-O(L)}$ is feasible if and only if $Ax \le b$ is. Proof sketch: if the relaxed system is feasible, its Farkas certificate for infeasibility of the original, which has bounded size, would still apply for tiny $\varepsilon$.

04 · ellipsoid

Useless in practice, indispensable in theory

why nobody solves LPs with it

The iteration count is close to its worst case on every instance. It never gets lucky the way simplex does.
It needs high numerical precision, and double precision breaks the volume argument.
Each iteration uses a dense $n \times n$ matrix, so sparsity doesn't help.

what it really needs

Look at step 3 again, or at the example: the method only ever asked "which constraint fails at $(1.18, 1.18)$?". It never reads the constraint matrix. It only asks, of one point: "is this feasible, and if not, give me a violated inequality."

That is a separation oracle. It works even when there are exponentially many constraints, provided violated ones can be found efficiently. This is unit 09.

section 3

Interior-point methods

Stay away from the boundary until the very end.

04 · interior point

The log barrier and the central path

Central path of $\max 2x + y$ on unit 01's polygon, computed by the lab's reference solution ($t = 0.01 \cdot 3^k$).

Replace the constraints by a penalty that is infinite on the boundary, the log barrier:

$$\min_x\; f_t(x) = -t\, c^\top x - \sum_{i=1}^m \log(b_i - a_i^\top x)$$

$t \to 0$: the objective hardly matters, and the minimizer is the analytic centre.
$t \to \infty$: the objective dominates, and the minimizer tends to an optimal point.

The minimizers $x^*(t)$ form a smooth curve, the central path. Finding one of them is centering, done by Newton's method.

The next slide shows all of this on one variable, where it can be solved by hand.

explore interactively

04 · interior point · example

The barrier on one variable

$\max x$ subject to $0 \le x \le 1$: two constraints, slacks $x$ and $1 - x$, optimum 1.

$$f_t(x) = -t\,x - \log x - \log(1 - x)$$

Set the derivative to zero, $-t - \tfrac1x + \tfrac{1}{1-x} = 0$, and solve the quadratic:

$$x^*(t) = \frac{t - 2 + \sqrt{t^2 + 4}}{2t}$$

The central path of this LP: from the analytic centre $\tfrac12$ toward the optimum 1, never touching the wall.

$t$	$x^*(t)$	$y = \frac{1}{ts}$	gap $= y_2 - x$	$1 - x^*$
→0	0.5	huge	$2/t$	0.5
1	0.618	(1.618, 2.618)	2	0.382
10	0.910	(0.110, 1.110)	0.2	0.090
100	0.990	(0.010, 1.010)	0.02	0.0099

what to read off

The centre. At $t \to 0$ the barrier alone is minimized at $\tfrac12$, as far from both walls as possible.

A free dual. Always $-y_1 + y_2 = 1 = c$: dual feasible. The gap $b^\top y - c^\top x = y_2 - x$ is exactly $2/t$, with $m = 2$.

The true error $1 - x^*$ is below the gap, so the gap is a safe stopping test.

explore interactively

Solving for $x^*(t)$. Multiply $-t - \frac1x + \frac{1}{1-x} = 0$ by $x(1-x)$: $-t x (1 - x) - (1 - x) + x = 0$, that is $t x^2 + (2 - t) x - 1 = 0$. The root in $(0, 1)$ is the formula on the slide. At $t = 1$ it is $(\sqrt5 - 1)/2 = 0.618$, the golden ratio's reciprocal. As $t \to 0$, $\sqrt{t^2 + 4} \approx 2 + t^2/4$, so the numerator is about $t$ and $x^* \to \tfrac12$.

The dual point. The constraints are $-x \le 0$ and $x \le 1$, with slacks $s_1 = x$, $s_2 = 1 - x$. So $y_1 = 1/(tx)$, $y_2 = 1/(t(1-x))$, and $A^\top y = -y_1 + y_2$. At $t = 1$: $1/0.618 = 1.618$ and $1/0.382 = 2.618$, and indeed $2.618 - 1.618 = 1$. The dual objective is $b^\top y = 0 \cdot y_1 + 1 \cdot y_2$.

Why the gap is exactly $2/t$: that is the next slide's computation, which here reads $y_1 s_1 + y_2 s_2 = \frac1t + \frac1t$. For large $t$ the true error $1 - x^*$ is about $1/t$, half the gap: the dual point is as far above the optimum as the primal point is below it.

The table's values come from the closed form and match the reference newton_center started at $x = 0.5$ to four decimals.

04 · interior point

Duality comes for free on the path

At the minimizer the gradient vanishes, and that equation says $A^\top y = c$ for the dual point $\;y = \dfrac{1}{t\, s} > 0$. So $y$ is dual feasible, and the duality gap is exactly

$$b^\top y - c^\top x = \sum_i y_i s_i = \sum_i \frac{1}{t} = \frac{m}{t}.$$

derivation in the notes (N)

Path following: centre at $t$, then set $t \leftarrow \mu t$ and re-centre from there. With $\mu = 10$ the number of centerings is $\lceil \log_{10}(m/(t_0\varepsilon)) \rceil$.

Theory: with $\mu = 1 + 1/\sqrt{m}$, each centering needs $O(1)$ Newton steps, for $O(\sqrt m \log(m/\varepsilon))$ in total.

consequence

$x^*(t)$ is within $m/t$ of optimal, with a certificate. To reach accuracy $\varepsilon$, run until $t \ge m/\varepsilon$.

on the polygon, one point of the path

At $t = 7.29$: $x = (2.316, 1.597)$, slacks $s = (2.316, 1.597, 0.087, 0.282)$

$y = 1/(ts) = (0.059, 0.086, 1.573, 0.487)$

$A^\top y$: $\ -0.059 + 1.573 + 0.487 = 2$, $\ -0.086 + 1.573 - 0.487 = 1$ ✓

Gap $4 \cdot 1.573 + 0.487 - 6.228 = 0.549 = 4/7.29$ ✓. True error: $6.5 - 6.228 = 0.272$.

explore interactively

Derivation

With $s = b - Ax$, the gradient of $f_t$ is $-t\,c + A^\top (1/s)$, where $1/s$ is taken entrywise. Setting it to zero and dividing by $t$ gives $A^\top \frac{1}{ts} = c$. The gap is $b^\top y - c^\top x = b^\top y - (A^\top y)^\top x = y^\top (b - Ax) = y^\top s$, and each $y_i s_i = \frac{1}{t s_i} s_i = \frac1t$.

The example box. The constraints of the polygon are $-x \le 0$, $-y \le 0$, $x + y \le 4$, $x - y \le 1$, so the columns of $A$ are $(-1, 0, 1, 1)$ for $x$ and $(0, -1, 1, -1)$ for $y$, and $b^\top y = 4y_3 + y_4$. Watch $y$ along the path: at $t = 5314$ it is $(0.0001, 0.0001, 1.5001, 0.5)$, heading for $(0, 0, \tfrac32, \tfrac12)$, unit 02's shadow prices of constraints 3 and 4. The path computes the dual optimum as it goes. All numbers are from the reference central_path with $t_0 = 0.01$, $\mu = 3$, started at $(0.5, 0.5)$.

Compare with complementary slackness (unit 03), which asks for $y_i s_i = 0$. On the central path $y_i s_i = 1/t$ for every $i$: approximate complementary slackness, with the same value for every constraint. Primal–dual interior-point methods are built exactly on this. They solve the perturbed KKT system $Ax + s = b$, $A^\top y = c$, $y_i s_i = \mu$ by Newton's method in all of $(x, y, s)$ together, and drive $\mu$ to zero.

The $O(\sqrt m \log(1/\varepsilon))$ bound comes from self-concordance theory (Nesterov & Nemirovski, 1994), which also explains why the log barrier is the right choice. Practical codes use long steps and a predictor–corrector scheme (Mehrotra 1992), and take 20–60 iterations almost regardless of problem size. The lab's then saw 10 to 25 as $n$ went from 100 to 20 000.

04 · interior point

Where floating point ends the pure barrier method

$t$	$m/t$	$b^\top y - c^\top x$	Newton steps
1e4	1.40e-3	1.40e-3	6
1e5	1.40e-4	1.40e-4	6
1e6	1.40e-5	1.42e-5	6
1e7	1.40e-6	1.31e-3	100

Lab then, part 2: 3 variables, 14 constraints.

On the polygon at $t = 5314$ the slacks are already $(2.5, 1.5, 0.0001, 0.0004)$: entries of $1/s^2$ range from $0.16$ to about $6 \times 10^7$. Each factor of 10 in $t$ adds a factor of 100 to that spread.

why

The Hessian $A^\top \mathrm{diag}(1/s^2) A$ has entries of size $1/s_i^2$. Near the optimum some slacks are $\approx 10^{-7}$ and others $\approx 1$, so its condition number reaches $10^{14}$ and beyond, and the Newton direction becomes numerical noise.

This is floating-point breakdown. Production IPMs (primal–dual) solve a better-behaved system with careful factorizations, and hand over to simplex (crossover) for an exact vertex.

At $t = 10^7$ the Newton iteration hit its 100-step cap without meeting the stopping test, and the resulting point is worse than at $t = 10^6$. At $t = 10^8$ it happened to recover. Erratic behaviour past a threshold is the signature of floating-point breakdown, not of a bug. The lab's tests stop checking centrality at $t = 10^5$ for this reason, and the default stopping gap is $10^{-5}$.

Crossover: an IPM returns a point in the relative interior of the optimal face. Simplex-based branch-and-bound and sensitivity analysis need a basis (a vertex). Crossover (Megiddo 1991) pushes the interior point to a vertex of the same objective value with a sequence of simplex-like steps. It is usually quick, and occasionally the slowest part of the whole solve.

section 4

Choosing between them

Cold solves favour interior point, and warm solves favour simplex.

04 · practice

Cold solves: IPM iterations barely grow

HiGHS, random sparse LPs with 4 nonzeros per column and $m = n/2$, presolve off:

$n$	simplex s	iterations	IPM s	iterations
300	0.006	283	0.006	13
1 000	0.069	1 359	0.025	17
3 000	0.95	5 268	0.14	18
20 000	199	66 941	12	25

A cold solve: simplex iterations grew faster than $n$ here. Each IPM iteration is more expensive (a sparse Cholesky factorization), but there are few of them.

Crossover on this family is at $n \approx 300$.

The $n = 20\,000$ row comes from then --big. Structured real LPs vary widely: some favour simplex even at large sizes, which is why solvers offer concurrent solving.

04 · practice

Warm solves: why branch-and-bound wants simplex

Same LP family, $n = 3000$. Solve, add one constraint that cuts off the optimum, re-solve (a warm solve):

	cold	re-solve
HiGHS simplex	0.832 s	0.033 s (22 it)
HiGHS IPM	0.131 s	0.136 s (18 it)

Simplex re-solves 25× faster than it solves cold. IPM re-solves no faster than cold.

why IPMs can't warm-start well

The old optimum sits on the boundary, with some slacks $\approx 0$. That is the worst starting point for a barrier method, which needs to begin well inside and on the central path. The old basis, by contrast, is one dual simplex pivot or so from optimal (unit 03).

04 · lab

Lab 04 — what you build

about 2½ hours · uv run co lab 04 · numpy, floats

The barrier: value, gradient and Hessian, checked by finite differences.
Newton centering with backtracking that also keeps you feasible.
Path following, with the dual point $y = 1/(ts)$ and the gap $m/t$ checked on every point.
The ellipsoid method: the update, and a search for a feasible point.
The iteration bound from the volume argument, which your method must respect.

uv run co test 04
uv run co then 04        # central path plot · float breakdown · IPM vs simplex cold and warm · ellipsoid

04 · checkpoint

Done when

checkpoint

All five lab steps are ✓.
You can choose between simplex and barrier for a given LP and defend it: size, sparsity, and whether a warm start is available.
You can explain why branch-and-bound still wants simplex at every node.

Reading

Bertsimas & Tsitsiklis, ch. 8 (§8.1–8.3) and ch. 9 (§9.1–9.5).
Boyd & Vandenberghe, Convex Optimization, §11.1–11.3 (free PDF).
Grötschel, Lovász & Schrijver, ch. 3, if you want the full ellipsoid proof.

05 formulations, where integer programming starts. 09 separation ⇔ optimization, where the ellipsoid method pays off, after 08.

Polynomial-time LP: ellipsoid and interior point

The route

Is LP in P?

What was known, and what wasn't

The ellipsoid method

The idea

Two ellipsoid steps, with numbers

The update, and why the volume shrinks

From feasibility to polynomial-time LP

Useless in practice, indispensable in theory

Interior-point methods

The log barrier and the central path

The barrier on one variable

Duality comes for free on the path

Where floating point ends the pure barrier method

Choosing between them

Cold solves: IPM iterations barely grow

Warm solves: why branch-and-bound wants simplex

Lab 04 — what you build

Done when

Reading

Next