OMP Algorithm¶
- class fastmat.algorithms.OMP¶
Bases:
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
- fmatAfm.Matrix
the system matrix
- arrBnp.ndarray
the measurement vector
- numMaxStepsint
the desired sparsity order
- Returns
- np.ndarray
solution array
- __init__(*args, **kwargs)¶
- arrA¶
- arrB¶
- arrC¶
- arrResidual¶
- arrSupport¶
- arrX¶
- arrXtmp¶
- fmatA¶
- fmatC¶
- matPinv¶
- newCols¶
- newIndex¶
- numL¶
- numM¶
- numMaxSteps¶
- numN¶
- numStep¶
- v2¶
- v2n¶
- v2y¶