Introduction

Every airline schedule, every factory production plan, every portfolio allocation — at their mathematical core they are all linear programs: maximize (or minimize) a linear objective subject to linear constraints. The simplex method had been solving these in practice since the 1940s, and it was blazingly fast. But nobody could prove it was theoretically efficient. In the worst case, simplex visits exponentially many vertices.

For decades the question sat open: is linear programming (LP) in P — solvable in polynomial time? Then in 1979 a Soviet mathematician, Leonid Khachiyan, answered yes, using an idea borrowed from geometry. Instead of walking along the edges of a polytope the way simplex does, his algorithm wraps an ellipsoid around the feasible region and shrinks it, step by step, until it collapses onto a solution.

The algorithm was slower than simplex in every practical test, but that was beside the point. It settled a fundamental theoretical question: LP is in P. And the ideas it introduced — separation oracles, volumetric arguments, and solving problems implicitly — went on to shape modern optimization.

Watch It Shrink

The demo below runs the ellipsoid algorithm on a small 2D linear program. The blue region is the feasible polytope — the set of points satisfying all constraints. The algorithm starts with a large ellipse guaranteed to contain the optimum and, at each step, cuts the ellipse in half with the most-violated constraint, then replaces it with the smallest ellipse containing the surviving half.

<p class="hint">A 2D linear program: maximize <b>c·x</b> subject to <b>Ax ≤ b</b>.
    The blue region is the feasible polytope. The ellipse shrinks each step by cutting
    at the most-violated constraint and taking the minimum enclosing ellipse of the half
    that could still contain the optimum.</p>
<canvas id="cv" width="340" height="280"></canvas>
<div class="info" id="info">Press <b>Step</b> or <b>Run</b> to start.</div>
<div class="btns">
  <button id="btnStep" type="button">Step</button>
  <button id="btnRun"  type="button">Run</button>
  <button id="btnReset" type="button" class="ghost">Reset</button>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .85rem; color: #444; margin: 0 0 .6rem; line-height: 1.45; }
