Theory: How and Why Compressed Sensing Works¶
1. The Problem: Underdetermined Systems¶
We are trying to solve a linear system $y = Ax$, where:
- $A$ is an $m \times n$ measurement matrix.
- $y$ is the vector of $m$ measurements (observations).
- $x$ is the vector of $n$ unknowns (the signal we want to recover).
In our example, $m = 50$ and $n = 500$. Because $m \ll n$, the system is heavily underdetermined. Geometrically, the equation $Ax = y$ defines a high-dimensional hyperplane, meaning there are infinitely many solutions for $x$ that perfectly satisfy the measurements.
2. The Sparsity Prior¶
To find the correct solution among infinitely many, we need a prior assumption. In Compressed Sensing, we assume the true signal $x$ is $k$-sparse. This means that out of $n = 500$ unknowns, only $k = 5$ are non-zero. We don't know which 5 indices are non-zero, but we know the count is very small.
3. Why $L_2$ Minimization (Least Squares) Fails¶
The standard approach to underdetermined systems is to find the solution with the minimum $L_2$ norm: $$ \min \|x\|_2 \quad \text{subject to} \quad Ax = y $$ Geometrically, the $L_2$ "unit ball" is a smooth, round sphere. When the affine subspace $Ax = y$ expands and touches this sphere, the point of contact will almost never lie exactly on a coordinate axis. As a result, the $L_2$ solution spreads the "energy" evenly across all 500 variables, resulting in a dense, noisy-looking vector that completely fails to recover the original 5 spikes.
4. Why $L_1$ Minimization (Compressed Sensing) Succeeds¶
Compressed sensing replaces the $L_2$ norm with the $L_1$ norm: $$ \min \|x\|_1 \quad \text{subject to} \quad Ax = y $$ Geometrically, the $L_1$ "unit ball" in high dimensions is a cross-polytope (like a diamond in 2D or an octahedron in 3D). Crucially, this shape has sharp corners (spikes) that lie exactly on the coordinate axes.
When the affine subspace $Ax = y$ expands and intersects this $L_1$ ball, it is highly probable that the first point of contact will be one of these sharp corners. A corner corresponds to a solution where most coordinates are exactly zero. This geometric property is what allows $L_1$ minimization to naturally promote sparsity and recover the true signal.
5. The Mathematical Guarantee: Restricted Isometry Property (RIP)¶
For $L_1$ minimization to work reliably, the measurement matrix $A$ must not "hide" sparse signals. This is formalized by the Restricted Isometry Property (RIP). A matrix satisfies the RIP if it acts almost like an orthogonal matrix when restricted to sparse vectors (i.e., it preserves the lengths of sparse vectors).
Random matrices with independent Gaussian entries (like np.random.randn in our code) satisfy the RIP with overwhelming probability. Theoretical guarantees state that exact recovery is possible provided the number of measurements $m$ satisfies:
$$ m \ge \mathcal{O}\left(k \log\left(\frac{n}{k}\right)\right) $$
For $k=5$ and $n=500$, the required $m$ is roughly $5 \log(100) \approx 23$. Our choice of $m=50$ is more than enough to guarantee perfect recovery!
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
# Set random seed for reproducibility
np.random.seed(17)
# Problem dimensions
n = 500 # Number of unknowns (signal length)
m = 50 # Number of equations (measurements)
k = 5 # Number of non-zero elements (sparsity)
print(f"System dimensions: {m} equations, {n} unknowns, {k} non-zeros")
System dimensions: 50 equations, 500 unknowns, 5 non-zeros
# 1. Generate a random Gaussian measurement matrix A
# Gaussian matrices satisfy the Restricted Isometry Property (RIP) with high probability
A = np.random.randn(m, n)
# 2. Generate the true k-sparse signal x_true
x_true = np.zeros(n)
# Randomly choose k indices to be non-zero
nonzero_indices = np.random.choice(n, k, replace=False)
# Assign random values to these indices
x_true[nonzero_indices] = np.random.randn(k) * 10
# 3. Generate the measurements y
y = A @ x_true
print(f"True non-zero indices: {nonzero_indices}")
print(f"True non-zero values: {x_true[nonzero_indices]}")
True non-zero indices: [336 178 43 135 316] True non-zero values: [-5.67326905 -2.72795721 8.78510681 -8.09743749 -6.63558559]
# To solve min ||x||_1 subject to Ax = y using linear programming,
# we split x into positive and negative parts: x = u - v, where u >= 0, v >= 0.
# Then ||x||_1 = sum(u) + sum(v).
# Objective function: minimize sum(u) + sum(v)
c = np.ones(2 * n)
# Equality constraint: A(u - v) = y => [A, -A] @ [u; v] = y
A_eq = np.hstack([A, -A])
b_eq = y
# Bounds: u >= 0, v >= 0
bounds = [(0, None) for _ in range(2 * n)]
# Solve the linear program using HiGHS (modern, robust SciPy solver)
res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs')
if res.success:
# Reconstruct x from u and v
u_v = res.x
x_l1 = u_v[:n] - u_v[n:]
print("L1 Minimization (Compressed Sensing) succeeded!")
else:
print("L1 Minimization failed:", res.message)
x_l1 = np.zeros(n)
L1 Minimization (Compressed Sensing) succeeded!
# L2 minimization (Minimum energy solution) using Moore-Penrose pseudo-inverse
x_l2 = np.linalg.pinv(A) @ y
print("L2 Minimization (Least Squares) completed.")
L2 Minimization (Least Squares) completed.
pwd
'/data/shared/examples/Compressed_Sensing'
# Calculate reconstruction errors
error_l1 = np.linalg.norm(x_true - x_l1, 2)
error_l2 = np.linalg.norm(x_true - x_l2, 2)
print(f"\nL2 Reconstruction Error: {error_l2:.4e}")
print(f"L1 Reconstruction Error: {error_l1:.4e}")
# Plotting
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
# Plot 1: True Signal
axes[0].stem(x_true, markerfmt='C0o', linefmt='C0-', basefmt=" ")
axes[0].set_title(f"True Signal (k={k} non-zeros)")
axes[0].set_xlabel("Index")
axes[0].set_ylabel("Amplitude")
axes[0].set_xlim(0, n)
# Plot 2: L1 Reconstruction (Compressed Sensing)
axes[1].stem(x_l1, markerfmt='C1o', linefmt='C1-', basefmt=" ")
axes[1].set_title(f"L1 Minimization (Error: {error_l1:.2e})")
axes[1].set_xlabel("Index")
axes[1].set_ylabel("Amplitude")
axes[1].set_xlim(0, n)
# Plot 3: L2 Reconstruction (Least Squares)
axes[2].stem(x_l2, markerfmt='C2o', linefmt='C2-', basefmt=" ")
axes[2].set_title(f"L2 Minimization (Error: {error_l2:.2e})")
axes[2].set_xlabel("Index")
axes[2].set_ylabel("Amplitude")
axes[2].set_xlim(0, n)
plt.tight_layout()
plt.show()
L2 Reconstruction Error: 1.4367e+01 L1 Reconstruction Error: 1.0436e-13
pwd
'/data/shared/examples/Compressed_Sensing'
%load_ext watermark
%watermark --python --date --iversions --machine
Python implementation: CPython Python version : 3.12.2 IPython version : 8.25.0 Compiler : GCC 12.3.0 OS : Linux Release : 6.17.0-35-generic Machine : x86_64 Processor : x86_64 CPU cores : 32 Architecture: 64bit numpy : 1.26.4 matplotlib: 3.8.4 sys : 3.12.2 | packaged by conda-forge | (main, Feb 16 2024, 20:50:58) [GCC 12.3.0] json : 2.0.9
!jupyter nbconvert --to HTML compressed_sensing.ipynb
[NbConvertApp] Converting notebook compressed_sensing.ipynb to HTML [NbConvertApp] WARNING | Alternative text is missing on 1 image(s). [NbConvertApp] Writing 399114 bytes to compressed_sensing.html