Orthogonal Matching Pursuit & Grouped Sparsity
Greedy sparse recovery, and extending sparsity to whole groups of coefficients.
On this page
Beginner: Lasso (2.3.3) gets sparse coefficients almost as a side effect — it minimizes squared error plus an L1 penalty, and coordinate descent's soft-thresholding happens to snap many coefficients to exactly zero along the way. Orthogonal Matching Pursuit (OMP) takes the opposite strategy: it builds a sparse model on purpose, one feature at a time, like assembling a small team by always recruiting whoever best covers whatever the current team is still missing, then re-auditioning the whole team together before deciding who to recruit next.
Intermediate: concretely, OMP starts with a residual equal to the full target y (nothing explained yet). At each round it asks: "of all the features not yet in my active set, which one is most correlated with what I still haven't explained?" It adds that single feature to the active set, then throws away its previous coefficients entirely and re-fits ordinary least squares (OLS) on every feature currently in the active set jointly — not just the new one. The residual is recomputed from that fresh joint fit, and the loop repeats. Stop after a fixed number of steps, or once the residual is small enough.
Advanced: that full re-fit at every step is what separates OMP from the plainer "Matching Pursuit" it's named after, which only nudges the coefficient along the newly chosen direction and never revisits earlier ones. The Derivation below proves exactly what that extra re-fit buys you: a structural orthogonality guarantee that plain Matching Pursuit doesn't have.
Advanced — beyond individual coefficients: Lasso and OMP both operate at the granularity of a single coefficient at a time being zero or not. Sometimes that's the wrong unit of sparsity. If a categorical variable like "region" is one-hot encoded into five dummy columns, it rarely makes sense for the model to keep two of those dummies and drop the other three — either "region" matters to the outcome or it doesn't, as a whole variable. Group Lasso extends the L1-sparsity idea from individual coefficients to whole groups of coefficients, so an entire group is switched on or off together. Fused Lasso applies a related structured-penalty idea to ordered coefficients (like time steps or spatial positions), encouraging neighboring coefficients to be equal rather than merely small. Multi-task Lasso ties several related regression problems together so they're encouraged to agree on which features matter at all, even if the coefficient values themselves differ per task.
OMP's greedy selection rule: among features not yet in the active set S, pick whichever column x_j has the largest-magnitude dot product with the current residual r — the feature most correlated with whatever is still unexplained.
The Group Lasso penalty: partition the coefficients into G predefined groups, take the (un-squared) L2 norm of each group's sub-vector, and sum — an "L1 of L2s." L1 across groups is what drives whole groups to exactly zero; L2 within a group is what keeps a non-zero group's coefficients from being individually sparse.
Fused Lasso adds this penalty on consecutive differences (on top of an ordinary L1 term) — it is soft-thresholding applied to differences between neighbors instead of to the coefficients themselves, which is what makes flat, piecewise-constant runs of equal coefficients cheap under the penalty.
Multi-task Lasso: T related regression tasks share the same feature index j, and the group in the penalty is "coefficient j across all T tasks" — the same Group Lasso mechanism as above, just with the groups defined across tasks instead of within one task's dummy variables.
Part 1 — the residual is orthogonal to every selected feature. Recall the normal equations from ordinary least squares (2.3.1): whenever you fit OLS on some design matrix X_S (the columns in active set S) and target y, the resulting residual r = y − X_S θ_S satisfies:
That's just the stationarity condition of the least-squares minimizer, stated as a fact about OLS in general — it holds for any design matrix you feed OLS, not anything specific to OMP. Read column-by-column, it says: after an OLS fit, the residual is orthogonal to every column of whatever matrix was just fit, i.e. x_j^Tr = 0 for every j in S.
Now apply this fact to OMP directly. At step k, OMP picks a new feature j* (the argmax-correlation rule above), forms the enlarged active set S' = S ∪ {j*}, and re-fits OLS on X_S' from scratch. The new residual r' is the residual of that OLS fit, so the same normal- equations fact applies immediately: X_S'^T r' = 0, meaning r' is orthogonal to every column in S' — every previously selected feature and the brand-new one, simultaneously.
This is exactly why the full re-fit matters, and exactly what the plainer "Matching Pursuit" algorithm gives up by skipping it. Plain Matching Pursuit only updates the coefficient along the newly chosen direction and leaves the older coefficients alone. If the features aren't mutually orthogonal (the typical case), that partial update generally leaves the residual with some leftover correlation against features selected several steps ago — so a later step can end up re-selecting an already-chosen feature, or something nearly redundant with one, wasting a step re-explaining variance that was already explained. OMP's full re-fit drives that correlation to exactly zero every single step, for every selected feature at once — so an already-selected feature can never again be the argmax of |x_j^Tr| (its score is pinned at zero) unless the whole residual is already zero, in which case there's nothing left to explain anyway.
Where this is used: this orthogonality property is the backbone of the exact-recovery theory behind OMP in compressed sensing — it's the mechanical reason greedy sparse-recovery algorithms can be proven, under the right conditions on X, to identify the true sparse support in a bounded number of steps rather than drifting or cycling.
Part 2 — why Group Lasso zeros out whole groups. This runs the same subgradient argument as the scalar L1 case (2.2.5), one level up: from scalars to vectors. Consider a single group's local sub-problem — the proximal step, minimizing ½‖θ_g − b_g‖² + λ‖θ_g‖₂ for a group's coefficient sub-vector θ_g, where b_g plays the same role the scalar correlation b played in the L1 derivation: everything the data-fit gradient contributes for that group, holding every other group fixed.
The L2 norm ‖θ_g‖₂ is not differentiable at θ_g = 0, exactly like |θ| wasn't at the scalar origin — but its subgradient set at that point is now the entire closed unit ball, every vector u with ‖u‖₂ ≤ 1, the direct multivariate generalization of the scalar interval [−1, 1].
Optimality at θ_g = 0 requires zero to be a valid subgradient of the whole objective there: 0 ∈ (0 − b_g) + λ·{u : ‖u‖≤1}, i.e. b_g itself must lie inside the ball of radius λ. That's a single condition on the entire group's gradient vector at once — ‖b_g‖₂ ≤ λ — not a separate condition per coordinate. So the whole group collapses to the zero vector together, or none of it does; there's no subgradient-consistent way for part of a group to sit at zero while the rest doesn't.
When ‖b_g‖₂ > λ instead, the quadratic term pulls θ_g toward b_g's direction, while the penalty term only cares about magnitude, not direction — so by symmetry the minimizer must point exactly along b_g. Writing θ_g = c · b_g/‖b_g‖ for a scalar c ≥ 0 reduces the whole sub-problem to a 1D minimization over c, whose stationarity condition gives c = ‖b_g‖ − λ. Combining both cases into one closed form:
— the group-wise soft-thresholding update, and the direct vector analogue of the scalar soft(θ, λ) = sign(θ)max(|θ|−λ, 0) derived in 2.2.5. Instead of shrinking a single number toward zero, it shrinks an entire group's vector uniformly toward the origin, all coordinates at once, and clamps the whole thing to the zero vector the moment the group's overall gradient magnitude drops at or below λ.
Where this is used: this exact group-wise soft-thresholding update is what packages like glmnet's grouped-lasso mode and Python's group-lasso solvers iterate to convergence — analogous to how plain coordinate descent iterates the scalar soft- thresholding update for ordinary Lasso.
Blue bars are candidate features competing this round; a bar flashes amber the instant it wins the argmax correlation, then locks violet. Notice locked bars collapse toward zero height in later rounds — that's the orthogonality proof above, visible. This toy dataset's only truly relevant features are x2 and x5.
Six coefficients arranged into three groups of two (imagine three one-hot-encoded region dummies). Top row: plain per-coefficient soft-thresholding, soft(θ,λ) — each bar shrinks and zeros independently. Bottom row: Group Lasso's soft_2(θ_g,λ) from Part 2 of the Derivation — every coefficient in a group is scaled by the same factor, so a group is either fully alive or fully zero, never split.
#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>
// Solve A x = b via Gaussian elimination with partial pivoting -- no library.
std::vector<double> solveLinearSystem(std::vector<std::vector<double>> A, std::vector<double> b) {
int n = static_cast<int>(b.size());
for (int col = 0; col < n; ++col) {
int pivot = col;
for (int r = col + 1; r < n; ++r)
if (std::fabs(A[r][col]) > std::fabs(A[pivot][col])) pivot = r;
std::swap(A[col], A[pivot]);
std::swap(b[col], b[pivot]);
for (int r = col + 1; r < n; ++r) {
double factor = A[r][col] / A[col][col];
for (int c = col; c < n; ++c) A[r][c] -= factor * A[col][c];
b[r] -= factor * b[col];
}
}
std::vector<double> x(n, 0.0);
for (int row = n - 1; row >= 0; --row) {
double sum = b[row];
for (int c = row + 1; c < n; ++c) sum -= A[row][c] * x[c];
x[row] = sum / A[row][row];
}
return x;
}
// OLS refit restricted to the currently selected columns: (X_S^T X_S) theta_S = X_S^T y.
std::vector<double> refitSelected(
const std::vector<std::vector<double>>& X,
const std::vector<double>& y,
const std::vector<int>& selected
) {
int n = static_cast<int>(X.size());
int k = static_cast<int>(selected.size());
std::vector<std::vector<double>> XtX(k, std::vector<double>(k, 0.0));
std::vector<double> Xty(k, 0.0);
for (int a = 0; a < k; ++a) {
for (int i = 0; i < n; ++i) Xty[a] += X[i][selected[a]] * y[i];
for (int c = 0; c < k; ++c)
for (int i = 0; i < n; ++i) XtX[a][c] += X[i][selected[a]] * X[i][selected[c]];
}
return solveLinearSystem(XtX, Xty);
}
int main() {
// Small synthetic dataset -- only features 1 and 3 (0-indexed) are truly relevant.
std::vector<std::vector<double>> X = {
{1.2, -0.4, 0.8, 0.3, -1.1},
{-0.7, 1.6, -0.2, 0.9, 0.4},
{0.5, 2.1, 0.1, -0.6, 1.3},
{1.8, -1.0, 0.6, 0.2, -0.5},
{-1.3, 0.9, -0.9, 1.4, 2.0},
{0.2, 1.4, 0.4, -1.2, -1.6},
{0.9, -0.6, 1.1, 0.7, -0.3},
{-0.4, 1.9, 0.3, -0.8, 0.6}
};
std::vector<double> trueTheta = {0.0, 3.0, 0.0, -2.0, 0.0};
int n = static_cast<int>(X.size());
int p = static_cast<int>(trueTheta.size());
std::vector<double> y(n, 0.0);
for (int i = 0; i < n; ++i)
for (int j = 0; j < p; ++j) y[i] += X[i][j] * trueTheta[j];
std::vector<double> residual = y;
std::vector<int> selected;
const int steps = 2;
for (int s = 0; s < steps; ++s) {
int bestJ = -1;
double bestScore = -1.0;
for (int j = 0; j < p; ++j) {
if (std::find(selected.begin(), selected.end(), j) != selected.end()) continue;
double corr = 0.0;
for (int i = 0; i < n; ++i) corr += X[i][j] * residual[i];
if (std::fabs(corr) > bestScore) { bestScore = std::fabs(corr); bestJ = j; }
}
selected.push_back(bestJ);
std::vector<double> thetaSel = refitSelected(X, y, selected);
residual = y;
for (size_t a = 0; a < selected.size(); ++a)
for (int i = 0; i < n; ++i) residual[i] -= X[i][selected[a]] * thetaSel[a];
std::printf("Step %d: selected feature x%d (|corr|=%.3f)\n", s + 1, bestJ, bestScore);
}
std::printf("Final selected features (in order): ");
for (int j : selected) std::printf("x%d ", j);
std::printf("\n");
return 0;
}- Compressed sensing and signal reconstruction — OMP's classical, foundational application. Reconstruct a signal known in advance to be sparse in some basis (e.g. a handful of nonzero frequency components) from far fewer linear measurements than the signal's raw dimension, well below what plain linear algebra alone would need.
- Sparse dictionary learning uses OMP (or close variants) as the inner "sparse coding" step: representing each data point as a sparse combination of a learned set of basis atoms. Covered in full in a later Dimensionality Reduction module lesson; OMP as introduced here is exactly the subroutine it relies on.
- Group Lasso for categorical variables — a categorical feature one-hot encoded into several dummy columns (e.g. five region indicators) should generally be included or excluded from the model as a unit. Grouping those dummies under one Group Lasso penalty term enforces exactly that.
- Fused Lasso for signals and time series — when neighboring coefficients represent adjacent time points or spatial positions and are expected to change slowly (piecewise-constant structure), Fused Lasso's penalty on consecutive differences encourages runs of equal coefficients rather than merely small ones.
- Multi-task Lasso for jointly predicting related outcomes — forecasting several correlated outputs at once (say, demand for several substitute products) with a shared design matrix, where all the tasks are expected to depend on roughly the same handful of input features even if the exact coefficient values differ per task.
- Genomics pathway analysis — grouping genetic markers by the biological pathway or gene they belong to and applying Group Lasso lets the model decide "this whole pathway matters" rather than sparsely picking individual markers with no biological grouping logic behind the selection.
- OMP's greedy nature is structural, not incidental: once a feature is locked into the active set, it is never reconsidered or removed, so an early suboptimal pick can never be undone later. Lasso, by contrast, solves one globally-optimized objective and can end up with a different, better support entirely. Greedy and globally-optimal sparse solutions are not guaranteed to agree.
- Choosing groups for Group Lasso that don't reflect genuine domain structure defeats the purpose. Grouping coefficients arbitrarily (say, just by their order in a spreadsheet) gives you the computational cost of structured sparsity with none of its benefit — the groups need to correspond to something meaningful: all levels of one categorical feature, all lags of one time series, all markers in one gene.
- Assuming OMP is guaranteed to recover the true sparse support in general. It provably does under specific conditions on how correlated the candidate features are with each other — in spirit similar to Lasso's "irrepresentable condition" from 2.3.3 — but with strongly correlated features, OMP can lock onto the wrong one early and never recover, since it has no mechanism to backtrack.
Going deeper
OMP is a genuine mechanical cousin of LARS (2.3.5), and it's worth being precise about exactly how they relate. Both build up an active set of features one at a time, both choose the next feature to add based on correlation with the current residual — but there the mechanical similarity ends. LARS, once a feature joins the active set, moves all active coefficients continuously in the "equiangular direction" — the one direction that keeps every active feature's correlation with the residual tied and decreasing together — stopping the instant some inactive feature's correlation catches up to join the active set. OMP skips all of that continuous movement: the moment a feature is selected, it jumps straight to the full OLS solution on the current active set, in one discrete step. LARS traces out a smooth, piecewise-linear path of solutions as its regularization strength varies (which is what makes it able to reproduce the entire Lasso solution path); OMP produces one discrete solution per fixed number of selected features, with no continuous path connecting them.
The theoretical conditions under which OMP is guaranteed to recover the exact sparse support are studied formally in compressed sensing under names like the restricted isometry property (RIP) and mutual-coherence bounds — precise, checkable statements about how close to orthogonal the candidate features need to be for the orthogonality argument in the Derivation above to actually translate into a correctness guarantee, rather than just a "never wastes a step re-picking the same feature" property.
After OMP selects a feature and refits OLS on the active set, why can that exact same feature never again be the argmax of |x_j^T r| in a later step (short of the residual already being zero)?
Because refitting OLS on the active set S makes the normal equations hold exactly: X_S^T r = 0. That means x_j^T r = 0 for every feature j already in S -- its correlation with the residual is pinned at exactly zero right after the refit. Later steps only ever grow the active set (features are never removed), so once a feature's column is orthogonal to the residual, it stays a fixed zero-score candidate unless the entire residual collapses to zero first. A feature with a guaranteed zero score can never be the argmax unless every other candidate also scores zero -- i.e., unless there's nothing left to explain.
Lasso gets sparsity from a globally-optimized penalty; OMP builds it explicitly, one greedy step at a time, backed by a hard guarantee — an OLS refit's normal equations make the residual exactly orthogonal to everything already selected, so the algorithm structurally cannot waste a step. Group Lasso, Fused Lasso, and multi-task Lasso all extend that same soft-thresholding machinery from individual coefficients to whatever grouping actually reflects the problem's real structure — categorical levels, neighboring time points, or shared support across related tasks. Together with Lasso and LARS, this closes out the Regression module's arc on sparse and structured regularization; OMP's role as the sparse- coding subroutine inside dictionary learning is picked back up in the Dimensionality Reduction module ahead.