Exceptional GitLab maintenance is scheduled for June 24, 2025. The service will be unavailable between 8:00 am and 9:30 am. Please do not work on the platform until an announcement indicates that maintenance is complete.
"<a name=\"[1]\">[1]</a> [Tropp, J. A. (2004). Greed is Good: Algorithmic Results for Sparse Approximation. IEEE Transactions on Information Theory, 50(10), 2231–2242](http://doi.org/10.1109/TIT.2004.834793)."
"<a name=\"[1]\">[1]</a> [Tropp, J. A. (2004). Greed is Good: Algorithmic Results for Sparse Approximation. IEEE Transactions on Information Theory, 50(10), 2231–2242](http://doi.org/10.1109/TIT.2004.834793)."
]
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Note**: this notebook was executed using the following pyfaust version:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pyfaust\n",
"pyfaust.version()"
]
}
}
],
],
"metadata": {
"metadata": {
...
@@ -421,7 +438,7 @@
...
@@ -421,7 +438,7 @@
"name": "python",
"name": "python",
"nbconvert_exporter": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"pygments_lexer": "ipython3",
"version": "3.7.6"
"version": "3.7.10"
}
}
},
},
"nbformat": 4,
"nbformat": 4,
...
...
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
# Using the FAµST API in Algorithms
# Using the FAµST API in Algorithms
After the little tour we've done in the previous notebooks, about the [creation of Faust objects](#creation_links), and their [manipulation](#manip_links), we shall see in this third notebook how the FAµST API can be deployed seamlessly in algorithms.
After the little tour we've done in the previous notebooks, about the [creation of Faust objects](#creation_links), and their [manipulation](#manip_links), we shall see in this third notebook how the FAµST API can be deployed seamlessly in algorithms.
Our example, already alluded in the [second notebook](#manip_links), will be the Orthogonal Matching Pursuit algorithm (OMP).
Our example, already alluded in the [second notebook](#manip_links), will 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 with it.
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 with it.
There is not so much to say so let's go straight to the code example.
There is not so much to say so let's go straight to the code example.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 1. The Toy OMP Algorithm Implementation
### 1. The Toy OMP Algorithm Implementation
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
frompyfaustimport*
frompyfaustimport*
fromnumpyimportzeros,copy,argmax
fromnumpyimportzeros,copy,argmax
deftomp(y,D,niter):
deftomp(y,D,niter):
nrows,ncols=D.shape
nrows,ncols=D.shape
x=zeros((ncols,1))
x=zeros((ncols,1))
supp=[]
supp=[]
res=copy(y)
res=copy(y)
i=0
i=0
K=min(nrows,ncols)
K=min(nrows,ncols)
while (len(supp)<Kandi<niter):
while (len(supp)<Kandi<niter):
j=argmax(abs(D.T@res))
j=argmax(abs(D.T@res))
supp+=[j]
supp+=[j]
x[supp,:]=pinv(D[:,supp])@y
x[supp,:]=pinv(D[:,supp])@y
res=y-D[:,supp]@x[supp,:]
res=y-D[:,supp]@x[supp,:]
i+=1
i+=1
returnx
returnx
```
```
%% Cell type:markdown id: tags:
%% 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.
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 (or rather as a ``numpy.matrix``) 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.matrix`` object.
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 (or rather as a ``numpy.matrix``) 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.matrix`` object.
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](#4.-An-OMP-Cholesky-Implementation) and in particular allows to define the algorithm stopping criterion according to the error tolerance the user wants.
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](#4.-An-OMP-Cholesky-Implementation) and in particular 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.
Next we will test this implementation in both cases. But first, let us define a test case.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 2. The Test Case Dictionary
### 2. The Test Case Dictionary
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
For convenience, we shall set up a dictionary which guarantees uniqueness of sufficiently sparse representations.
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.
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:
Below is the block matrix of our dictionary:
$
$
D =
D =
\left[
\left[
\begin{array}{c|c}
\begin{array}{c|c}
I_n & H_n \\
I_n & H_n \\
\end{array}
\end{array}
\right]
\right]
$
$
$I_n$ is the identity or Dirac matrix and $H_n$ the orthonormal Hadamard matrix, with n being a power of two.
$I_n$ is the identity or Dirac 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:
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 case of our specially crafted dictionary $\mu = {1 \over \sqrt n}$.
$ \| 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 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.
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]](#[1]).
Note that, for a better view and understanding you might consult this article [[1]](#[1]).
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
frompyfaustimporteye,wht
frompyfaustimporteye,wht
n=128
n=128
FD=hstack((eye(n),wht(n)))
FD=hstack((eye(n),wht(n)))
D=FD.toarray()
D=FD.toarray()
print(D)
print(D)
```
```
%% Cell type:markdown id: tags:
%% 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$.
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:
%% Cell type:code id: tags:
``` python
``` python
fromnumpyimportzeros,count_nonzero
fromnumpyimportzeros,count_nonzero
fromnumpy.randomimportrandn,permutationasrandperm
fromnumpy.randomimportrandn,permutationasrandperm
frommathimportfloor,sqrt
frommathimportfloor,sqrt
x0=zeros((2*n,1))# NB: FD.shape[1] == 2*n
x0=zeros((2*n,1))# NB: FD.shape[1] == 2*n
nnz=floor(.5*(1+sqrt(n)))
nnz=floor(.5*(1+sqrt(n)))
nonzero_inds=randperm(2*n)[:nnz]
nonzero_inds=randperm(2*n)[:nnz]
# we got nnz indices, now build the vector x0
# we got nnz indices, now build the vector x0
x0[nonzero_inds,0]=randn()
x0[nonzero_inds,0]=randn()
print("l0 norm of x0: ",count_nonzero(x0))
print("l0 norm of x0: ",count_nonzero(x0))
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
It remains to compute $y$.
It remains to compute $y$.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
y=FD@x0
y=FD@x0
```
```
%% Cell type:markdown id: tags:
%% 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!
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:
%% Cell type:markdown id: tags:
### 3. Running the Algorithm
### 3. Running the Algorithm
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
x=tomp(y,FD,nnz)
x=tomp(y,FD,nnz)
fromnumpyimportallclose
fromnumpyimportallclose
assert(allclose(x,x0))
assert(allclose(x,x0))
print("We succeeded to retrieve x0, OMP works!")
print("We succeeded to retrieve x0, OMP works!")
```
```
%% Cell type:markdown id: tags:
%% 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``?
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:
%% Cell type:code id: tags:
``` python
``` python
x=tomp(y,D,nnz)
x=tomp(y,D,nnz)
assert(allclose(x,x0))
assert(allclose(x,x0))
print("We succeeded to retrieve x0, OMP works!")
print("We succeeded to retrieve x0, OMP works!")
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
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 if needed in the future.
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 if needed in the future.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 4. An OMP-Cholesky Implementation
### 4. An OMP-Cholesky Implementation
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.
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 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.
This algorithm is implemented into the ``tools`` module of pyfaust.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
frompyfaust.toolsimportomp
frompyfaust.toolsimportomp
help(omp)
help(omp)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
This implementation is integrated into pyfaust as a tool for the Brain Source Localization (BSL) demo which is documented [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1demo_1_1bsl.html).
This implementation is integrated into pyfaust as a tool for the Brain Source Localization (BSL) demo which is documented [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1demo_1_1bsl.html).
To show off a little, let's run this demo.
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``.
**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``.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%%capture
%%capture
# It will take some time (sorry), many Faust-s are compared to the original MEG matrix
# 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
# However running this cell is not mandatory, it'll be skipped if you switch use_precomputed_data to True in the next code cell
frompyfaust.demoimportbsl
frompyfaust.demoimportbsl
bsl.run()
bsl.run()
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%%capture--no-display
%%capture--no-display
%matplotlibinline
%matplotlibinline
frompyfaust.demoimportbsl
frompyfaust.demoimportbsl
bsl.fig_time_cmp(use_precomputed_data=False)
bsl.fig_time_cmp(use_precomputed_data=False)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
What we see in this figure is that it takes a few dozens of milliseconds (the median time) to compute the BSL experiment on the dense matrix M. This is well above the time it takes with Faust approximates $\hat M_6$ to $\hat M_{23}$ in which the numbers 6 and 23 denote the Faust [RCG](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1Faust.html#a6a51a05c20041504a0b8f2a73dd8d05a). The greater the RCG the better the computation time is, as we already saw in the [notebook about Faust manipulations](#manip_links).
What we see in this figure is that it takes a few dozens of milliseconds (the median time) to compute the BSL experiment on the dense matrix M. This is well above the time it takes with Faust approximates $\hat M_6$ to $\hat M_{23}$ in which the numbers 6 and 23 denote the Faust [RCG](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1Faust.html#a6a51a05c20041504a0b8f2a73dd8d05a). The greater the RCG the better the computation time is, as we already saw in the [notebook about Faust manipulations](#manip_links).
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
As a complementary test, let's verify that the two runs of ``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 $
As a complementary test, let's verify that the two runs of ``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 $
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
nnz=98
nnz=98
x1=zeros((2*n,1))
x1=zeros((2*n,1))
nnz_inds=randperm(2*n)[:nnz]
nnz_inds=randperm(2*n)[:nnz]
x1[nnz_inds,0]=[randn()foriinrange(0,nnz)]
x1[nnz_inds,0]=[randn()foriinrange(0,nnz)]
y=FD@x1
y=FD@x1
x2=omp(y,D,maxiter=nnz)
x2=omp(y,D,maxiter=nnz)
x3=omp(y,FD,maxiter=nnz)
x3=omp(y,FD,maxiter=nnz)
# verify if the solutions differ
# verify if the solutions differ
print("Are x2 and x3 solutions almost equal?",norm(x2-x3)/norm(x3)<10**-12)
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 x2?",allclose(x1,x2))
print("Is x1 retrieved into x3?",allclose(x1,x3))
print("Is x1 retrieved into x3?",allclose(x1,x3))
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
As expected, we didn't retrieve our starting x1 (the reason is the condition already discussed in [2.](#2.-The-Test-Case-Dictionary)). 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.
As expected, we didn't retrieve our starting x1 (the reason is the condition already discussed in [2.](#2.-The-Test-Case-Dictionary)). 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.
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.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
fromnumpy.randomimportrandn
fromnumpy.randomimportrandn
defdirac_hadamard(n):
defdirac_hadamard(n):
FD=hstack((eye(n),wht(n)))
FD=hstack((eye(n),wht(n)))
returnFD
returnFD
n=1024
n=1024
FD=dirac_hadamard(n)
FD=dirac_hadamard(n)
D=FD.toarray()
D=FD.toarray()
x=randn(2*n,1)
x=randn(2*n,1)
%timeitD@x
%timeitD@x
%timeitFD@x
%timeitFD@x
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
For a dimension as smaller than 128, it's possible on particular machines to obtain a slower FD multiplication comparatively to the D multiplication.
For a dimension as smaller than 128, it's possible on particular machines to obtain a slower FD multiplication comparatively to the D multiplication.
This is essentially because the speedup offered by Faust appears rather for higher matrix dimensions.
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.
Let us illustrate the speedup more generally by repeating the experiment for various dimensions n.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%matplotlibinline
%matplotlibinline
fromnumpy.randomimportrandn
fromnumpy.randomimportrandn
fromnumpyimportarray
fromnumpyimportarray
frompyfaustimporteye,wht
frompyfaustimporteye,wht
fromtimeimportprocess_time
fromtimeimportprocess_time
frommatplotlib.pyplotimport(plot,show,
frommatplotlib.pyplotimport(plot,show,
legend,semilogy,
legend,semilogy,
xlabel,ylabel,
xlabel,ylabel,
title,scatter)
title,scatter)
d_times=[]
d_times=[]
fd_times=[]
fd_times=[]
dims=[]
dims=[]
num_muls=100
num_muls=100
fornin[2**iforiinrange(8,13)]:
fornin[2**iforiinrange(8,13)]:
FD=dirac_hadamard(n)
FD=dirac_hadamard(n)
D=FD.toarray()
D=FD.toarray()
x=randn(2*n,1)
x=randn(2*n,1)
dims+=[n]
dims+=[n]
t=process_time()
t=process_time()
[D@xforiinrange(0,num_muls)]
[D@xforiinrange(0,num_muls)]
d_times+=[process_time()-t]
d_times+=[process_time()-t]
t=process_time()
t=process_time()
[FD@xforiinrange(0,num_muls)]
[FD@xforiinrange(0,num_muls)]
fd_times+=[process_time()-t]
fd_times+=[process_time()-t]
p1=semilogy(dims,fd_times,label="FD*x times")
p1=semilogy(dims,fd_times,label="FD*x times")
p2=semilogy(dims,d_times,label="D*x times")
p2=semilogy(dims,d_times,label="D*x times")
xlabel("dimension")
xlabel("dimension")
ylabel("time")
ylabel("time")
legend()
legend()
title("D*x and FD*x time comparison")
title("D*x and FD*x time comparison")
show()
show()
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
As shown for dimensions n above 1024 an actual speedup occurs, the speedup figure below confirms this result.
As shown for dimensions n above 1024 an actual speedup occurs, the speedup figure below confirms this result.
Improving such a speedup and decreasing the dimensions where it occurs is part of the roadmap for future developments of pyfaust.
Improving such a speedup and decreasing the dimensions where it occurs is part of the roadmap for future developments of pyfaust.
***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!
***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!
<aname="[1]">[1]</a> [Tropp, J. A. (2004). Greed is Good: Algorithmic Results for Sparse Approximation. IEEE Transactions on Information Theory, 50(10), 2231–2242](http://doi.org/10.1109/TIT.2004.834793).
<aname="[1]">[1]</a> [Tropp, J. A. (2004). Greed is Good: Algorithmic Results for Sparse Approximation. IEEE Transactions on Information Theory, 50(10), 2231–2242](http://doi.org/10.1109/TIT.2004.834793).
%% Cell type:markdown id: tags:
**Note**: this notebook was executed using the following pyfaust version:
