8.1 Sparsity and compressibility

Sparsity and Compressibility

Sparsity and compressibility describe how signals can be represented using only a few significant components. These properties are the foundation of compressive sensing: if a signal is sparse (or close to it), you can acquire, store, and reconstruct it far more efficiently than classical sampling theory would suggest. This topic covers the formal definitions, the representations that exploit sparsity, and the algorithms and applications that follow from it.

Sparsity in Signal Processing

A signal is sparse when most of its coefficients are zero in some representation. This single property unlocks efficient compression, denoising, and reconstruction from incomplete measurements.

Sparse Signals

A signal $x \in \mathbb{R}^N$ is called $k$-sparse if it has at most $k$ non-zero coefficients in a particular basis or dictionary, where $k \ll N$. Because only $k$ values matter, you can represent the entire signal compactly using just those few basis functions.

• An audio signal with only a handful of active frequency components is sparse in the Fourier basis.
• A piecewise-constant image has sparse gradients (few non-zero entries in a finite-difference representation).
• Sparsity is what makes compression possible: store the $k$ significant coefficients and discard the rest.

Compressible Signals

Most real-world signals aren't exactly sparse, but they're compressible: their sorted coefficient magnitudes decay rapidly. Formally, a signal is compressible if its sorted coefficients $|c_{(1)}| \geq |c_{(2)}| \geq \cdots$ obey a power-law decay:

$|c_{(i)}| \leq C \cdot i^{-q}$

for some constants $C > 0$ and $q > 1$. The faster the decay (larger $q$), the better the signal can be approximated by keeping only its largest coefficients.

Natural images and audio recordings are classic examples. They aren't strictly sparse, but truncating to the top $k$ coefficients gives an excellent approximation.

Sparsity vs. Compressibility

The distinction matters because theoretical guarantees for sparse recovery (exact reconstruction) are stronger than those for compressible signals (bounded approximation error). Both can be exploited, but you should know which regime your signal falls into.

Sparse Representations

A sparse representation expresses a signal as a linear combination of a small number of elements drawn from a basis or dictionary:

$x = \sum_{i=1}^{k} \alpha_i \phi_i$

where $\alpha_i$ are the non-zero coefficients and $\phi_i$ are the corresponding basis functions (or atoms). The goal is to capture the signal's essential information with as few terms as possible.

Basis Functions

Basis functions are the building blocks of a representation. The choice of basis determines whether a given signal looks sparse or not.

• Fourier basis (complex sinusoids): ideal for signals with a few dominant frequencies.
• Wavelet basis (localized, multi-scale functions): well-suited for signals with transient features or piecewise-smooth structure.
• Gabor functions: provide joint time-frequency localization, useful for analyzing non-stationary signals.

A signal that looks dense in one basis may be sparse in another. Choosing the right basis is critical.

Overcomplete Dictionaries

An overcomplete dictionary $D \in \mathbb{R}^{N \times M}$ with $M > N$ contains more atoms than the signal dimension. The atoms are not linearly independent, which means the representation of any signal is no longer unique.

This redundancy is actually useful: it gives you more flexibility to find a sparse representation. Examples include concatenations of multiple bases (e.g., Fourier + wavelet), random dictionaries, and dictionaries learned directly from training data (e.g., via K-SVD).

The tradeoff is that finding the sparsest representation in an overcomplete dictionary is computationally harder than in an orthonormal basis.

Sparse Coding

Sparse coding is the problem of finding the sparsest coefficient vector $\alpha$ such that $x \approx D\alpha$. The ideal formulation minimizes the $\ell_0$ "norm" (the count of non-zero entries):

$\min \|\alpha\|_0 \quad \text{subject to} \quad \|x - D\alpha\|_2 \leq \epsilon$

This problem is NP-hard in general, so practical approaches use tractable relaxations or greedy methods:

• Basis pursuit relaxes $\ell_0$ to $\ell_1$ (convex optimization).
• Orthogonal matching pursuit greedily selects atoms one at a time.
• Iterative thresholding alternates gradient steps with coefficient truncation.

