SMILE — Compressed Sensing User Guide

Overview
Quick Start
Core Concepts
Measurement Matrix
Recovery Algorithms
Hyperparameter Tuning
Choosing a Sensing Matrix
Noisy Measurements
Wavelet-Domain Sparsity
Complete Examples
Mathematical Background
References

Overview

The smile.cs package implements Compressed Sensing (CS), a signal-acquisition paradigm that exploits signal sparsity to recover a high-dimensional signal from far fewer measurements than the classical Nyquist–Shannon theorem requires.

Given a signal x ∈ ℝⁿ that is k-sparse (has at most k non-zero entries), and a measurement matrix A ∈ ℝ^{m×n} with m ≪ n, CS recovers x from y = Ax exactly (or approximately in the noisy case) when m = O(k log(n/k)).

Class	Algorithm	Type	Best for
`MeasurementMatrix`	—	Sensing operator	Building `y = Φ x`
`OMP`	Orthogonal Matching Pursuit	Greedy	Fast, known sparsity `k`
`CoSaMP`	Compressive Sampling Matching Pursuit	Greedy	Provable guarantees, near-optimal
`BasisPursuit`	L1-minimization via interior point	Convex	Best recovery quality, noisy data

Quick Start

java

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

MathEx.setSeed(42);

// 1. Build a k-sparse ground-truth signal
int n = 128, m = 50, k = 5;
double[] xTrue = new double[n];
int[] locs = MathEx.permutate(n);
for (int i = 0; i < k; i++) xTrue[locs[i]] = 1.0 + MathEx.random();

// 2. Create a Gaussian sensing matrix and take measurements
MeasurementMatrix mm = MeasurementMatrix.gaussian(m, n);
double[] y = mm.measure(xTrue);

// 3. Recover with OMP (fastest, requires known k)
OMP omp = OMP.fit(mm.phi(), y, k);
System.out.printf("OMP support: %s%n", java.util.Arrays.toString(omp.support()));

// 4. Recover with CoSaMP (provable guarantees)
CoSaMP cosamp = CoSaMP.fit(mm.phi(), y, k);

// 5. Recover with Basis Pursuit (best quality, convex)
BasisPursuit bp = BasisPursuit.fit(mm.phi(), y);

// 6. Check relative L2 error
double err = relError(omp.x(), xTrue);
System.out.printf("OMP relative error: %.2e%n", err);

java

static double relError(double[] xRec, double[] xTrue) {
    double num = 0, den = 0;
    for (int i = 0; i < xTrue.length; i++) {
        double d = xRec[i] - xTrue[i];
        num += d * d;
        den += xTrue[i] * xTrue[i];
    }
    return Math.sqrt(num / den);
}

Core Concepts

Sparsity

A signal x ∈ ℝⁿ is k-sparse if it has at most k non-zero entries:

‖x‖₀  =  |{ i : x[i] ≠ 0 }|  ≤  k

The signal need not be sparse in the original domain — it suffices that it is sparse in some known basis Ψ. For example, natural images are sparse in a wavelet basis; audio signals are sparse in the Fourier basis.

The Measurement Model

The forward model is:

y  =  A x  =  Φ Ψ s

where:

Symbol	Meaning
`x ∈ ℝⁿ`	Original signal
`y ∈ ℝᵐ`	Measurement vector (`m ≪ n`)
`Φ ∈ ℝ^{m×n}`	Sensing (measurement) matrix
`Ψ ∈ ℝ^{n×n}`	Sparsifying basis (e.g. wavelet transform)
`s ∈ ℝⁿ`	Sparse coefficients in Ψ-domain
`A = Φ Ψ`	Compound measurement matrix

When the signal is already sparse in the canonical basis, Ψ = I and A = Φ.

Restricted Isometry Property (RIP)

The RIP is the key theoretical condition guaranteeing recovery. A matrix A satisfies the RIP of order k with constant δ_k if, for all k-sparse vectors x:

(1 − δ_k) ‖x‖₂²  ≤  ‖Ax‖₂²  ≤  (1 + δ_k) ‖x‖₂²

Gaussian and Bernoulli matrices satisfy the RIP with high probability when m ≥ C · k · log(n/k) for a universal constant C ≈ 1–2. This means the sensing matrix approximately preserves the geometry of sparse signals.

Sample Complexity

The minimum number of measurements for reliable recovery is:

