Introduction

Suppose you want to find a sparse solution to an underdetermined system of equations — fewer unknowns on are nonzero than the total number of variables. This is the central challenge in compressed sensing, medical imaging and signal processing. The natural objective adds a smooth data-fit term and a nonsmooth L1 penalty that promotes sparsity:

F(x)=12∥Ax−b∥2⏟f(x) smooth+λ∥x∥1⏟g(x) nonsmoothF(x) = \underbrace{\tfrac{1}{2}\|Ax - b\|^2}_{f(x)\text{ smooth}} + \underbrace{\lambda \|x\|_1}_{g(x)\text{ nonsmooth}}

Gradient descent works beautifully on the smooth part ff, but it cannot handle gg because the L1 norm has corners — points where no gradient exists. Proximal gradient methods split the problem cleanly: take a gradient step in ff, then solve a tiny standalone minimization problem for gg. That second step is called the proximal operator, and for the L1 norm it has a simple closed form: soft-thresholding.

The basic algorithm is ISTA (Iterative Shrinkage-Thresholding Algorithm). Its convergence rate is O(1/k)O(1/k). In 2009, Beck and Teboulle showed that adding a single momentum term — extrapolating the current iterate using the previous one — yields FISTA (Fast ISTA) with rate O(1/k2)O(1/k^2), at essentially zero extra cost.

Try It: Sparse Recovery

Below is a 1-D sparse recovery problem: the algorithm must recover a sparse vector xx from noisy measurements. Use the sliders to change the L1 penalty λ\lambda and the step size tt, then press Run to watch ISTA and FISTA converge. The left panel tracks the objective value; the right panel shows the current estimate vs. the true sparse signal.

<!-- Proximal Gradient demo: ISTA vs FISTA on a sparse 1-D recovery problem -->
<div class="controls">
  <label>λ (L1 penalty): <span id="lval">0.10</span>
    <input type="range" id="lambda" min="0.01" max="0.5" step="0.01" value="0.10">
  </label>
  <label>t (step size): <span id="tval">0.10</span>
    <input type="range" id="step" min="0.01" max="0.19" step="0.01" value="0.10">
  </label>
  <label>Iterations: <span id="kval">60</span>
    <input type="range" id="iters" min="10" max="150" step="5" value="60">
  </label>
</div>
<div class="btn-row">
  <button id="run">Run</button>
  <button id="reset" class="ghost">Reset</button>
</div>
<div class="charts">
  <canvas id="objChart" width="280" height="200" title="Objective value vs iterations"></canvas>
  <canvas id="sigChart" width="280" height="200" title="Estimated vs true signal"></canvas>
</div>
<div id="status" class="status"></div>