These are covered in more detail in the recovery algorithms section below.

Compressed Sensing

Compressed sensing (CS) is a framework for acquiring and reconstructing sparse or compressible signals from far fewer measurements than the Nyquist rate requires. If a signal is $k$-sparse in some basis, CS theory shows you can recover it from roughly $\mathcal{O}(k \log(N/k))$ random linear measurements.

Sampling Sparse Signals

In CS, you don't sample the signal directly. Instead, you acquire $m$ linear measurements:

$y = \Phi x$

where $\Phi \in \mathbb{R}^{m \times N}$ is the measurement matrix and $m \ll N$. The system is underdetermined, but sparsity makes the solution unique (under the right conditions on $\Phi$).

The number of measurements $m$ needed scales with the sparsity level $k$, not the ambient dimension $N$. This is what makes sub-Nyquist acquisition possible.

Incoherent Sampling

For CS to work, the measurement matrix $\Phi$ must be incoherent with the sparsifying basis $\Psi$. Incoherence means the rows of $\Phi$ cannot sparsely represent the columns of $\Psi$ (and vice versa). In practical terms, each measurement captures a little information about every sparse coefficient rather than concentrating on a few.

• Random Gaussian and Bernoulli matrices are highly incoherent with virtually any fixed basis, which is why random measurements are so popular in CS.
• Partial Fourier matrices (randomly selected rows of the DFT matrix) are incoherent with the canonical (time-domain) basis and are useful in applications like MRI.

Restricted Isometry Property

The restricted isometry property (RIP) is the key theoretical condition guaranteeing stable recovery. A matrix $\Phi$ satisfies the RIP of order $k$ with constant $\delta_k \in (0,1)$ if, for all $k$-sparse vectors $x$:

$(1 - \delta_k)\|x\|_2^2 \leq \|\Phi x\|_2^2 \leq (1 + \delta_k)\|x\|_2^2$

This says $\Phi$ approximately preserves the $\ell_2$ norm of every $k$-sparse signal. If $\delta_{2k}$ is small enough, then $\ell_1$ minimization (basis pursuit) recovers every $k$-sparse signal exactly and every compressible signal with bounded error.

Random Gaussian matrices satisfy the RIP with high probability when $m = \mathcal{O}(k \log(N/k))$. Verifying the RIP for a specific matrix is itself NP-hard, so in practice you rely on probabilistic guarantees.

Sparse Recovery Algorithms

These algorithms reconstruct a sparse signal from compressed or incomplete measurements. They differ in their approach (convex optimization vs. greedy vs. iterative), their computational cost, and the conditions under which they succeed.

Basis Pursuit

Basis pursuit replaces the intractable $\ell_0$ problem with a convex $\ell_1$ minimization:

$\min \|\alpha\|_1 \quad \text{subject to} \quad y = \Phi \alpha$

When noise is present, the constraint is relaxed to $\|y - \Phi \alpha\|_2 \leq \epsilon$ (this variant is called basis pursuit denoising or BPDN, also known as LASSO in statistics).

Steps for solving via standard methods:

1. Formulate the problem as a linear program (LP) or second-order cone program (SOCP).
2. Solve using an LP solver, or use specialized algorithms like ADMM or proximal gradient methods.
3. The solution $\hat{\alpha}$ is the recovered sparse coefficient vector.

Basis pursuit provides strong theoretical guarantees (exact recovery under RIP) and is stable to noise, but it can be slower than greedy methods for very large problems.

Orthogonal Matching Pursuit

Orthogonal matching pursuit (OMP) is a greedy algorithm that builds up the sparse representation one atom at a time:

1. Initialize the residual $r_0 = y$ and the support set $S_0 = \emptyset$.

2. Find the dictionary atom most correlated with the current residual: $j^* = \arg\max_j |\langle \phi_j, r_t \rangle|$.

3. Add $j^*$ to the support set: $S_{t+1} = S_t \cup \{j^*\}$.

