Dictionary Learning & Sparse Coding

Chapter 1: The Sparse Coding Problem

Formally, given a signal x ∈ ℝ^m and a dictionary D ∈ ℝ^m×n (with n > m, overcomplete), find the sparsest coefficient vector α such that x = Dα:

The L0 "norm" ||α||₀ counts the number of non-zero entries. We want the representation with fewest active atoms.

Why Overcomplete?

If D has exactly m columns (a square basis like the DFT or wavelets), there's only one representation. With n > m columns, there are infinitely many solutions to x = Dα. The sparsity constraint selects the "simplest" one.

The Optimization Landscape

The L0 minimization is NP-hard in general. We need algorithms that find (approximately) sparse solutions efficiently. The next chapters cover the two main approaches: greedy pursuit (Matching Pursuit, OMP) and convex relaxation (LASSO/Basis Pursuit).

Chapter 2: Matching Pursuit — Greedy Atom Selection

Matching Pursuit (MP) is the simplest sparse coding algorithm. It's greedy: at each step, find the atom most correlated with the current residual, subtract its contribution, repeat.

python
def matching_pursuit(x, D, n_atoms=5, tol=1e-4):
    """Greedy atom selection. D: [m, n], x: [m]."""
    r = x.copy()           # residual
    alpha = np.zeros(D.shape[1])
    selected = []

    for _ in range(n_atoms):
        # Find most correlated atom
        corrs = D.T @ r     # [n] inner products
        k = np.argmax(np.abs(corrs))
        selected.append(k)

        # Update
        alpha[k] += corrs[k]
        r = r - corrs[k] * D[:, k]

        if np.linalg.norm(r) < tol:
            break

    return alpha, selected

Chapter 3: OMP — Orthogonal Matching Pursuit

Orthogonal Matching Pursuit (OMP) improves on MP with one key change: after selecting each new atom, it re-solves the least-squares problem over ALL selected atoms. This ensures the residual is orthogonal to the span of selected atoms — hence "orthogonal."

The Key Difference

This is just the least-squares projection of x onto the columns of D indexed by S. The residual is r = x - D_Sα_S, which is guaranteed orthogonal to span(D_S).

python
def omp(x, D, n_atoms=5, tol=1e-4):
    """OMP: re-solve LS after each atom selection."""
    r = x.copy()
    support = []
    alpha = np.zeros(D.shape[1])

    for _ in range(n_atoms):
        # Select: same as MP
        corrs = D.T @ r
        k = np.argmax(np.abs(corrs))
        support.append(k)

        # KEY DIFFERENCE: re-solve LS over entire support
        Ds = D[:, support]
        alpha_s = np.linalg.lstsq(Ds, x, rcond=None)[0]
        for i, idx in enumerate(support):
            alpha[idx] = alpha_s[i]

        # Residual orthogonal to span(Ds)
        r = x - Ds @ alpha_s

        if np.linalg.norm(r) < tol:
            break

    return alpha, support

Chapter 4: LASSO — L1 Convex Relaxation

The L0 problem is combinatorial. Greedy methods (MP, OMP) are fast but may not find the globally sparsest solution. Is there a convex relaxation that's both solvable and provably correct?

Yes: replace ||α||₀ with ||α||₁ = ∑|α_i|. This is the LASSO (Least Absolute Shrinkage and Selection Operator):

The L1 norm is the tightest convex relaxation of L0. It promotes sparsity because the L1 ball has corners on the axes — the optimal solution tends to land at a corner where some coordinates are exactly zero.

Why L1 Gives Sparsity

Geometrically: the L1 ball is a diamond (in 2D) or cross-polytope. The level sets of the quadratic loss are ellipsoids. When an ellipsoid first touches the diamond, it typically hits a corner (a sparse point). The L2 ball is a sphere — it touches ellipsoids at arbitrary points, giving non-sparse solutions.

Solving LASSO: ISTA