Regime	Required `m`
Exact recovery (noiseless)	`m ≥ C k log(n/k)`
Stable recovery (noisy)	`m ≥ C k log(n/k)` (same order)
Practical rule of thumb	`m ≥ 4k` to `m ≥ 5k`

For example, recovering a k=5-sparse signal in n=128 dimensions requires approximately m ≥ 20–25 measurements in practice.

Measurement Matrix

MeasurementMatrix is an immutable record that encapsulates:

The sensing matrix Φ (a DenseMatrix), and
An optional wavelet sparsifying basis Ψ.

All factory methods produce sensing matrices whose entries are scaled so that E[‖Φ x‖₂²] ≈ ‖x‖₂² (isometry in expectation).

Gaussian

java

MeasurementMatrix mm = MeasurementMatrix.gaussian(m, n);

Entries are i.i.d. N(0, 1/m). Gaussian matrices satisfy the RIP with optimal sample complexity and are incoherent with every fixed basis. They are the standard choice for general-purpose compressed sensing.

When to use: default choice; strongest theoretical guarantees.

Bernoulli

java

MeasurementMatrix mm = MeasurementMatrix.bernoulli(m, n);

Entries are i.i.d. ±1/√m with equal probability. Equivalent RIP guarantees to Gaussian, but entries are binary which can be advantageous for hardware implementations.

When to use: hardware-friendly sensing, or when Gaussian random numbers are expensive to generate.

Partial Identity

A random row-subsampling of the n × n identity — equivalently, measuring m randomly selected components of x:

java

// Automatic random selection of m rows
MeasurementMatrix mm = MeasurementMatrix.partial(m, n);

// Explicit row indices (must be in [0, n))
int[] rows = {0, 7, 15, 22, 40};
MeasurementMatrix mm2 = MeasurementMatrix.partial(rows, n);

When to use: when only a random subset of signal samples can be acquired (e.g. missing data, random sub-Nyquist sampling).
Requirement: the signal must be sparse in a basis incoherent with the identity (e.g. the DFT basis); it will not work if the signal is sparse in the canonical (time/space) domain itself.

Wavelet Sparsifying Basis

When the signal is sparse in a wavelet domain, chain the sensing matrix with a wavelet basis using withWavelet(Wavelet):

java

import smile.wavelet.*;

// Compound operator A = Φ · IDWT  (m × n, implicit)
MeasurementMatrix wm = MeasurementMatrix.gaussian(m, n)
                                        .withWavelet(new DaubechiesWavelet(8));

// Take measurements: y = Φ · IDWT(s)  where s is the sparse wavelet vector
double[] ySparse = wm.measureSparse(s);

// Back-project to wavelet domain: v = DWT(Φ^T y)
double[] v = wm.backProjectSparse(y);

// Get the implicit compound matrix for use with solvers
Matrix A = wm.toMatrix();

The wavelet operator is applied implicitly (never forming the full n × n matrix), making this memory-efficient for large n.

Supported wavelets: any subclass of smile.wavelet.Wavelet — HaarWavelet, D4Wavelet, DaubechiesWavelet, CoifletWavelet, SymletWavelet, BestLocalizedWavelet.

Adjoint Operations

Method	Computes	Notes
`mm.measure(x)`	`y = Φ x`	Forward measurement
`mm.measureSparse(s)`	`y = Φ · IDWT(s)`	Forward in wavelet domain
`mm.backProject(y)`	`v = Φᵀ y`	Back-projection (adjoint)
`mm.backProjectSparse(y)`	`v = DWT(Φᵀ y)`	Adjoint in wavelet domain

The adjoint identity holds exactly (to floating-point precision):

⟨Φ x, y⟩  =  ⟨x, Φᵀ y⟩

Implicit Matrix Interface

mm.toMatrix() returns a smile.tensor.Matrix that supports both mv (forward) and tv (adjoint) products without materializing the full compound matrix. Pass this directly to the solvers:

java

Matrix A = mm.toMatrix();
OMP result = OMP.fit(A, y, k);

Recovery Algorithms

All three algorithms share the same call convention:

java

SolverResult result = Solver.fit(A, y, /* options */);
double[] xRecovered = result.x();       // recovered signal, length n
int[]    support    = result.support(); // non-zero index set (OMP, CoSaMP)
int      iters      = result.iter();    // iterations performed

