Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 5884d087 authored by hhakim's avatar hhakim
Browse files

Update/Review MHTP algorithms jupyter notebook.

parent 9efc183e
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
gen_doc/MHTP_algorithms.html
+ 5
− 5
View file @ 5884d087
......@@ -14259,24 +14259,24 @@ 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-PALM4MSA-MHTP-Algorithm">Using the PALM4MSA-MHTP Algorithm<a class="anchor-link" href="#Using-the-PALM4MSA-MHTP-Algorithm">&#182;</a></h1><p>In this notebook we shall see how to use the PALM4MSA-MHTP algorithm. A notebook has already been written on the Hierarchical PALM4MSA algorithm and its wrappers and is a prerequisite to the reading of this notebook.</p>
<h1 id="Using-the-PALM4MSA-MHTP-Algorithm">Using the PALM4MSA-MHTP Algorithm<a class="anchor-link" href="#Using-the-PALM4MSA-MHTP-Algorithm">&#182;</a></h1><p>In this notebook we shall see how to use the PALM4MSA-MHTP algorithm. A <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/Faust_factorization.html">notebook</a> has already been written on the Hierarchical PALM4MSA algorithm and its wrappers and is a prerequisite to the reading of this notebook.</p>
<p>The PALM4MSA-MHTP is a variant of PALM4MSA in which intervenes the Multilinear Hard Tresholdhing Pursuit algorithm (MHTP).</p>
<p>The interest of this variant is to avoid the situation in which PALM4MSA tends to stuck on some matrix supports without any way out. MHTP allows to explore the support more freely and hopefully find a more accurate factorization at the cost of just a few more dozens iterations of the gradient descent algorithm.</p>
<p>For more information on the theory, you can read the following paper in which is treated the particular case of the BHTP (Bilinear HTP, that is running the MHTP on only two factors).</p>
<p><a name="[1]">[1]</a> Quoc-Tung Le, Rémi Gribonval. Structured Support Exploration For Multilayer Sparse Matrix Fac- torization. ICASSP 2021 - IEEE International Conference on Acoustics, Speech and Signal Processing, Jun 2021, Toronto, Ontario, Canada. pp.1-5. <a href="https://hal.inria.fr/hal-03132013/document">hal-03132013</a>.</p>
<h2 id="Configuring-and-Running-PALM4MSA-MHTP">Configuring and Running PALM4MSA-MHTP<a class="anchor-link" href="#Configuring-and-Running-PALM4MSA-MHTP">&#182;</a></h2><p>This variant works very similarly to a classic run of PALM4MSA, that is with at least the same set of parameters. The main difference is that periodically (in term of PALM4MSA iterations) the MHTP algorithm is launched to renew each layer of the Faust being refined.</p>
<h2 id="Configuring-and-Running-PALM4MSA-MHTP">Configuring and Running PALM4MSA-MHTP<a class="anchor-link" href="#Configuring-and-Running-PALM4MSA-MHTP">&#182;</a></h2><p>This variant works very similarly to a classic run of PALM4MSA, that is with at least the same set of parameters. The main difference is that periodically (in term of PALM4MSA number of iterations) the MHTP algorithm is launched to renew each layer of the Faust being refined.</p>
<p>Hence running the PALM4MSA-MHTP needs two sets of parameters: <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ParamsPalm4MSA.html">ParamsPalm4MSA</a> and <a href="https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1MHTPParams.html">MHTPParams</a> objects. The former should not be really new if you are used to PALM4MSA, the latter is dedicated to the configuartion of the MHTP part of PALM4MSA-MHTP.</p>
<p>The arguments to configure <code>MHTPParams</code> are basically:</p>
<ul>
<li><code>num_its</code>: the number of iterations MHTP runs on each layer of the Faust. Remember that this number of iterations is for each factor. If you have two factors the overall number of iterations is <code>2 x num_its</code> (exactly as it is for PALM4MSA).</li>
<li><code>constant_step_size</code> and <code>step_size</code>: that determines if the MHTP gradient descent will be ran according to a constant step size, and in that case how long is the step size. By default, the step size is not constant and recomputed dynamically with the Lipschitz coefficient as in PALM4MSA. In most case, it is recommended to not use a constant step size to achieve a better loss function.</li>
<li><code>constant_step_size</code> and <code>step_size</code>: that determines if the MHTP gradient descent will be ran according to a constant step size, and in that case how long is the step size. By default, the step size is not constant and recomputed dynamically with the Lipschitz coefficient as in PALM4MSA. In most cases, it is recommended to not use a constant step size to achieve a better loss function.</li>
<li><code>palm4msa_period</code>: which governs how many times to evenly run the MHTP algorithm inside PALM4MSA itself. By default, the value is 50. It means that for example if PALM4MSA is running for 200 iterations, MHTP will run 4 times: at iterations 0, 49, 99, 149 and 199 of PALM4MSA. Every time it runs MHTP will run for <code>num_its</code> iterations.</li>
<li><code>updating_lambda</code>: this boolean when set to <code>True</code> allows to update the scale factor of the Faust (the same one that is used in PALM4MSA) in the end of each iteration of MHTP.</li>
</ul>
<p>So let's run PALM4MSA-MHTP on a small example: we propose to factorize a 500x32 matrix into two factors.</p>
<b>First</b> we configure PALM4MSA as usual:
<ul>
<li>The number of iterations of PALM4SA with the StoppingCriterion (here 200 iterations).</li>
<li>The number of iterations of PALM4MSA with the <code>StoppingCriterion</code> (here 200 iterations).</li>
<li>Then we define the constraints / projectors to use, here the SPLIN projector for the first factor of size 500x32 into which we want to count 5 nonzeros per row and the NORMCOL projector for the second factor in which each column must be normalized.</li>
</ul>
 
