Compressed Sensing

Compressed sensing recovers a high-dimensional signal from far fewer linear measurements than the Nyquist–Shannon rate when the signal is sparse (or sparse in a known basis). SMILE’s smile.cs package provides sensing matrices and three recovery algorithms:

With a k-sparse x ∈ ℝⁿ and measurements y = A x where A ∈ ℝm×n and m ≪ n, recovery is typically reliable when m = O(k log(n/k)) (a practical rule of thumb is m ≥ 4k–5k).

Quick Start

Build a sparse ground truth, take Gaussian measurements, and recover with OMP and CoSaMP. Relative L2 error is machine-epsilon on this noiseless instance.

import java.util.Arrays
import smile.cs.*
import smile.math.MathEx

MathEx.setSeed(19650218)
int n = 128, m = 50, k = 5

double[] xTrue = new double[n]
int[] locs = MathEx.permutate(n)
for (int i = 0; i < k; i++) {
    double amp = 1.0 + MathEx.random()
    xTrue[locs[i]] = MathEx.random() < 0.5 ? amp : -amp
}

var mm = MeasurementMatrix.gaussian(m, n)
double[] y = mm.measure(xTrue)

var omp = OMP.fit(mm.phi(), y, k)
Arrays.toString(omp.support())

double[] d = omp.x().clone()
for (int i = 0; i < n; i++) d[i] -= xTrue[i]
MathEx.norm(d) / MathEx.norm(xTrue)

var cosamp = CoSaMP.fit(mm.phi(), y, k)
d = cosamp.x().clone()
for (int i = 0; i < n; i++) d[i] -= xTrue[i]
MathEx.norm(d) / MathEx.norm(xTrue)
cosamp.iter()

[56, 110, 20, 34, 114]
2.59e-16
2.59e-16
1

Measurement Matrix

MeasurementMatrix wraps the sensing operator Φ (and an optional sparsifying wavelet Ψ). Factory methods scale entries so that E[‖Φ x‖₂²] ≈ ‖x‖₂². Use measure / backProject for forward and adjoint maps, and toMatrix() for an implicit Matrix that supports both mv and tv.

import smile.cs.*
import smile.math.MathEx

MathEx.setSeed(0)
var g = MeasurementMatrix.gaussian(30, 100)
var b = MeasurementMatrix.bernoulli(20, 64)
var p = MeasurementMatrix.partial(20, 64)

g.nrow() + " x " + g.ncol()
b.nrow() + " x " + b.ncol()
p.nrow()

// Adjoint identity: ⟨Φx, y⟩ = ⟨x, Φᵀ y⟩
double[] x = new double[100]
double[] yv = new double[30]
for (int i = 0; i < 100; i++) x[i] = MathEx.random() - 0.5
for (int i = 0; i < 30; i++) yv[i] = MathEx.random() - 0.5
double lhs = MathEx.dot(g.measure(x), yv)
double rhs = MathEx.dot(x, g.backProject(yv))
lhs
rhs
Math.abs(lhs - rhs)

30 x 100
20 x 64
20
-0.3969560959
-0.3969560959
3.33e-16

Orthogonal Matching Pursuit

OMP greedily grows a support of size k: each step picks the column of A most correlated with the residual, then re-solves least squares on the active set (incremental QR). With a Gaussian matrix and m ≥ C k log n, exact recovery holds when the RIP constant satisfies δ2k < √2 − 1.

Use OMP.fit(A, y, k) or new OMP.Options(k, tol) for early stopping on ‖r‖₂ / ‖y‖₂. Prefer OMP when k is known and speed matters.

CoSaMP

CoSaMP maintains a 2k-candidate support, prunes to k, and re-solves least squares—yielding near-optimal guarantees and better noise tolerance than plain OMP. It is robust to moderate over-estimates of k (requires 2k < m). Call CoSaMP.fit(A, y, k) or pass CoSaMP.Options(k, maxIter, tol).