Orthogonal Matching Pursuit (OMP)

OMP is a greedy algorithm that builds the support set one column at a time. At each iteration it selects the column of A most correlated with the current residual, augments the support, and solves a small least-squares problem on the active support.

Convergence guarantee: for a Gaussian matrix with m ≥ C k log n, OMP recovers any k-sparse x exactly in exactly k iterations when δ_{2k} < √2 − 1 ≈ 0.41.

java

// With known sparsity k
OMP result = OMP.fit(A, y, k);

// With explicit options (custom tolerance)
var opts = new OMP.Options(k, 1e-8);
OMP result = OMP.fit(A, y, opts);

// Inspect results
System.out.println("Support: " + Arrays.toString(result.support()));
System.out.println("Iterations: " + result.iter());

Algorithm sketch:

Initialize residual r = y, support S = ∅.
Find j* = argmax_j |⟨A_j, r⟩| / ‖A_j‖₂.
Add j* to S; solve min ‖y − A_S c‖₂ via QR (Gram–Schmidt).
Update residual r = y − A_S c_S.
Repeat until |S| = k or ‖r‖₂ / ‖y‖₂ < tol.

Complexity per iteration: O(mn) correlation + O(mk²) least-squares.
Total: O(k(mn + mk²)).

When to use:

Sparsity k is known or can be estimated.
Speed is the primary concern.
Signal-to-noise ratio is high.

CoSaMP

CoSaMP (Compressive Sampling Matching Pursuit) is a greedy algorithm with provably near-optimal recovery guarantees — it achieves essentially the same error as the best k-term approximation.

Unlike OMP, CoSaMP maintains a candidate support of size 2k, prunes it back to k at each step, and re-solves least-squares on the pruned support.

java

// With known sparsity k
CoSaMP result = CoSaMP.fit(A, y, k);

// Over-estimated sparsity (2k) still recovers correctly
CoSaMP result = CoSaMP.fit(A, y, 2 * k);

// With explicit options
var opts = new CoSaMP.Options(k, 100, 1e-8);
CoSaMP result = CoSaMP.fit(A, y, opts);

Algorithm sketch (one iteration):

Proxy: e = Aᵀ r (gradient of data-fit term).
Identify: find the 2k largest entries of |e|; form union Ω with current support.
Least-squares: b = argmin ‖y − A_Ω b‖₂ on the union support Ω.
Prune: retain only the k largest entries of b.
Update: x_new = least-squares on the pruned support; r = y − A x_new.

Convergence: O(log(‖x‖₂ / ε)) iterations to reach error ε.

When to use:

Provable guarantees are needed.
Sparsity may be slightly over-estimated.
Moderate noise tolerance is required.

Basis Pursuit / BPDN

Basis Pursuit solves the convex L1-minimization problem:

Basis Pursuit (noiseless):
  minimise  ‖x‖₁
  subject to  Ax = y

Basis Pursuit DeNoising — BPDN (noisy):
  minimise  ‖x‖₁
  subject to  ‖Ax − y‖₂ ≤ ε

