OMP Algorithm

class fastmat.algorithms.OMP

Bases: fastmat.algorithms.Algorithm.Algorithm

Orthogonal Matching Pursuit

Definition and Interface: For a given matrix \(A \in \mathbb{C}^{m \times N}\) with \(m \ll N\) and a vector \(b \in \mathbb{C}^m\) we approximately solve

\[\min\limits_{ x \in \mathbb{C}^N} \Vert x \Vert_0 \quad\mathrm{s.t.}\quad A \cdot x = x.\]

If it holds that \(b = A \cdot x_0\) for some \(k\)-sparse \(x_0\) and \(k\) is low enough, we can recover \(x_0\) via OMP [1].

This type of problem as the one described above occurs in Compressed Sensing and Sparse Signal Recovery, where signals are approximated by sparse representations.

>>> # import the packages
>>> import numpy.linalg as npl
>>> import numpy as np
>>> import fastmat as fm
>>> import fastmat.algorithms as fma
>>> # define the dimensions
>>> # and the sparsity
>>> n, k = 512, 3
>>> # define the sampling positions
>>> t = np.linspace(0, 20 * np.pi, n)
>>> # construct the convolution matrix
>>> c = np.cos(2 * t)
>>> C = fm.Circulant(c)
>>> # create the ground truth
>>> x = np.zeros(n)
>>> x[np.random.choice(range(n), k, replace=0)] = 1
>>> b = C * x
>>> # reconstruct it
>>> omp = fma.OMP(C, numMaxSteps=100)
>>> y = omp.process(b)
>>> # test if they are close in the
>>> # domain of C
>>> print(npl.norm(C * y - b))

We describe a sparse deconvolution problem, where the signal is in \(\mathbb{R}^{512}\) and consists of \(3\) windowed cosine pulses of the form \(c\) with circulant displacement. Then we take the convolution and try to recover the location of the pulses using the OMP algorithm.

Note

The algorithm exploits two mathematical shortcuts. First it obviously uses the fast transform of the involved system matrix during the correlation step and second it uses a method to calculate the pseudo inverse after a rank-\(1\) update of the matrix.

Todo

  • optimize einsum-stuff
[1]S. G. Mallat, Z. Zhang, “Matching pursuits with time-frequency dictionaries”, IEEE Transactions on Signal Processing, vol. 41, no. 12, pp. 3397-3415, Dec 1993
Parameters:
fmatA : fm.Matrix

the system matrix

arrB : np.ndarray

the measurement vector

numMaxSteps : int

the desired sparsity order

Returns:
np.ndarray

solution array

__init__(*args, **kwargs)

Initialize self. See help(type(self)) for accurate signature.

arrA
arrB
arrC
arrResidual
arrSupport
arrX
arrXtmp
fmatA
fmatC
matPinv
newCols
newIndex
numL
numM
numMaxSteps
numN
numStep
v2
v2n
v2y