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.
"Uhat.factors(range(0,min(10,Uhat.numfactors()))).imshow(name='Uhat0-9') # For the sake of visibility \n",
"Uhat.factors(range(0,min(10,Uhat.numfactors()))).imshow(name='Uhat0-9') # For the sake of readabilityk \n",
" # only the first nine factors of Uhat \n",
" # only the first nine factors of Uhat \n",
" # are displayed, as you can see they are\n",
" # are displayed, as you can see they are\n",
" # sparse and Dhat is a diagonal matrix.\n",
" # sparse and Dhat is a diagonal matrix.\n",
"_a = plt.matshow(Dhat.toarray())\n",
"Dhat = diag(Dhat)\n",
"_a = plt.matshow(Dhat)\n",
"title = _a.axes.set_title('Dhat\\n')\n",
"title = _a.axes.set_title('Dhat\\n')\n",
"plt.show()"
"plt.show()"
]
]
...
@@ -447,7 +449,7 @@
...
@@ -447,7 +449,7 @@
"metadata": {},
"metadata": {},
"outputs": [],
"outputs": [],
"source": [
"source": [
"norm(Uhat*Dhat.toarray()*Uhat.T-L)/norm(L)"
"norm(Uhat*Dhat*Uhat.T-L)/norm(L)"
]
]
},
},
{
{
...
...
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
# Using the FAµST Factorization Wrappers
# Using the FAµST Factorization Wrappers
Previous notebooks already addressed a part of the [pyfaust's](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html) functionalities however ``pyfaust.fact`` which is a central module was not covered. It is mainly dedicated to the algorithms that (approximately) factorize a dense matrix into a [Faust](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1Faust.html) object.
Previous notebooks already addressed a part of the [pyfaust's](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html) functionalities however ``pyfaust.fact`` which is a central module was not covered. It is mainly dedicated to the algorithms that (approximately) factorize a dense matrix into a [Faust](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1Faust.html) object.
This is the subject covered (partly) in this notebook.
This is the subject covered (partly) in this notebook.
You will see how to launch the algorithms easily, how to feed them with your own complex set of parameters, how to run them for only few steps when possible, how to use a variant or another, how to enable debugging information, etc.
You will see how to launch the algorithms easily, how to feed them with your own complex set of parameters, how to run them for only few steps when possible, how to use a variant or another, how to enable debugging information, etc.
**NOTE**: the notebook is made to be executed sequentially, otherwise, skipping some cells, you would end up on an import error.
**NOTE**: the notebook is made to be executed sequentially, otherwise, skipping some cells, you would end up on an import error.
### Table of Contents
### Table of Contents
[**1. The Hierarchical PALM4MSA Algorithm**](#1.-The-Hierarchical-PALM4MSA-Algorithm)<br/>
[**1. The Hierarchical PALM4MSA Algorithm**](#1.-The-Hierarchical-PALM4MSA-Algorithm)<br/>
[1.1 Generating a Hadamard Matrix](#1.1-Generating-a-Hadamard-Matrix)<br/>
[1.1 Generating a Hadamard Matrix](#1.1-Generating-a-Hadamard-Matrix)<br/>
[1.2 Factorizing a Hadamard Matrix](#1.2-Factorizing-a-Hadamard-Matrix)<br/>
[1.2 Factorizing a Hadamard Matrix](#1.2-Factorizing-a-Hadamard-Matrix)<br/>
[1.2.1 Defining the Constraints](#1.2.1-Defining-the-Constraints)<br/>
[1.2.1 Defining the Constraints](#1.2.1-Defining-the-Constraints)<br/>
[1.2.2 Setting the Rest of the Parameters and Running the Algorithm](#1.2.2-Setting-the-Rest-of-the-Parameters-and-Running-the-Algorithm)<br/>
[1.2.2 Setting the Rest of the Parameters and Running the Algorithm](#1.2.2-Setting-the-Rest-of-the-Parameters-and-Running-the-Algorithm)<br/>
[**2. The FGFT Algorithms**](#2.-The-FGFT-Algorithms)<br/>
[**2. The FGFT Algorithms**](#2.-The-FGFT-Algorithms)<br/>
[2.1 The Givens / Jacobi Truncated Algorithms](#2.1-The-Givens-/-Jacobi-Truncated-Algorithms)<br/>
[2.1 The Givens / Jacobi Truncated Algorithms](#2.1-The-Givens-/-Jacobi-Truncated-Algorithms)<br/>
[2.2 The Hierarchical PALM4SMA based FGFT Algorithms](#2.2-The-Hierarchical-PALM4SMA-based-FGFT-Algorithms)<br/>
[2.2 The Hierarchical PALM4SMA based FGFT Algorithms](#2.2-The-Hierarchical-PALM4SMA-based-FGFT-Algorithms)<br/>
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## 1. The Hierarchical PALM4MSA Algorithm
## 1. The Hierarchical PALM4MSA Algorithm
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 1.1 Generating a Hadamard Matrix
### 1.1 Generating a Hadamard Matrix
Before to tackle Hadamard matrices factorization, let's introduce one pyfaust's function which is actually directly related: [wht](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html#a7dca1d7342d67faaf765f6c80ce16008).
Before to tackle Hadamard matrices factorization, let's introduce one pyfaust's function which is actually directly related: [wht](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust.html#a7dca1d7342d67faaf765f6c80ce16008).
This method allows you to generate sparse Hadamard transforms.
This method allows you to generate sparse Hadamard transforms.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%matplotlibinline
%matplotlibinline
importmatplotlib.pyplotasplt
importmatplotlib.pyplotasplt
frompyfaustimportwht
frompyfaustimportwht
fromnumpy.linalgimportnorm
fromnumpy.linalgimportnorm
fromnumpy.randomimportrand
fromnumpy.randomimportrand
# generate a Hadamard Faust of size 32x32
# generate a Hadamard Faust of size 32x32
FH=wht(32,normed=False)# normed=False is to avoid column normalization
FH=wht(32,normed=False)# normed=False is to avoid column normalization
H=FH.toarray()# the dense matrix version
H=FH.toarray()# the dense matrix version
FH.imshow()
FH.imshow()
plt.show()
plt.show()
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
fromnumpyimportcount_nonzero
fromnumpyimportcount_nonzero
count_nonzero(H)
count_nonzero(H)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
All the factors are the same, let's print the first one.
All the factors are the same, let's print the first one.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
plt.imshow(FH.factors(0).toarray())
plt.imshow(FH.factors(0).toarray())
plt.show()
plt.show()
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
You might want to verify if $H$ is a proper Hadamard matrix.
You might want to verify if $H$ is a proper Hadamard matrix.
The ugly code above basically verifies that all column or row vectors are mutually orthogonal and made only of -1, 0 and 1 coefficients, that is exactly the Hadamard matrix definition. The response is yes, so we can go ahead.
The ugly code above basically verifies that all column or row vectors are mutually orthogonal and made only of -1, 0 and 1 coefficients, that is exactly the Hadamard matrix definition. The response is yes, so we can go ahead.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 1.2 Factorizing a Hadamard Matrix
### 1.2 Factorizing a Hadamard Matrix
Let's begin by factorizing $H$ the easy way with the automatic parametrization provided by pyfaust.
Let's begin by factorizing $H$ the easy way with the automatic parametrization provided by pyfaust.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
frompyfaustimportfaust_fact
frompyfaustimportfaust_fact
FH2=faust_fact(H,'hadamard')
FH2=faust_fact(H,'hadamard')
print(FH2)
print(FH2)
FH2.imshow(name='FH2')
FH2.imshow(name='FH2')
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Interestingly, the FH2's factors are not the same as FH's but the resulting numpy array is exactly the same.
Interestingly, the FH2's factors are not the same as FH's but the resulting numpy array is exactly the same.
Just to be sure let's verify the relative error against H.
Just to be sure let's verify the relative error against H.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
print('err=',(FH2-H).norm()/norm(H))
print('err=',(FH2-H).norm()/norm(H))
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Good! The argument ``'hadamard'`` of ``faust_fact`` is a shorthand to a set of parameters fitting the factorization of Hadamard matrices.
Good! The argument ``'hadamard'`` of ``faust_fact`` is a shorthand to a set of parameters fitting the factorization of Hadamard matrices.
Speaking about shorthand, ```faust_fact``` is an alias of [pyfaust.fact.hierarchical](file:///home/hinria/faust/build/doc/html/namespacepyfaust_1_1fact.html#ad2108e97e89a17da27ce480d65155eb4).
Speaking about shorthand, ```faust_fact``` is an alias of [pyfaust.fact.hierarchical](file:///home/hinria/faust/build/doc/html/namespacepyfaust_1_1fact.html#ad2108e97e89a17da27ce480d65155eb4).
This function implements the hiearchical factorization based on PALM4MSA. For details about the theory behind these algorithms you can read the thesis of [Luc Le Magoarou](https://tel.archives-ouvertes.fr/tel-01412558).
This function implements the hiearchical factorization based on PALM4MSA. For details about the theory behind these algorithms you can read the thesis of [Luc Le Magoarou](https://tel.archives-ouvertes.fr/tel-01412558).
We won't go into further details nevertheless it's important to see how to define the parameters manually in order to proceed to your own custom factorizations. For that purpose we need to describe approximately the algorithm but please keep in mind this is just an insight (for instance, the norms below are voluntarily not defined).
We won't go into further details nevertheless it's important to see how to define the parameters manually in order to proceed to your own custom factorizations. For that purpose we need to describe approximately the algorithm but please keep in mind this is just an insight (for instance, the norms below are voluntarily not defined).
Basically, the hierarchical algorithm works like this:
Basically, the hierarchical algorithm works like this:
- If you want to decompose a matrix $M$ in $J$ factors, the algorithm will iterate $(J-1)$ times. Each iteration follows two steps:<br/>
- If you want to decompose a matrix $M$ in $J$ factors, the algorithm will iterate $(J-1)$ times. Each iteration follows two steps:<br/>
**1.** (_local optimization_) At each iteration $i$, the currently considered matrix $R_{i-1}$ (with $R_0=M$) is decomposed in two factors by the PALM4MSA algorithm as respect to the minimization $\| R_{i-1} - S_{i}R_{i} \|$. <br/>$S_i$ is the resulting factor while $R_{i}$ is the _residual_ of our factorization.<br/>
**1.** (_local optimization_) At each iteration $i$, the currently considered matrix $R_{i-1}$ (with $R_0=M$) is decomposed in two factors by the PALM4MSA algorithm as respect to the minimization $\| R_{i-1} - S_{i}R_{i} \|$. <br/>$S_i$ is the resulting factor while $R_{i}$ is the _residual_ of our factorization.<br/>
**2.** (_global optimization_) At the end of each iteration $i$, PALM4MSA is called again to ideally compute the $argmin_{\{S_1, ..., S_i\}, R_i}$ $\|M - (\prod_{j=1}^i S_j) R_i\|$
**2.** (_global optimization_) At the end of each iteration $i$, PALM4MSA is called again to ideally compute the $argmin_{\{S_1, ..., S_i\}, R_i}$ $\|M - (\prod_{j=1}^i S_j) R_i\|$
- So ideally at the end of the iteration (J-1) you'll get something like $M \approx \prod_{i=1}^J S_i$ taking $S_J = R_{J-1} $.
- So ideally at the end of the iteration (J-1) you'll get something like $M \approx \prod_{i=1}^J S_i$ taking $S_J = R_{J-1} $.
### 1.2.1 Defining the Constraints
### 1.2.1 Defining the Constraints
The explanation above is eluding something essential: the sparsity constraints. Indeed, the purpose of the algorithm is not only to decompose a matrix but also to enhance its sparsity, as hinted by the FAµST acronym: Flexible Approximate Multi-layer Sparse Transform.
The explanation above is eluding something essential: the sparsity constraints. Indeed, the purpose of the algorithm is not only to decompose a matrix but also to enhance its sparsity, as hinted by the FAµST acronym: Flexible Approximate Multi-layer Sparse Transform.
So you'll need to feed the algorithm with sparsity constraints. In fact, you'll define one pair of constraints per iteration, the first is for the factor $S_{i}$ and the second for $R_i$ (the _residual_).
So you'll need to feed the algorithm with sparsity constraints. In fact, you'll define one pair of constraints per iteration, the first is for the factor $S_{i}$ and the second for $R_i$ (the _residual_).
The pyfaust API is here to help you define the constraints in one shot but you can if you want define constraints one by one as we'll see later.
The pyfaust API is here to help you define the constraints in one shot but you can if you want define constraints one by one as we'll see later.
Let's unveil the factorization constraints to decompose the Hadamard matrix $H \in \mathbb R^{n \times n}$.
Let's unveil the factorization constraints to decompose the Hadamard matrix $H \in \mathbb R^{n \times n}$.
Let's go back to the code and define our constraints as lists (provided by the module ```pyfaust.factparams```):
Let's go back to the code and define our constraints as lists (provided by the module ```pyfaust.factparams```):
The ```'splincol'``` constraint used here is for defining the maximum number of nonzeros elements to respect for any row or column of the considered matrix (here $S_i$ or $R_i$).
The ```'splincol'``` constraint used here is for defining the maximum number of nonzeros elements to respect for any row or column of the considered matrix (here $S_i$ or $R_i$).
Looking at ```S_constraints``` initialization, we see in the first line ```'splincol', 2, n, n```; the value 2 means we target 2 nonzeros per-column and per-row, the next arguments define the size of the matrix to constrain (its number of rows and columns).
Looking at ```S_constraints``` initialization, we see in the first line ```'splincol', 2, n, n```; the value 2 means we target 2 nonzeros per-column and per-row, the next arguments define the size of the matrix to constrain (its number of rows and columns).
So in the example of ```S_constraints``` all the targeted matrices have a 0-norm equal to 2n.
So in the example of ```S_constraints``` all the targeted matrices have a 0-norm equal to 2n.
More generally in pyfaust, the constraints are defined in terms of norms, in our case of Hadamard factorization we used the following 0-norms:
More generally in pyfaust, the constraints are defined in terms of norms, in our case of Hadamard factorization we used the following 0-norms:
Enough of explanations! I let you decrypt why ```R_constraints``` correspond likewise to (C2).
Enough of explanations! I let you decrypt why ```R_constraints``` correspond likewise to (C2).
But wait a minute, what would happen if we have a set of one hundred constraints? Don't worry! The pyfaust API allows to alternatively set the constraints one by one:
But wait a minute, what would happen if we have a set of one hundred constraints? Don't worry! The pyfaust API allows to alternatively set the constraints one by one:
from pyfaust.factparams import ConstraintInt
from pyfaust.factparams import ConstraintInt
S_constraints = [ConstraintInt('splincol', n, n, 2) for i in range(0,4)]
S_constraints = [ConstraintInt('splincol', n, n, 2) for i in range(0,4)]
R_constraints = [ConstraintInt('splincol', n, n, int(n/2**i)) for i in range(1,5)]
R_constraints = [ConstraintInt('splincol', n, n, int(n/2**i)) for i in range(1,5)]
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
```S_constraints``` and ```R_constraints``` are still the exact same set of constraints as before.
```S_constraints``` and ```R_constraints``` are still the exact same set of constraints as before.
[ConstraintInt](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintInt.html) is a family of integer constraints. More globally, there is a hierarchy of classes whose parent class is [ConstraintGeneric](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintGeneric.html) where you'll find other kind of constraints ; [ConstraintMat](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintMat.html), [ConstraintReal](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintReal.html).
[ConstraintInt](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintInt.html) is a family of integer constraints. More globally, there is a hierarchy of classes whose parent class is [ConstraintGeneric](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintGeneric.html) where you'll find other kind of constraints ; [ConstraintMat](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintMat.html), [ConstraintReal](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintReal.html).
The table [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/constraint.png) can give you an idea of each constraint definition. Most constraints are associated to a proximal operator which projects the matrix into the set of matrices defined by the constraint or at least it tries to (because it's not necessarily possible).
The table [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/constraint.png) can give you an idea of each constraint definition. Most constraints are associated to a proximal operator which projects the matrix into the set of matrices defined by the constraint or at least it tries to (because it's not necessarily possible).
In the latest version of FAµST it's possible to ```project``` matrices from the wrappers themselves, let's try one ConstraintInt:
In the latest version of FAµST it's possible to ```project``` matrices from the wrappers themselves, let's try one ConstraintInt:
We asked for a sparsity of 18 nonzeros and that's what we got after calling ```project```.
We asked for a sparsity of 18 nonzeros and that's what we got after calling ```project```.
The function ```project``` can help you debug a factorization or even understand exactly how a _prox_ works for a specific constraint and matrix.
The function ```project``` can help you debug a factorization or even understand exactly how a _prox_ works for a specific constraint and matrix.
Another API doc link completes the definition of constraints : [ConstraintName](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintName.html)
Another API doc link completes the definition of constraints : [ConstraintName](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ConstraintName.html)
If you find it complicated to parameterize your own set of constraints, well, you are not the only one! You may be happy to know that it is one of the priorities of the development team to provide a simplified API for the factorization algorithms in an upcoming release
If you find it complicated to parameterize your own set of constraints, well, you are not the only one! You may be happy to know that it is one of the priorities of the development team to provide a simplified API for the factorization algorithms in an upcoming release
### 1.2.2 Setting the Rest of the Parameters and Running the Algorithm
### 1.2.2 Setting the Rest of the Parameters and Running the Algorithm
OK, let's continue defining the algorithm parameters, the constraints are the most part of it. You have optional parameters too, but one key parameter is the stopping criterion of the PALM4MSA algorithm.
OK, let's continue defining the algorithm parameters, the constraints are the most part of it. You have optional parameters too, but one key parameter is the stopping criterion of the PALM4MSA algorithm.
You have to define two stopping criteria, one for the _local optimization_ and another for the _global optimization_.
You have to define two stopping criteria, one for the _local optimization_ and another for the _global optimization_.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
from pyfaust.factparams import StoppingCriterion
from pyfaust.factparams import StoppingCriterion
loc_stop = StoppingCriterion(num_its=30)
loc_stop = StoppingCriterion(num_its=30)
glob_stop = StoppingCriterion(num_its=30)
glob_stop = StoppingCriterion(num_its=30)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
The type of stopping criterion used here is the number of iterations of PALM4MSA, that is the number of times all the current factors ($R_i$ included) are updated either for the local or the global optimization. So, with the initialization of ```loc_stop``` and ```glob_stop``` above, if you ask J factors (through the number of constraints you set) you'd count a total of $(30+30)*(J-1)$ iterations of PALM4MSA.
The type of stopping criterion used here is the number of iterations of PALM4MSA, that is the number of times all the current factors ($R_i$ included) are updated either for the local or the global optimization. So, with the initialization of ```loc_stop``` and ```glob_stop``` above, if you ask J factors (through the number of constraints you set) you'd count a total of $(30+30)*(J-1)$ iterations of PALM4MSA.
Ok, I think we're done with the algorithm parameters, we can pack them into one object and give it to ```faust_fact```.
Ok, I think we're done with the algorithm parameters, we can pack them into one object and give it to ```faust_fact```.
is_update_way_R2L=True) # the argument order matters!
is_update_way_R2L=True) # the argument order matters!
# launch the factorization
# launch the factorization
FH3 = faust_fact(H, params)
FH3 = faust_fact(H, params)
FH3
FH3
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
You might wonder what is the boolean ```is_update_way_R2L```. Its role is to define if the PALM4MSA algorithm will update the $S_i$ and last $R_i$ factors from the left to right (if ```False```) or toward the opposite direction (if ```True```). It does change the results of factorization!
You might wonder what is the boolean ```is_update_way_R2L```. Its role is to define if the PALM4MSA algorithm will update the $S_i$ and last $R_i$ factors from the left to right (if ```False```) or toward the opposite direction (if ```True```). It does change the results of factorization!
I must mention that there are several other optional arguments you'd want to play with for configuring, just dive into the [documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html).
I must mention that there are several other optional arguments you'd want to play with for configuring, just dive into the [documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html).
You also might call the PALM4MSA algorithm directly (not the hierarchical algorithm which makes use of it), the function is [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a62842439447dad2e2e696f2e695d6399). You can theoretically reproduce the hierarchical algorithm step by step calling ```fact_palm4msa()``` by yourself.
You also might call the PALM4MSA algorithm directly (not the hierarchical algorithm which makes use of it), the function is [here](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a62842439447dad2e2e696f2e695d6399). You can theoretically reproduce the hierarchical algorithm step by step calling ```fact_palm4msa()``` by yourself.
Well, there would be a lot to say and show about these two pyfaust algorithms and their parameters but at least with this example, I hope, you got an insight of how it works.
Well, there would be a lot to say and show about these two pyfaust algorithms and their parameters but at least with this example, I hope, you got an insight of how it works.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## 2. The FGFT Algorithms
## 2. The FGFT Algorithms
Since version 2.4, FAµST is able to compute Fast Graph Fourier Transforms.
Since version 2.4, FAµST is able to compute Fast Graph Fourier Transforms.
FAµST includes two algorithms with this goal:
FAµST includes two algorithms with this goal:
- The first one is a variant of the hierarchical algorithm discussed in [1](#1.-The-Hierarchical-PALM4MSA-Algorithm).
- The first one is a variant of the hierarchical algorithm discussed in [1](#1.-The-Hierarchical-PALM4MSA-Algorithm).
- The second one is a truncated version of the Jacobi algorithm for matrix diagonalization.
- The second one is a truncated version of the Jacobi algorithm for matrix diagonalization.
The two of them are implemented in C++ and Python wrappers make them available from pyfaust.
The two of them are implemented in C++ and Python wrappers make them available from pyfaust.
Here again the theory isn't what matters but feel free to take a look at the following papers:
Here again the theory isn't what matters but feel free to take a look at the following papers:
- Le Magoarou L., Gribonval R., Tremblay N., [“Approximate fast graph Fourier transforms via multi-layer sparse approximation“,Transactions on Signal and Information Processing over Networks](https://hal.inria.fr/hal-01416110)
- Le Magoarou L., Gribonval R., Tremblay N., [“Approximate fast graph Fourier transforms via multi-layer sparse approximation“,Transactions on Signal and Information Processing over Networks](https://hal.inria.fr/hal-01416110)
- Le Magoarou L. and Gribonval R., ["Are there approximate Fast Fourier Transforms on graphs ?", ICASSP, 2016](https://hal.inria.fr/hal-01254108)
- Le Magoarou L. and Gribonval R., ["Are there approximate Fast Fourier Transforms on graphs ?", ICASSP, 2016](https://hal.inria.fr/hal-01254108)
### 2.0 The pygsp Toolbox
### 2.0 The pygsp Toolbox
Since we're going to work with graph Laplacians it would be a pity not to plot things. The pygsp toolbox provides what you need to easily obtain Laplacians of different graph families. Please install the [package](https://pypi.org/project/PyGSP/) with pip (or maybe conda).
Since we're going to work with graph Laplacians it would be a pity not to plot things. The pygsp toolbox provides what you need to easily obtain Laplacians of different graph families. Please install the [package](https://pypi.org/project/PyGSP/) with pip (or maybe conda).
Let's plot the Minnesota graph (on the right) and its associated Laplacian matrix (on the left).
Let's plot the Minnesota graph (on the right) and its associated Laplacian matrix (on the left).
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
from pygsp.graphs import minnesota
from pygsp.graphs import minnesota
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
G = minnesota.Minnesota()
G = minnesota.Minnesota()
fig, axes = plt.subplots(1, 2)
fig, axes = plt.subplots(1, 2)
_ = axes[0].spy(G.L, markersize=2)
_ = axes[0].spy(G.L, markersize=2)
_ = G.plot(ax=axes[1])
_ = G.plot(ax=axes[1])
plt.show()
plt.show()
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
The goal of FGFT algorithms is to obtain an orthonormal Faust that nearly diagonalizes a symmetric matrix. In terms of graph signal processing, the considered symmetric matrix is typically a graph Laplacian $L$, and the Faust is expected to approximate the graph Fourier basis, i.e., the basis of eigenvectors of $L$. Given $L$, the algorithm outputs an orthonormal Faust $\hat U$ and a diagonal factor $\hat D$ such that: $L \approx \hat U \hat D \hat U^T$.
The goal of FGFT algorithms is to obtain an orthonormal Faust that nearly diagonalizes a symmetric matrix. In terms of graph signal processing, the considered symmetric matrix is typically a graph Laplacian $L$, and the Faust is expected to approximate the graph Fourier basis, i.e., the basis of eigenvectors of $L$. Given $L$, the algorithm outputs an orthonormal Faust $\hat U$ and a diagonal factor $\hat D$ such that: $L \approx \hat U \hat D \hat U^T$.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 2.1 The Truncated Jacobi Algorithms
### 2.1 The Truncated Jacobi Algorithms
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Let's run the first algorithm (truncated Jacobi for eigen decomposition, hence the name ```eigtj```) on the Minnesota Laplacian example. You'll see it's really straightforward!
Let's run the first algorithm (truncated Jacobi for eigen decomposition, hence the name ```eigtj```) on the Minnesota Laplacian example. You'll see it's really straightforward!
Uhat.factors(range(0,min(10,Uhat.numfactors()))).imshow(name='Uhat0-9') # For the sake of visibility
Uhat.factors(range(0,min(10,Uhat.numfactors()))).imshow(name='Uhat0-9') # For the sake of readabilityk
# only the first nine factors of Uhat
# only the first nine factors of Uhat
# are displayed, as you can see they are
# are displayed, as you can see they are
# sparse and Dhat is a diagonal matrix.
# sparse and Dhat is a diagonal matrix.
_a = plt.matshow(Dhat.toarray())
Dhat = diag(Dhat)
_a = plt.matshow(Dhat)
title = _a.axes.set_title('Dhat\n')
title = _a.axes.set_title('Dhat\n')
plt.show()
plt.show()
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Write some code to check the first factor of $\hat U $ is orthogonal and then that $\hat U$ itself is orthogonal (with a tolerance near to the machine espilon).
Write some code to check the first factor of $\hat U $ is orthogonal and then that $\hat U$ itself is orthogonal (with a tolerance near to the machine espilon).
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
from numpy import eye
from numpy import eye
Uhat0 = Uhat.factors(0)
Uhat0 = Uhat.factors(0)
n = Uhat0.shape[0]
n = Uhat0.shape[0]
((Uhat0*Uhat0.T).toarray() == eye(n)).all()
((Uhat0*Uhat0.T).toarray() == eye(n)).all()
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
((Uhat*Uhat.T)-eye(n)).norm()
((Uhat*Uhat.T)-eye(n)).norm()
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Ignoring the fact that $\hat U$ is approximate and maybe some slight numerical error, yes it is orthogonal.
Ignoring the fact that $\hat U$ is approximate and maybe some slight numerical error, yes it is orthogonal.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
norm(Uhat*Dhat.toarray()*Uhat.T-L)/norm(L)
norm(Uhat*Dhat*Uhat.T-L)/norm(L)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
Good! The Laplacian approximation defined by $\hat U$ and $\hat D$ seems not so bad according to the relative error. Do you want a smaller error ? It's time to read the [documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a82967b7913a66c945f32b7bd1ec6625a), the ```eigtj()```'s keyword argument ```J``` can help you in that goal.
Good! The Laplacian approximation defined by $\hat U$ and $\hat D$ seems not so bad according to the relative error. Do you want a smaller error ? It's time to read the [documentation](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a82967b7913a66c945f32b7bd1ec6625a), the ```eigtj()```'s keyword argument ```J``` can help you in that goal.
Is that all? No, we can evaluate the memory saving $\hat U$ brings as a Faust relatively to the matrix U obtained by ```numpy.linalg.eig()```.
Is that all? No, we can evaluate the memory saving $\hat U$ brings as a Faust relatively to the matrix U obtained by ```numpy.linalg.eig()```.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
from numpy.linalg import eig
from numpy.linalg import eig
from numpy import diag, count_nonzero
from numpy import diag, count_nonzero
D, U = eig(L)
D, U = eig(L)
D = diag(D)
D = diag(D)
print(norm(U@D@U.T-L)/norm(L))
print(norm(U@D@U.T-L)/norm(L))
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
Uhat.nnz_sum()
Uhat.nnz_sum()
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
count_nonzero(U)
count_nonzero(U)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
The memory space to store $\hat U$ is smaller since its number of nonzeros elements is less than that of $U$.
The memory space to store $\hat U$ is smaller since its number of nonzeros elements is less than that of $U$.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
x = rand(n,1)
x = rand(n,1)
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
timeit Uhat@x # note: @ and * are both the matrix product when a Faust is one of the operands
timeit Uhat@x # note: @ and * are both the matrix product when a Faust is one of the operands
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
timeit U@x # note: * isn't equivalent to matrix product when using numpy arrays
timeit U@x # note: * isn't equivalent to matrix product when using numpy arrays
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
The execution time for Faust-vector multiplication is really better too (compared to the matrix-vector multiplication) !
The execution time for Faust-vector multiplication is really better too (compared to the matrix-vector multiplication) !
Applying the Graph Fourier Transform $\hat U$ on x is faster.
Applying the Graph Fourier Transform $\hat U$ on x is faster.
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
### 2.2 The Hierarchical PALM4MSA based FGFT Algorithms
### 2.2 The Hierarchical PALM4MSA based FGFT Algorithms
To finish this notebook briefly, we won't go into further details about the second algorithm for computing approximate FGFTs because it's an extension of the algorithms discussed in [1](#1.-The-Hierarchical-PALM4MSA-Algorithm).
To finish this notebook briefly, we won't go into further details about the second algorithm for computing approximate FGFTs because it's an extension of the algorithms discussed in [1](#1.-The-Hierarchical-PALM4MSA-Algorithm).
However, you'll find below a copy of the main example of [fgft_palm()](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a82967b7913a66c945f32b7bd1ec6625a). This function is the entry point of the algorithm and again it's very similar to what you've seen above.
However, you'll find below a copy of the main example of [fgft_palm()](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/namespacepyfaust_1_1fact.html#a82967b7913a66c945f32b7bd1ec6625a). This function is the entry point of the algorithm and again it's very similar to what you've seen above.
You can try to unroll it progressively and understand how it tries to do the same thing as the eigtj() on another Laplacian matrix.
You can try to unroll it progressively and understand how it tries to do the same thing as the eigtj() on another Laplacian matrix.
It's much more complicated to configure but starting from an example you'll be able to write your own FGFT calculation.
It's much more complicated to configure but starting from an example you'll be able to write your own FGFT calculation.
About the Laplacian matrix used in the code, note that pyfaust embeds some data examples (matrices in matlab v7 format) but depending on how you installed the package you may need to download them. It's not a big deal normally because ```get_data_dirpath()``` should handle it automatically (otherwise please read the [help doc](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/FAQ.html)). Note that depending on your connection the download can take some time.
About the Laplacian matrix used in the code, note that pyfaust embeds some data examples (matrices in matlab v7 format) but depending on how you installed the package you may need to download them. It's not a big deal normally because ```get_data_dirpath()``` should handle it automatically (otherwise please read the [help doc](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/FAQ.html)). Note that depending on your connection the download can take some time.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
from pyfaust.fact import fgft_palm
from pyfaust.fact import fgft_palm
from pyfaust.factparams import *
from pyfaust.factparams import *
from pyfaust.demo import get_data_dirpath
from pyfaust.demo import get_data_dirpath
from os.path import sep
from os.path import sep
from scipy.io import loadmat, savemat
from scipy.io import loadmat, savemat
from numpy.linalg import eig, eigh, norm
from numpy.linalg import eig, eigh, norm
from numpy import sort, argsort, log2, size, copy, diag
from numpy import sort, argsort, log2, size, copy, diag
d = loadmat(sep.join((get_data_dirpath(),'Laplacian_256_ring.mat')))
d = loadmat(sep.join((get_data_dirpath(),'Laplacian_256_ring.mat')))
Lap = d['Lap'].astype('float')
Lap = d['Lap'].astype('float')
D, U = eig(Lap)
D, U = eig(Lap)
indices = argsort(D)
indices = argsort(D)
D = D[indices].astype('float')
D = D[indices].astype('float')
U = U[:,indices].astype('float')
U = U[:,indices].astype('float')
print(D.shape, type(D))
print(D.shape, type(D))
print("eig(Lap), U error:", norm(Lap.dot(U)-U.dot(diag(D))))
print("eig(Lap), U error:", norm(Lap.dot(U)-U.dot(diag(D))))
***The third notebook is ending here***, I hope you'll be interested to dig into pyfaust API and maybe even give some feedback later. The API has made a lot of progress lastly but it remains a work-in-progress and some bugs might appear here or there. Besides, the documentation and the API can always receive some enhancement so don't hesitate, any suggestion is welcome!
***The third notebook is ending here***, I hope you'll be interested to dig into pyfaust API and maybe even give some feedback later. The API has made a lot of progress lastly but it remains a work-in-progress and some bugs might appear here or there. Besides, the documentation and the API can always receive some enhancement so don't hesitate, any suggestion is welcome!