4. Update the coefficient estimate by solving the least-squares problem over the current support: $\hat{\alpha}_{S_{t+1}} = \arg\min \|y - \Phi_{S_{t+1}} \alpha\|_2$.

5. Update the residual: $r_{t+1} = y - \Phi_{S_{t+1}} \hat{\alpha}_{S_{t+1}}$.

6. Repeat steps 2-5 until $k$ atoms are selected or the residual is below a tolerance.

OMP is fast and simple to implement. It works well when the signal is highly sparse and the dictionary atoms are incoherent, but it can fail if atoms are too correlated.

Iterative Thresholding

Iterative thresholding algorithms alternate between a gradient step and a sparsity-promoting thresholding step:

1. Perform a gradient update: $z_{t+1} = \hat{x}_t + \mu \Phi^T(y - \Phi \hat{x}_t)$, where $\mu$ is a step size.

2. Apply a thresholding operator to enforce sparsity: $\hat{x}_{t+1} = \mathcal{T}_\lambda(z_{t+1})$.

3. Repeat until convergence.

The two main variants differ in the thresholding operator $\mathcal{T}_\lambda$:

• Iterative hard thresholding (IHT): keeps the $k$ largest-magnitude entries, sets the rest to zero.
• Iterative soft thresholding (IST): shrinks all coefficients toward zero by $\lambda$, i.e., $\mathcal{T}_\lambda(z_i) = \text{sign}(z_i)\max(|z_i| - \lambda, 0)$. This is equivalent to proximal gradient descent for the $\ell_1$-penalized problem.

These methods are memory-efficient and easy to implement, making them attractive for large-scale problems.

Applications of Sparsity

Image Compression

Image compression standards exploit the fact that natural images are sparse (or compressible) in transform domains. JPEG uses the DCT; JPEG2000 uses wavelets. In both cases, the encoder transforms the image, quantizes the (mostly near-zero) coefficients, and stores only the significant ones. Compressed sensing extends this idea to the acquisition stage: rather than sampling densely and then compressing, you acquire compressed measurements directly.

Audio Compression

Audio signals are compressible in time-frequency representations. MP3 uses the MDCT, and AAC uses similar transforms. The encoder discards coefficients below a perceptual threshold (exploiting both sparsity and psychoacoustic masking). The result is compression ratios of 10:1 or more with minimal perceptible quality loss.

Biomedical Signal Processing

Sparsity-based methods are widely used in biomedical applications:

• MRI acceleration: CS allows reconstruction of MR images from undersampled k-space data, reducing scan times significantly (often by a factor of 4-8x). This is one of the most successful real-world deployments of compressed sensing.
• EEG/ECG compression: Sparse representations enable efficient telemetry for wearable health monitors with limited bandwidth and battery life.
• Feature extraction: Sparse coding can identify diagnostically relevant patterns in physiological signals.

Sparsity-Based Denoising

Denoising exploits a simple idea: the signal is sparse in some transform domain, but noise is not. By transforming the noisy signal, suppressing small coefficients (which are mostly noise), and inverting the transform, you recover a cleaner version of the signal.

Sparse Signal Denoising

Given a noisy observation $y = x + n$, where $x$ is sparse in basis $\Psi$:

1. Transform: compute coefficients $c = \Psi^T y$.
2. Threshold: apply a thresholding operator to suppress noise-dominated coefficients.
3. Reconstruct: $\hat{x} = \Psi \hat{c}$.

This approach (wavelet shrinkage, introduced by Donoho and Johnstone) is simple and effective. Dictionary learning-based methods (e.g., K-SVD denoising) adapt the basis to the data for better performance.

Sparse Image Denoising

Image denoising methods like BM3D and K-SVD denoising go beyond global transforms by exploiting local sparsity and self-similarity:

• BM3D groups similar patches, applies a collaborative 3D transform, thresholds the coefficients, and averages the results. It remains a benchmark method.
• K-SVD learns an overcomplete dictionary from the noisy image patches themselves, then uses sparse coding to denoise each patch.