The implementation uses a log-barrier interior-point method (primal–dual Newton's method). Each Newton step requires solving a linear system (A D A^T) ν = b where D is a positive-definite diagonal matrix; this is solved with a preconditioned conjugate-gradient (CG) method.

java

// Exact Basis Pursuit (epsilon = 0)
BasisPursuit bp = BasisPursuit.fit(A, y);

// BPDN with explicit noise bound epsilon
double epsilon = 0.01 * MathEx.norm(y);
var opts = new BasisPursuit.Options(epsilon, 10.0, 1e-8, 200, 50, 1e-3);
BasisPursuit bp = BasisPursuit.fit(A, y, opts);

// Check feasibility
double[] Ax = BasisPursuit.matvec(A, bp.x(), m, n, false);

Convergence: O(log(1/ε)) outer Newton iterations; each inner CG solve takes O(mn) per iteration.

When to use:

Best recovery quality is needed, especially in moderate noise.
k is unknown (L1-minimization automatically promotes sparsity).
The problem is small-to-medium sized (inner CG per Newton step is expensive).

Algorithm Comparison

Property	OMP	CoSaMP	Basis Pursuit
Type	Greedy	Greedy	Convex
Recovery guarantee	Exact (RIP δ_{2k} < √2−1)	Near-optimal (RIP δ_{4k} < 0.1)	Near-optimal (δ_{2k} < √2−1)
Requires known `k`	Yes	Yes (robust to over-estimate)	No
Noisy robustness	Moderate	Good	Best
Speed	Fastest	Fast	Slowest
Memory	O(mk)	O(mn)	O(mn)
Iterations	Exactly `k`	O(log(1/ε))	O(log(1/ε)) outer + CG inner
Output	x, support, iter	x, support, iter	x, iter

Hyperparameter Tuning

OMP Options

java

OMP.Options opts = new OMP.Options(
    k,     // sparsity:  target number of non-zero components (> 0)
    1e-6   // tol:       stop early when ‖r‖₂/‖y‖₂ < tol (≥ 0)
);

Parameter	Default	Effect
`sparsity`	required	Maximum support size; algorithm terminates after exactly `k` iterations unless `tol` is reached first
`tol`	`1e-6`	Early stopping on relative residual; set to `0` to always run exactly `k` iterations

Guidelines:

Set sparsity to your best estimate of k. An over-estimate increases the number of iterations but usually does not hurt accuracy (the extra selected atoms will have near-zero coefficients).
Tighten tol (e.g. 1e-10) for high-precision noiseless recovery.
Loosen tol (e.g. 1e-4) for noisy data to avoid over-fitting the noise.

CoSaMP Options

java

CoSaMP.Options opts = new CoSaMP.Options(
    k,     // sparsity:  target sparsity level (> 0)
    50,    // maxIter:   maximum number of iterations (> 0)
    1e-6   // tol:       relative residual tolerance (≥ 0)
);

Parameter	Default	Effect
`sparsity`	required	Target support size `k`; CoSaMP is robust to moderate over-estimation
`maxIter`	`50`	Iteration budget; CoSaMP typically converges in 10–30 iterations
`tol`	`1e-6`	Stop when `‖r‖₂/‖y‖₂ < tol`

Guidelines:

CoSaMP requires 2k < m for the internal 2k-proxy step to be well-defined. If 2k ≥ m, a warning is logged and the result may be unreliable.
For noisy data, set tol to a value slightly larger than the expected noise-to-signal ratio (e.g. 1e-3 for 0.1% noise).

BasisPursuit Options

java

BasisPursuit.Options opts = new BasisPursuit.Options(
    epsilon,   // L2 constraint bound (0.0 for exact BP)
    10.0,      // mu: barrier update factor (> 1)
    1e-8,      // cgtol: inner CG convergence tolerance (> 0)
    200,       // cgMaxIter: maximum CG iterations per Newton step (> 0)
    50,        // maxIter: maximum outer Newton iterations (> 0)
    1e-3       // tol: duality-gap stopping criterion (> 0)
);

Parameter	Default	Effect
`epsilon`	`0.0`	Noise bound `ε`; set to `0` for exact BP (noiseless), or to an estimate of `‖noise‖₂` for BPDN
`mu`	`10.0`	Barrier parameter growth rate per outer iteration; larger → fewer outer iters but less stable
`cgtol`	`1e-8`	Inner CG solver tolerance; tighten for higher precision
`cgMaxIter`	`200`	CG budget per Newton step; increase for ill-conditioned problems
`maxIter`	`50`	Outer Newton iteration budget; typically 10–30 suffice
`tol`	`1e-3`	Duality-gap stopping criterion; tighten for higher precision

Guidelines:

For exact BP (epsilon = 0), use tight cgtol = 1e-10 and a generous cgMaxIter = 300 to ensure the Newton steps are accurate enough.
For BPDN (epsilon > 0), a natural choice for epsilon is epsilon = σ √m where σ is the standard deviation of the measurement noise. Alternatively use epsilon = 0.05 * ‖y‖₂ as a conservative estimate.
Increase mu to 20 or 50 if the algorithm is converging too slowly; decrease to 5 if it is oscillating.

Persisting Options

All Options records support serialization to/from java.util.Properties:

java

// Serialize
OMP.Options opts = new OMP.Options(5, 1e-6);
Properties props = opts.toProperties();
props.store(new FileWriter("omp.properties"), "OMP options");

// Deserialize
Properties loaded = new Properties();
loaded.load(new FileReader("omp.properties"));
OMP.Options opts2 = OMP.Options.of(loaded);

The same pattern applies to CoSaMP.Options and BasisPursuit.Options.

Choosing a Sensing Matrix

Signal type	Sensing matrix	Notes
General `k`-sparse in canonical basis	`gaussian(m, n)`	Best theoretical guarantees
Binary/integer signal values	`bernoulli(m, n)`	Hardware-friendly
Signal sparse in frequency (DFT)	`partial(m, n)`	Random sub-sampling in time
Signal sparse in wavelet domain	`gaussian(m, n).withWavelet(w)`	Memory-efficient implicit operator
Very large `n` (implicit operator)	Any `mm.toMatrix()`	DWT applied without forming Ψ

Rule of thumb for m: aim for m ≥ 4k at minimum. For reliable recovery with high probability use m ≥ 4k log(n/k). Check the recovery error and increase m if it is unsatisfactory.

Noisy Measurements

In practice, measurements always contain noise: y = Ax + ε. All three algorithms handle noise, but with different mechanisms:

OMP (noisy)

Use the tol option to stop early before the algorithm fits the noise:

java

// Stop when relative residual < noise level
double noiseLevel = sigma * Math.sqrt(m) / MathEx.norm(y);
OMP result = OMP.fit(A, y, new OMP.Options(k, noiseLevel));

CoSaMP (noisy)

CoSaMP has built-in noise stability. The recovery error satisfies:

‖x_rec − x‖₂  ≤  C · ‖ε‖₂ / √m

Use the default options and set tol to match the noise level:

java

double sigma = 0.01;  // per-measurement noise std
double[] yNoisy = y.clone();
for (int i = 0; i < m; i++) yNoisy[i] += sigma * rng.nextGaussian();

CoSaMP result = CoSaMP.fit(A, yNoisy, new CoSaMP.Options(k, 100, sigma));

Basis Pursuit Denoising (BPDN)

Set epsilon to bound the expected noise norm ‖ε‖₂ ≈ σ √m:

java

double sigma = 0.01;
double epsilon = sigma * Math.sqrt(m);  // expected ‖ε‖₂

var opts = new BasisPursuit.Options(epsilon, 10.0, 1e-8, 200, 50, 1e-3);
BasisPursuit result = BasisPursuit.fit(A, yNoisy, opts);

Wavelet-Domain Sparsity

Natural signals such as audio and images are often not sparse in the canonical basis, but become sparse after a wavelet transform. The MeasurementMatrix.withWavelet() API handles this transparently.

Setup

java

import smile.wavelet.*;

int n = 256;  // signal length (must be power of 2)
int m = 60;   // number of measurements
int k = 8;    // sparsity in wavelet domain

// True sparse wavelet coefficients
double[] sSparse = new double[n];
sSparse[3] = 2.5; sSparse[17] = -1.8; // ... k non-zeros

// Synthesise the signal from its wavelet coefficients
double[] xSignal = sSparse.clone();
new DaubechiesWavelet(8).inverse(xSignal);   // x = IDWT(s)

// Build compound measurement matrix A = Φ · IDWT
MeasurementMatrix wm = MeasurementMatrix.gaussian(m, n)
                                        .withWavelet(new DaubechiesWavelet(8));

// Measurements: y = Φ · IDWT(s) = Φ x
double[] y = wm.measureSparse(sSparse);

Recovery

The recovery algorithms operate on the wavelet coefficient domain — they recover the sparse s, not x directly:

java

// Get the implicit compound matrix A = Φ Ψ
Matrix A = wm.toMatrix();   // m × n, implicit

// Recover sparse coefficients s via OMP
OMP result = OMP.fit(A, y, k);
double[] sRec = result.x();

// Synthesise the recovered signal: x_rec = IDWT(s_rec)
double[] xRec = sRec.clone();
new DaubechiesWavelet(8).inverse(xRec);

Choosing a Wavelet for CS

Wavelet	Notes
`HaarWavelet`	Fastest; good for piecewise-constant signals
`DaubechiesWavelet(8)`	Good general choice for smooth signals
`DaubechiesWavelet(20)`	Best sparsity for very smooth signals
`CoifletWavelet(12)`	Good when scaling coefficients must match samples
`SymletWavelet(8)`	Near-symmetric; reduced boundary effects

Complete Examples

Example 1 — Basic Sparse Recovery with All Three Algorithms

java

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

public class BasicRecoveryExample {
    public static void main(String[] args) {
        MathEx.setSeed(2024);
        int n = 256, m = 60, k = 8;

        // Build k-sparse ground truth
        double[] xTrue = new double[n];
        int[] locs = MathEx.permutate(n);
        for (int i = 0; i < k; i++) {
            xTrue[locs[i]] = (i % 2 == 0 ? 1 : -1) * (1.0 + MathEx.random());
        }

        // Gaussian sensing matrix
        MeasurementMatrix mm = MeasurementMatrix.gaussian(m, n);
        double[] y = mm.measure(xTrue);

        // OMP
        OMP omp = OMP.fit(mm.phi(), y, k);
        System.out.printf("OMP:   rel error = %.2e, support = %s%n",
                relError(omp.x(), xTrue), Arrays.toString(omp.support()));

        // CoSaMP
        CoSaMP cosamp = CoSaMP.fit(mm.phi(), y, k);
        System.out.printf("CoSaMP: rel error = %.2e, iters = %d%n",
                relError(cosamp.x(), xTrue), cosamp.iter());

        // Basis Pursuit (exact)
        var bpOpts = new BasisPursuit.Options(0.0, 10.0, 1e-10, 300, 50, 1e-4);
        BasisPursuit bp = BasisPursuit.fit(mm.phi(), y, bpOpts);
        System.out.printf("BP:    rel error = %.2e, iters = %d%n",
                relError(bp.x(), xTrue), bp.iter());
    }

    static double relError(double[] xRec, double[] xTrue) {
        double num = 0, den = 0;
        for (int i = 0; i < xTrue.length; i++) {
            double d = xRec[i] - xTrue[i];
            num += d * d; den += xTrue[i] * xTrue[i];
        }
        return Math.sqrt(num / den);
    }
}

Example 2 — Noisy Recovery with BPDN

java

import smile.cs.*;
import smile.math.MathEx;
import java.util.Random;

public class NoisyRecoveryExample {
    public static void main(String[] args) {
        MathEx.setSeed(42);
        Random rng = new Random(42);
        int n = 128, m = 50, k = 5;
        double sigma = 0.05;   // per-measurement noise std

        // Ground truth
        double[] xTrue = new double[n];
        int[] locs = MathEx.permutate(n);
        for (int i = 0; i < k; i++) xTrue[locs[i]] = 1.0 + MathEx.random();

        // Noisy measurements
        MeasurementMatrix mm = MeasurementMatrix.gaussian(m, n);
        double[] yClean = mm.measure(xTrue);
        double[] yNoisy = yClean.clone();
        for (int i = 0; i < m; i++) yNoisy[i] += sigma * rng.nextGaussian();

        System.out.printf("Noise level: ‖ε‖₂ = %.4f%n",
                MathEx.norm(subtract(yNoisy, yClean)));

        // BPDN: set epsilon = sigma * sqrt(m)
        double epsilon = sigma * Math.sqrt(m);
        var opts = new BasisPursuit.Options(epsilon, 10.0, 1e-8, 200, 50, 1e-3);
        BasisPursuit bp = BasisPursuit.fit(mm.phi(), yNoisy, opts);

        // Verify residual constraint
        double[] Ax = BasisPursuit.matvec(mm.phi(), bp.x(), m, n, false);
        double resNorm = MathEx.norm(subtract(Ax, yNoisy));
        System.out.printf("BPDN: ‖Ax−y‖₂ = %.4f, epsilon = %.4f%n", resNorm, epsilon);
        System.out.printf("BPDN: rel error = %.4f%n",
                MathEx.norm(subtract(bp.x(), xTrue)) / MathEx.norm(xTrue));
    }

    static double[] subtract(double[] a, double[] b) {
        double[] r = a.clone();
        for (int i = 0; i < r.length; i++) r[i] -= b[i];
        return r;
    }
}

Example 3 — Wavelet-Domain Sparse Recovery

java

import smile.cs.*;
import smile.math.MathEx;
import smile.wavelet.*;

public class WaveletRecoveryExample {
    public static void main(String[] args) {
        MathEx.setSeed(100);
        int n = 128, m = 35, k = 6;
        Wavelet psi = new DaubechiesWavelet(8);

        // k-sparse signal in wavelet domain
        double[] sSparse = new double[n];
        int[] locs = MathEx.permutate(n);
        for (int i = 0; i < k; i++) sSparse[locs[i]] = 1.5 + MathEx.random();

        // Synthesise signal x = IDWT(s)
        double[] xSignal = sSparse.clone();
        psi.inverse(xSignal);

        // Compound measurement matrix A = Φ IDWT
        MeasurementMatrix wm = MeasurementMatrix.gaussian(m, n).withWavelet(psi);
        double[] y = wm.measureSparse(sSparse);

        // Recover sparse coefficients via OMP on A
        OMP result = OMP.fit(wm.toMatrix(), y, k);
        double[] sRec = result.x();

        // Reconstruct signal
        double[] xRec = sRec.clone();
        psi.inverse(xRec);

        // Report
        double coeffErr = MathEx.norm(subtract(sRec, sSparse)) / MathEx.norm(sSparse);
        double signalErr = MathEx.norm(subtract(xRec, xSignal)) / MathEx.norm(xSignal);
        System.out.printf("Wavelet coeff error: %.2e%n", coeffErr);
        System.out.printf("Signal error:         %.2e%n", signalErr);
        System.out.printf("Support: %s%n", java.util.Arrays.toString(result.support()));
    }

    static double[] subtract(double[] a, double[] b) {
        double[] r = a.clone(); for (int i = 0; i < r.length; i++) r[i] -= b[i]; return r;
    }
}

Example 4 — Persisting and Loading Options

java

import smile.cs.*;
import java.util.Properties;
import java.io.*;

public class OptionsSerializationExample {
    public static void main(String[] args) throws IOException {
        // Create and save OMP options
        OMP.Options ompOpts = new OMP.Options(8, 1e-8);
        Properties p = ompOpts.toProperties();
        try (var w = new FileWriter("omp.properties")) {
            p.store(w, "OMP recovery options");
        }

        // Load back
        Properties loaded = new Properties();
        try (var r = new FileReader("omp.properties")) {
            loaded.load(r);
        }
        OMP.Options restored = OMP.Options.of(loaded);
        System.out.printf("Restored: sparsity=%d, tol=%.2e%n",
                restored.sparsity(), restored.tol());

        // Same pattern for BasisPursuit
        BasisPursuit.Options bpOpts = new BasisPursuit.Options(0.02, 10.0, 1e-8, 200, 50, 1e-3);
        Properties bp = bpOpts.toProperties();
        BasisPursuit.Options bpRestored = BasisPursuit.Options.of(bp);
        System.out.printf("BP epsilon=%.4f, mu=%.1f%n",
                bpRestored.epsilon(), bpRestored.mu());
    }
}

Example 5 — Adjoint / Incoherence Verification

java

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

public class AdjointCheckExample {
    public static void main(String[] args) {
        MathEx.setSeed(0);
        int m = 30, n = 100;
        MeasurementMatrix mm = MeasurementMatrix.gaussian(m, n);

        double[] x = new double[n];
        double[] yv = new double[m];
        for (int i = 0; i < n; i++) x[i] = MathEx.random() - 0.5;
        for (int i = 0; i < m; i++) yv[i] = MathEx.random() - 0.5;

        // Verify ⟨Φx, y⟩ = ⟨x, Φᵀy⟩
        double[] Ax  = mm.measure(x);
        double[] Aty = mm.backProject(yv);
        double lhs = MathEx.dot(Ax, yv);
        double rhs = MathEx.dot(x, Aty);
        System.out.printf("⟨Φx,y⟩  = %.10f%n", lhs);
        System.out.printf("⟨x,Φᵀy⟩ = %.10f%n", rhs);
        System.out.printf("Error:   = %.2e%n", Math.abs(lhs - rhs));
    }
}

Mathematical Background

Sparse Recovery as L0 Minimization

The ideal (but NP-hard) formulation is:

minimise  ‖x‖₀   subject to  Ax = y

where ‖x‖₀ counts the number of non-zeros. CS theory shows that under the RIP condition, the convex relaxation

minimise  ‖x‖₁   subject to  Ax = y

achieves the same solution. This is the Basis Pursuit problem.

Why L1 Promotes Sparsity

The unit ball of ‖·‖₁ in ℝⁿ is a cross-polytope (diamond in 2D). When a hyperplane {x : Ax = y} is tangent to this ball, it touches a vertex or edge — corresponding to a sparse vector. This geometric intuition explains why L1 minimization recovers sparse solutions.

RIP and Incoherence

Two equivalent conditions guarantee recovery:

RIP: A has δ_{2k} < √2 − 1. Gaussian and Bernoulli matrices satisfy this with m = O(k log(n/k)) rows.
Mutual incoherence: The coherence μ(A) = max_{i≠j} |⟨A_i, A_j⟩| between columns of the normalized A should be small. OMP and BP provably recover any k-sparse signal when k < (1 + 1/μ) / 2.

Universal Threshold for `m`

For Gaussian Φ and a fixed sparsifying basis Ψ, exact recovery holds with probability ≥ 1 − δ when:

m  ≥  C · k · log(n/k) · log(1/δ)

For k = 5, n = 128, δ = 0.01:

m  ≥  2 · 5 · log(25.6) · log(100)  ≈  75

In practice the constant C ≈ 1–2 and m ≥ 4k often suffices empirically.

Gram–Schmidt Least-Squares in OMP

OMP solves the active-support least-squares problem via incremental QR (Gram–Schmidt orthogonalization). At step t, the QR factors of A_{S_t} are updated in O(mt) time by orthogonalizing the new column against the existing Q-vectors. The coefficient estimate is then c = Qᵀ y (projection) followed by back-substitution through R.

Interior-Point Method for Basis Pursuit

Basis Pursuit is solved by lifting to:

minimise  1ᵀu   subject to  Ax = y,  −u ≤ x ≤ u

and applying a log-barrier (primal–dual interior-point) method:

minimise  t · 1ᵀu − Σ log(u_i + x_i) − Σ log(u_i − x_i)
subject to  Ax = y

The barrier parameter t is multiplied by μ > 1 at each outer iteration. Each Newton step requires solving a system of the form (A D A^T) ν = b where D = diag(1/w) is computed from the barrier Hessian; this is solved with CG.

References

E. J. Candès and M. B. Wakin, "An introduction to compressive sampling," IEEE Signal Process. Mag., 25(2):21–30, 2008.
E. J. Candès and T. Tao, "Near-optimal signal recovery from random projections: Universal encoding strategies?" IEEE Trans. Inf. Theory, 52(12):5406–5425, 2006.
D. L. Donoho, "Compressed sensing," IEEE Trans. Inf. Theory, 52(4):1289–1306, 2006.
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.
D. Needell and J. A. Tropp, "CoSaMP: Iterative signal recovery from incomplete and inaccurate samples," Appl. Comput. Harmon. Anal., 26(3):301–321, 2009.
S. S. Chen, D. L. Donoho and M. A. Saunders, "Atomic decomposition by basis pursuit," SIAM Rev., 43(1):129–159, 2001.
E. J. Candès, J. Romberg and T. Tao, "Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information," IEEE Trans. Inf. Theory, 52(2):489–509, 2006.

SMILE — Compressed Sensing User Guide

SMILE — Compressed Sensing User Guide

Table of Contents

Overview

Quick Start

Core Concepts

Sparsity

The Measurement Model

Restricted Isometry Property (RIP)

Sample Complexity

Measurement Matrix

Gaussian

Bernoulli

Partial Identity

Wavelet Sparsifying Basis

Adjoint Operations

Implicit Matrix Interface

Recovery Algorithms

Orthogonal Matching Pursuit (OMP)

CoSaMP

Basis Pursuit / BPDN

Algorithm Comparison

Hyperparameter Tuning

OMP Options

CoSaMP Options

BasisPursuit Options

Persisting Options

Choosing a Sensing Matrix

Noisy Measurements

OMP (noisy)

CoSaMP (noisy)

Basis Pursuit Denoising (BPDN)

Wavelet-Domain Sparsity

Setup

Recovery

Choosing a Wavelet for CS

Complete Examples

Example 1 — Basic Sparse Recovery with All Three Algorithms

Example 2 — Noisy Recovery with BPDN

Example 3 — Wavelet-Domain Sparse Recovery

Example 4 — Persisting and Loading Options

Example 5 — Adjoint / Incoherence Verification

Mathematical Background

Sparse Recovery as L0 Minimization

Why L1 Promotes Sparsity

RIP and Incoherence

Universal Threshold for m

Gram–Schmidt Least-Squares in OMP

Interior-Point Method for Basis Pursuit

References

Universal Threshold for `m`