After the little tour we've done in the previous notebooks, about the creation of Faust objects, and their manipulation, we shall see in this third notebook how the FAµST API can be deployed seamlessly in algorithms. Our example, already mentioned in the second notebook, will be the Orthogonal Matching Pursuit algorithm (OMP).

This algorithm comes up in the dictionary learning problem. Assuming that the reader is already familiar with this problem we will not treat the theory behind. There is not so much to say so let's go straight to the code example.

In [1]:

```
from pyfaust import *
from numpy import zeros, copy, argmax
def tomp(y, D, niter):
nrows, ncols = D.shape
x = zeros((ncols,1))
supp = []
res = copy(y)
i = 0
K = min(nrows, ncols)
while (len(supp) < K and i < niter):
j = argmax(abs(D.T@res))
supp += [j]
x[supp,:] = pinv(D[:,supp])@y
res = y-D[:,supp]@x[supp,:]
i += 1
return x
```

The most important point to notice in this code is that except the import part in the header, all the code seems to be a natural numpy implementation of OMP.

This is in fact the core philosophy of the FAµST API, as explained in previous notebooks and also in the API documentation, we made sure that a Faust can be seen as a numpy array. Hence this code is in fact totally compatible with the two APIs: the function argument D, which is the dictionary, can be indifferently a `pyfaust.Faust`

object or a numpy array.

A secondary point is that this implementation is more like a toy concept (as indicated by the "t" in the function name). A more advanced and optimized version is introduced in the last part of this notebook. In particular this version allows to define the algorithm stopping criterion according to the error tolerance the user wants.

Next we will test this implementation in both cases. But first, let us define a test case.

For convenience, we shall set up a dictionary which guarantees uniqueness of sufficiently sparse representations. The dictionary is the concatenation of an identity matrix and a Hadamard matrix, and because we work with Faust objects, this concatenation will be a Faust object.

Below is the block matrix of our dictionary:

$
D =
\left[
\begin{array}{c|c}
I_n & H_n \\
\end{array}
\right]
$

$I_n$ is the identity matrix and $H_n$ the orthonormal Hadamard matrix, with n being a power of two.

The condition on which the uniqueness of the sparse representation $x$ of a vector $y$ is ensured is defined by the following inequality:

$ \| x \|_0 < (1 + 1/\mu)/2 $ where $\mu$ denotes the coherence of the dictionary and in the case of our specially crafted dictionary $\mu = {1 \over \sqrt n}$.

So let's construct the Faust of D, compute y for a sparse enough x and test our OMP implementation to find out if we effectively retrieve this unique x as we should according to this theorem.

Note that, for a better view and understanding you might consult this article [1].

In [2]:

```
from pyfaust import eye, wht
n = 128
FD = hstack((eye(n), wht(n)))
D = FD.toarray()
print(D)
```

In [3]:

```
from numpy import zeros, count_nonzero
from numpy.random import randn, permutation as randperm
from math import floor, sqrt
x0 = zeros((2*n, 1)) # NB: FD.shape[1] == 2*n
nnz = floor(.5*(1+sqrt(n)))
nonzero_inds = randperm(2*n)[:nnz]
# we got nnz indices, now build the vector x0
x0[nonzero_inds,0] = randn()
print("l0 norm of x0: ", count_nonzero(x0))
```

l0 norm of x0: 6

It remains to compute $y$.

In [4]:

```
y = FD@x0
```

In [5]:

```
x = tomp(y, FD, nnz)
from numpy import allclose
assert(allclose(x,x0))
print("We succeeded to retrieve x0, OMP works!")
```

We succeeded to retrieve x0, OMP works!

In [6]:

```
x = tomp(y, D, nnz)
assert(allclose(x,x0))
print("We succeeded to retrieve x0, OMP works!")
```

We succeeded to retrieve x0, OMP works!

Speaking of the OMP algorithm and the possibility to implement other optimization algorithms with FAµST, it would be a pity not to mention that the library is delivered with another implementation of OMP.

This implementation is actually an optimized version which takes advantage of the Cholesky factorization to simplify the least-square problem to solve at each iteration.
This algorithm is implemented into the `tools`

module of pyfaust.

In [7]:

```
from pyfaust.tools import omp
help(omp)
```