......@@ -14502,7 +14502,7 @@ PALM4MSA error 0.36661565702911475
<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">
<p>The MEG (for magnetoencephalography) matrix is also used in <a href="#[1]">[1]</a> to compare PALM4MSA and PALM4MSA-MHTP performance.<br>
The goal is to factorize the MEG matrix as $M_{MEG} \approx A \times B$ with $M_{MEG} \in \mathbb{R}^{8293 \times 204}, A \in \mathbb{R}^{8193 \times 204}$ and $B \in \mathbb{R}^{204 \times 204}$. A and B are subject to sparsity constraints. Here we'll test only one sparsity configuration of the two factors ($k_0$ = 100 and $k_1 = 25$ being respectively the per-row number of nonzeros of A and B).</p>
The goal is to factorize the MEG matrix as $M_{MEG} \approx A \times B$ with $M_{MEG} \in \mathbb{R}^{8193 \times 204}, A \in \mathbb{R}^{8193 \times 204}$ and $B \in \mathbb{R}^{204 \times 204}$. A and B are subject to sparsity constraints. Here we'll test only one sparsity configuration of the two factors ($k_0$ = 100 and $k_1 = 25$ being respectively the per-row number of nonzeros of A and B).</p>
<p>Let's load the MEG matrix which is embedded in FAµST data package (which should be downloaded automatically).</p>
 
