Greedy Sparse Recovery¶
Algorithm versions
Several algorithms are available in multiple versions.
The library allows a dictionary or a sensing process to be represented as either:
A matrix of real/complex values
A linear operator with fast implementation (see
cr.sparse.lopmodule)
E.g. a partial sensing sensing matrix can be implemented more efficiently using a linear operator consisting of applying the fourier transform followed by selecting a subset of fourier measurements.
A prefix
matrix_indicates an implementation which accepts matrices as dictionaries or compressive sensors.A prefix
operator_indicates an implementation which accepts linear operators described incr.sparse.lopmodule as dictionaries or compressive sensors.A suffix
_jitmeans that it is the JIT (Just In Time ) compiled version of the original implementation.A suffix
_multimeans that it is the version of implementation which can process multiple signals/measurement vectors simultaneously. The recovery problem \(y = \Phi x + e\) is extended to \(Y = \Phi X + E\) such that:Each column of \(Y\) represents one signal/measurement vector
Each column of \(X\) represents one representation vector to be recovered
Each column of \(E\) representation corresponding measurement error/noise.
Conditions on dictionaries/sensing matrices
Different algorithms have different requirements on the dictionaries or sensing matrices:
Some algorithms accept overcomplete dictionaries/sensing matrices with unit norm columns
Some algorithms accept overcomplete dictionaries/sensing matrices with orthogonal rows
Some algorithms accept any overcomplete dictionary
Basic Matching Pursuit Based Algorithms¶
Matching Pursuit
|
Solves the sparse recovery problem \(y = \Phi x + e\) using matching pursuit algorithm |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using matching pursuit algorithm |
Orthogonal Matching Pursuit
|
Solves the recovery/approximation problem \(y = \Phi x + e\) using Orthogonal Matching Pursuit |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Orthogonal Matching Pursuit for linear operators |
|
Solves the recovery/approximation problem \(y = \Phi x + e\) using Orthogonal Matching Pursuit |
|
Solves the recovery/approximation problem \(y = \Phi x + e\) using Orthogonal Matching Pursuit |
|
Solves the MMV recovery/approximation problem \(Y = \Phi X + E\) using Orthogonal Matching Pursuit |
Compressive Sensing Matching Pursuit (CSMP) Algorithms
Compressive Sampling Matching Pursuit
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Compressive Sampling Matching Pursuit for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Compressive Sampling Matching Pursuit for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Compressive Sampling Matching Pursuit for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Compressive Sampling Matching Pursuit for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Compressive Sampling Matching Pursuit for linear operators |
Subspace Pursuit
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Subspace Pursuit for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Subspace Pursuit for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Subspace Pursuit for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Subspace Pursuit for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Subspace Pursuit for linear operators |
Hard Thresholding Based Algorithms¶
Iterative Hard Thresholding
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Iterative Hard Thresholding for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Iterative Hard Thresholding for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Iterative Hard Thresholding for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Iterative Hard Thresholding for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Iterative Hard Thresholding for linear operators |
Hard Thresholding Pursuit
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Hard Thresholding Pursuit for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Hard Thresholding Pursuit for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Hard Thresholding Pursuit for matrices |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Hard Thresholding Pursuit for linear operators |
|
Solves the sparse recovery problem \(y = \Phi x + e\) using Hard Thresholding Pursuit for linear operators |
Data Types¶
Represents the solution of a sparse recovery problem |
Utilities¶
|
Returns the index of entry with highest magnitude |
|
Incrementally updates the Cholesky factorization \(G = L L^T\) where \(G = \Phi^T \Phi\) |
Using the greedy algorithms¶
These algorithms solve the inverse problem \(y = \Phi x + e\) where \(\Phi\) and \(y\) are known and \(x\) is desired with \(\Phi\) being an overcomplete dictionary or a sensing matrix.
For sparse approximation problems, we require the following to invoke any of these algorithms:
A sparsifying dictionary \(\Phi\).
A signal \(y\) which is expected to have a sparse or compressible representation \(x\) in \(\Phi\).
For sparse recovery problems, we require the following to invoke any of these algorithms:
A sensing matrix \(\Phi\) with suitable RIP or other properties.
A measurement vector \(y\) generated by applying \(\Phi\) to a sparse signal \(x\)
A synthetic example
Build a Gaussian dictionary/sensing matrix:
from jax import random
import cr.sparse.dict as crdict
M = 128
N = 256
key = random.PRNGKey(0)
Phi = crdict.gaussian_mtx(key, M,N)
Build a K-sparse signal with Gaussian non-zero entries:
import cr.sparse.data as crdata
import jax.numpy as jnp
K = 16
key, subkey = random.split(key)
x, omega = crdata.sparse_normal_representations(key, N, K, 1)
x = jnp.squeeze(x)
Build the measurement vector:
y = Phi @ x
We have built the necessary inputs for a sparse recovery problem. It is time to run the solver.
Import a sparse recovery solver:
from cr.sparse.pursuit import cosamp
Solve the recovery problem:
solution = cosamp.matrix_solve(Phi, y, K)
You can choose any other solver.
The support for the non-zero entries in the solution is given by solution.I and
the values for non-zero entries are given by solution.x_I. You can build the
sparse representation as follows:
from cr.sparse import build_signal_from_indices_and_values
x_hat = build_signal_from_indices_and_values(N, solution.I, solution.x_I)
Finally, you can use the utility to evaluate the quality of reconstruction:
from cr.sparse.ef import RecoveryPerformance
rp = RecoveryPerformance(Phi, y, x, x_hat)
rp.print()
This would output something like:
M: 128, N: 256, K: 16
x_norm: 3.817, y_norm: 3.922
x_hat_norm: 3.817, h_norm: 1.55e-06, r_norm: 1.72e-06
recovery_snr: 127.83 dB, measurement_snr: 127.16 dB
T0: [ 27 63 79 85 88 111 112 124 131 137 160 200 230 234 235 250]
R0: [ 27 63 79 85 88 111 112 124 131 137 160 200 230 234 235 250]
Overlap: [ 27 63 79 85 88 111 112 124 131 137 160 200 230 234 235 250], Correct: 16
success: True