Basis Pursuit / BPDN

Basis Pursuit solves the convex program minimize ‖x‖₁ subject to A x = y (or ‖A x − y‖₂ ≤ ε for BPDN). SMILE uses a log-barrier interior-point method with CG Newton steps. Unlike the greedy methods, k need not be known. Prefer BP/BPDN for best quality under noise on small-to-medium problems.

Exact BP is more reliable when m is not too small relative to n. The example below uses a compact instance and tight options; the residual ‖A x − y‖₂ / ‖y‖₂ shows feasibility.

import smile.cs.*
import smile.math.MathEx
import smile.tensor.Vector

MathEx.setSeed(7777)
int n = 32, m = 20, k = 3
double[] xTrue = new double[n]
int[] locs = MathEx.permutate(n)
for (int i = 0; i < k; i++)
    xTrue[locs[i]] = (i + 1) * (MathEx.random() < 0.5 ? 1 : -1)

var mm = MeasurementMatrix.gaussian(m, n)
double[] y = mm.measure(xTrue)

var opts = new BasisPursuit.Options(0.0, 10.0, 1e-10, 300, 50, 1e-4)
var bp = BasisPursuit.fit(mm.phi(), y, opts)
bp.iter()

var Ax = mm.phi().mv(Vector.column(bp.x()))
double[] res = new double[m]
for (int i = 0; i < m; i++) res[i] = Ax.get(i) - y[i]
MathEx.norm(res) / MathEx.norm(y)

5
2.10e-12

For noisy data set epsilon ≈ σ √m (expected ‖ε‖₂) in BasisPursuit.Options. OMP/CoSaMP instead loosen tol to about the noise-to-signal ratio.

Choosing an Algorithm

Aim for m ≥ 4k at least; for high-probability recovery use roughly m ≥ 4k log(n/k). Increase m if the relative error stays large.

Wavelet-Domain Sparsity

Images and audio are often dense in the time/space domain but sparse after a wavelet transform. Chain MeasurementMatrix.gaussian(m, n).withWavelet(psi) so solvers operate on the coefficient domain (A = Φ · IDWT). Recover sparse coefficients with OMP, then synthesize the signal via psi.inverse.

import java.util.Arrays
import smile.cs.*
import smile.math.MathEx
import smile.wavelet.HaarWavelet

MathEx.setSeed(19650218)
int n = 128, m = 50, k = 5

double[] s = new double[n]
int[] locs = MathEx.permutate(n)
for (int i = 0; i < k; i++) {
    double amp = 1.0 + MathEx.random()
    s[locs[i]] = MathEx.random() < 0.5 ? amp : -amp
}

var wm = MeasurementMatrix.gaussian(m, n).withWavelet(new HaarWavelet())
double[] y = wm.measureSparse(s)

var omp = OMP.fit(wm.toMatrix(), y, k)
Arrays.toString(omp.support())

double[] d = omp.x().clone()
for (int i = 0; i < n; i++) d[i] -= s[i]
MathEx.norm(d) / MathEx.norm(s)

[34, 110, 56, 20, 114]
2.50e-16

Haar is a good default for piecewise-constant structure; DaubechiesWavelet(8) or higher suits smoother signals. Signal length must be a power of two for the DWT.

Options Persistence

All Options records serialize through java.util.Properties:

import smile.cs.*

var opts = new OMP.Options(8, 1e-8)
var props = opts.toProperties()
var restored = OMP.Options.of(props)
restored.sparsity()
restored.tol()

8
1.0E-8

References

  1. E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag., 25(2):21–30, 2008.
  2. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, 52(4):1289–1306, 2006.
  3. J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Trans. Inf. Theory, 53(12):4655–4666, 2007.
  4. D. Needell and J. A. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples,” Appl. Comput. Harmon. Anal., 26(3):301–321, 2009.
  5. S. S. Chen, D. L. Donoho and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Rev., 43(1):129–159, 2001.
←   Wavelet Graph  →