Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 34d00676 authored by hhakim's avatar hhakim
Browse files

Update/Review Faust projectors jupyter notebook.

parent cc820941
No related branches found
No related tags found
No related merge requests found
Pipeline #834203 skipped
......@@ -14260,7 +14260,7 @@ a.anchor-link {
<div class="jp-Cell-inputWrapper"><div class="jp-InputPrompt jp-InputArea-prompt">
</div><div class="jp-RenderedHTMLCommon jp-RenderedMarkdown jp-MarkdownOutput " data-mime-type="text/markdown">
<h1 id="Using-The-FA&#181;ST-Projectors-API">Using The FA&#181;ST Projectors API<a class="anchor-link" href="#Using-The-FA&#181;ST-Projectors-API">&#182;</a></h1><p>This notebook put the focus on the <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1proj.html"><code>proj</code></a> module of the pyfaust API. This module provides a bunch of projectors which are for instance necessary to implement proximal operators in <a href="https://link.springer.com/article/10.1007/s10107-013-0701-9">PALM</a> algorithms.</p>
<p>Indeed these projectors matter in the parametrization of the <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a686e523273cf3e38b1b614a69f4b48af">PALM4MSA</a> and <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a7ff9e21a4f0b4acd2107629d788c441c">hierarchical factorization</a> algorithms, so let's maintain their configuration as simple as possible by using projectors!</p>
<p>Indeed these projectors matter in the parameterization of the <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a686e523273cf3e38b1b614a69f4b48af">PALM4MSA</a> and <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a7ff9e21a4f0b4acd2107629d788c441c">hierarchical factorization</a> algorithms, so let's maintain their configuration as simple as possible by using projectors!</p>
<p>First let's explain some generalities about projectors:</p>
<ul>
<li>They are all functor objects (objects that you can call as a function).</li>
......@@ -14843,7 +14843,7 @@ For detailed definitions of these kinds of matrices you can refer to wikipedia:<
</div>
<div class="jp-Cell-inputWrapper"><div class="jp-InputPrompt jp-InputArea-prompt">
</div><div class="jp-RenderedHTMLCommon jp-RenderedMarkdown jp-MarkdownOutput " data-mime-type="text/markdown">
<h3 id="The-NORMLIN-and-NORMCOL-projectors">The NORMLIN and NORMCOL projectors<a class="anchor-link" href="#The-NORMLIN-and-NORMCOL-projectors">&#182;</a></h3><p>The <code>pyfaust.proj</code> module provides two projectors, <code>normlin</code> (resp. <code>normcol</code>) that project a matrix onto the closest matrix with rows (resp. columns) of a prescribed 2-norm.</p>
<h3 id="The-NORMLIN-and-NORMCOL-projectors">The NORMLIN and NORMCOL projectors<a class="anchor-link" href="#The-NORMLIN-and-NORMCOL-projectors">&#182;</a></h3><p>The <code>pyfaust.proj</code> module provides two projectors, <code>normlin</code> (resp. <code>normcol</code>) that projects a matrix onto the closest matrix with rows (resp. columns) of a prescribed 2-norm.</p>
<p>Let's try them:</p>
 
</div>
......
%% Cell type:markdown id: tags:
# Using The FAµST Projectors API
This notebook put the focus on the [``proj``](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1proj.html) module of the pyfaust API. This module provides a bunch of projectors which are for instance necessary to implement proximal operators in [PALM](https://link.springer.com/article/10.1007/s10107-013-0701-9) algorithms.
Indeed these projectors matter in the parametrization of the [PALM4MSA](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a686e523273cf3e38b1b614a69f4b48af) and [hierarchical factorization](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a7ff9e21a4f0b4acd2107629d788c441c) algorithms, so let's maintain their configuration as simple as possible by using projectors!
Indeed these projectors matter in the parameterization of the [PALM4MSA](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a686e523273cf3e38b1b614a69f4b48af) and [hierarchical factorization](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a7ff9e21a4f0b4acd2107629d788c441c) algorithms, so let's maintain their configuration as simple as possible by using projectors!
First let's explain some generalities about projectors:
- They are all functor objects (objects that you can call as a function).
- They are all types defined by child classes of the parent abstract class [``proj_gen``](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1proj__gen.html).
The general pattern to use a projector unfolds in two steps:
1. Instantiate the projector passing the proper arguments.
2. Call this projector (again, as a function) on the matrix you're working on. This step is optional or for test purposes, because generally it is the algorithm implementation responsibility to call the projectors. You just need to feed the algorithms (PALM4MSA for example) with them.
Let's see how to define and use the projectors in the code. For the brief math definitions, I let you consult this [document](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/constraint.png).
Remember that the projector API is documented too, you'll find the link for each projector below. Last, if you're looking for a reference about proximal operators here it is: [proximity operators](http://proximity-operator.net/).
%% Cell type:markdown id: tags:
### The SP projector (projection onto matrices with a prescribed sparsity)
%% Cell type:markdown id: tags:
This projector performs a projection onto matrices with a prescribed sparsity. It governs the global sparsity of a matrix given an integer k.
The matrix $ A \in \mathbb{R}^{m \times n }, A = (a_{ij}), {0<= i < m, 0<= j < n}$, is projected to the closest matrix $ B = (b_{ij})$ such that $ \|B\|_0 = \#\{(i,j): b_{ij} \neq 0 \} \leq k$ which implies , if $k < mn$, that some entries of $A$ are kept in $B$ and others are set to zero. The projector keeps the $k$ most significant values (in term of absolute values or magnitude).
Let's try on an example, here a random matrix.
%% Cell type:code id: tags:
``` python
from numpy.random import rand
A = rand(5,5)*100
print("A=\n", A)
```
%% Cell type:code id: tags:
``` python
import pyfaust
from pyfaust.proj import sp
# 1. instantiate the projector
k = 2
p = sp(A.shape, k, normalized=False)
# 2. project the matrix through it
B = p(A)
print("B=\n", B)
```
%% Cell type:markdown id: tags:
The projector is simply defined by the input matrix shape and the integer k to specify the targeted sparsity.
**Optional normalization**:
As you noticed, the argument ``normalized`` is set to ``False`` in the projector definition. This is the default behaviour. When normalized is ``True``, the result $B$ is normalized according to its Frobenius norm.
The next example gives you a concrete view of what happens when ``normalized`` is True.
%% Cell type:code id: tags:
``` python
from numpy.linalg import norm
from numpy import allclose
pnorm = sp(A.shape, k, normalized=True)
C = pnorm(A)
print("B/norm(B, 'fro') == C: ", allclose(B/norm(B,'fro'),C))
```
%% Cell type:markdown id: tags:
**Sparsity and optional positivity**:
It is also possible to "filter" the negative entries of A by setting the [``pos``](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1sp.html) argument of ``sp`` to ``True``.
You can see the projector as a pipeline, the first stage is to filter out the negative values, then the sp projector is applied and finally the resulting image is normalized if ``normalized==True``.
The following example shows how the projector operates depending on combinations of ``pos`` and ``normalized`` values.
%% Cell type:code id: tags:
``` python
p_pos = sp(A.shape, k, normalized=False, pos=True)
print(p_pos(A))
```
%% Cell type:markdown id: tags:
Well, it's exactly the same as the ``In [2]`` output. The reason is quite obvious, it's because A doesn't contain any negative value, so let's try on a copy of A where we set the ``p_pos(A)`` nonzeros to negative values.
%% Cell type:code id: tags:
``` python
from numpy import nonzero
D = A.copy()
D[nonzero(p_pos(A))] *= -1
print(p_pos(D))
```
%% Cell type:markdown id: tags:
The entries selected when ``p_pos(A)`` is applied are now skipped because they are negative in D, so ``p_pos(D)`` selects the new greatest two values of A in term of magnitude.
What happens now if all values of the matrix are negative ? Let's see it in the next example.
%% Cell type:code id: tags:
``` python
E = - A.copy()
print(p_pos(E))
```
%% Cell type:markdown id: tags:
That should not be surprising that the resulting matrix is a zero matrix, indeed E contains only negative values which are all filtered by setting ``pos=True`` in the ``p_pos`` definition.
%% Cell type:markdown id: tags:
A last question remains: what would happen if we normalize the output matrix when ``pos==True`` and the input matrix is full of negative values ?
The response is simple: a division by zero error would be raised because the norm of a zero matrix is zero, hence it's not possible to normalize.
%% Cell type:markdown id: tags:
### The SPLIN and SPCOL projectors
They are very similar to the ``sp`` projector except that ``splin`` governs the integer sparsity on a row basis and ``spcol`` does it by columns as indicated by the suffix name.
Look at the two short examples, just to be sure.
%% Cell type:code id: tags:
``` python
from pyfaust.proj import splin, spcol
pl = splin(A.shape, k) # reminder: k == 2
pc = spcol(A.shape, k)
B1 = pl(A)
B2 = pc(A)
print("B1=\n", B1)
print("\bB2=\n", B2)
```
%% Cell type:markdown id: tags:
Here the k most significant values are chosen (by rows for splin or by columns for spcol) and the image normalization is disabled.
As for the SP projector, it is possible to incorporate a normalization and/or positivity constraint passing ``normalized=True`` and ``pos=True`` to the functor constructor.
%% Cell type:markdown id: tags:
### The SPLINCOL projector
%% Cell type:markdown id: tags:
The projector [``splincol``](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1splincol.html) tries to constrain the sparsity both by columns and by rows and I wrote it "tries" because there is not always a solution. The use is again the same.
%% Cell type:code id: tags:
``` python
from pyfaust.proj import splincol
plc = splincol(A.shape, k)
print(plc(A))
```
%% Cell type:markdown id: tags:
The image matrix support is in fact the union set of the supports obtained through ``splin`` and ``spcol`` projectors (that's the reason why there is not always a solution). You can refer to this [documentation page](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1splincol.html) which demonstrates in an example how this union is defined.
Another projector for the same purpose but more precise is available for square matrices only. It is named skperm, you'll find its API doc [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1skperm.html). In brief, it is based on a derivate of the hungarian algorithm.
%% Cell type:markdown id: tags:
### The BLOCKDIAG projector
Another matrix projector is the ``blockdiag`` projector. As its name suggests it projects onto the closest block-diagonal matrix with a prescribed structure.
The block-diagonal structure can be defined by the list of the shapes of the block diagonal submatrices you want to keep from the input matrix into the output matrix.
An example will show how it works easily:
%% Cell type:code id: tags:
``` python
%matplotlib inline
from pyfaust.proj import blockdiag
from pyfaust import Faust
R = rand(15,25)
pbd = blockdiag(R.shape, [(1,1), (1,12) ,(R.shape[0]-2, R.shape[1]-13)]) # shapes of blocks in second argument
# show it as a Faust composed of a single factor
Faust([pbd(R)]).imshow()
```
%% Cell type:markdown id: tags:
This blockdiag projector above is defined in order to keep three blocks of the input matrix A, from the upper-left to the lower-right: the first block is the singleton block composed only of the entry (0,0), the second block is a bit of the next row starting from entry (1,1) and finishing to the entry (1, 12) (its shape is (1,12)) and the final block starts from the element (2,13) to finish on the element (R.shape[0]-1, R.shape[1]-1). It's important that the list of blocks covers the whole matrix from its entry (0,0) to its entry (R.shape[0]-1, R.shape[1]-1) or the projector will end up on an error. In other words, if you sum together the first (resp. the second) coordinates of all pair of shapes you must find R.shape[0] (resp. R.shape[1]).
%% Cell type:markdown id: tags:
### The CIRC, TOEPLITZ and HANKEL projectors
These projectors all return the closest corresponding structured matrix (circ being the short name for circulant).
For detailed definitions of these kinds of matrices you can refer to wikipedia:
- [circulant matrix](https://en.wikipedia.org/wiki/Circulant_matrix)
- [toeplitz matrix](https://en.wikipedia.org/wiki/Toeplitz_matrix)
- [hankel matrix](https://en.wikipedia.org/wiki/Hankel_matrix)
The output is constant along each 'diagonal' (resp. 'anti-diagonal'), where the corresponding constant value is the mean of the values of the input matrix along the same diagonal (resp. anti-diagonal).
In the following example, a Faust is constructed with in first factor a circulant matrix, in second position a toeplitz matrix and at the end a hankel matrix.
%% Cell type:code id: tags:
``` python
from pyfaust.proj import circ, toeplitz, hankel
CI = rand(10,10) # circ proj input
TI = rand(10,15) # toeplitz proj input
HI = rand(15,10) # hankel proj input
cp = circ(CI.shape)
tp = toeplitz(TI.shape)
hp = hankel(HI.shape)
F = Faust([cp(CI), tp(TI), hp(HI)])
F.imshow()
```
%% Cell type:markdown id: tags:
You should recognize clearly the structures of a circulant matrix, a toeplitz matrix and a hankel matrix.
Note that these projectors are also capable to receive the ``normalized`` and ``pos`` keyword arguments we've seen before.
The API documentation will give you other examples:
- [circ(culant)](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1circ.html),
- [toeplitz](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1toeplitz.html),
- [hankel](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1hankel.html)
%% Cell type:markdown id: tags:
### The NORMLIN and NORMCOL projectors
The ``pyfaust.proj`` module provides two projectors, ``normlin`` (resp. ``normcol``) that project a matrix onto the closest matrix with rows (resp. columns) of a prescribed 2-norm.
The ``pyfaust.proj`` module provides two projectors, ``normlin`` (resp. ``normcol``) that projects a matrix onto the closest matrix with rows (resp. columns) of a prescribed 2-norm.
Let's try them:
%% Cell type:code id: tags:
``` python
from pyfaust.proj import normcol, normlin
from numpy.linalg import norm
pnl = normlin(A.shape, .2)
pnc = normcol(A.shape, .2)
# let's verify the norm for one column obtained by normlin
print(norm(pnl(A)[2,:]))
# and the same for normcol
print(norm(pnc(A)[:,2]))
```
%% Cell type:markdown id: tags:
Something that is important to notice is the particular case of zero columns or rows. When the NORMLIN (resp. NORMCOL) projector encounters a zero row (resp. a zero column) it simply ignores it.
Let's try:
%% Cell type:code id: tags:
``` python
B = A.copy()
B[:,2] = 0
print(pnc(B)[:,2])
```
%% Cell type:markdown id: tags:
The column 2 is set to zero in B and stays to zero in the projector image matrix.
%% Cell type:markdown id: tags:
### The SUPP projector
The ``supp`` projector projects a matrix onto the closest matrix with a prescribed support. In other words, it preserves the matrix entries lying on this support. The others are set to zero.
The support must be defined by a binary matrix (the ``dtype`` must be the same as the input matrix though).
%% Cell type:code id: tags:
``` python
from pyfaust.proj import supp
from numpy import eye
# keep only the diagonal of A
ps = supp(eye(*A.shape)) # by default normalized=False and pos=False
print(ps(A))
```
%% Cell type:markdown id: tags:
This projector is also capable to receive the ``normalized`` and ``pos`` keyword arguments.
%% Cell type:markdown id: tags:
### The CONST projector
This last projector is really simple, it returns a constant matrix whatever is the input matrix. The way it's instantiated is very similar to the SUPP projector.
Look at its documentation to get an example: [const(ant) projector](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1proj_1_1const.html).
----------------------
%% Cell type:markdown id: tags:
**Thanks** for reading this notebook, you'll find others on the [FAµST website](faust.inria.fr). Any feedback is welcome, all contact information is as well available on the website.
%% Cell type:markdown id: tags:
**Note**: this notebook was executed using the following pyfaust version:
%% Cell type:code id: tags:
``` python
import pyfaust
pyfaust.version()
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment