Introduction to Megiddo's LP Algorithm

Introduction

Linear programming (LP) is one of the most useful tools in all of optimization: given a set of linear inequalities and a linear objective, find the point that satisfies all the constraints and maximizes (or minimizes) the goal. It models everything from factory schedules to diet plans.

The general LP problem with n constraints and d variables was long known to be solvable in polynomial time — the simplex method is fast in practice, and the ellipsoid method proved polynomial complexity in 1979. But could small-dimensional LP be dramatically faster?

In 1984, Nimrod Megiddo answered yes. His prune-and-search algorithm solves LP in the plane — n linear constraints, objective along a line — in $O(n)$ time. No other approach had cracked the linear barrier for LP. The result was one of the first to show that the hard cases of polynomial-time problems can sometimes collapse to linear time when the dimension is held fixed.

The idea is elegant: instead of testing the optimum directly, pair up constraints, use each pair to prune half of the remaining candidates in constant time, and recurse on the surviving half. Each level of recursion halves the problem, so only $O(\log n)$ levels are needed — each costing $O(n)$ — giving $O(n)$ total after a careful implementation.

Try It: Prune-and-Search

Below is a 2D linear programming instance: maximize cx (i.e. move as far right as possible) subject to n half-plane constraints of the form ax + by ≤ c. The feasible region is the intersection of all half-planes; the optimum is the rightmost feasible point.

<div id="controls">
  <button id="btnStep" type="button">Step</button>
  <button id="btnSolve" type="button">Solve instantly</button>
  <button id="btnReset" type="button" class="ghost">Reset</button>
  <span id="info">n = <b id="nCount">0</b> constraints active</span>