canvas { display: block; border: 1px solid #cdd9e3; border-radius: 8px;
         background: #f8fafc; max-width: 100%; }
.info  { font-size: .9rem; font-weight: 600; margin: .5rem 0; min-height: 1.4em; color: #1d3557; }
.info.done { color: #0a7d33; }
.btns  { display: flex; gap: .5rem; flex-wrap: wrap; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem;
         border: 1px solid #1d3557; background: #1d3557; color: #fff;
         border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Ellipsoid method in 2D (simplified but correct).
// Problem: maximize c·x subject to A x <= b (x in R^2).
// We find the minimum by minimizing -c·x.
// Constraints define a feasible polygon rendered in blue.
//
// State: ellipsoid E = { x : (x-c)^T P^{-1} (x-c) <= 1 }
// represented by center vector 'ctr' and 2×2 positive-definite matrix 'M'
// (M = P, the shape matrix).
//
// Update rule (Khachiyan):  given a violated constraint a·x > b at center c:
//   g = P a / sqrt(a^T P a)            (normalised gradient in ellipsoid metric)
//   c_new = c - (1/(n+1)) g
//   P_new = (n^2/(n^2-1)) ( P - (2/(n+1)) g g^T )
// Volume reduction factor: e^{-1/(2n)} per step (n=2 => ~0.78 per step).

const W = 340, H = 280;
const canvas = document.getElementById('cv');
const ctx = canvas.getContext('2d');

// World coordinates: x in [-1, 10], y in [-1, 10]
const XMIN = -1, XMAX = 10, YMIN = -1, YMAX = 10;
function toScreen(wx, wy) {
  return [
    (wx - XMIN) / (XMAX - XMIN) * W,
    H - (wy - YMIN) / (YMAX - YMIN) * H
  ];
}
function toWorld(sx, sy) {
  return [
    sx / W * (XMAX - XMIN) + XMIN,
    (H - sy) / H * (YMAX - YMIN) + YMIN
  ];
}

// Feasible region: A x <= b  (each row of A, each entry of b)
// x >= 0, y >= 0, x + y <= 8, 2x + y <= 12, x + 2y <= 12
const A = [
  [-1,  0],   // -x <= 0  => x >= 0
  [ 0, -1],   // -y <= 0  => y >= 0
  [ 1,  1],   // x + y <= 8
  [ 2,  1],   // 2x + y <= 12
  [ 1,  2],   // x + 2y <= 12
];
const b = [0, 0, 8, 12, 12];

// Objective: maximize c·x = maximize x + y  =>  minimise -(x+y)
const objC = [1, 1];          // objective direction (maximise)
const gradMin = [-1, -1];     // gradient for minimisation (-objC)

// True optimum: x=4, y=4 (corner of x+y=8, x+2y=12 => y=4, x=4)
const OPT = [4, 4];

// 2×2 matrix helpers
function mScale(s, M) { return [[s*M[0][0], s*M[0][1]], [s*M[1][0], s*M[1][1]]]; }
function mAdd(A2, B2) { return [[A2[0][0]+B2[0][0], A2[0][1]+B2[0][1]], [A2[1][0]+B2[1][0], A2[1][1]+B2[1][1]]]; }
function mSub(A2, B2) { return [[A2[0][0]-B2[0][0], A2[0][1]-B2[0][1]], [A2[1][0]-B2[1][0], A2[1][1]-B2[1][1]]]; }
function mMulVec(M, v) { return [M[0][0]*v[0]+M[0][1]*v[1], M[1][0]*v[0]+M[1][1]*v[1]]; }
function dot(u, v) { return u[0]*v[0] + u[1]*v[1]; }
function outerProduct(u, v) { return [[u[0]*v[0], u[0]*v[1]], [u[1]*v[0], u[1]*v[1]]]; }
function vScale(s, v) { return [s*v[0], s*v[1]]; }
function vSub(u, v) { return [u[0]-v[0], u[1]-v[1]]; }
function vAdd(u, v) { return [u[0]+v[0], u[1]+v[1]]; }

// Find a violated constraint at point x; returns index or -1 if feasible.
function findViolated(x) {
  let worst = -1, worstSlack = 0;
  for (let i = 0; i < A.length; i++) {
    const slack = dot(A[i], x) - b[i];
    if (slack > worstSlack) { worstSlack = slack; worst = i; }
  }
  return worst;
}

// Ellipsoid update: minimising gradMin direction, or cutting at constraint i.
// If x is feasible, we cut in the objective gradient direction to improve.
// Returns { ctr, M } for new ellipsoid.
function ellipsoidStep(ctr, M) {
  const n = 2;
  let a; // the separating hyperplane normal

  const vi = findViolated(ctr);
  if (vi >= 0) {
    a = A[vi];            // violated constraint gradient
  } else {
    a = gradMin;          // feasible but can improve (cut towards optimum)
  }

  const Pa = mMulVec(M, a);
  const aPa = dot(a, Pa);
  if (aPa <= 0) return { ctr, M }; // degenerate
  const sqrtAPa = Math.sqrt(aPa);

  // g = Pa / sqrt(a^T P a)
  const g = vScale(1 / sqrtAPa, Pa);

  // c_new = c - (1/(n+1)) g
  const ctrNew = vSub(ctr, vScale(1 / (n + 1), g));

  // P_new = (n^2/(n^2-1)) * (P - (2/(n+1)) * g g^T)
  const factor = (n * n) / (n * n - 1);
  const ggT = outerProduct(g, g);
  const Mnew = mScale(factor, mSub(M, mScale(2 / (n + 1), ggT)));

  return { ctr: ctrNew, M: Mnew };
}

// Initial ellipsoid: large ball centred at (5, 5), radius 9 (covers feasible region)
const R0 = 9;
function initialState() {
  return {
    ctr: [5, 5],
    M: [[R0*R0, 0], [0, R0*R0]],
    step: 0,
    done: false,
    history: []
  };
}
let state = initialState();

// Compute vertices of feasible polygon for drawing (clipped to world bounds).
function feasiblePolygon() {
  // Grid-sample to approximate polygon boundary (good enough for display).
  const pts = [];
  const eps = 0.05;
  for (let xi = XMIN; xi <= XMAX; xi += eps) {
    for (let yi = YMIN; yi <= YMAX; yi += eps) {
      const x = [xi, yi];
      let ok = true;
      for (let i = 0; i < A.length; i++) if (dot(A[i], x) > b[i] + 1e-6) { ok = false; break; }
      if (ok) pts.push([xi, yi]);
    }
  }
  return pts;
}
// Compute actual polygon corners (intersections of constraint pairs).
function feasibleCorners() {
  const corners = [];
  for (let i = 0; i < A.length; i++) {
    for (let j = i+1; j < A.length; j++) {
      const det = A[i][0]*A[j][1] - A[i][1]*A[j][0];
      if (Math.abs(det) < 1e-9) continue;
      const x = (b[i]*A[j][1] - b[j]*A[i][1]) / det;
      const y = (A[i][0]*b[j] - A[j][0]*b[i]) / det;
      let feasible = true;
      for (let k = 0; k < A.length; k++) if (dot(A[k], [x,y]) > b[k] + 1e-6) { feasible = false; break; }
      if (feasible && x >= XMIN-1 && x <= XMAX+1 && y >= YMIN-1 && y <= YMAX+1) {
        corners.push([x, y]);
      }
    }
  }
  // Sort corners by angle for polygon drawing
  const cx = corners.reduce((s,p)=>s+p[0],0)/corners.length;
  const cy = corners.reduce((s,p)=>s+p[1],0)/corners.length;
  corners.sort((a,c) => Math.atan2(a[1]-cy,a[0]-cx) - Math.atan2(c[1]-cy,c[0]-cx));
  return corners;
}
const corners = feasibleCorners();

function drawEllipse(ctr, M) {
  // Eigendecompose M to get axes and angle.
  // For 2x2 symmetric [[a,b],[b,d]]:
  const a = M[0][0], bv = M[0][1], d = M[1][1];
  const tr = a + d, det = a*d - bv*bv;
  const disc = Math.sqrt(Math.max(0, tr*tr/4 - det));
  const lam1 = tr/2 + disc;
  const lam2 = tr/2 - disc;
  if (lam1 <= 0 || lam2 <= 0) return;
  const rx = Math.sqrt(lam1), ry = Math.sqrt(Math.max(lam2, 1e-9));
  // Eigenvector for lam1: (M - lam1 I) v = 0
  let angle = 0;
  if (Math.abs(bv) > 1e-10) {
    angle = Math.atan2(lam1 - a, bv);
  } else {
    angle = (a >= d) ? 0 : Math.PI/2;
  }
  // Convert axes from world to screen scale (approximate — aspect ratio)
  const scaleX = W / (XMAX - XMIN);
  const scaleY = H / (YMAX - YMIN);
  // Centre in screen
  const [cx, cy] = toScreen(ctr[0], ctr[1]);
  ctx.save();
  ctx.translate(cx, cy);
  ctx.rotate(-angle); // screen y flipped
  ctx.scale(rx * scaleX, ry * scaleY);
  ctx.beginPath();
  ctx.arc(0, 0, 1, 0, 2*Math.PI);
  ctx.restore();
}

function draw() {
  ctx.clearRect(0, 0, W, H);

  // Draw feasible polygon
  if (corners.length > 2) {
    ctx.beginPath();
    const [sx, sy] = toScreen(corners[0][0], corners[0][1]);
    ctx.moveTo(sx, sy);
    for (let i = 1; i < corners.length; i++) {
      const [px, py] = toScreen(corners[i][0], corners[i][1]);
      ctx.lineTo(px, py);
    }
    ctx.closePath();
    ctx.fillStyle = 'rgba(29, 53, 87, 0.15)';
    ctx.fill();
    ctx.strokeStyle = '#1d3557';
    ctx.lineWidth = 1.5;
    ctx.stroke();
  }

  // Draw grid axes
  ctx.strokeStyle = '#cdd9e3';
  ctx.lineWidth = 0.8;
  const [ox, oy] = toScreen(0, 0);
  ctx.beginPath(); ctx.moveTo(ox, 0); ctx.lineTo(ox, H); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(0, oy); ctx.lineTo(W, oy); ctx.stroke();

  // Draw ellipse (current)
  drawEllipse(state.ctr, state.M);
  ctx.strokeStyle = 'rgba(230, 57, 70, 0.7)';
  ctx.lineWidth = 1.8;
  ctx.stroke();
  ctx.fillStyle = 'rgba(230, 57, 70, 0.08)';
  ctx.fill();

  // Draw history ellipses (faded)
  for (let h of state.history) {
    drawEllipse(h.ctr, h.M);
    ctx.strokeStyle = 'rgba(230, 57, 70, 0.2)';
    ctx.lineWidth = 1;
    ctx.stroke();
  }

  // Draw optimum star
  const [ox2, oy2] = toScreen(OPT[0], OPT[1]);
  ctx.fillStyle = '#0a7d33';
  ctx.beginPath(); ctx.arc(ox2, oy2, 5, 0, 2*Math.PI); ctx.fill();
  ctx.fillStyle = '#fff';
  ctx.beginPath(); ctx.arc(ox2, oy2, 2, 0, 2*Math.PI); ctx.fill();

  // Draw current center
  const [ccx, ccy] = toScreen(state.ctr[0], state.ctr[1]);
  ctx.fillStyle = '#e63946';
  ctx.beginPath(); ctx.arc(ccx, ccy, 4, 0, 2*Math.PI); ctx.fill();
}

function isFeasible(x) {
  for (let i = 0; i < A.length; i++) if (dot(A[i], x) > b[i] + 1e-4) return false;
  return true;
}

function step() {
  if (state.done) return;
  state.history.push({ ctr: state.ctr, M: state.M });
  if (state.history.length > 30) state.history.shift();
  const result = ellipsoidStep(state.ctr, state.M);
  state.ctr = result.ctr;
  state.M = result.M;
  state.step++;

  // Check if volume is negligibly small (converged)
  const vol = Math.sqrt(Math.max(0, state.M[0][0]*state.M[1][1] - state.M[0][1]*state.M[1][0]));

  const feasible = isFeasible(state.ctr);
  const info = document.getElementById('info');
  if (vol < 0.01 || state.step >= 120) {
    state.done = true;
    const dist = Math.sqrt((state.ctr[0]-OPT[0])**2 + (state.ctr[1]-OPT[1])**2);
    info.textContent = fmt('Converged in %step% steps — centre (%v%, %v2%), dist to opt: %dist%', { step: state.step, v: state.ctr[0].toFixed(2), v2: state.ctr[1].toFixed(2), dist: dist.toFixed(3) });
    info.className = 'info done';
  } else {
    info.textContent = fmt('Step %step% — centre (%v%, %v2%)%v3% | vol ∝ %vol%', { step: state.step, v: state.ctr[0].toFixed(2), v2: state.ctr[1].toFixed(2), v3: feasible ? ' [feasible]' : ' [infeasible]', vol: vol.toFixed(3) });
    info.className = 'info';
  }
  draw();
}

let timer = null;
document.getElementById('btnStep').addEventListener('click', () => {
  clearInterval(timer); timer = null;
  step();
});
document.getElementById('btnRun').addEventListener('click', () => {
  if (timer) { clearInterval(timer); timer = null; return; }
  timer = setInterval(() => { step(); if (state.done) { clearInterval(timer); timer = null; } }, 120);
});
document.getElementById('btnReset').addEventListener('click', () => {
  clearInterval(timer); timer = null;
  state = initialState();
  document.getElementById('info').textContent = 'Press Step or Run to start.';
  document.getElementById('info').className = 'info';
  draw();
});

draw();
console.log('Ellipsoid method demo loaded. True optimum:', OPT, 'obj =', dot(objC, OPT));
Result
Console

Notice two things. First, the ellipse always contains at least one optimal point — the algorithm never misses the target. Second, each iteration shrinks the volume by a fixed factor, so after polynomially many steps the ellipse is too small to hold anything except (approximately) the optimum. The red dot marks the current center; when it lands inside the feasible region it becomes the candidate solution.

The Real Complexity

What did Khachiyan actually prove, and why does it matter?

  • LP is in P. The ellipsoid method runs in time polynomial in the number of variables, constraints, and the bit-length of the input — specifically O(n6L2)O(n^{6} L^{2}) bit operations, where n is the number of variables and L is the encoding length. This is the first polynomial-time LP algorithm.
  • Volume shrinks by a constant factor each step. The key lemma: when you cut an ellipsoid in half with a hyperplane and take the smallest enclosing ellipsoid of the surviving half, its volume drops by a factor of at least e−1/2ne^{-1/2n}. After O(n2L)O(n^{2} L) iterations, the ellipse is too small to contain a feasible point of any meaningful size.
  • Separation oracle model. The method only needs to test whether a point is feasible (and if not, find a violated constraint). This abstraction let it solve problems with exponentially many constraints efficiently — as long as violations can be found quickly.
  • Slower than simplex in practice. Despite its theoretical superiority, the ellipsoid method is rarely used directly. The constant factors are large and the algorithm is numerically sensitive. It was superseded by interior-point methods (Karmarkar, 1984), which are both polynomial and fast in practice.
  • Theoretical legacy. The ellipsoid method remains central to complexity theory. It proves linear programming is in P, solves problems in combinatorial optimization (like the ellipsoid-plus-separation-oracle approach for max-weight matching), and shows that many problems previously thought hard are actually polynomial.

Where It Matters

The ellipsoid method's theoretical influence extends far beyond linear programs:

  • Combinatorial optimization: problems with exponentially many constraints (like the max-weight perfect matching polytope or the stable-set polytope of perfect graphs) can be solved in polynomial time whenever a separation oracle exists — even though you could never list all the constraints explicitly.
  • Semidefinite programming (SDP): SDPs generalize LP by replacing linear inequalities with matrix positive-semidefiniteness. The ellipsoid method (and its successor, interior-point SDP solvers) made SDP tractable and powered decades of advances in approximation algorithms.
  • Game theory: the ellipsoid method can compute Nash equilibria in zero-sum games and other structured games in polynomial time, connecting optimization to Nash equilibrium theory.
  • Separation = optimization: the deep insight is that for convex problems, testing feasibility and optimizing are equivalent in polynomial time. If you can separate, you can optimize. This principle drives modern convex optimization.
  • Interior-point methods: Karmarkar's 1984 algorithm shares the ellipsoid method's polynomial guarantee but is practical. Today's commercial LP solvers use interior-point methods for large instances — the direct descendant of Khachiyan's theoretical breakthrough.

Conclusion

The ellipsoid method answered a question that had haunted optimization for decades: linear programming is in P. It did so by replacing the combinatorial edge-walking of simplex with a purely geometric idea — shrink an ellipsoid by a fixed volumetric factor at every step, and after polynomially many steps nothing remains except the answer.

In practice, engineers kept using simplex (and later Karmarkar's interior-point method). But in theory, the ellipsoid method was a watershed: it introduced the separation oracle model, showed that optimization and feasibility testing are equivalent for convex problems, and planted the seeds for semidefinite programming and a generation of approximation algorithms.

Sometimes the best proof is not the fastest algorithm, but the one that finally shows something is possible at all. For LP, that proof was an ellipsoid, shrinking.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/ellipsoid-method/Content licensed under CC BY-NC 4.0.