Table Of Contents

Previous topic

PyMP Tutorial 2. Audio compression

Next topic

PyMP Tutorial 4. Similarity Measures

This Page

PyMP Tutorial 3. Greedy decomposition variants

There are many variants of the MP algorithm. Only a few of them are directly implemented in PyMP. In the first tutorial, we saw some variants of the selection step (e.g. RSSMP). Here we’re interested in the update step. In order to use this variants, we use the mp.greedy() function of the mp module.

Here’s a example of use with a simple multiscale MDCT dictionary and the standard MP algorithm:

>>> import numpy as np
>>> from PyMP.mdct import Dico
>>> from PyMP import mp, Signal
>>> sig = Signal('data/ClocheB.wav', mono=True)  # Load Signal
>>> sig.crop(0, 4.0 * sig.fs)     # Keep only 4 seconds
>>> # atom of scales 8, 64 and 512 ms
>>> scales = [(s * sig.fs / 1000) for s in (8, 64, 512)]
>>> sig.pad(scales[-1])
>>> mp_dico = Dico(scales)
>>> # Launching decomposition, stops either at 20 dB of SRR or 2000 iterations
>>> mp_approx, mp_decay = mp.greedy(sig, mp_dico, 20, 2000, pad=False, update='mp')
>>> mp_approx.atom_number

565

Note

This is equivalent to using the mp.mp() method. However, the additionnal update keyword argument allows to change the nature of the update step.

Orthogonal Matching Pursuit

The principle of OMP is to ensure that at each iteration, chosen atom weights are the best possible ones (in a least squares sense). An orthogobnal projection of the residual onto the subspace spanned by all previously selected atoms is thus performed. Original paper is [1]. The keyword to use OMP is omp.

Warning

OMP is extremely resource consuming. You should not try to use it on signals with lengths greater than a few thousand samples. Try the local version instead or GP.

Gradient Pursuit

Gradient Pursuits have been introduced by Blumensath and Davis in [2]. The idea is to replace the costly orthogonal projection by a much more affordable directionnal update step.

Warning

Orthogonal and directionnal pursuits implementation are not compatible (yet) with FullDico is implemented.

Local versions of OMP/GP

Mailhé et al [3] have presented variants of OMP and GP that are much more efficient since they only perform the update in a neighborhood of the last selected atom. For the kind of dictionaries that are implemented in PyMP (i.e. time-localized MDCT and Wavelets), this is a great news. The keyword arguments to use the local OMP is locomp and local GP is locgp

Cyclic Matching Pursuit

Coming Soon!

example

The following example build a k-sparse signal and compares decompositions using MP, locOMP and locGP:

import numpy as np
from PyMP import Signal, mp
from PyMP.mdct.dico import Dico, LODico
from PyMP.mdct.atom import Atom
print "Running MP, OMP and local versions on synthetic k-sparse"
scales = [16, 64, 256]
dico = Dico(scales)
M = len(scales)
L = 256 * 4
k = 0.2*L
# create a k-sparse signal
sp_vec = np.zeros(M*L,)
from PyMP.tools import mdct
random_indexes = np.arange(M*L)
np.random.shuffle(random_indexes)
random_weights = np.random.randn(M*L)
sp_vec[random_indexes[0:k]] = random_weights[0:k]
sparse_data = np.zeros(L,)
for m in range(M):
    sparse_data += mdct.imdct(sp_vec[m*L:(m+1)*L], scales[m])
signal_original = Signal(sparse_data, Fs=8000, mono=True, normalize=False)
signal_original.data += 0.01 * np.random.random(signal_original.length,)
n_atoms = k 
signal_original.pad(dico.get_pad())
app_2, dec2 = mp.greedy(signal_original, dico, 100,
                        n_atoms, debug=0, pad=False, update='locgp')
app_1, dec1 = mp.greedy(signal_original, dico, 100,
                        n_atoms, debug=0, pad=False, update='mp')
app_3, dec3 = mp.greedy(signal_original, dico, 100,
                        n_atoms, debug=0, pad=False, update='locomp')

Then you can plot the decays

(Source code, png, hires.png, pdf)

_images/OMP_example.png

Bibliography

1 : Pati, Y. C. C., Rezaiifar, R., & Krishnaprasad, P. S. S. Orthogonal Matching Pursuit: Recursive Function Approximation with Applications to Wavelet Decomposition (pdfomp).
Proceedings of the 27 th Annual Asilomar Conference on Signals, Systems, and Computers (pp. 40–44).
2 : Blumensath, T., & Davies, M. . (2008). Gradient Pursuits.
IEEE Transactions on Signal Processing, 56(6), 2370–2382 (pdfgp).
3 : Mailhé, B., Gribonval, R., Vandergheynst, P., & Bimbot, F. (2011). Fast orthogonal sparse approximation algorithms over local dictionaries.
Signal Processing vol 91(12) (2822-2835) (pdfloc).