After the little tour we've done in the previous notebooks, about the [creation of Faust objects](#creation_links), their [manipulation](#manip_links) and their creation using [factorization](#facto_links) algorithms, we shall see in this fourth notebook how the FAµST API can be deployed seamlessly in algorithms.
Our example, already alluded in the second notebook, we'll be the Orthogonal Matching Pursuit algorithm (OMP).
This algorithm intervenes in the dictionary learning problem. Of course, I will not treat the theory behind but I assume the reader is already familiar to.
There is not so much to say so let's go straight to the code example.
%% Cell type:markdown id: tags:
### 1. The OMP Algorithm Implementation
%% Cell type:code id: tags:
``` python
frompyfaustimport*
fromnumpyimportzeros,copy,argmax
defomp(y,D,niter):
m,n=D.shape
x=zeros((n,1))
supp=[]
res=copy(y)
i=0
K=min(m,n)
while (len(supp)<Kandi<niter):
j=argmax(abs(D.T*res))
supp+=[j]
x[supp,:]=pinv(D[:,supp])*y
res=y-D[:,supp]*x[supp,:]
i+=1
returnx
```
%% Cell type:markdown id: tags:
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 sure that a Faust can be seen as a numpy array (or rather as a ``numpy.matrix``) hence this code is in fact totally compatible with the two APIs: the D function argument, which is the dictionary, can be indifferently a ``pyfaust.Faust`` object or a ``numpy.matrix`` object.
Next we'll test this implementation in the both cases. But first, let us define a test case.
%% Cell type:markdown id: tags:
### 2. The Test Case Dictionary
%% Cell type:markdown id: tags:
For convenience, we shall set up a dictionary which garantees uniqueness of inefficiently 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 or Dirac matrix and $H_n$ the 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 insured is defined by the following inequality:
$ \| x \|_0 < (1 + 1/\mu)/2 $ where $\mu$ denotes the coherence of the dictionary and in case of our specially crafted dictionary $\mu = {1 \over \sqrt n}$.
So let's construct the Faust of D, compute y for a unique sparse x and test our OMP implementation to find out if we effectively retrieve the unique x as we should according to this theorem.
Note that, for a better view and understanding you might consult this document: [Rémi Gribonval HDR](http://videos.rennes.inria.fr/hdr/gribonval/HDR2007.pdf).
%% Cell type:code id: tags:
``` python
frompyfaustimportFaustFactoryasFF
fromnumpyimportlog2
n=128
log2n=log2(n)
FD=FF.eye(n).concatenate(FF.wht(log2n),axis=1)
FD=hstack((FF.eye(n),FF.wht(log2n)))
D=FD.todense()
print(D)
```
%% Cell type:markdown id: tags:
Now that we have our dictionary both defined as a Faust (FD) and as a matrix (D), let's construct our reference sparse vector x, we'll call it $x_0$.
%% Cell type:code id: tags:
``` python
fromnumpyimportzeros,count_nonzero
fromnumpy.randomimportrandn
fromrandomimportrandint
frommathimportfloor,sqrt
fromitertoolsimportcombinations
fromscipy.specialimportbinom
x0=zeros((FD.shape[1],1))
nnz=max(1,randint(5,floor(.5*(1+sqrt(n)))))# min 0-norm is 5
c=combinations(list(range(0,n)),nnz)# all possible combinations of nnz indices
nc=binom(n,nnz)
# pick random nnz indices
nonzero_inds=[next(c)foriinrange(0,randint(0,min(nc,20)))][-1]# limiting to 20 iterations to accelerate
It remains to compute $y$, normalizing FD and D before.
%% Cell type:code id: tags:
``` python
FD=FD.normalize()
D=FD.todense()
y=FD*x0
```
%% Cell type:markdown id: tags:
Our test case is complete, we are fully prepared to run the OMP algorithm using a well-defined dictionary as a Faust or as numpy array, this should retrieve our $x_0$ from the vector y. Let's try!
%% Cell type:markdown id: tags:
### 3. Running the Algorithm
%% Cell type:code id: tags:
``` python
x=omp(y,FD,nnz)
fromnumpyimportallclose
assert(allclose(x,x0))
print("We succeeded to retrieve x0, OMP works!")
```
%% Cell type:markdown id: tags:
We tested OMP on a Faust, go ahead and verify what I was aiming at in the first part of the notebook: is this OMP implementation really working identically on a Faust and a ``numpy.matrix``?
%% Cell type:code id: tags:
``` python
x=omp(y,D,nnz)
assert(allclose(x,x0))
print("We succeeded to retrieve x0, OMP works!")
```
%% Cell type:markdown id: tags:
As a complementary test, let's verify that the two runs of ``omp()`` on FD and D give the same results even if the vector to retrieve is not sparse. Here for instance, $ \| x_1 \|_0 = 98 $
We can conclude that the algorithm is indeed available to both numpy and Faust worlds, and we can imagine surely that other algorithms are reachable through the FAµST API. That's anyway in that purpose that the FAµST library will be extended in the future (starting from the latest version which is 2.5.1).
%% Cell type:markdown id: tags:
***The fourth 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!