</div>
......
gen_doc/MHTP_algorithms.ipynb
+ 32
− 7
View file @ 5884d087
%% Cell type:markdown id: tags:
%% Cell type:markdown id:813f8f3e tags:
# Using the PALM4MSA-MHTP Algorithm
In this notebook we shall see how to use the PALM4MSA-MHTP algorithm. A notebook has already been written on the Hierarchical PALM4MSA algorithm and its wrappers and is a prerequisite to the reading of this notebook.
In this notebook we shall see how to use the PALM4MSA-MHTP algorithm. A [notebook](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/Faust_factorization.html) has already been written on the Hierarchical PALM4MSA algorithm and its wrappers and is a prerequisite to the reading of this notebook.
The PALM4MSA-MHTP is a variant of PALM4MSA in which intervenes the Multilinear Hard Tresholdhing Pursuit algorithm (MHTP).
The interest of this variant is to avoid the situation in which PALM4MSA tends to stuck on some matrix supports without any way out. MHTP allows to explore the support more freely and hopefully find a more accurate factorization at the cost of just a few more dozens iterations of the gradient descent algorithm.
For more information on the theory, you can read the following paper in which is treated the particular case of the BHTP (Bilinear HTP, that is running the MHTP on only two factors).
<a name="[1]">[1]</a> Quoc-Tung Le, Rémi Gribonval. Structured Support Exploration For Multilayer Sparse Matrix Fac- torization. ICASSP 2021 - IEEE International Conference on Acoustics, Speech and Signal Processing, Jun 2021, Toronto, Ontario, Canada. pp.1-5. [hal-03132013](https://hal.inria.fr/hal-03132013/document).
## Configuring and Running PALM4MSA-MHTP
This variant works very similarly to a classic run of PALM4MSA, that is with at least the same set of parameters. The main difference is that periodically (in term of PALM4MSA iterations) the MHTP algorithm is launched to renew each layer of the Faust being refined.
This variant works very similarly to a classic run of PALM4MSA, that is with at least the same set of parameters. The main difference is that periodically (in term of PALM4MSA number of iterations) the MHTP algorithm is launched to renew each layer of the Faust being refined.
Hence running the PALM4MSA-MHTP needs two sets of parameters: [ParamsPalm4MSA](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1ParamsPalm4MSA.html) and [MHTPParams](https://faustgrp.gitlabpages.inria.fr/faust/last-doc/html/classpyfaust_1_1factparams_1_1MHTPParams.html) objects. The former should not be really new if you are used to PALM4MSA, the latter is dedicated to the configuartion of the MHTP part of PALM4MSA-MHTP.
The arguments to configure ``MHTPParams`` are basically:
- ``num_its``: the number of iterations MHTP runs on each layer of the Faust. Remember that this number of iterations is for each factor. If you have two factors the overall number of iterations is ``2 x num_its`` (exactly as it is for PALM4MSA).
- ``constant_step_size`` and ``step_size``: that determines if the MHTP gradient descent will be ran according to a constant step size, and in that case how long is the step size. By default, the step size is not constant and recomputed dynamically with the Lipschitz coefficient as in PALM4MSA. In most case, it is recommended to not use a constant step size to achieve a better loss function.
- ``constant_step_size`` and ``step_size``: that determines if the MHTP gradient descent will be ran according to a constant step size, and in that case how long is the step size. By default, the step size is not constant and recomputed dynamically with the Lipschitz coefficient as in PALM4MSA. In most cases, it is recommended to not use a constant step size to achieve a better loss function.
- ``palm4msa_period``: which governs how many times to evenly run the MHTP algorithm inside PALM4MSA itself. By default, the value is 50. It means that for example if PALM4MSA is running for 200 iterations, MHTP will run 4 times: at iterations 0, 49, 99, 149 and 199 of PALM4MSA. Every time it runs MHTP will run for ``num_its`` iterations.
- ``updating_lambda``: this boolean when set to ``True`` allows to update the scale factor of the Faust (the same one that is used in PALM4MSA) in the end of each iteration of MHTP.
So let's run PALM4MSA-MHTP on a small example: we propose to factorize a 500x32 matrix into two factors.
**First** we configure PALM4MSA as usual:
- The number of iterations of PALM4SA with the StoppingCriterion (here 200 iterations).
- The number of iterations of PALM4MSA with the ``StoppingCriterion`` (here 200 iterations).
- Then we define the constraints / projectors to use, here the SPLIN projector for the first factor of size 500x32 into which we want to count 5 nonzeros per row and the NORMCOL projector for the second factor in which each column must be normalized.
%% Cell type:code id: tags:
%% Cell type:code id:f9a11d48 tags:
``` python
from pyfaust.factparams import ParamsPalm4MSA, StoppingCriterion
from pyfaust.proj import splin, normcol
projs = [ splin((500,32), 5), normcol((32,32), 1.0)]
stop_crit = StoppingCriterion(num_its=200)
param = ParamsPalm4MSA(projs, stop_crit)
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:66835080 tags:
**Second** we define the ``MHTPParams`` structure to configure the MHTP pass of PALM4MSA-MHTP
Hence we define the arguments described above: ``num_its``, etc. We let all of them to their default values so there is not much to do.
%% Cell type:code id: tags:
%% Cell type:code id:c0cee9d0 tags:
``` python
from pyfaust.factparams import MHTPParams
mhtp_param = MHTPParams()
print(mhtp_param)
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:d3ab06f1 tags:
It's now time to run the PALM4MSA-MHTP algorithm passing the two structures of parameters. Before we generate a random matrix ``M`` to factorize.
%% Cell type:code id: tags:
%% Cell type:code id:e42eac67 tags:
``` python
import numpy as np
from pyfaust.fact import palm4msa_mhtp
M = np.random.rand(500, 32)
F = palm4msa_mhtp(M, param, mhtp_param)
F
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:c97afd11 tags:
As you see it's pretty similar to running PALM4MSA, which we could have done with the following code.
%% Cell type:code id: tags:
%% Cell type:code id:2b447263 tags:
``` python
from pyfaust.fact import palm4msa
G = palm4msa(M, param)
G
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:c34d7821 tags:
We can verify that the results are however not the same:
%% Cell type:code id: tags:
%% Cell type:code id:928a748c tags:
``` python
print("PALM4MSA-MHTP error:", (F-M).norm()/np.linalg.norm(M))
print("PALM4MSA error", (G-M).norm()/np.linalg.norm(M))
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:9377f592 tags:
They are very close though! In the next part of this notebook we'll demonstrate how PALM4MSA-MHTP can really enhance the accuracy of the Faust approximate and will do that on the MEG matrix (this matrix is also discussed and factorized in <a href="#[1]">[1]</a>).
%% Cell type:markdown id: tags:
%% Cell type:markdown id:7a7b0a83 tags:
## Factorizing the MEG matrix using the PALM4MSA-MHTP algorithm
%% Cell type:markdown id: tags:
%% Cell type:markdown id:47489974 tags:
The MEG (for magnetoencephalography) matrix is also used in <a href="#[1]">[1]</a> to compare PALM4MSA and PALM4MSA-MHTP performance.
The goal is to factorize the MEG matrix as $M_{MEG} \approx A \times B$ with $M_{MEG} \in \mathbb{R}^{8293 \times 204}, A \in \mathbb{R}^{8193 \times 204}$ and $B \in \mathbb{R}^{204 \times 204}$. A and B are subject to sparsity constraints. Here we'll test only one sparsity configuration of the two factors ($k_0$ = 100 and $k_1 = 25$ being respectively the per-row number of nonzeros of A and B).
The goal is to factorize the MEG matrix as $M_{MEG} \approx A \times B$ with $M_{MEG} \in \mathbb{R}^{8193 \times 204}, A \in \mathbb{R}^{8193 \times 204}$ and $B \in \mathbb{R}^{204 \times 204}$. A and B are subject to sparsity constraints. Here we'll test only one sparsity configuration of the two factors ($k_0$ = 100 and $k_1 = 25$ being respectively the per-row number of nonzeros of A and B).
Let's load the MEG matrix which is embedded in FAµST data package (which should be downloaded automatically).
%% Cell type:code id: tags:
%% Cell type:code id:71a15b5f tags:
``` python
from pyfaust.demo import get_data_dirpath
from scipy.io import loadmat
MEG = loadmat(get_data_dirpath()+"/matrix_MEG.mat")['matrix']
MEG.shape
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:cf94923d tags:
Going ahead we set the PALM4MSA parameters:
%% Cell type:code id: tags:
%% Cell type:code id:31fdc250 tags:
``` python
from pyfaust.proj import *
from pyfaust.factparams import ParamsPalm4MSA
k0 = 100
k1 = 25
# the relevant projectors for our sparsity constraints
projs = [splin((8193, 204), k0), splin((204, 204), k1)]
# the number of iterations of PALM4MSA
stop_crit = StoppingCriterion(num_its=2000)
param = ParamsPalm4MSA(projs, stop_crit, is_verbose=True)
param.factor_format = 'dense' # working with dense matrices
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:64a6f5ec tags:
It remains the ``MHTPParams`` configuration (it's easy, we use the default parameters) :
%% Cell type:code id: tags:
%% Cell type:code id:d9f5b222 tags:
``` python
from pyfaust.factparams import MHTPParams
mhtp_p = MHTPParams()
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:f8313943 tags:
Now we are able to launch PALM4MSA and PALM4MSA-MHTP and compare the errors: the computation takes some time, it can last about 30 minutes.
%% Cell type:code id: tags:
%% Cell type:code id:e2d4a744 tags:
``` python
from pyfaust.fact import palm4msa, palm4msa_mhtp
F_palm4msa = palm4msa(MEG, param)
F_mhtp = palm4msa_mhtp(MEG, param, mhtp_p)
err_palm4msa = (F_palm4msa-MEG).norm()/np.linalg.norm(MEG)
err_palm4msa_mhtp = (F_mhtp-MEG).norm()/np.linalg.norm(MEG)
print("PALM4MSA error:", err_palm4msa)
print("PALM4MSA-MHTP error:", err_palm4msa_mhtp)
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:a8aa5b45 tags:
As you see the MHTP variant is twice accurate than PALM4MSA on this configuration.
## Using the Hierarchical PALM4MSA-MHTP algorithm
Exactly the same way you can use the hierarchical factorization with PALM4MSA, it is possible to use the function ``pyfaust.fact.hierachical_mhtp`` to run a hierarchical factorization based on PALM4MSA-MHTP instead of simply PALM4MSA.
The launch of the algorithm function is very similar, you just need to add a ``MHTPParams`` instance to the argument list.
Let's show how with this simple example:
%% Cell type:code id: tags:
%% Cell type:code id:720c6d2b tags:
``` python
from pyfaust.fact import hierarchical_mhtp
from pyfaust.factparams import ParamsHierarchical, StoppingCriterion
from pyfaust.factparams import MHTPParams
from pyfaust.proj import sp, normcol, splin
import numpy as np
M = np.random.rand(500, 32)
fact_cons = [splin((500, 32), 5), sp((32,32), 96), sp((32,32), 96)]
res_cons = [normcol((32,32), 1), sp((32,32), 666), sp((32,32), 333)]
stop_crit1 = StoppingCriterion(num_its=200)
stop_crit2 = StoppingCriterion(num_its=200)
# 50 iterations of MHTP will run every 100 iterations of PALM4MSA (each time PALM4MSA is called by the hierarchical algorithm)
mhtp_param = MHTPParams(num_its=50, palm4msa_period=100)
param = ParamsHierarchical(fact_cons, res_cons, stop_crit1, stop_crit2)
param.is_verbose = True
F = hierarchical_mhtp(M, param, mhtp_param)
F
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id:19011532 tags:
This notebook is ending here. Please note that although the article <a href="#[1]">[1]</a> tackles the optimization problem of approximately factorizing a matrix in two sparse factors with the Bilinear Hard Tresholding Pursuit (BHTP) algorithm, the MHTP is a generalization to N factors that needs further experiences to be mature. Hence the function palm4msa_mhtp and moreover the function hierarchical_mhtp should be considered as experimental code and might evolve significantly in the future.
%% Cell type:markdown id: tags:
%% Cell type:markdown id:eb047dfe tags:
Thanks for reading this notebook! Many other are available at [faust.inria.fr](https://faust.inria.fr).
**Note:** this notebook was executed using the following pyfaust version:
%% Cell type:code id: tags:
%% Cell type:code id:b80572ec tags:
``` python
import pyfaust
pyfaust.version()
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment