Signed Spikes, Gaussian Measurements

In this example we have

1. A signal \(\bx\) of length \(n=2560\) with \(k=20\) signed spikes. Each spike has a magnitude of 1. The sign for each spike is randomly assigned. The locations of spikes in the signal are also randomly chosen.

2. A Gaussian sensing matrix \(\Phi\) of shape \(m \times n = 600 \times 2560\) making 600 random measurements in a vector \(\bb\) given by the sensing equation \(\bb = \Phi \bx\). The columns of the sensing matrix are unit normalized.

The signal is sparse in standard basis. This is a relatively easy sparse recovery problem and focuses on the compressive sensing process modeled as:

(1)\[\bb = \bA \bx = \Phi \bx.\]

See also:

• Test Problems

• Linear Operators

# Configure JAX to work with 64-bit floating point precision.
from jax.config import config
config.update("jax_enable_x64", True)

import jax.numpy as jnp
import cr.nimble as crn
import cr.sparse.plots as crplot

Setup

We shall construct our test signal, measurements and sensing matrix using our test problems module.

from cr.sparse import problems
k=20
m=600
n=2560
prob = problems.generate('signed-spikes:dirac:gaussian', k=k, m=m, n=n)
fig, ax = problems.plot(prob)

Let us access the relevant parts of our test problem

# The Gaussian sensing matrix operator
A = prob.A
# The measurements
b0 = prob.b
# The sparse signal
x0 = prob.x

Sparse Recovery using Subspace Pursuit

We shall use subspace pursuit to reconstruct the signal.

import cr.sparse.pursuit.sp as sp
# We will try to estimate a k-sparse representation
sol = sp.solve(A, b0, k)

This utility function helps us quickly analyze the quality of reconstruction

problems.analyze_solution(prob, sol)

m: 600, n: 2560
b_norm: original: 4.441 reconstruction: 4.441 SNR: 297.61 dB
x_norm: original: 4.472 reconstruction: 4.472 SNR: 298.08 dB
Sparsity: original: 20, reconstructed: 20, overlap: 20, ratio: 1.000
Iterations: 2

The estimated sparse signal

x = sol.x

Let us reconstruct the measurements from this signal

b = A.times(x)

Let us visualize the original and reconstructed signal

def plot_signals(x0, x):
    ax = crplot.h_plots(2)
    ax[0].stem(x0, markerfmt='.')
    ax[0].set_title('Original signal')
    ax[1].stem(x, markerfmt='.')
    ax[1].set_title('Reconstructed signal')
plot_signals(x0, x)

Let us visualize the original and reconstructed measurements

def plot_measurments(b0, b):
    ax = crplot.h_plots(2)
    ax[0].plot(b0)
    ax[0].set_title('Original measurements')
    ax[1].plot(b)
    ax[1].set_title('Reconstructed measurements')
plot_measurments(b0, b)

Sparse Recovery using Compressive Sampling Matching Pursuit

We shall now use compressive sampling matching pursuit to reconstruct the signal.

import cr.sparse.pursuit.cosamp as cosamp
# We will try to estimate a k-sparse representation
sol = cosamp.solve(A, b0, k)
problems.analyze_solution(prob, sol)

m: 600, n: 2560
b_norm: original: 4.441 reconstruction: 4.441 SNR: 293.05 dB
x_norm: original: 4.472 reconstruction: 4.472 SNR: 292.62 dB
Sparsity: original: 20, reconstructed: 20, overlap: 20, ratio: 1.000
Iterations: 2

The estimated sparse signal

x = sol.x

Let us reconstruct the measurements from this signal

b = A.times(x)

Let us visualize the original and reconstructed signals

plot_signals(x0, x)

Let us visualize the original and reconstructed measurements

plot_measurments(b0, b)

Sparse Recovery using SPGL1

import cr.sparse.cvx.spgl1 as crspgl1
options = crspgl1.SPGL1Options()
sol = crspgl1.solve_bp_jit(A, b0, options=options)
problems.analyze_solution(prob, sol)

m: 600, n: 2560
b_norm: original: 4.441 reconstruction: 4.441 SNR: 93.60 dB
x_norm: original: 4.472 reconstruction: 4.472 SNR: 89.45 dB
Sparsity: original: 20, reconstructed: 20, overlap: 20, ratio: 1.000
Iterations: 36 n_times: 39, n_trans: 37

The estimated sparse signal

x = sol.x

Let us reconstruct the measurements from this signal

b = A.times(x)

Let us visualize the original and reconstructed signals

plot_signals(x0, x)

Let us visualize the original and reconstructed measurements

plot_measurments(b0, b)