/* Layout and chart container styles */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: 4px; }
.controls { display: flex; flex-direction: column; gap: .4rem; margin-bottom: .6rem; }
label { font-size: .85rem; display: flex; align-items: center; gap: .5rem; }
input[type=range] { flex: 1; accent-color: #1d3557; }
.btn-row { display: flex; gap: .5rem; margin-bottom: .7rem; }
button { font: 600 14px system-ui; padding: .4rem .9rem;
         border: 1px solid #1d3557; background: #1d3557; color: #fff;
         border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
.charts { display: flex; gap: 8px; flex-wrap: wrap; }
canvas { border: 1px solid #d0d7de; border-radius: 6px; background: #fafbfc; }
.status { font-size: .88rem; font-weight: 600; margin-top: .5rem; min-height: 1.3em; color: #0a7d33; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Sparse 1-D recovery: minimize (1/2)||Ax-b||^2 + lambda*||x||_1
// True signal TRUE_X is sparse; A is a random measurement matrix with M < N rows.

// Ground-truth sparse signal (10 coordinates, 3 nonzero)
const TRUE_X = [0, 0, 1.5, 0, 0, -2.0, 0, 0, 0.8, 0];
const N = TRUE_X.length;

// Random measurement matrix A (M x N) with M < N (underdetermined).
// We use a simple LCG for reproducibility (no external library).
const SEED = 42;
function lcg(s){ return (s * 1664525 + 1013904223) & 0xffffffff; }
let rng = SEED;
function rand(){ rng = lcg(rng); return (rng >>> 0) / 0xffffffff - 0.5; }

const M = 7; // 7 measurements for 10 unknowns: underdetermined system
const A = [];
for(let i = 0; i < M; i++){ const row = []; for(let j = 0; j < N; j++) row.push(rand()); A.push(row); }

// Normalize columns of A to unit L2 norm for conditioning
for(let j = 0; j < N; j++){
  let norm = 0; for(let i = 0; i < M; i++) norm += A[i][j]*A[i][j]; norm = Math.sqrt(norm);
  for(let i = 0; i < M; i++) A[i][j] /= norm;
}

// Compute noisy measurements b = A * TRUE_X + small noise
const B = [];
for(let i = 0; i < M; i++){
  let s = 0; for(let j = 0; j < N; j++) s += A[i][j] * TRUE_X[j];
  B.push(s + 0.02 * rand());
}

// Matrix-vector product helper
function matvec(mat, x){ return mat.map(row => row.reduce((s,v,j) => s + v*x[j], 0)); }
// Residual r = A*x - b
function residual(x){ return matvec(A, x).map((v,i) => v - B[i]); }
// Gradient of smooth part f(x) = (1/2)||Ax-b||^2
function grad_f(x){ const r = residual(x); const g = new Array(N).fill(0);
  for(let i=0; i<M; i++) for(let j=0; j<N; j++) g[j] += A[i][j]*r[i]; return g; }
// L1 norm (sum of absolute values)
function l1(x){ return x.reduce((s,v) => s + Math.abs(v), 0); }
// Full composite objective F(x) = f(x) + lambda*||x||_1
function objective(x, lam){ const r = residual(x); const sq = r.reduce((s,v) => s+v*v, 0);
  return 0.5*sq + lam*l1(x); }
// Soft-thresholding: proximal operator for lambda*||.||_1
function softThreshold(x, thresh){ return x.map(v => Math.sign(v) * Math.max(Math.abs(v) - thresh, 0)); }

// One ISTA step: gradient step on f, then soft-threshold for g
function istaStep(x, t, lam){
  const g = grad_f(x);
  const y = x.map((v,j) => v - t*g[j]);
  return softThreshold(y, lam*t);
}

// Run ISTA and FISTA for K iterations, return objective histories and final iterates
function runFISTA(lam, t, K){
  let x = new Array(N).fill(0);
  let xPrev = new Array(N).fill(0);
  let theta = 1;
  const objISTA = [], objFISTA = [];
  let xISTA = new Array(N).fill(0);

  for(let k = 0; k < K; k++){
    // ISTA plain gradient step
    xISTA = istaStep(xISTA, t, lam);
    objISTA.push(objective(xISTA, lam));

    // FISTA momentum update: extrapolate with previous iterate
    const xNew = istaStep(x, t, lam);
    const thetaNew = (1 + Math.sqrt(1 + 4*theta*theta)) / 2;
    const beta = (theta - 1) / thetaNew;
    const y = xNew.map((v,j) => v + beta*(v - xPrev[j]));
    xPrev = xNew; x = y; theta = thetaNew;
    objFISTA.push(objective(xNew, lam));
  }
  return { objISTA, objFISTA, xISTA, xFISTA: xPrev };
}

let currentResult = null;

function draw(lam, t, K){
  const res = runFISTA(lam, t, K);
  currentResult = res;
  drawObj(res.objISTA, res.objFISTA, K);
  drawSig(res.xFISTA);
  const finalISTA = res.objISTA[K-1].toFixed(4);
  const finalFISTA = res.objFISTA[K-1].toFixed(4);
  document.getElementById('status').textContent =
    fmt('After %K% iterations — objective:: ISTA=%finalISTA%  FISTA=%finalFISTA%', { K: K, finalISTA: finalISTA, finalFISTA: finalFISTA });
}

function drawObj(istaVals, fistaVals, K){
  const canvas = document.getElementById('objChart');
  const ctx = canvas.getContext('2d');
  const W = canvas.width, H = canvas.height;
  ctx.clearRect(0,0,W,H);
  const pad = { l:38, r:8, t:18, b:28 };
  const pw = W - pad.l - pad.r, ph = H - pad.t - pad.b;

  const allVals = [...istaVals, ...fistaVals];
  const minV = Math.min(...allVals), maxV = Math.max(...allVals);
  const range = maxV - minV || 1;

  // Draw coordinate axes
  ctx.strokeStyle = '#888'; ctx.lineWidth = 1;
  ctx.beginPath(); ctx.moveTo(pad.l, pad.t); ctx.lineTo(pad.l, pad.t+ph); ctx.lineTo(pad.l+pw, pad.t+ph); ctx.stroke();

  // Draw ISTA convergence curve (blue)
  ctx.strokeStyle = '#1d6aaa'; ctx.lineWidth = 1.5;
  ctx.beginPath();
  istaVals.forEach((v,i) => { const x = pad.l + (i/(K-1))*pw; const y = pad.t + ph - ((v-minV)/range)*ph;
    i===0 ? ctx.moveTo(x,y) : ctx.lineTo(x,y); });
  ctx.stroke();

  // Draw FISTA convergence curve (orange)
  ctx.strokeStyle = '#e07b39'; ctx.lineWidth = 1.5;
  ctx.beginPath();
  fistaVals.forEach((v,i) => { const x = pad.l + (i/(K-1))*pw; const y = pad.t + ph - ((v-minV)/range)*ph;
    i===0 ? ctx.moveTo(x,y) : ctx.lineTo(x,y); });
  ctx.stroke();

  // Draw axis labels and legend
  ctx.fillStyle = '#444'; ctx.font = '10px system-ui'; ctx.textAlign = 'center';
  ctx.fillText('Iteration', pad.l + pw/2, H - 4);
  ctx.save(); ctx.translate(10, pad.t + ph/2); ctx.rotate(-Math.PI/2);
  ctx.fillText('F(x)', 0, 0); ctx.restore();
  ctx.fillStyle = '#1d6aaa'; ctx.fillText('ISTA', pad.l + pw*0.25, pad.t+12);
  ctx.fillStyle = '#e07b39'; ctx.fillText('FISTA', pad.l + pw*0.65, pad.t+12);
}

function drawSig(xFISTA){
  const canvas = document.getElementById('sigChart');
  const ctx = canvas.getContext('2d');
  const W = canvas.width, H = canvas.height;
  ctx.clearRect(0,0,W,H);
  const pad = { l:38, r:8, t:18, b:28 };
  const pw = W - pad.l - pad.r, ph = H - pad.t - pad.b;

  const allVals = [...TRUE_X, ...xFISTA];
  const minV = Math.min(...allVals) - 0.3, maxV = Math.max(...allVals) + 0.3;
  const range = maxV - minV || 1;
  const bar = pw / (N * 2.5);

  // Draw horizontal zero line
  ctx.strokeStyle = '#ccc'; ctx.lineWidth = 1;
  const y0 = pad.t + ph - ((-minV)/range)*ph;
  ctx.beginPath(); ctx.moveTo(pad.l, y0); ctx.lineTo(pad.l+pw, y0); ctx.stroke();

  for(let j = 0; j < N; j++){
    const cx = pad.l + (j+0.5)/(N) * pw;
    // True signal bar (dark blue, semi-transparent)
    ctx.fillStyle = 'rgba(29,53,87,0.35)';
    const yt = pad.t + ph - ((TRUE_X[j]-minV)/range)*ph;
    ctx.fillRect(cx - bar - 1, Math.min(yt, y0), bar, Math.abs(yt - y0));
    // FISTA estimate bar (orange)
    ctx.fillStyle = 'rgba(224,123,57,0.85)';
    const ye = pad.t + ph - ((xFISTA[j]-minV)/range)*ph;
    ctx.fillRect(cx + 1, Math.min(ye, y0), bar, Math.abs(ye - y0));
  }

  // Draw coordinate axes for signal chart
  ctx.strokeStyle = '#888'; ctx.lineWidth = 1;
  ctx.beginPath(); ctx.moveTo(pad.l, pad.t); ctx.lineTo(pad.l, pad.t+ph); ctx.lineTo(pad.l+pw, pad.t+ph); ctx.stroke();

  ctx.fillStyle = '#444'; ctx.font = '10px system-ui'; ctx.textAlign = 'center';
  ctx.fillText('Coordinate index', pad.l + pw/2, H - 4);
  ctx.fillStyle = 'rgba(29,53,87,0.7)'; ctx.fillText('True', pad.l + pw*0.25, pad.t+12);
  ctx.fillStyle = 'rgba(224,123,57,0.9)'; ctx.fillText('FISTA est.', pad.l + pw*0.7, pad.t+12);
}

function readParams(){ return {
  lam: parseFloat(document.getElementById('lambda').value),
  t:   parseFloat(document.getElementById('step').value),
  K:   parseInt(document.getElementById('iters').value)
}; }

document.getElementById('lambda').addEventListener('input', e => {
  document.getElementById('lval').textContent = parseFloat(e.target.value).toFixed(2); });
document.getElementById('step').addEventListener('input', e => {
  document.getElementById('tval').textContent = parseFloat(e.target.value).toFixed(2); });
document.getElementById('iters').addEventListener('input', e => {
  document.getElementById('kval').textContent = e.target.value; });

document.getElementById('run').addEventListener('click', () => {
  const p = readParams(); draw(p.lam, p.t, p.K); });
document.getElementById('reset').addEventListener('click', () => {
  document.getElementById('lambda').value = 0.10;
  document.getElementById('step').value   = 0.10;
  document.getElementById('iters').value  = 60;
  document.getElementById('lval').textContent = '0.10';
  document.getElementById('tval').textContent = '0.10';
  document.getElementById('kval').textContent = '60';
  document.getElementById('status').textContent = '';
  const c1 = document.getElementById('objChart'), c2 = document.getElementById('sigChart');
  c1.getContext('2d').clearRect(0,0,c1.width,c1.height);
  c2.getContext('2d').clearRect(0,0,c2.width,c2.height);
});

// Run automatically on page load with default parameters
draw(0.10, 0.10, 60);
Result
Console

Notice how FISTA (orange) consistently reaches low objective values in far fewer steps than ISTA (blue). The gap widens when the problem is poorly conditioned. Soft-thresholding is what drives the sparsity: components smaller than λt\lambda t are zeroed out at every iteration, and the rest are shrunk toward zero by the same amount.

The Real Complexity

How fast do these algorithms converge?

  • ISTA performs, at each iteration kk, two cheap operations: a gradient step x←y−t∇f(y)x \leftarrow y - t \nabla f(y), then a soft-threshold xi←sign(xi)max⁡(∣xi∣−λt,0)x_i \leftarrow \text{sign}(x_i)\max(|x_i| - \lambda t, 0). After kk steps the objective gap satisfies F(xk)−F∗≤O(1/k)F(x_k) - F^* \leq O(1/k).
  • FISTA adds a momentum coefficient θk\theta_k and an extrapolated point yk+1=xk+θk−1θk+1(xk−xk−1)y_{k+1} = x_k + \tfrac{\theta_k - 1}{\theta_{k+1}}(x_k - x_{k-1}), where θk=1+1+4θk−122\theta_k = \tfrac{1 + \sqrt{1 + 4\theta_{k-1}^2}}{2}. This lifts the rate to O(1/k2)O(1/k^2) with no extra gradient evaluations.
  • Optimal in the first-order oracle model. Nesterov's lower bound (1983) shows that no algorithm using only gradient information can do better than O(1/k2)O(1/k^2) for smooth convex objectives. FISTA matches this bound even though gg is nonsmooth, because the proximal step counts as exact minimization of gg.
  • Step size matters. The step size tt must satisfy t≤1/Lt \leq 1/L, where LL is the Lipschitz constant of ∇f\nabla f. Too large a step causes divergence; too small slows convergence. Line-search variants remove the need to estimate LL.

The elegance of the method is that O(1/k2)O(1/k^2) acceleration comes entirely from reusing the previous iterate — not from second-order information or a more expensive oracle. Compare to interior-point methods for linear programming, which need O(n)O(\sqrt{n}) iterations but each costs O(n3)O(n^3) per step.

Where It Matters

Proximal gradient methods are the workhorses of modern large-scale optimization wherever sparse structure is expected or desired:

  • LASSO regression: the canonical example. Minimizing 12∥Ax−b∥2+λ∥x∥1\tfrac{1}{2}\|Ax-b\|^2 + \lambda\|x\|_1 selects a sparse subset of features automatically. ISTA/FISTA solve it efficiently even for millions of variables.
  • Compressed sensing and MRI: accelerated MRI acquires far fewer measurements than the Nyquist rate, then reconstructs the image by promoting sparsity in a wavelet domain — solved by FISTA.
  • Image denoising (TV regularization): total-variation denoising is a proximal-gradient problem; the proximal step for the TV norm is a projection onto an ℓ∞\ell_\infty ball, solved by Chambolle's algorithm.
  • Neural network pruning: weight sparsity during training can be encouraged by proximal-gradient updates on the L1 penalty, yielding smaller and faster deployed models.
  • Portfolio optimization: sparse portfolio selection — holding few assets — maps directly to an L1-penalized quadratic, solvable by the same two-step recipe.

The method extends naturally to matrix variables (nuclear-norm minimization for low-rank matrix completion, used in recommendation systems) and to block-separable penalties where each coordinate group has its own proximal operator. Related algorithms like ADMM and the Douglas–Rachford splitting address the same composite structure with different decompositions — see non-convex optimization for what happens when ff loses convexity.

Conclusion

Proximal gradient methods reduce a hard composite problem to a simple two-step recipe: slide down the smooth slope, then apply a closed-form operator that handles the nonsmooth penalty in one shot. For L1 regularization that operator is soft-thresholding — among the simplest formulas in all of optimization.

FISTA shows that momentum is not a heuristic trick but a mathematically certified accelerator: one extra subtraction per step squares the convergence rate from O(1/k)O(1/k) to O(1/k2)O(1/k^2), matching the theoretical lower bound for first-order methods. The result is that compressed sensing, LASSO regression, and image reconstruction problems with millions of variables can now be solved in seconds on a laptop.

The deeper lesson is structural: many practical objectives are smooth-plus-nonsmooth composites, and splitting them along that seam — rather than trying to treat the whole thing with one tool — unlocks both theoretical guarantees and practical speed. That splitting idea ripples through all of modern non-convex optimization and beyond.

Share this article

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

Comments

Loading comments...

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