</div>
<canvas id="canvas" width="560" height="310"></canvas>
<div id="log"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; background: #f7f9fb; }
#controls { display: flex; gap: .5rem; align-items: center; flex-wrap: wrap; padding: .5rem .5rem .3rem; }
button { font: 600 13px system-ui; padding: .38rem .8rem; border-radius: 7px; cursor: pointer;
         border: 1px solid #1d3557; background: #1d3557; color: #fff; }
button.ghost { background: #fff; color: #1d3557; }
#info { font-size: .85rem; color: #555; margin-left: .3rem; }
canvas { display: block; width: 100%; max-width: 560px; height: auto; border: 1px solid #cdd5de;
         border-radius: 8px; background: #fff; margin: 0 auto; }
#log { font-size: .8rem; color: #555; min-height: 1.6em; padding: .3rem .5rem; text-align: center; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Megiddo's prune-and-search (pedagogical): maximize x over half-planes.
// Each constraint is   y <= (c - a*x)/b   (upper)  or   y >= (c - a*x)/b  (lower).
// The feasible region is the x-range where lowerBound(x) <= upperBound(x).
// We find the rightmost feasible x using prune-and-search:
//   pair each upper with a lower, compute their intersection x (pivot), find median pivot, test feasibility on one side, discard non-binding constraints.

const W = 560, H = 310, PAD = 36;
const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');

// -- coordinate transform --
const xLo=-3.5, xHi=3.8, yLo=-4, yHi=4;
function cx(x){ return PAD + (x-xLo)/(xHi-xLo)*(W-2*PAD); }
function cy(y){ return H - PAD - (y-yLo)/(yHi-yLo)*(H-2*PAD); }

// -- constraint evaluators --
function evalConstraint(c, x){ return (c.c - c.a*x) / c.b; }
function upperAt(active, cs, x){
  let u=Infinity; for(const i of active) if(cs[i].type==='upper') u=Math.min(u,evalConstraint(cs[i],x)); return u;
}
function lowerAt(active, cs, x){
  let l=-Infinity; for(const i of active) if(cs[i].type==='lower') l=Math.max(l,evalConstraint(cs[i],x)); return l;
}
function feasible(active, cs, x){ return lowerAt(active,cs,x) <= upperAt(active,cs,x)+1e-9; }

// -- state --
let CS=[], active=[], xLeft, xRight, optX=null, optY=null, solved=false, stepN=0;

// -- pivot: x where upper constraint i equals lower constraint j --
function pivot(cs, i, j){
  const u=cs[i], v=cs[j];
  // evalU(x) = evalV(x)  =>  (u.c-u.a*x)/u.b = (v.c-v.a*x)/v.b
  // v.b*(u.c-u.a*x) = u.b*(v.c-v.a*x)
  // x*(u.b*v.a - v.b*u.a) = u.b*v.c - v.b*u.c
  const denom = u.b*v.a - v.b*u.a;
  if(Math.abs(denom)<1e-9) return null;
  return (u.b*v.c - v.b*u.c)/denom;
}

function median(arr){
  const s=arr.slice().sort((a,b)=>a-b); return s[Math.floor(s.length/2)];
}

// -- generate instance guaranteed feasible in [-3, optXTarget] --
function generate(){
  // pick a random optimal x in [0.5, 2.5] so the instance is non-trivial
  const optTarget = 0.5 + Math.random()*2;
  const optY0 = -1 + Math.random()*2; // y-value at optimum

  CS = [];
  // Build upper and lower constraints that cross at (optTarget, optY0)
  // plus random ones that are non-binding at the optimum
  const N = 6; // pairs
  for(let i=0;i<N;i++){
    const a = (Math.random()-0.5)*2;
    const b = 0.5 + Math.random()*1.5;
    // upper: y <= (c-a*x)/b  at (optTarget,optY0): optY0 = (c-a*optTarget)/b => c = optY0*b + a*optTarget
    // add a small positive slack so it's not exactly tight (except one)
    const slack = (i===0) ? 0 : 0.5+Math.random()*1.5; // first upper is tight
    const c = optY0*b + a*optTarget + slack;
    CS.push({a, b, c, type:'upper'});
  }
  for(let i=0;i<N;i++){
    const a = (Math.random()-0.5)*2;
    const b = 0.5 + Math.random()*1.5;
    // lower: y >= (c-a*x)/b; must be satisfied at optimum
    const slack = (i===0) ? 0 : 0.5+Math.random()*1.5; // first lower is tight
    const c = optY0*b + a*optTarget - slack;
    CS.push({a, b, c, type:'lower'});
  }

  active = CS.map((_,i)=>i);
  xLeft = -3; xRight = 3;
  optX = null; optY = null; solved = false; stepN = 0;
  updateInfo();
  render();
  log(fmt('Press Step to watch prune-and-search eliminate constraints. Expected opt x ≈ %optTarget%.', { optTarget: optTarget.toFixed(2) }));
}

// -- one prune step --
function step(){
  if(solved){ log('Already solved — press Reset for a new instance.'); return; }
  const ups  = active.filter(i=>CS[i].type==='upper');
  const lows = active.filter(i=>CS[i].type==='lower');
  if(ups.length===0 || lows.length===0){ finalize(); return; }

  // pair ups with lows and collect pivots inside current x-range
  const pairs = Math.min(ups.length, lows.length);
  const pivots = [];
  for(let k=0;k<pairs;k++){
    const p = pivot(CS, ups[k], lows[k]);
    if(p!==null && p>xLeft-1e-9 && p<xRight+1e-9) pivots.push({x:p, ui:ups[k], li:lows[k]});
  }
  if(!pivots.length){ finalize(); return; }

  const medX = median(pivots.map(p=>p.x));
  const goRight = feasible(active, CS, xRight - 1e-9);

  const toRemove = new Set();
  if(goRight){
    // optimum is in [medX, xRight]; for each pair, if lower < upper on [medX, xRight] => lower non-binding
    for(const {x, ui, li} of pivots){
      if(x <= medX){
        // lower constraint li was already overtaken by upper before medX => non-binding for x >= medX
        toRemove.add(li);
      }
    }
    xLeft = Math.max(xLeft, medX);
  } else {
    // optimum is in [xLeft, medX]; for each pair, if upper < lower on [xLeft, medX] => upper non-binding
    for(const {x, ui, li} of pivots){
      if(x >= medX){
        toRemove.add(ui);
      }
    }
    xRight = Math.min(xRight, medX);
  }

  const removed = toRemove.size;
  active = active.filter(i=>!toRemove.has(i));
  stepN++;
  const dir = goRight ? '→ right' : '← left';
  log(fmt('Step %stepN%: median pivot x=%medX%, searching %dir%, removed %removed% → %length% constraints remain', { stepN: stepN, medX: medX.toFixed(2), dir: dir, removed: removed, length: active.length }));
  updateInfo();
  if(active.length <= 2 || removed===0){ finalize(); return; }
  render();
}

function finalize(){
  // binary search for the rightmost feasible x in [xLeft, xRight]
  let lo=xLeft, hi=xRight;
  for(let i=0;i<80;i++){
    const mid=(lo+hi)/2;
    if(feasible(active,CS,mid)) lo=mid; else hi=mid;
  }
  optX = lo;
  const u=upperAt(active,CS,optX), l=lowerAt(active,CS,optX);
  optY = (u+l)/2;
  solved = true;
  updateInfo();
  render();
  log(fmt('✓ Optimal x ≈ %optX%, y ≈ %optY% — found in %stepN% prune steps.', { optX: optX.toFixed(3), optY: optY.toFixed(3), stepN: stepN }));
  // console check
  console.log(fmt('megiddo-lp demo: optX=%optX% feasible=%feasible% stepN=%stepN%', { optX: optX.toFixed(4), feasible: feasible(active,CS,optX), stepN: stepN }));
}

function updateInfo(){ document.getElementById('nCount').textContent = active.length; }

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

  // axes
  ctx.strokeStyle='#ccd5de'; ctx.lineWidth=1; ctx.setLineDash([]);
  ctx.beginPath(); ctx.moveTo(cx(0),PAD); ctx.lineTo(cx(0),H-PAD); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(PAD,cy(0)); ctx.lineTo(W-PAD,cy(0)); ctx.stroke();

  // axis labels
  ctx.fillStyle='#aab'; ctx.font='10px system-ui';
  ctx.fillText('x',W-PAD+4,cy(0)+4); ctx.fillText('y',cx(0)+4,PAD-4);

  // shade feasible band over current x-range
  const STEPS=300;
  ctx.fillStyle='rgba(29,53,87,0.07)';
  ctx.beginPath();
  let first=true;
  for(let s=0;s<=STEPS;s++){
    const x=xLeft+(xRight-xLeft)*s/STEPS;
    const u=Math.min(upperAt(active,CS,x),yHi), l=Math.max(lowerAt(active,CS,x),yLo);
    if(l>u) continue;
    if(first){ ctx.moveTo(cx(x),cy(u)); first=false; } else ctx.lineTo(cx(x),cy(u));
  }
  for(let s=STEPS;s>=0;s--){
    const x=xLeft+(xRight-xLeft)*s/STEPS;
    const u=Math.min(upperAt(active,CS,x),yHi), l=Math.max(lowerAt(active,CS,x),yLo);
    if(l>u) continue;
    ctx.lineTo(cx(x),cy(l));
  }
  ctx.closePath(); ctx.fill();

  // draw constraints
  CS.forEach((c,i)=>{
    const isActive=active.includes(i);
    ctx.strokeStyle=isActive?(c.type==='upper'?'#1d6fa5':'#c0392b'):'#dde3ea';
    ctx.lineWidth=isActive?1.5:0.7;
    ctx.setLineDash(isActive?[]:[4,3]);
    ctx.beginPath();
    ctx.moveTo(cx(xLo), cy(evalConstraint(c,xLo)));
    ctx.lineTo(cx(xHi), cy(evalConstraint(c,xHi)));
    ctx.stroke();
  });
  ctx.setLineDash([]);

  // x-range brackets
  ctx.strokeStyle='#e67e22'; ctx.lineWidth=2; ctx.setLineDash([5,3]);
  ctx.beginPath(); ctx.moveTo(cx(xLeft),PAD); ctx.lineTo(cx(xLeft),H-PAD); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(cx(xRight),PAD); ctx.lineTo(cx(xRight),H-PAD); ctx.stroke();
  ctx.setLineDash([]);

  // optimum dot
  if(optX!==null && isFinite(optY)){
    ctx.fillStyle='#27ae60'; ctx.strokeStyle='#fff'; ctx.lineWidth=2;
    ctx.beginPath(); ctx.arc(cx(optX),cy(optY),7,0,Math.PI*2); ctx.fill(); ctx.stroke();
    ctx.fillStyle='#27ae60'; ctx.font='bold 11px system-ui';
    ctx.fillText(fmt('opt ≈ (%optX%, %optY%)', { optX: optX.toFixed(2), optY: optY.toFixed(2) }), cx(optX)+10, cy(optY)+4);
  }

  // legend
  ctx.font='10px system-ui';
  ctx.fillStyle='#1d6fa5'; ctx.fillText('— upper (y ≤ ...)', PAD, H-8);
  ctx.fillStyle='#c0392b'; ctx.fillText('— lower (y ≥ ...)', PAD+120, H-8);
  ctx.fillStyle='#e67e22'; ctx.fillText('| search range', PAD+240, H-8);
}

function log(msg){ document.getElementById('log').textContent=msg; }

document.getElementById('btnStep').addEventListener('click', step);
document.getElementById('btnSolve').addEventListener('click', ()=>{ for(let i=0;i<30&&!solved;i++) step(); });
document.getElementById('btnReset').addEventListener('click', generate);

generate();

<div id="controls">
  <button id="btnStep" type="button">Step</button>
  <button id="btnSolve" type="button">Solve instantly</button>
  <button id="btnReset" type="button" class="ghost">Reset</button>
  <span id="info">n = <b id="nCount">0</b> constraints active</span>
</div>
<canvas id="canvas" width="560" height="310"></canvas>
<div id="log"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; background: #f7f9fb; }
#controls { display: flex; gap: .5rem; align-items: center; flex-wrap: wrap; padding: .5rem .5rem .3rem; }
button { font: 600 13px system-ui; padding: .38rem .8rem; border-radius: 7px; cursor: pointer;
         border: 1px solid #1d3557; background: #1d3557; color: #fff; }
button.ghost { background: #fff; color: #1d3557; }
#info { font-size: .85rem; color: #555; margin-left: .3rem; }
canvas { display: block; width: 100%; max-width: 560px; height: auto; border: 1px solid #cdd5de;
         border-radius: 8px; background: #fff; margin: 0 auto; }
#log { font-size: .8rem; color: #555; min-height: 1.6em; padding: .3rem .5rem; text-align: center; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Megiddo's prune-and-search (pedagogical): maximize x over half-planes.
// Each constraint is   y <= (c - a*x)/b   (upper)  or   y >= (c - a*x)/b  (lower).
// The feasible region is the x-range where lowerBound(x) <= upperBound(x).
// We find the rightmost feasible x using prune-and-search:
//   pair each upper with a lower, compute their intersection x (pivot), find median pivot, test feasibility on one side, discard non-binding constraints.

const W = 560, H = 310, PAD = 36;
const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');

// -- coordinate transform --
const xLo=-3.5, xHi=3.8, yLo=-4, yHi=4;
function cx(x){ return PAD + (x-xLo)/(xHi-xLo)*(W-2*PAD); }
function cy(y){ return H - PAD - (y-yLo)/(yHi-yLo)*(H-2*PAD); }

// -- constraint evaluators --
function evalConstraint(c, x){ return (c.c - c.a*x) / c.b; }
function upperAt(active, cs, x){
  let u=Infinity; for(const i of active) if(cs[i].type==='upper') u=Math.min(u,evalConstraint(cs[i],x)); return u;
}
function lowerAt(active, cs, x){
  let l=-Infinity; for(const i of active) if(cs[i].type==='lower') l=Math.max(l,evalConstraint(cs[i],x)); return l;
}
function feasible(active, cs, x){ return lowerAt(active,cs,x) <= upperAt(active,cs,x)+1e-9; }

// -- state --
let CS=[], active=[], xLeft, xRight, optX=null, optY=null, solved=false, stepN=0;

// -- pivot: x where upper constraint i equals lower constraint j --
function pivot(cs, i, j){
  const u=cs[i], v=cs[j];
  // evalU(x) = evalV(x)  =>  (u.c-u.a*x)/u.b = (v.c-v.a*x)/v.b
  // v.b*(u.c-u.a*x) = u.b*(v.c-v.a*x)
  // x*(u.b*v.a - v.b*u.a) = u.b*v.c - v.b*u.c
  const denom = u.b*v.a - v.b*u.a;
  if(Math.abs(denom)<1e-9) return null;
  return (u.b*v.c - v.b*u.c)/denom;
}

function median(arr){
  const s=arr.slice().sort((a,b)=>a-b); return s[Math.floor(s.length/2)];
}

// -- generate instance guaranteed feasible in [-3, optXTarget] --
function generate(){
  // pick a random optimal x in [0.5, 2.5] so the instance is non-trivial
  const optTarget = 0.5 + Math.random()*2;
  const optY0 = -1 + Math.random()*2; // y-value at optimum

  CS = [];
  // Build upper and lower constraints that cross at (optTarget, optY0)
  // plus random ones that are non-binding at the optimum
  const N = 6; // pairs
  for(let i=0;i<N;i++){
    const a = (Math.random()-0.5)*2;
    const b = 0.5 + Math.random()*1.5;
    // upper: y <= (c-a*x)/b  at (optTarget,optY0): optY0 = (c-a*optTarget)/b => c = optY0*b + a*optTarget
    // add a small positive slack so it's not exactly tight (except one)
    const slack = (i===0) ? 0 : 0.5+Math.random()*1.5; // first upper is tight
    const c = optY0*b + a*optTarget + slack;
    CS.push({a, b, c, type:'upper'});
  }
  for(let i=0;i<N;i++){
    const a = (Math.random()-0.5)*2;
    const b = 0.5 + Math.random()*1.5;
    // lower: y >= (c-a*x)/b; must be satisfied at optimum
    const slack = (i===0) ? 0 : 0.5+Math.random()*1.5; // first lower is tight
    const c = optY0*b + a*optTarget - slack;
    CS.push({a, b, c, type:'lower'});
  }

  active = CS.map((_,i)=>i);
  xLeft = -3; xRight = 3;
  optX = null; optY = null; solved = false; stepN = 0;
  updateInfo();
  render();
  log(fmt('Press Step to watch prune-and-search eliminate constraints. Expected opt x ≈ %optTarget%.', { optTarget: optTarget.toFixed(2) }));
}

// -- one prune step --
function step(){
  if(solved){ log('Already solved — press Reset for a new instance.'); return; }
  const ups  = active.filter(i=>CS[i].type==='upper');
  const lows = active.filter(i=>CS[i].type==='lower');
  if(ups.length===0 || lows.length===0){ finalize(); return; }

  // pair ups with lows and collect pivots inside current x-range
  const pairs = Math.min(ups.length, lows.length);
  const pivots = [];
  for(let k=0;k<pairs;k++){
    const p = pivot(CS, ups[k], lows[k]);
    if(p!==null && p>xLeft-1e-9 && p<xRight+1e-9) pivots.push({x:p, ui:ups[k], li:lows[k]});
  }
  if(!pivots.length){ finalize(); return; }

  const medX = median(pivots.map(p=>p.x));
  const goRight = feasible(active, CS, xRight - 1e-9);

  const toRemove = new Set();
  if(goRight){
    // optimum is in [medX, xRight]; for each pair, if lower < upper on [medX, xRight] => lower non-binding
    for(const {x, ui, li} of pivots){
      if(x <= medX){
        // lower constraint li was already overtaken by upper before medX => non-binding for x >= medX
        toRemove.add(li);
      }
    }
    xLeft = Math.max(xLeft, medX);
  } else {
    // optimum is in [xLeft, medX]; for each pair, if upper < lower on [xLeft, medX] => upper non-binding
    for(const {x, ui, li} of pivots){
      if(x >= medX){
        toRemove.add(ui);
      }
    }
    xRight = Math.min(xRight, medX);
  }

  const removed = toRemove.size;
  active = active.filter(i=>!toRemove.has(i));
  stepN++;
  const dir = goRight ? '→ right' : '← left';
  log(fmt('Step %stepN%: median pivot x=%medX%, searching %dir%, removed %removed% → %length% constraints remain', { stepN: stepN, medX: medX.toFixed(2), dir: dir, removed: removed, length: active.length }));
  updateInfo();
  if(active.length <= 2 || removed===0){ finalize(); return; }
  render();
}

function finalize(){
  // binary search for the rightmost feasible x in [xLeft, xRight]
  let lo=xLeft, hi=xRight;
  for(let i=0;i<80;i++){
    const mid=(lo+hi)/2;
    if(feasible(active,CS,mid)) lo=mid; else hi=mid;
  }
  optX = lo;
  const u=upperAt(active,CS,optX), l=lowerAt(active,CS,optX);
  optY = (u+l)/2;
  solved = true;
  updateInfo();
  render();
  log(fmt('✓ Optimal x ≈ %optX%, y ≈ %optY% — found in %stepN% prune steps.', { optX: optX.toFixed(3), optY: optY.toFixed(3), stepN: stepN }));
  // console check
  console.log(fmt('megiddo-lp demo: optX=%optX% feasible=%feasible% stepN=%stepN%', { optX: optX.toFixed(4), feasible: feasible(active,CS,optX), stepN: stepN }));
}

function updateInfo(){ document.getElementById('nCount').textContent = active.length; }

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

  // axes
  ctx.strokeStyle='#ccd5de'; ctx.lineWidth=1; ctx.setLineDash([]);
  ctx.beginPath(); ctx.moveTo(cx(0),PAD); ctx.lineTo(cx(0),H-PAD); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(PAD,cy(0)); ctx.lineTo(W-PAD,cy(0)); ctx.stroke();

  // axis labels
  ctx.fillStyle='#aab'; ctx.font='10px system-ui';
  ctx.fillText('x',W-PAD+4,cy(0)+4); ctx.fillText('y',cx(0)+4,PAD-4);

  // shade feasible band over current x-range
  const STEPS=300;
  ctx.fillStyle='rgba(29,53,87,0.07)';
  ctx.beginPath();
  let first=true;
  for(let s=0;s<=STEPS;s++){
    const x=xLeft+(xRight-xLeft)*s/STEPS;
    const u=Math.min(upperAt(active,CS,x),yHi), l=Math.max(lowerAt(active,CS,x),yLo);
    if(l>u) continue;
    if(first){ ctx.moveTo(cx(x),cy(u)); first=false; } else ctx.lineTo(cx(x),cy(u));
  }
  for(let s=STEPS;s>=0;s--){
    const x=xLeft+(xRight-xLeft)*s/STEPS;
    const u=Math.min(upperAt(active,CS,x),yHi), l=Math.max(lowerAt(active,CS,x),yLo);
    if(l>u) continue;
    ctx.lineTo(cx(x),cy(l));
  }
  ctx.closePath(); ctx.fill();

  // draw constraints
  CS.forEach((c,i)=>{
    const isActive=active.includes(i);
    ctx.strokeStyle=isActive?(c.type==='upper'?'#1d6fa5':'#c0392b'):'#dde3ea';
    ctx.lineWidth=isActive?1.5:0.7;
    ctx.setLineDash(isActive?[]:[4,3]);
    ctx.beginPath();
    ctx.moveTo(cx(xLo), cy(evalConstraint(c,xLo)));
    ctx.lineTo(cx(xHi), cy(evalConstraint(c,xHi)));
    ctx.stroke();
  });
  ctx.setLineDash([]);

  // x-range brackets
  ctx.strokeStyle='#e67e22'; ctx.lineWidth=2; ctx.setLineDash([5,3]);
  ctx.beginPath(); ctx.moveTo(cx(xLeft),PAD); ctx.lineTo(cx(xLeft),H-PAD); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(cx(xRight),PAD); ctx.lineTo(cx(xRight),H-PAD); ctx.stroke();
  ctx.setLineDash([]);

  // optimum dot
  if(optX!==null && isFinite(optY)){
    ctx.fillStyle='#27ae60'; ctx.strokeStyle='#fff'; ctx.lineWidth=2;
    ctx.beginPath(); ctx.arc(cx(optX),cy(optY),7,0,Math.PI*2); ctx.fill(); ctx.stroke();
    ctx.fillStyle='#27ae60'; ctx.font='bold 11px system-ui';
    ctx.fillText(fmt('opt ≈ (%optX%, %optY%)', { optX: optX.toFixed(2), optY: optY.toFixed(2) }), cx(optX)+10, cy(optY)+4);
  }

  // legend
  ctx.font='10px system-ui';
  ctx.fillStyle='#1d6fa5'; ctx.fillText('— upper (y ≤ ...)', PAD, H-8);
  ctx.fillStyle='#c0392b'; ctx.fillText('— lower (y ≥ ...)', PAD+120, H-8);
  ctx.fillStyle='#e67e22'; ctx.fillText('| search range', PAD+240, H-8);
}

function log(msg){ document.getElementById('log').textContent=msg; }

document.getElementById('btnStep').addEventListener('click', step);
document.getElementById('btnSolve').addEventListener('click', ()=>{ for(let i=0;i<30&&!solved;i++) step(); });
document.getElementById('btnReset').addEventListener('click', generate);

generate();

Result

Console

Press Step to watch the prune-and-search strategy eliminate constraints one batch at a time. Notice how each step cuts the active constraint count roughly in half — that halving is the heart of the $O(n)$ guarantee. Press Solve instantly to jump straight to the optimum, or Reset to start over with a new random instance.

The Real Complexity

How fast can you solve LP in two dimensions? Megiddo's answer is surprising even today.

The naive bound is $O(n^{2})$ : try all $O(n^{2})$ vertices of the constraint polytope and pick the best. Faster approaches reduced this, but not to linear.
Megiddo's key insight (1984): pair up the n constraints arbitrarily. For each pair, compute a pivot value — the x-coordinate at which one constraint "takes over" from the other in limiting the objective. These $O(n)$ pivot values can be median-found in $O(n)$ time. Then determine (in $O(n)$ time) whether the optimum lies to the left or right of the median pivot; this single decision lets you discard at least n/4 constraints as permanently non-binding in any direction. Recurse.
Recurrence: T(n) = T(3n/4) + $O(n)$ , which solves to T(n) = $O(n)$ — linear time.
Extension to d fixed dimensions: the same prune-and-search idea extends to any fixed d, giving $O(n)$ with the constant depending exponentially on d. For d = 3 Megiddo proved $O(n)$ as well; Seidel later gave a randomized $O(d! \cdot n)$ algorithm, and Clarkson an expected $O(d^{2}n)$ one.
Status: LP in fixed dimension is solved — $O(n)$ is a lower bound too (you must read all constraints), so Megiddo's result is optimal in dimension 2 and 3. In high dimension the best deterministic bound is still super-linear and the precise complexity remains open.

The result also inspired the broader prune-and-search paradigm, used today across computational geometry for problems like smallest enclosing balls, ham-sandwich cuts, and weighted median finding. See also linear programming and selection and median finding.

Where It Matters

Linear time LP in low dimensions is not just a theoretical prize — it changes what is feasible in practice:

Computational geometry: smallest enclosing disk/ball, halfspace intersection, linear separability, and ham-sandwich cuts all reduce to low-dimensional LP and benefit directly from $O(n)$ algorithms.
Robot motion planning: collision avoidance with convex obstacles in 2D or 3D reduces to linear programs with a bounded number of variables, making real-time guarantees achievable.
Embedded and real-time systems: when constraints arrive in a stream and the dimension is fixed (e.g. two resource limits), $O(n)$ LP means you can certify feasibility without expensive solvers.
Randomized LP algorithms: Megiddo's framework motivated Seidel's and Clarkson's randomized $O(n)$ algorithms, which are simpler to implement and widely used in practice for 2D/3D problems.
Teaching algorithm design: the prune-and-search pattern — pair, pivot, discard, recurse — is a clean paradigm that sits alongside divide-and-conquer and the master theorem in any advanced algorithms course.

Understand Megiddo's algorithm and you hold the key to a wide family of efficient geometry algorithms that underpin graphics engines, GIS tools, and planning software worldwide.

Conclusion

Megiddo's 1984 algorithm carries a clean message: polynomial time is not the end of the story. General LP is polynomial, but by fixing the dimension to 2 or 3 the same problem collapses to linear time — the best possible, since you cannot avoid reading every constraint.

The prune-and-search paradigm behind the result has since become a cornerstone of computational geometry, appearing wherever a divide-and-discard strategy can exploit the low dimension of the problem. It is a reminder that structure in a problem's geometry can be as powerful as any algebraic insight.

If you want to dig deeper, explore how linear programming scales in higher dimensions, or how selection and median finding — the linear-time subroutine at the heart of Megiddo's method — works in its own right.

Comments

Loading comments...

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

Cookies & privacy

Linear Programming in Linear Time

Megiddo's prune-and-search: an $O(n)$ algorithm for 2D LP