0
$\begingroup$

I have recently attempted to implement a simple version of the Matching Pursuit algorithm, following the description detailed on Wikipedia. Although the algorithm itself seems dead simple, my implementation produces very suboptimal results in some cases (i.e. the inner product between the represented signal and the sparse reconstruction is very small).

Here is the Python function I have written, x is the input signal and D is the dictionary matrix with normalized rows.

import numpy as np
import numpy.linalg as la

def matching_pursuit(x: np.ndarray,
                     D: np.ndarray,
                     eps_min: np.float64 = 1e-3,
                     iter_max: int = 1000) -> np.ndarray:

    s = np.zeros(D.shape[0], dtype=np.float64)

    r = x.astype(np.float64)
    for _ in range(iter_max):
        dot = np.dot(D, r)
        i = np.argmax(np.abs(dot))

        s[i] = dot[i]
        r -= s[i] * D[i, :]

        if la.norm(r) < eps_min:
            break

    return s

As an example, for a randomly generated matrix D equal to:

[[0.09365858 0.65561007 0.74926865]
 [0.24806947 0.62017367 0.74420841]
 [0.08804509 0.70436073 0.70436073]
 [0.54545455 0.54545455 0.63636364]
 [0.62469505 0.46852129 0.62469505]
 [0.68100522 0.59587957 0.42562827]
 [0.64465837 0.72524067 0.24174689]
 [0.24077171 0.48154341 0.84270097]
 [0.58834841 0.19611614 0.78446454]
 [0.80178373 0.26726124 0.53452248]
 [0.45584231 0.56980288 0.68376346]
 [0.36214298 0.81482171 0.45267873]
 [0.19802951 0.69310328 0.69310328]
 [0.74535599 0.59628479 0.2981424 ]
 [0.76249285 0.45749571 0.45749571]
 [0.19802951 0.69310328 0.69310328]
 [0.46829291 0.74926865 0.46829291]
 [0.70014004 0.14002801 0.70014004]
 [0.76626103 0.28734789 0.57469577]
 [0.93704257 0.31234752 0.15617376]]

and an input signal x = [1, 2, 3] the function returns

[ 0.          0.          0.          0.          0.          0.
  0.00119178 -0.00095663 -0.00167583  0.          0.          0.
  0.          0.          0.          0.          0.          0.
  0.          0.        ]

which results, when multiplied with the dictionary matrix again, in x = [-4.48011162e-04, 7.50118419e-05, -1.83267255e-03].

Now what I'm wondering is whether I have made some mistake in my implementation or whether this is an expected outcome due to the nature of the algorithm and/or properties of these specific sample inputs.

$\endgroup$
  • $\begingroup$ Where did your data come from? did you just randomize it? $\endgroup$ – havakok Dec 26 '18 at 13:00
  • $\begingroup$ Yes, I just used np.random.randint to create D. $\endgroup$ – Peter Dec 26 '18 at 13:01
  • $\begingroup$ It is difficult to folllow your program without having the main part. However, the approximation of $x$ is equal to $\sum_i s_i D_i$. It does not seem to be provided by the program $\endgroup$ – Damien Dec 26 '18 at 14:15
1
$\begingroup$

In MP it is very important which dictionary is used. For example, if the atoms are orthogonal The problem will converge faster.

Another interesting example, where two entries of the dictionary are the same, you can wither find both corresponding values in the sparse signal as 1 and zeros or as 0.5 and 0.5. You can find any number of representation for such an example. Which one will the algorithm converge to?

In general, the picking of the dictionary $D$ is a very important point in the MP algorithm. I recommend you take a look in the following link as well as researching a bit more yourself. The Fourier matrix, for example, is a common choice for the dictionary.

To answer your question, while I could not find a problem with your implementation of the algorithm, the data you are using might be problematic.

$\endgroup$
  • $\begingroup$ Okay, I had just assumed that even with a random dictionary the algorithm would converge to an okay-ish solution but I'll take a look at the article you posted, thanks. $\endgroup$ – Peter Dec 26 '18 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.