Both methods preserve fine details and textures better than simple global thresholding.

Thresholding Techniques

Thresholding is the core operation in transform-domain denoising. The two standard choices:

• Hard thresholding: $\hat{c}_i = c_i$ if $|c_i| \geq \lambda$, else $\hat{c}_i = 0$. Preserves coefficient magnitudes but can introduce artifacts at the threshold boundary.
• Soft thresholding: $\hat{c}_i = \text{sign}(c_i)\max(|c_i| - \lambda, 0)$. Shrinks all coefficients toward zero, producing smoother results but introducing a systematic bias.

The threshold $\lambda$ is typically set proportional to the noise standard deviation $\sigma$. A common choice is $\lambda = \sigma\sqrt{2\log N}$ (the universal threshold), which provides near-optimal risk for soft thresholding in the wavelet domain.

Sparsity in Machine Learning

Sparsity is a powerful tool for controlling model complexity. In high-dimensional settings where the number of features exceeds the number of samples, sparsity-inducing methods prevent overfitting and produce interpretable models.

Sparse Feature Selection

LASSO (Least Absolute Shrinkage and Selection Operator) adds an $\ell_1$ penalty to the regression objective:

$\min_\beta \|y - X\beta\|_2^2 + \lambda\|\beta\|_1$

The $\ell_1$ penalty drives many coefficients to exactly zero, effectively selecting a subset of features. Elastic Net combines $\ell_1$ and $\ell_2$ penalties, which helps when features are correlated (LASSO tends to arbitrarily pick one from a group of correlated features; Elastic Net selects the group).

Sparse Regularization

Beyond feature selection, $\ell_1$ regularization appears throughout machine learning:

• Group sparsity (group LASSO) sets entire groups of coefficients to zero, useful when features have a natural grouping structure.
• Sparse + low-rank decompositions separate a data matrix into sparse and low-rank components (e.g., robust PCA for separating foreground from background in video).
• The regularization parameter $\lambda$ controls the sparsity level: larger $\lambda$ means fewer non-zero coefficients. Cross-validation is the standard way to select it.

Sparse Deep Learning

Modern neural networks are heavily overparameterized, and sparsity techniques can reduce their size without significant accuracy loss:

• Magnitude pruning: remove weights with the smallest absolute values, then fine-tune the remaining network. Pruning rates of 90%+ are common for large models.
• Structured pruning: remove entire filters or neurons rather than individual weights, yielding actual speedups on standard hardware.
• Weight quantization: reduce the precision of weights (e.g., from 32-bit to 4-bit), which is a form of approximate sparsity in the representation of each weight.

These techniques are essential for deploying large models on edge devices with limited memory and compute.

Challenges and Limitations

Computational Complexity

Sparse recovery algorithms range from cheap (OMP scales roughly as $\mathcal{O}(kmN)$) to expensive (basis pursuit via general-purpose LP solvers can scale as $\mathcal{O}(N^3)$ or worse). For large-scale problems, first-order methods like ISTA/FISTA and ADMM are preferred because they scale more gracefully, but convergence can be slow for ill-conditioned problems.

Robustness to Noise

All sparse recovery methods assume the signal is (approximately) sparse. When this assumption is violated, or when measurement noise is high, recovery quality degrades. Basis pursuit denoising and its variants provide bounded error guarantees under the RIP, but the bounds can be loose in practice. Developing algorithms that are robust to model mismatch (the signal isn't as sparse as assumed) remains an active research area.

Selection of Sparsifying Basis

The performance of any sparsity-based method depends heavily on choosing a basis or dictionary in which the signal is actually sparse. For some signal classes (e.g., piecewise-smooth signals in a wavelet basis), good choices are well-established. For others, the right basis isn't obvious, and you need data-driven approaches:

• Dictionary learning (e.g., K-SVD) adapts the dictionary to training data.
• Deep unrolling methods learn both the dictionary and the recovery algorithm end-to-end.

The computational cost of dictionary learning and the risk of overfitting to training data are practical concerns, especially when the signal statistics are non-stationary.