The Iterative Shrinkage-Thresholding Algorithm (ISTA) alternates gradient step + soft thresholding:

where L = largest eigenvalue of D^TD, and Shrink_τ(z) = sign(z)max(|z|-τ, 0).

python
def ista(x, D, lam, n_iter=500):
    """ISTA for LASSO: min 0.5*||x-D@alpha||^2 + lam*||alpha||_1"""
    m, n = D.shape
    L = np.linalg.norm(D.T @ D, ord=2)  # Lipschitz constant
    alpha = np.zeros(n)

    for _ in range(n_iter):
        # Gradient step
        grad = D.T @ (D @ alpha - x)
        z = alpha - grad / L

        # Soft thresholding (proximal of L1)
        alpha = np.sign(z) * np.maximum(np.abs(z) - lam/L, 0)

    return alpha

Chapter 5: K-SVD — Learning the Dictionary

So far we assumed the dictionary D is given. But the real power comes from learning D from data. The K-SVD algorithm (Aharon, Elad, Bruckstein, 2006) alternates between sparse coding and dictionary update:

where X = [x₁ | ... | x_N] are training signals and Α = [α₁ | ... | α_N] are their sparse codes.

The Two-Stage Algorithm

Atom Update (the "SVD" in K-SVD)

To update atom k, collect all signals that use atom k (call this set ω_k). The error matrix restricted to these signals, with all other atoms' contributions removed, is:

We want d_k and the k-th row of Α to best approximate E_k. This is a rank-1 approximation problem — solved by the largest singular vector of E_k:

python
def ksvd(X, n_atoms, sparsity, n_iter=50):
    """K-SVD dictionary learning."""
    m, N = X.shape
    # Initialize dictionary (random normalized columns)
    D = np.random.randn(m, n_atoms)
    D /= np.linalg.norm(D, axis=0)

    for _ in range(n_iter):
        # Stage 1: Sparse code each signal with OMP
        Alpha = np.zeros((n_atoms, N))
        for i in range(N):
            Alpha[:, i], _ = omp(X[:, i], D, n_atoms=sparsity)

        # Stage 2: Update each atom
        for k in range(n_atoms):
            # Signals using atom k
            omega = np.where(Alpha[k, :] != 0)[0]
            if len(omega) == 0:
                continue

            # Error without atom k
            E = X[:, omega] - D @ Alpha[:, omega] + \
                np.outer(D[:, k], Alpha[k, omega])

            # Rank-1 SVD
            U, s, Vt = np.linalg.svd(E, full_matrices=False)
            D[:, k] = U[:, 0]
            Alpha[k, omega] = s[0] * Vt[0, :]

    return D, Alpha

Chapter 6: Showcase — Interactive K-SVD

Watch K-SVD learn a dictionary from synthetic signal patches. You control the sparsity level, number of atoms, and noise. The dictionary starts random and converges to meaningful structure.

Things to Try

• Increase noise to 0.4: The atoms become noisier but still emerge. K-SVD is robust to moderate noise.
• Set sparsity to 1: Each signal uses only one atom. K-SVD converges fastest here (essentially k-means on signal directions).
• Set dictionary larger than needed: Extra atoms become "dead" (rarely used). K-SVD naturally handles overparameterization.

Chapter 7: Beyond — Connections

Dictionary learning sits at the crossroads of signal processing, machine learning, and neuroscience.

The Neuroscience Connection

Olshausen & Field (1996) showed that sparse coding on natural image patches produces atoms that look like oriented Gabor filters — exactly what V1 neurons in the visual cortex compute. The brain may literally be doing dictionary learning. Sparsity isn't just a computational convenience; it may be a fundamental principle of biological representation.

Related Lessons

Previous: NMF & Clustering — non-negative factorization as an alternative to signed sparse coding.

Previous: Autoencoders & Robust PCA — linear AE = PCA, the dense counterpart to dictionary learning.