This implementation is integrated into pyfaust as a tool for the Brain Source Localization (BSL) demo which is documented here.

To show off a little, let's run this demo.

**Hint**: the demo takes a few minutes, it you find it too annoying to wait all this time you can jump to the next cell and render the precomputed data by changing the boolean `use_precomputed_data`

to `True`

.

In [8]:

```
%%capture
# It will take some time (sorry), many Faust-s are compared to the original MEG matrix
# However running this cell is not mandatory, it'll be skipped if you switch use_precomputed_data to True in the next code cell
from pyfaust.demo import bsl
bsl.run()
```

In [9]:

```
%%capture --no-display
%matplotlib inline
from pyfaust.demo import bsl
bsl.fig_time_cmp(use_precomputed_data=False)
```

`omp()`

on FD and D we constructed before for the toy OMP give the same results even if the vector to retrieve is not very sparse. Here for instance, $ \| x_1 \|_0 = 98 $

In [10]:

```
nnz = 98
x1 = zeros((2*n, 1))
nnz_inds = randperm(2*n)[:nnz]
x1[nnz_inds, 0] = [randn() for i in range(0,nnz)]
y = FD@x1
x2 = omp(y, D, maxiter=nnz)
x3 = omp(y, FD, maxiter=nnz)
# verify if the solutions differ
print("Are x2 and x3 solutions almost equal?", norm(x2-x3)/norm(x3) < 10**-12)
print("Is x1 retrieved into x2?", allclose(x1,x2))
print("Is x1 retrieved into x3?", allclose(x1,x3))
```

/home/hinria/tmp_notebook_bsl/pyfaust_venv/lib64/python3.10/site-packages/pyfaust/__init__.py:1165: UserWarning: The numpy array to multiply is not F_CONTIGUOUS, costly conversion on the fly. Please use np.asfortranarray by yourself. w = lambda: warnings.warn("The numpy array to multiply is not "

Are x2 and x3 solutions almost equal? True Is x1 retrieved into x2? False Is x1 retrieved into x3? False

As expected, we didn't retrieve our starting x1 (the reason is the condition already discussed in 2.). However let's mention that here again (like it was with the toy OMP) it works the same with the Faust API or with numpy.

Finally, let's check the computation time for applying our dictionary to a vector both for the numpy and Faust versions. Although, in order to avoid major differences in results calculated on distinct computer configurations the comparison is performed on a larger dimension than before.

In [11]:

```
from numpy.random import randn
def dirac_hadamard(n):
FD = hstack((eye(n), wht(n)))
return FD
n = 1024
FD = dirac_hadamard(n)
D = FD.toarray()
x = randn(2*n,1)
%timeit D@x
%timeit FD@x
```

This is essentially because the speedup offered by Faust appears rather for higher matrix dimensions.

Let us illustrate the speedup more generally by repeating the experiment for various dimensions n.

In [12]:

```
%matplotlib inline
from numpy.random import randn
from numpy import array
from pyfaust import eye, wht
from time import process_time
from matplotlib.pyplot import (plot, show,
legend, semilogy,
xlabel, ylabel,
title, scatter)
d_times = []
fd_times = []
dims = []
num_muls = 100
for n in [2**i for i in range(8,13)]:
FD = dirac_hadamard(n)
D = FD.toarray()
x = randn(2*n,1)
dims += [n]
t = process_time()
[D@x for i in range(0,num_muls)]
d_times += [process_time()-t]
t = process_time()
[FD@x for i in range(0,num_muls)]
fd_times += [process_time()-t]
p1 = semilogy(dims, fd_times, label="FD*x times")
p2 = semilogy(dims, d_times, label="D*x times")
xlabel("dimension")
ylabel("time")
legend()
title("D*x and FD*x time comparison")
show()
```

In [13]:

```
scatter(dims, array(d_times)/array(fd_times), label="speedup")
title("Speedup factor of FD*x relatively to D*x")
show()
```

** The third notebook is ending here**, I hope you'll be interested in trying yourself to write another algorithm with the FAµST API and maybe discovering any current limitation. Don't hesitate to contact us in that case, we'll appreciate any feedback!

**Note**: this notebook was executed using the following pyfaust version:

In [14]:

```
import pyfaust
pyfaust.version()
```

Out[14]:

'3.27.2'