# Sur l'approximation d'opérateurs linéaires par une somme de produits de matrices creuses

Moujahid Bou-Laouz

<figure><img src="image_1.png" /></figure>

## To cite this version:

Moujahid Bou-Laouz. Sur l'approximation d'opérateurs linéaires par une somme de produits de matrices creuses. Aix-Marseille Université. 2023. hal-04276822

HAL Id: hal-04276822

https://hal.science/hal-04276822v1

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
Master Data science
Mémoire de stage de fin d'étude

Sur l'approximation d'opérateurs linéaires
par une somme de produits de matrices
creuses

Rédigé par :

Moujahid BOU-LAOUZ

Sous la direction de :

Valentin EMIYA

Caroline CHAUX

Année universitaire : 2020-2021
À ma chère mère
Je n'arriverai jamais à t'exprimer mon amour profond
Et surtout mon respect infini ...
REMERCIEMENTS

Je remercie mes superviseurs Monsieur Valentin EMIYA et Madame Caroline CHAUX de
m'avoir proposé ce stage de recherche qui était en totale adéquation avec mes attentes. Votre
rigueur, votre disponibilité et vos qualités humaines m'ont profondément touché...

Mes remerciements vont aussi à Monsieur Hakim HADJ-DJILANI pour son aide précieuse, pour
sa disponibilité et pour sa gentillesse.

Je remercie également Monsieur Pierre PUDLO pour son implication dans mon parcours et pour
avoir répondu à mes nombreuses demandes tout au long de l'année.

Finalement, je tiens à remercier Faiz EL ARJI, Ama Marina KREME, Clothilde MELOT et
Alexandre GAUDILLIERE pour leurs supports et leurs conseils.
Résumé

Les opérateurs linéaires sont les principaux responsables de la complexité temporelle et spatiale de plusieurs méthodes en apprentissage automatique et en traitement du signal. Une façon de réduire le coût du stockage et de calcul d'un opérateur est de l'approximer par un produit de matrices creuses. On propose dans cette étude une généralisation de ce principe. On approxime l'opérateur par une somme de produits de matrices creuses. On présente d'abord les motivations de cette généralisation, on montre ensuite l'efficacité de cette méthode sur plusieurs opérateurs linéaires.

Table des matières

<table><tr><td>1 Introduction et définition</td><td>3</td></tr><tr><td>2 L'algorithme PALM</td><td>5</td></tr><tr><td>2.1 Problème d'optimisation</td><td>5</td></tr><tr><td>2.2 Algorithme</td><td>5</td></tr><tr><td>3 Le modèle FAμST</td><td>6</td></tr><tr><td>3.1 Fonction objectif</td><td>7</td></tr><tr><td>3.2 Une précision sur la projection</td><td>7</td></tr><tr><td>3.2.1 Projection sur E<sub>s</sub></td><td>8</td></tr><tr><td>3.2.2 Projection sur L<sub>s</sub> et C<sub>s</sub></td><td>8</td></tr><tr><td>3.2.3 Projection sur R<sub>s</sub></td><td>8</td></tr><tr><td>3.3 L'algorithme palm4MSA</td><td>8</td></tr><tr><td>3.4 palm4MSA hiérarchique</td><td>9</td></tr><tr><td>3.5 Remarques sur les algorithmes</td><td>10</td></tr><tr><td>4 Motivation du modèle Σ-FAμST</td><td>10</td></tr><tr><td>4.1 Expressivité</td><td>11</td></tr><tr><td>4.2 Flexibilité</td><td>11</td></tr><tr><td>5 Algorithmes proposés pour le modèle Σ-FAμST</td><td>13</td></tr><tr><td>5.1 Formulation du problème</td><td>13</td></tr><tr><td>5.2 Gradient et constante de Lipschitz</td><td>13</td></tr><tr><td>5.3 L'algorithme palm4SMSA</td><td>15</td></tr><tr><td>5.4 Discussion sur l'initialisation</td><td>15</td></tr><tr><td>5.5 Palm4SMSA hiérarchique</td><td>18</td></tr><tr><td>6 Expériences</td><td>19</td></tr><tr><td>6.1 Vers une meilleure distribution des facteurs</td><td>20</td></tr><tr><td>6.2 Sur la difficulté d'apprendre des transformées rapides</td><td>21</td></tr><tr><td>6.3 Sur l'approximation des opérateurs de grande taille</td><td>22</td></tr><tr><td>6.4 Sur la Réduction de coût du stockage d'une base de donnée</td><td>24</td></tr><tr><td>7 Discussion</td><td>24</td></tr><tr><td>Appendices</td><td>27</td></tr><tr><td>A Preuve de la proposition 1</td><td>27</td></tr><tr><td>B Sur la sensibilité de palm4MSA hiérarchique</td><td>27</td></tr></table>
# 1 Introduction et définition

Les matrices jouent un rôle crucial en apprentissage automatique et en traitement du signal. On peut les retrouver sous la forme de tableaux de données, d'images, ou encore, sous la forme d'opérateurs linéaires dans des modèles d'apprentissage automatique. Quand ces matrices sont grandes, ce qui est de plus en plus le cas dans une ère du "Big Data", le coût du stockage et de calcul lors d'une opération de multiplication devient exorbitant, cela engendre des problèmes majeurs en terme d'apprentissage, de prédiction, d'inférence et de sécurité. La réduction du coût de stockage et de calcul est alors un enjeu décisif en apprentissage automatique et en traitement de signal.

La transformée de Fourier rapide ouvre un chemin prometteur vers la réduction du coût de stockage et calcul des opérateurs linéaires. En effet, la représentation d'un signal **x** de taille *N* dans le domaine de Fourier n'est qu'une transformation linéaire **F**x avec **F** la matrice de Fourier discrète de taille *N* × *N*. Cette transformation nécessite un coût total de calcul de l'ordre $O(N^2)$. Cet ordre indique que, lorsque *N* devient très grand, le coût de calcul devient important. Heureusement, on peut écrire **F** sous la forme d'un produit de log(*N*) matrices, chacune contenant exactement deux éléments non nuls par ligne et par colonne [1]. Cette nouvelle écriture nous permet de faire le produit **F**x sous un coût total de l'ordre $O(N \log(N))$ (le produit matrice/vecteur se fait de droite à gauche).

Écrire une matrice sous la forme d'un produit de matrices creuses s'applique parfaitement à la plupart des transformées usuelles [1], telles que la transformée de Hadamard, la transformée en sinus discrète (DST), la transformée en cosinus discrète (DCT), la transformée de Legendre discrète (DLT), et la convolution discrète.
On appelle toute matrice qui peut s'écrire sous la forme d'une combinaison de matrices creuses une **transformée rapide**.

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – la matrice de Hadamard (à gauche) de taille $N \times N$ avec $N = 16$, le coût du stockage est de l'ordre de $O(N^2 = 256)$. La factorisation par 4 matrices ayant exactement deux éléments non nuls par ligne et par colonne (à droite), le coût du stockage est de l'ordre de $O(2 \cdot N \log_2(N)) = 128)$. Les éléments en vert sont de valeurs nulles, les éléments en violet sont de valeurs -1 et les éléments en jaune sont de valeurs 1.</figcaption></figure>

Les auteurs de [2] proposent une généralisation du principe des transformées rapides à un plus grand nombre de matrices, cette généralisation concerne aussi le masque de parcimonie des facteurs : on généralise la contrainte de deux éléments non nuls par ligne et par colonne à une famille plus grande de matrices parcimonieuses. [2] propose le modèle appelé **FAμST** (Flexible Approximate
Multi-layer Sparse Transforms). Il est formulé comme suit :

$$A = \prod_{j=1}^{J} S_j,$$

avec **A** une matrice de taille $M \times N$ (entrée) et $S_j$ des matrices parcimonieuses (à apprendre). [2] propose une méthode d'apprentissage de ces matrices parcimonieuses, notamment l'algorithme **palm4MSA hiérarchique** [2], et démontre son efficacité sur plusieurs opérateurs. Cependant, il est évident qu'écrire tout opérateur sous la forme d'un produit de matrices creuses est impossible, l'introduction d'un nouveau modèle plus expressif nous semble alors nécessaire.

On propose dans ce rapport l'étude du modèle suivant :

$$A = \sum_{k=1}^{K} \prod_{j=1}^{J_k} S_{kj}.$$

avec $S_{kj}$ des matrices parcimonieuses à apprendre. Ce modèle peut être vu comme une extension du modèle **FAμST**: pour $K = 1$, on retrouve le modèle **FAμST**. On espère, avec ce nouveau modèle, obtenir de meilleures approximations des opérateurs et retrouver de nouvelles transformées rapides. On appelle ce nouveau modèle **Σ-FAμST**<fnref n="1" /> (Flexible Approximate Sum of Multi-layers Sparse Transforms).

L'approximation d'une matrice dense **A** de taille $M \times N$ par un **Σ-FAμST** vise à réduire le coût du stockage et de calcul lors d'une multiplication matrice/vecteur. On utilise la complexité relative (**CR**) comme une mesure du gain :

**Definition 1** — La complexité relative du stockage d'un **Σ-FAμST** (**CRS**) est le rapport entre le nombre d'éléments non nuls dans le **Σ-FAμST** et le nombre d'éléments non nuls de l'entrée **A** :

$$CRS = \frac{s_{tot}}{\|A\|_0},$$

$$\text{avec } s_{tot} = \sum_{k=1}^{K} \sum_{j=1}^{J} \|S_{kj}\|_0.$$

— La complexité relative de calcul d'un **Σ-FAμST** (**CRC**) est une mesure du gain lors d'une multiplication **Σ-FAμST**/vecteur.

$$CRC^2 = \frac{s_{tot} + K.N}{\|A\|_0}.$$

Quand les dimensions de **A** sont grandes et $K$ est petit, $\frac{KN}{\|A\|_0}$ tends vers 0, on retrouve alors **CRC** ≈ **CRS**, on utilise dans ce rapport alors une seule mesure du gain de calcul et de stockage appelée complexité relative (**CR**) :

$$CR = \frac{s_{tot}}{\|A\|_0}. \qquad (1)$$

<footnote n="1">1. Pour une meilleure expérience scientifique, on garde les mêmes notations que celles utilisées dans [2]. On rajoute des "Σ" ou "s" afin de référer à la somme pour toute nouvelle appellation</footnote>

<footnote n="2">2. La multiplication **Σ-FAμST**/vecteur se fait d'abord par la multiplication de chaque terme/vecteur, cela nécessite $s_{tot}$ opérations. Additionner ensuite les $K$ vecteurs obtenus se fait en $K.N$ opérations</footnote>
En plus de la réduction du coût de stockage et calcul, le modèle Σ-FAμST bénéfice de tous les avancages des modèles parcimonieux lors d'un apprentissage automatique : une meilleure interprétation statistique des paramètres et une meilleure inférence (évitement de sur-apprentissage).

Les algorithmes d'apprentissage proposés dans ce rapport se basent principalement sur l'algorithme d'optimisation **PALM** [3], on commence alors dans la section 2 par une brève présentation de cet algorithme. Ensuite on présente dans la section 3 le modèle FAμST ainsi que les algorithmes d'apprentissage proposés dans [2]. Dans la suite on présente dans la section 4 les motivations théoriques et expérimentales qui nous ont menés à proposer le modèle Σ-FAμST, puis des algorithmes d'apprentissage dans la section 5. On finira par présenter des résultats expérimentaux dans la section 6.

# 2 L'algorithme PALM

On présente dans cette partie l'algorithme **PALM** (Proximal alternating linearized minimization) [3]. Cet algorithme est au cœur de tous les algorithmes d'apprentissage présentés dans ce rapport, que ce soit les algorithmes du modèle FAμST présenté dans [2], ou les algorithmes du modèle Σ-FAμST présenté dans la section 5. On renvoie le lecteur à [3] pour plus de détails.

## 2.1 Problème d'optimisation

Soit le problème d'optimisation non convexe/non différentiable suivant :

$$ \underset{x \in \mathbb{R}^n, y \in \mathbb{R}^m}{\text{Minimiser}} \quad \psi(x, y) = H(x, y) + f(x) + h(y), \qquad (2) $$

avec $H$ une fonction différentiable, $\nabla_x H$ et $\nabla_y H$ sont globalement lipchitziennes, $f$ et $g$ sont propres, non nécessairement convexes, et semi-continues inférieurement et $\psi$ est semi-algébriques dans le sens de la définition 5 de [3]. On note $L_x(\text{resp. } L_y)$ la constante de Lipschitz de $\nabla_x H$ (resp. $\nabla_y H$). Le choix de deux blocs de variables est fait dans le but de simplifier l'explication. Tous les résultats sont valables pour un nombre fini de blocs de variables.

Le problème est non-convexe, et on ne dispose actuellement d'aucun algorithme qui garantit la convergence vers un minimum global. Cependant, l'algorithme **PALM** garantit la convergence vers un minimum local.

## 2.2 Algorithme

La présentation de l'algorithme nécessite d'abord l'introduction de la notion de l'opérateur proximal.

On définit l'opérateur proximal d'une fonction $J: \mathbb{R}^n \to \mathbb{R}$ par :

$$ \text{prox}_{J,\lambda}(u) = \underset{v \in \mathbb{R}^n}{\text{Minimiser}} \left( J(v) + \frac{1}{\lambda} \|v - u\|_2^2 \right) $$

Dans le cas de la fonction nulle : $J: u \mapsto 0, \forall u$, on a :

$$ \text{prox}_{J,\lambda}(u) = u. $$
Dans le cas de la fonction indicatrice d'un ensemble : $J(u) = \delta_{\mathcal{E}}(u)$, on a :

$$ \mathrm{prox}_{J,\lambda}(u) = \mathcal{P}_{\mathcal{E}}(u), $$

avec $\mathcal{P}_{\mathcal{E}}$ l'opérateur de projection orthogonale sur l'ensemble $\mathcal{E}$.

L'algorithme 1 présente l'algorithme **PALM**. L'algorithme suit les étapes suivantes :

- sur $x$, effectuer la descente de gradient sur la partie différentiable de $\psi$

- sur $x$, effectuer une mise à jour proximale sur la partie non différentiable de $\psi$

- sur $y$, effectuer la descente de gradient sur la partie différentiable de $\psi$

- sur $y$, effectuer une mise à jour proximale sur la partie non différentiable de $\psi$

### Algorithme 1 : PALM

**Entrée** : Initialisation $x^0 \in \mathbb{R}^n$, $y^0 \in \mathbb{R}^m$

1. **pour** $i = 0$ à $n_{max} - 1$ **faire**

2. Prendre $c_i > L_1(y^i)$

3. $x^{i+1} \in \mathrm{prox}_{f,c_i}(x^k - \frac{1}{c_i^2} \nabla_x H(x^k, y^k))$

4. Prendre $d_i > L_2(x^{i+1})$

5. $y^{i+1} \in \mathrm{prox}_{g,d_i}(y^k - \frac{1}{d_i^2} \nabla_y H(x^{k+1}, y^k))$

6. **fin**

**Sortie** : $x^{n_{max}}, y^{n_{max}}$

Dans la suite **PALM** sera exploité pour résoudre le problème suivant :

$$ \text{Minimiser } H(x_1, \dots, x_J) + \sum_{j=1}^{J} \delta_{\varepsilon_j}(x_j) \qquad (3) $$

L'ordre de mises à jour des variables dans **PALM** n'a pas d'importance pour la convergence de l'algorithme. Néanmoins, l'algorithme converge pour chaque ordre vers un minimum local différent. Dans ce rapport, et pour des raisons de cohérence et de clarté, la mise à jour des variables sera faite suivant l'ordre écrit dans la fonction objectif, par exemple dans (3) on met à jour $x_1$ puis $x_2$ et ainsi de suite.

## 3 Le modèle FAμST

On présente dans cette partie une méthode d'approximation d'une matrice $\mathbf{A} \in \mathbb{K}^{N \times M}$ avec $\mathbb{K} = \mathbb{R}$ ou $\mathbb{C}$ par un produit de matrices creuses :

$$ \mathbf{A} \approx \prod_{j=1}^{J} \mathbf{S}_j, $$

avec $\mathbf{S}_j \in \mathbb{K}^{a_{j+1} \times a_j}$, $a_1 = M$ et $a_{J+1} = N$. $\mathbf{S}_1$ représente le dernier terme, et $\mathbf{S}_J$ le premier. Cette notation "non usuelle" est liée au fait que le modèle **FAμST** vise la réduction du coût d'une transformation linéaire, sachant que le produit matrice/vecteur se fait de droite à gauche. Une telle approximation est nommée **FAμST**. La qualité de l'approximation est évaluée avec l'erreur relative :
$$ \frac{\|\mathbf{A} - \mathbf{F}\mathbf{A}_{\mu}\mathbf{S}\mathbf{T}\|_{F}}{\|\mathbf{A}\|_{F}} $$

Sans perdre de généralité, et pour simplifier le langage, on note dans toute la suite $\mathbb{K} = \mathbb{R}$. Le lecteur peut se référer à [2] pour plus de détails.

## 3.1 Fonction objectif

Apprendre une telle décomposition passe naturellement par le problème suivant :

$$ \underset{\mathbf{S}_1, \dots, \mathbf{S}_J}{\text{Minimiser}} \left\| \mathbf{A} - \prod_{j=1}^{J} \mathbf{S}_j \right\|_{F}^{2} + \sum_{j=1}^{J} \delta_{\mathbf{S}_j}(\mathbf{S}_j) \quad (4) $$

avec $\delta_{\mathcal{E}_j}$ la fonction indicatrice qui vaut 0 quand $\mathbf{S}_j \in \mathcal{S}_j$ et $+\infty$ sinon; $\mathcal{S}_j$ représente un ensemble de contraintes du facteur $\mathbf{S}_j$. Par exemple, $\mathcal{S}_j$ peut être l'ensemble des matrices ayant exactement $s_j$ éléments non nuls, $\mathcal{S}_j = \{\mathbf{S} \in \mathbb{R}^{a_{j+1} \times a_j}, \|\mathbf{S}\|_0 = s_j\}$. Le 1er terme de (4) garantit la fidélité aux données, le 2ème terme garantit la parcimonie des facteurs (l'appartenance à l'espace de contraintes). Pour éviter les ambiguïtés d'échelle qui apparaissent naturellement lorsque les ensembles de contraintes sont (positivement) homogènes, il est courant [4], [5] de normaliser les facteurs et d'introduire un scalaire $\lambda$ dans le terme de fidélité des données. Pour cela, on introduit l'ensemble des facteurs unitaires $\mathcal{N}_j = \{\mathbf{S} \in \mathbb{R}^{a_{j+1} \times a_j}, \|\mathbf{S}\|_F = 1\}$. L'ensemble de contraintes devient alors $\mathcal{E}_j = \mathcal{S}_j \cap \mathcal{N}_j$, et la fonction objectif devient alors :

$$ \underset{\mathbf{S}_1, \dots, \mathbf{S}_J, \lambda}{\text{Minimiser}} \left\| \mathbf{A} - \lambda \prod_{j=1}^{J} \mathbf{S}_j \right\|_{F}^{2} + \sum_{j=1}^{J} \delta_{\mathcal{E}_j}(\mathbf{S}_j) \quad (5) $$

Il est immédiat de voir une similitude entre (5) et (3). On utilise, donc, **PALM** pour résoudre ce problème. (5) vérifie les conditions de convergence de **PALM** [2] et donc on a la garantie d'obtenir un minimum local.

## 3.2 Une précision sur la projection

La projection d'une matrice $\mathbf{A}$ de taille $M \times N$ sur un espace de contraintes $\mathcal{E}$ constitue une étape cruciale dans l'algorithme **PALM**. Cette projection doit être non seulement connue théoriquement mais aussi simple à calculer.

Dans ce rapport, on projette principalement sur les espaces de contraintes suivants :

$$ \mathcal{E}_{s}^{M \times N} = \{ \mathbf{S} \in \mathbb{R}^{M \times N}, \|\mathbf{S}\|_{0} = s, \|\mathbf{S}\|_{F} = 1 \}, \quad (6) $$

$\mathcal{E}_s^{M \times N}$ est l'ensemble des matrices unitaires de taille $M \times N$ ayant exactement $s$ éléments non nuls.

$$ \mathcal{L}_{s}^{M \times N} = \{ \mathbf{S} \in \mathbb{R}^{M \times N}, \forall i \in \{1, \dots, M\} \|\mathbf{S}_{i}\|_{0} = s, \|\mathbf{S}\|_{F} = 1 \}, \quad (7) $$

$\mathcal{L}_s^{M \times N}$ est l'ensemble des matrices unitaires de taille $M \times N$ ayant exactement $s$ éléments non nuls par ligne.

$$ \mathcal{C}_{s}^{M \times N} = \{ \mathbf{S} \in \mathbb{R}^{M \times N}, \forall j \in \{1, \dots, N\} \|\mathbf{S}_{j}\|_{0} = s, \|\mathbf{S}\|_{F} = 1 \}, \quad (8) $$
$\mathcal{C}_s^{M \times N}$ est l'ensemble des matrices unitaires de taille $M \times N$ ayant exactement $s$ éléments non nuls par colonne.

$$ \mathcal{R}_s^{M \times N} = \{ \mathbf{S} \in \mathbb{R}^{M \times N}, \forall i, j \; \|\mathbf{S}_i\|_0 = \|\mathbf{S}_j\|_0 = s, \|\mathbf{S}\|_F = 1 \}, \quad (9) $$

$\mathcal{R}_s^{M \times N}$ est l'ensemble des matrices unitaires de taille $M \times N$ ayant exactement $s$ éléments non nuls par ligne et par colonne, on appelle ces matrices des matrices s-régulières. C'est cet ensemble qui sera utilisé en générale dans la suite, car ça correspond aux transformées rapides connues [1]. Sauf mention contraire, on utilise l'abus de notation : $\mathcal{E}_s$, $\mathcal{L}_s$, $\mathcal{C}_s$, $\mathcal{R}_s$.

### 3.2.1 Projection sur $\mathcal{E}_s$

La projection de $\mathbf{A}$ sur l'espace $\mathcal{E}_s$ se fait en 2 étapes [2] :

- Normaliser $\mathbf{A}$ : $\mathbf{A}_u = \frac{\mathbf{A}}{\|\mathbf{A}\|_F}$

- Garder les $s$ plus grand éléments de $\mathbf{A}_u$ ayant les plus petits indices<fnref n="3" />(on compare le positionnement des éléments sur la ligne, et après le positionnement sur la colonne), tous les autres éléments sont remis à 0.

### 3.2.2 Projection sur $\mathcal{L}_s$ et $\mathcal{C}_s$

La projection sur les espaces $\mathcal{L}_s$ (resp. $\mathcal{C}_s$) se fait de la même façon, sauf qu'on garde les $s$ plus grands éléments par ligne (resp. par colonne) ayant les plus petits indices.

### 3.2.3 Projection sur $\mathcal{R}_s$

On définit la matrice union des deux matrices $\mathbf{A}$ et $\mathbf{B}$ comme suit :

$$ \mathbf{A}, \mathbf{B} \in \mathbb{R}^{N \times M} \quad (\mathbf{A} \cup \mathbf{B})_{ij} := \begin{cases} a_{ij} & \text{si } a_{ij} \ge b_{ij} \\ b_{ij} & \text{sinon} \end{cases} $$

La projection de $\mathbf{A}$ sur $\mathcal{R}_s$ est l'union de la projection de $\mathbf{A}$ sur $\mathcal{L}_s$ et sur $\mathcal{C}_s$.

Cette projection n'est pas exacte, la matrice projetée n'est pas forcément s-régulière. [6] et [7] proposent une méthode de projection exacte. Néanmoins, cette méthode simple donne une matrice "presque" s-régulière dans peu de temps comparée aux autres approches exactes. En effet les approches exactes repose sur des algorithmes plus complexes, et donc le temps de calcul devient très grand quand la taille de la matrice est grande.

## 3.3 L'algorithmes palm4MSA

PALM appliquée à (5) s'appelle **palm4MSA** (PALM for Multi-layer Sparse Approximations).

On peut décrire une seule itération de cet algorithme comme suit :

- une mise à jour alternée de tous les facteurs $\mathbf{S}_j$ par une descente de gradient suivie d'une projection $\mathcal{P}_{\mathcal{E}_j}$ sur l'espace $\mathcal{E}_j$;

- une mise à jour du facteur d'échelle par une descente de gradient (aucune contrainte n'est posée sur $\lambda$, la fonction indicatrice est nulle, et l'opérateur proximal est alors l'identité).

Ces principales étapes sont illustrées dans l'algorithme 2. En pratique, on utilise souvent un seuil comme critère d'arrêt : dès que la différence de deux erreurs relatives issus de deux itérations

<footnote n="3">3. Le choix des éléments ayant le plus petits indices garantit l'unicité de la projection</footnote>
## Algorithm 2 : palm4MSA

**Entrée** : **A** : Opérateur.

**J** : Nombre de facteurs.

`{E_j}_{j=1}^J` : L'ensemble des contraintes des facteurs.

`{S_j^0}_{j=1}^J` : Facteurs initiaux.

`λ^0` : valeurs initiales du facteur d'échelle.

`n_max` : critère d'arrêt, par exemple, nombres d'itérations.

1 **pour** i = 0 à n_max - 1 **itérations de la descente de gradient faire**

2 **pour** j = 1 à J facteurs faire

3 **L** ← $\prod_{l=j+1}^{J} S_{l}^{i}$

4 **R** ← $\prod_{l=1}^{j-1} S_{l}^{i+1}$

5 Prendre $c_j^i > (\lambda^i \|L_2\|_2 \cdot \|R_2\|_2)^2$

6 **Z** ← $S_j^i - \frac{1}{c_j^i} \lambda^i \mathbf{L}^T (\lambda^i \mathbf{L} S_j^i \mathbf{R} - \mathbf{A}) \mathbf{R}^T$

7 $S_j^{i+1} \leftarrow \mathcal{P}_{\mathcal{E}_j}(\mathbf{Z})$

8 **fin**

9 **Â** ← $\prod_{j=1}^{J} S_{j}^{i+1}$

10 $\lambda^{i+1} \leftarrow \frac{\mathrm{Tr}((\mathbf{A})^T \mathbf{A})}{\mathrm{Tr}(\mathbf{A}^T \mathbf{A})}$

11 **fin**

**Sortie** : La factorisation estimée : $\lambda, \{S_j\}_{j=1}^J$

consécutives est inférieur à un seuil, l'algorithme s'arrête. Ce critère permet d'éviter des itérations qui n'ont pas un grand impact sur l'erreur relative de l'approximation.

L'algorithme **palm4MSA** affiche de mauvais résultats lors de l'approximation des opérateurs linéaires notamment lors de l'approximation des transformées rapides : Fourier et Hadamard. Les auteurs de [2] proposent alors une méthode hiérarchique d'approximation.

## 3.4 palm4MSA hiérarchique

Au lieu de chercher une factorisation de **A** en **J** facteurs d'un seul coup, cet algorithme hiérarchique consiste à chercher une factorisation de l'entrée **A** en deux facteurs : le premier est parcimonieux, le deuxième est moins parcimonieux (on l'appelle résidu). On décompose le résidu en deux facteurs de la même façon et ainsi de suite. Le dernier résidu est parcimonieux. Après chaque étape, une optimisation globale (ligne 6 de l'algorithme 3) est introduite pour une meilleure fidélité aux données (l'opérateur **A**). On utilise alors l'algorithme **PALM4MSA** $J-1$ fois pour décomposer **A** en **J** facteurs parcimonieux. Cette décomposition hiérarchique permet d'obtenir une factorisation exacte de la matrice de Fourier et Hadamard. L'algorithme 3 présente **palm4MSA hiérarchique**. "default" dans la ligne 3 fait référence à une initialisation "par défaut" proposé par les auteurs de [2] : $\lambda_0 = 1$, $S_1 = I$, $\forall j \ge 2 S_j = 0$. avec **I** une matrice avec des 1 sur la diagonale et des 0 partout ailleurs (la matrice identité dans le cas des matrices carrées).
### Algorithm 3 : palm4MSA hiérarchique

**Entrée** : **A** : Opérateur.

**J** : Nombre de facteurs.

`{ε_j}_{j=1}^{J-1}` : L'ensemble des contraintes des facteurs.

`{ε̃_j}_{j=1}^{J-1}` : L'ensemble des contraintes des résidus.

**n_max** : Nombres d'itérations.

1. $T_0 \leftarrow A$

2. **pour** $l = 0 \text{ à } J - 1$ **itérations** **faire**

3. $$ \lambda', \{\mathbf{F}_2\mathbf{F}_1\} \leftarrow \text{palm4MSA}(T_{l-1}, 2, \{\tilde{\mathcal{E}}_l, \mathcal{E}_l\}, \text{init}=\text{default}, n_{\text{max}}) $$

4. $\mathbf{S}_l \leftarrow \mathbf{F}_1$

5. $\mathbf{T}_l \leftarrow \lambda'\mathbf{F}_2$

6. $$ \lambda, \{\mathbf{T}_l, \{\mathbf{S}_j\}_{j=1}^l\} \leftarrow \text{palm4MSA}(\mathbf{A}, l+1, \{\tilde{\mathcal{E}}_l, \{\mathcal{E}_j\}_{j=1}^l\}, \text{init}=\text{current}, n_{\text{max}}) $$

7. **fin**

8. $\mathbf{S}_J \leftarrow \mathbf{T}_{J-1}$

**Sortie** : La factorisation estimée : $\lambda, \{\mathbf{S}_j\}_{j=1}^J$

## 3.5 Remarques sur les algorithmes

Durant ce stage, on a effectué plusieurs expériences avec les 2 algorithmes ci-dessus, voici les résultats les plus pertinents :

- **palm4MSA hiérarchique** est très sensible à l'opérateur à factoriser : en effet, cet algorithme permet d'obtenir une factorisation exacte de la matrice de Hadamard **H**, Néanmoins, en ajoutant un petit bruit : $\mathbf{H}+10^{-8}\mathbf{M}$ avec **M** une matrice aléatoire, l'algorithme n'est plus capable de retrouver la factorisation, de plus, l'erreur relative est très grande. L'expérience complète est disponible dans l'annexe B.

- Il n'y a aucune famille de matrices où ces 2 algorithmes garantissent des meilleurs résultats systématiquement. On a essayé les classes suivantes : Matrices positives, négatives, orthogonales, symétriques.

## 4 Motivation du modèle Σ-FAμST

L'approximation parcimonieuse d'un opérateur dense peut se faire suivant plusieurs modèles. Au lieu d'apprendre la décomposition d'une matrice **A** sous le modèle **FAμST**, on propose de l'apprendre sous le modèle suivant baptisé **Σ-FAμST** :

$$ \mathbf{A} = \sum_{k=1}^{K} \prod_{j=1}^{J_k} \mathbf{s}_{kj}. $$

Dans tout ce qui suit, on considère $\forall k J_k = J$ pour alléger les notations.

Si le produit est un choix naturel (Transformées rapides, garantie de densité), la somme reste un choix à motiver. Cette section est consacrée à motiver ce choix de modèle avec des appuis théoriques et expérimentaux.
## 4.1 Expressivité

Soit

$$ \mathbf{A} = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 2 \end{bmatrix} $$

La proposition 1 ci-dessous montre qu'il est impossible d'écrire **A** sous la forme d'un produit de deux matrices 2-régulières.

**Proposition 1** Soit $\mathbf{A} \in \mathbb{R}^{4 \times 4}$ tel que $a_{ij} \neq 0, 1 \le i, j \le 4$, et $\mathbf{A} = \mathbf{BC}$ avec $\mathbf{B}, \mathbf{C} \in \mathbb{R}_2$ alors :

$$ \forall (i,j) \in [1,4]^2, \exists (i',j') \in [1,4]^2, i' \neq i, j' \neq j, \frac{a_{ij}}{a_{i'j}} = \frac{a_{ij'}}{a_{i'j'}} $$

*Preuve :* Voir annexe A

En gardant la même complexité relative **CR** (16 éléments non nuls), il est évident qu'on peut écrire **A** sous la forme d'une somme de deux matrices 2-régulières.

$$ \mathbf{A} = \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 2 \end{bmatrix} + \begin{bmatrix} 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \end{bmatrix} $$

Cet exemple nous laisse penser qu'introduire la somme va nous permettre d'avoir un modèle plus expressif. Les théorèmes suivant démontrés sur [8] et [9] soutiennent ces propos.

**Théorème 1** Soit $n \ge 2$ un entier. Toute matrice $\mathbf{A} \in M_n(\mathbb{C})$ peut s'écrire sous la forme d'une somme de matrices appartenant à $O_n(\mathbb{C})$.

**Théorème 2** Soit $n \ge 2$ un entier. Toute matrice $\mathbf{A} \in M_n(\mathbb{R})$ peut s'écrire sous la forme d'une somme de matrices appartenant à $O_n(\mathbb{R})$.

**Théorème 3** Soit $n \ge 2$ un entier, alors $\mathbf{A} \in M_n(\mathbb{C})$ peut s'écrire sous la forme d'une somme de $m$ matrices stables$^4$ si et seulement si $|\mathrm{Tr}(\mathbf{A})| \le nm$

## 4.2 Flexibilité

Motivé par la décomposition des transformées rapides, peut-on trouver une "bonne" approximation d'une matrice quelconque en se basant uniquement sur des facteurs tous $s$-réguliers ?
L'expérience suivante est conçue pour répondre à cette question. On propose l'approximation d'une matrice aléatoire par un **FAμST** contenant uniquement des facteurs $s$-réguliers générés par l'algorithme **palm4MSA hiérarchique**. On commence par une approximation par 2 facteurs, on augmente ensuite le nombre de facteurs et ceci jusqu'à ce que la complexité relative du **FAμST** soit supérieure à $1^5$.

<footnote n="4">4. La partie réelle de toutes valeurs propres est négative</footnote>

<footnote n="5">5. Aucune utilité de l'approximation quand $CR > 1$</footnote>
## Paramètres de l'expérience

- **A** une matrice de taille $N \times N$, tous les éléments suivent une loi normale $N(0, 1)$

- $N = 64$

- Tous les facteurs sont s-réguliers : $s \in \{2, \dots, 10\}$

- Le nombre des facteurs : $J \in \{2, \dots, C\}$ avec $C \times s \times 64 < 64^2$ (CR > 1)

- **Paramètres de l'algorithm e palm4MSA hiérarchique :**

- L'ensemble des contraintes sur les facteurs : $\forall j \quad \varepsilon_j = \mathcal{R}_2$

- L'ensemble des contraintes sur les résidus : $\forall j < J - 2 \quad \tilde{\varepsilon}_j = \mathcal{R}_{\frac{N^2}{2+j}}$, $\tilde{\varepsilon}_{J-1} = \mathcal{R}_2$

<figure><img src="image_4.png" /><figcaption>FIGURE 2 – Approximation d'une matrice aléatoire par un produit de matrices toutes s-régulières</figcaption></figure>

La figure 2 montre les résultats, l'augmentation de nombre de facteurs permet de retrouver une meilleure précision, cependant l'erreur relative reste assez grande, ceci nous incite à chercher à mieux distribuer les facteurs, par exemple au lieu de chercher une approximation suivant **FA**μ**ST** avec *J* facteurs, on peut chercher une approximation suivant le modèle **Σ-FA**μ**ST** avec 2 termes chacun contenant *J*/2 facteurs.

Le nouveau modèle dispose de plus de flexibilité : choix du masque de parcimonie, choix du nombre de termes, choix du nombre de facteurs sur chaque terme. Plusieurs possibilités alors surgissent lors de l'approximation d'un opérateur. ceci promet de meilleurs résultats en terme d'approximation.
# 5 Algorithmes proposés pour le modèle Σ-FAμST

On présente dans cette partie une méthode d'approximation d'une matrice $\mathbf{A} \in \mathbb{R}^{N \times M}$ par une somme de produit de matrices creuses :

$$ \mathbf{A} \approx \sum_{k=1}^{K} \prod_{j=1}^{J} \mathbf{S}_{kj}, $$

avec $\mathbf{S}_{kj} \in \mathbb{R}^{a_{k(j+1)} \times a_{kj}}$, $a_{k1} = M$ et $a_{k(J+1)} = N$ pour tout $k \in \{1, \dots, K\}$. Contrairement à la multiplication, l'addition est une opération commutative, on garde alors la notation usuelle des indices $k$ (l'indice augmente de gauche à droite). On évalue l'approximation avec l'erreur relative :

$$ \frac{\|\mathbf{A} - \Sigma\mathbf{FA}\mu\mathbf{ST}\|_F}{\|\mathbf{A}\|_F}. $$

## 5.1 Formulation du problème

Pour l'apprentissage du modèle Σ-FAμST, le problème d'optimisation à résoudre est le suivant :

$$ \underset{\mathbf{S}_{11}, \dots, \mathbf{S}_{1J}, \lambda_1, \dots, \mathbf{S}_{K1}, \dots, \lambda_K}{\text{Minimiser}} \quad \frac{1}{2} \left\| \mathbf{A} - \sum_{k=1}^{K} \lambda_k \prod_{j=1}^{J} \mathbf{S}_{kj} \right\|_{F}^{2} + \sum_{k=1}^{K} \sum_{j=1}^{J} \delta_{\varepsilon_{kj}}(\mathbf{S}_{kj}). \quad (10) $$

Il est direct de voir une ressemblance entre (10) et (3). En effet, en prenant $K(J+1)$ blocs de variables avec $x_{(k-1)(J+1)+j} = \mathbf{S}_{kj}$ pour $j \in \{1, \dots, J\}$, $k \in \{1, \dots, K\}$, $x_{k(J+1)} = \lambda_k$. $\varepsilon_{k(J+1)} = \mathbb{R}$ pour $k \in \{1, \dots, K\}$ et donc $\delta_{\varepsilon_{k(J+1)}} = 0$ (pas de contrainte sur les scalaires $\lambda_k$).

Pour appliquer **PALM** au problème (10), il suffit alors de vérifier les conditions de convergence : la fonction $H$ est différentiable (car quadratique), les fonctions caractéristiques sont propres, semi-continues inférieurement (par définition), la fonction objectif est semi-algébrique [3, Définition 5], $\nabla_{x_j} H$ est globalement lipschitzienne pour tout $j \in \{1, \dots, K(J+1)\}$ (Démonstration ci-dessous).

En appliquant **PALM** au problème (10), on a une garantie d'atteindre un minimum local. On rappelle que les blocs seront mis à jour suivant leurs ordres d'écriture dans la fonction objectif.

## 5.2 Gradient et constante de Lipschitz

Pour spécifier les itérations de **PALM** spécialisées au problème (10), on fixe l'itération $i$ et le facteur $j$ du terme $k$, $\mathbf{S}_{kj}^i$ est le facteur qu'on veut mettre à jour, on note **F** (Fixed) les termes constants lors de la mise à jours du terme $k$ :

$$ \mathbf{F} = \sum_{l=1}^{k-1} \lambda_l^{i+1} \prod_{j=1}^{J} \mathbf{S}_{lj}^{i+1} + \sum_{l=k+1}^{K} \lambda_l^{i} \prod_{j=1}^{J} \mathbf{S}_{lj}^{i}, $$
avec la convention $\sum_{l \in \emptyset} T_l = 0$. On note respectivement $\mathbf{L}$ la partie gauche (Left) et $\mathbf{R}$ la partie droite (Right) du facteur $\mathbf{S}_{kj}$ du terme $k$ :

$$ \mathbf{L} = \lambda_k^i \prod_{t=j+1}^{J} \mathbf{S}_{kt}^i, \quad \mathbf{R} = \prod_{t=1}^{j-1} \mathbf{S}_{kt}^{i+1}, $$

avec la convention $\prod_{t \in \emptyset} S_t$. Ces notations donnent, lors de la mise à jour du $kj$ème facteur $\mathbf{S}_{kj}^i$ :

$$ \begin{aligned} H(\mathbf{S}_{11}^{i+1}, \dots, \mathbf{S}_{k(j-1)}^{i+1}, \mathbf{S}_{kj}^{i}, \dots, \lambda_K) &= H(\mathbf{F}, \mathbf{L}, \mathbf{S}_{kj}^{i}, \mathbf{R}) \\ &= \frac{1}{2} \|\mathbf{A} - \mathbf{F} - \mathbf{L}\mathbf{S}_{kj}^{i}\mathbf{R}\|_{F}^{2} \end{aligned} $$

Calculons Le gradient de $H$ par rapport au $kj$ème facteur :

$$ \begin{aligned} \nabla_{\mathbf{S}_{kj}^i} H(\mathbf{F}, \mathbf{L}, \mathbf{S}_{kj}^i, \mathbf{R}) &= \frac{1}{2} \nabla_{\mathbf{S}_{kj}^i} \mathrm{Tr}((\mathbf{A} - \mathbf{F} - \mathbf{L}\mathbf{S}_{kj}^i\mathbf{R})^T(\mathbf{A} - \mathbf{F} - \mathbf{L}\mathbf{S}_{kj}^i\mathbf{R})) \\ &= \frac{1}{2} \nabla_{\mathbf{S}_{kj}^i} \mathrm{Tr}((\mathbf{F} - \mathbf{A})^T \mathbf{L} \mathbf{S}_{kj}^i \mathbf{R} + (\mathbf{F} - \mathbf{A})(\mathbf{L} \mathbf{S}_{kj}^i \mathbf{R})^T + \mathbf{R}^T \mathbf{S}_{kj}^i \mathbf{T} \mathbf{L} \mathbf{S}_{kj}^i \mathbf{R}) \\ &= \mathbf{L}^T (\mathbf{F} - \mathbf{A}) \mathbf{R}^T + \mathbf{L}^T \mathbf{L} \mathbf{S}_{kj}^i \mathbf{R} \mathbf{R}^T \\ &= \mathbf{L}^T (\mathbf{F} - \mathbf{A} + \mathbf{L} \mathbf{S}_{kj}^i \mathbf{R}) \mathbf{R}^T \end{aligned} $$

Pour montrer que $\nabla_{\mathbf{S}_{kj}^i} H$ est lipschitzienne, on commence par démontrer la proposition suivante :

**Proposition 2** Pour toute matrice $\mathbf{A} \in \mathbb{K}^{n \times p}$ et $\mathbf{B} \in \mathbb{K}^{p \times m}$, on a :

$$ \|\mathbf{A}\mathbf{B}\|_F \le \|\mathbf{A}\|_2 \|\mathbf{B}\|_F $$

**Preuve** : Par définition de la norme spectrale on a : $\|\mathbf{A}\mathbf{x}\|_2 \le \|\mathbf{A}\|_2 \|\mathbf{x}\|_2$, on note $\mathbf{B} = [\mathbf{b}_1, \dots, \mathbf{b}_m]$ alors, $\|\mathbf{B}\|_F^2 = \sum_{i=1}^m \|\mathbf{b}_i\|_2^2$ et on a :

$$ \|\mathbf{A}\mathbf{B}\|_F^2 = \sum_{i=1}^{n} \|\mathbf{A}\mathbf{b}_i\|_2^2 \le \|\mathbf{A}\|_2^2 \sum_{i=1}^{n} \|\mathbf{b}_i\|_2^2 = \|\mathbf{A}\|_2^2 \|\mathbf{B}\|_F^2 $$

Montrons que $\nabla_{\mathbf{S}_{kj}^i} H$ est lipschitzienne :

$$ \begin{aligned} \left\| \nabla_{\mathbf{S}_{kj}^i} H(\mathbf{F}, \mathbf{L}, \mathbf{S}_1, \mathbf{R}) - \nabla_{\mathbf{S}_{kj}^i} H(\mathbf{F}, \mathbf{L}, \mathbf{S}_2, \mathbf{R}) \right\|_F &= \| \mathbf{L}^T \mathbf{L} (\mathbf{S}_1 - \mathbf{S}_2) \mathbf{R} \mathbf{R}^T \|_F \\ &\le \| \mathbf{L}^T \mathbf{L} \|_2 \| \mathbf{R} \mathbf{R}^T \|_2 \| \mathbf{S}_1 - \mathbf{S}_2 \|_F \\ &= \| \mathbf{L} \|_2^2 \cdot \| \mathbf{R} \|_2^2 \| \mathbf{S}_1 - \mathbf{S}_2 \|_F \end{aligned} $$

Une fois que tous les facteurs $j$ du terme $k$ sont mis à jour, on met à jour le $k$ème facteur d'échelle $\lambda_k$. On pose :

$$ \hat{\mathbf{A}} = \prod_{j=1}^{J} S_{kj}^{i+1} $$
Ces notations donnent, lors de la mise à jour de $\lambda_k^i$ :

$$ 
\begin{aligned}
H(\mathbf{S}_{11}^{i+1}, \dots, \mathbf{S}_{k1}^{i+1}, \dots, \mathbf{S}_{kJ}^{i+1}, \lambda_k^i, \mathbf{S}_{(k+1)1}^i, \dots, \lambda_K^i) &= H(\mathbf{F}, \lambda_k^i, \hat{\mathbf{A}}) \\
&= \frac{1}{2} \|\mathbf{A} - \mathbf{F} - \lambda_k^i \hat{\mathbf{A}}\|_F^2
\end{aligned}
$$

Calculons Le gradient de $H$ par rapport à $\lambda_k^i$

$$ 
\begin{aligned}
\nabla_{\lambda_k^i} H(\mathbf{F}, \lambda_k^i, \hat{\mathbf{A}}) &= \frac{1}{2} \nabla_{\lambda_k^i} \mathrm{Tr}((\mathbf{A} - \mathbf{F} - \lambda_k^i \hat{\mathbf{A}})^T (\mathbf{A} - \mathbf{F} - \lambda_k^i \hat{\mathbf{A}})) \\
&= \frac{1}{2} \nabla_{\lambda_k^i} \mathrm{Tr}(\lambda_k^i (\mathbf{F} - \mathbf{A})^T \hat{\mathbf{A}} + \lambda_k^i \hat{\mathbf{A}}^T (\mathbf{F} - \mathbf{A}) + \lambda_k^i^2 \hat{\mathbf{A}}^T \hat{\mathbf{A}}) \\
&= \mathrm{Tr}((\mathbf{F} - \mathbf{A})^T \hat{\mathbf{A}}) + \lambda_k^i \mathrm{Tr}(\hat{\mathbf{A}}^T \hat{\mathbf{A}})
\end{aligned}
$$

La constante de lipschitz de $\nabla_{\lambda_k^i} H$ est $\mathrm{Tr}(\hat{\mathbf{A}}^T \hat{\mathbf{A}})$.

## 5.3 L'algorithmme palm4SMSA

PALM appliquée à (10) s'appelle **palm4SMSA** (PALM for Sum of Multi-layers Sparse Ap- proximations). Cette algorithme met à jour les termes $k$ un par un de la façon suivante :

- Une mise à jour alternée de tous les facteurs $\mathbf{S}_{kj}$ du terme $k$ par une descente de gradient suivie d'une projection $\mathcal{P}_{\mathcal{E}_j}$ sur l'espace $\mathcal{E}_j$

- Une mise à jour du facteur d'échelle du terme $k$ par une descente de gradient.

On note $\mathbf{T}_k$ (Term) l'ensemble des facteurs du terme $k$ :

$$ \mathbf{T}_k = \{\mathbf{S}_{kj}\}_{j=1}^J. \qquad (11) $$

Ces principales étapes sont illustrées dans l'algorithmme 4. Les lignes 4 à 12 de l'algorithmme **palm4SMSA** n'est que l'algorithme **palm4MSA** appliquée à $\mathbf{A} - \mathbf{F}$. L'algorithme 5 présente une version courte de l'algorithme **palm4SMSA**.

Sur la figure 3, on remarque que la courbe de convergence de la fonction objectif générée par l'algorithme **palm4SMSA** converge dès les premières itérations et vers un mauvais minimum local dans la plus part des cas, ceci est du à la complexité de l'espace où la fonction objectif évolue et aux choix d'initialisations.

## 5.4 Discussion sur l'initialisation

Le minimum local obtenu par l'algorithme **palm4SMSA** dépend directement de l'initialisation choisie, malheureusement, on ne dispose actuellement d'aucune méthode pour choisir une "bonne" initialisation. Soient les initialisations suivantes :

$$ \{\lambda_k^0\}_{k=1}^K, \{\mathbf{T}_k^0\}_{k=1}^K = \{\mathbf{1}\}_{k=1}^K, \{\mathbf{S}_{k1} = \mathbf{I}, \forall j > 2, \mathbf{S}_{kj} = \mathbf{0}\}_{k=1}^K \qquad (12) $$

(12) est une extension de l'initialisation par défaut proposé dans [2], tous les scalaires sont initiés avec 1, le dernier facteur de chaque terme est initié avec $\mathbf{I}$, tous les autres facteurs sont initiés avec $\mathbf{0}$.

$$ \{\lambda_k^0\}_{k=1}^K, \{\mathbf{T}_k^0\}_{k=1}^K = \{\mathbf{1}\}_{k=1}^K, \{\forall j, \mathbf{S}_{kj} = \mathbf{I}\}_{k=1}^K \qquad (13) $$
## Algorithme 4 : palm4SMSA version longue

Entrée : **A** : Opérateur.

K : Nombre de termes.

J : Nombre de facteurs.

$\{\mathcal{E}_k\}_{k=1}^K : \text{L'ensemble des contraintes, avec } \mathcal{E}_k = \{\mathcal{E}_{kj}\}_{j=1}^J$.

$\{\mathbf{T}_k^0\}_{k=1}^K : \text{Facteurs initiaux, avec } \mathbf{T}_k^0 = \{\mathbf{S}_{kj}^0\}_{j=1}^J$.

$\{\lambda_k^0\}_{k=1}^K : \text{valeurs initiales.}$

$n_{\max}$ : Nombres d'itérations.

1 pour $i = 0 \text{ à } n_{\max} - 1$ itérations de la descente de gradient faire

2 pour $k = 1 \text{ à } K$ faire

3 $\mathbf{F} \leftarrow \sum_{l=1}^{k-1} \lambda_l^{i+1} \prod_{j=1}^{J} \mathbf{S}_{lj}^{i+1} + \sum_{l=k+1}^{K} \lambda_l^i \prod_{j=1}^{J} \mathbf{S}_{lj}^i$

4 pour $j = 1 \text{ à } J$ faire

5 $\mathbf{L} \leftarrow \lambda_k^i \prod_{t=j+1}^{J} \mathbf{S}_{kt}^i$

6 $\mathbf{R} \leftarrow \prod_{t=1}^{j-1} \mathbf{S}_{kt}^{i+1}$

7 Prendre $c_{kj}^i > \|L_2\|_2^2 \cdot \|R_2\|_2^2$

8 $\mathbf{Z} \leftarrow \mathbf{S}_{kj}^i - \frac{1}{c_{jk}^i} \mathbf{L}^T (\mathbf{F} - \mathbf{A} + \mathbf{L} \mathbf{S}_{kj}^i \mathbf{R}) \mathbf{R}^T$

9 $\mathbf{S}_{jk}^{i+1} \leftarrow \mathcal{P}_{\mathcal{E}_{kj}}(\mathbf{Z})$

10 fin

11 $\hat{\mathbf{A}} \leftarrow \prod_{j=1}^{J} \mathbf{S}_{kj}^{i+1}$

12 $\lambda_k^{i+1} \leftarrow \frac{\text{Tr}((\mathbf{A}-\mathbf{F})^T \hat{\mathbf{A}})}{\text{Tr}(\hat{\mathbf{A}}^T \hat{\mathbf{A}})}$

13 fin

14 fin

Sortie : $\{\lambda_k\}_{k=1}^K, \{\mathbf{T}_k\}_{k=1}^K$

## Algorithme 5 : palm4SMSA version courte

1 pour $i = 0 \text{ à } n_{\max} - 1$ itérations faire

2 pour $k = 1 \text{ à } K$ faire

3 $\mathbf{F} = \sum_{l=1}^{k-1} \lambda_l^{i+1} \prod_{j=1}^{J} \mathbf{S}_{lj}^{i+1} + \sum_{l=k+1}^{K} \lambda_l^i \prod_{j=1}^{J} \mathbf{S}_{lj}^i$

4 fin

5 fin

Dans (13), tous les scalaires sont initiés avec 1, le dernier facteur de chaque terme est initié avec **I**, tous les autres facteurs sont initiés avec **0**.

$$ \{\lambda_k^0\}_{k=1}^K, \{\mathbf{T}_k^0\}_{k=1}^K = \{1\}_{k=1}^K, \{\mathbf{S}_{k1} = \mathbf{A}, \forall j > 2, \mathbf{S}_{kj} = \mathbf{0}\}_{k=1}^K \quad (14) $$

Dans (14), tous les scalaires sont initiés avec 1, le dernier facteur de chaque terme est initié avec **A** l'opérateur à approximer, tous les autres facteurs sont initiés avec **0**.

Ces initialisations donnent les meilleurs approximations en terme d'erreur relative sur toutes les
<figure><img src="image_5.png" /><figcaption>FIGURE 3 – La courbe de convergence (Erreur relative/itération) générée par l’algorithme **Palm4SMSA** associé au modèle $\Sigma$-FA$\mu$ST. La courbe provient de l’approximation d’une matrice aléatoire de taille 100 × 100 par le $\Sigma$-FA$\mu$ST $\sum_{k=1}^{3} \prod_{j=1}^{2} S_{kj}$, avec $S_{kj}$ sont tous 2-régulières expériences effectuées.</figcaption></figure>

La figure 4 montre l'erreur relative de l'approximation d'une matrice aléatoire qui est issue d'une initialisation aléatoire (tous les facteurs sont aléatoires), les erreurs relative sont très grandes et l'approximation est alors très pauvre, les box-plots montrent aussi une grande variance de l'erreur relative, ce qui indique une sensibilité de l'algorithme **palm4SMSA** à l'initialisation.

On garde l'initialisation (12) comme une initialisation par "défaut" de l'algorithme **palm4SMSA**.

<figure><img src="image_6.png" /><figcaption>FIGURE 4 – Boxplot représentant l’erreur relative de l’approximation d’une matrice aléatoire A par le FASμST $\sum_{k=1}^{4} \prod_{j=1}^{3} S_{kj}$, avec $S_{kj}$ sont tous 2-régulières, toutes les initialisations de l’algorithme **palm4SMSA** sont des matrices aléatoires dont les éléments suivent une loi normale $N(0, 1)$ (à gauche) et une loi uniforme $U([0, 1])$ (à droite).</figcaption></figure>
Étant donné le manque d'un protocole théorique ou expérimental de désignation d'une initialisation offrant des "bons" minimas locaux systématiquement, on propose d'autres algorithmes d'apprentissage basés sur une approche hiérarchique.

## 5.5 Palm4SMSA hiérarchique

On propose pour résoudre le problème (10) une méthode hiérarchique : Au lieu de retrouver une approximation de A en K termes, on commence par une approximation en 1 terme, ce terme sera ensuite approximé par un autre terme et ainsi de suite. On réduit ainsi la complexité de l'espace où la fonction objectif évolue et on espère retrouver alors des meilleurs minima locaux.

Notre méthode consiste à factoriser une matrice A en J facteurs avec l'algorithme **palm4MSA** ou **palm4MSA hiérarchique**$^6$. Le résidu $A - \text{palm4MSA}(A)$ sera encore factoriser et ainsi de suite. Dans chaque étape, on introduit une optimisation globale à l'aide de l'algorithme **palm4SMSA** pour assurer une meilleure fidélité aux données.

Cette méthode bénéficie de la propriété suivante : quand on augmente le nombre de termes, on augmente la précision. Cette propriété fera le sujet de l'expérience 6.3.

On appelle l'algorithme utilisant **palm4MSA** par **palm4SMSA hiérarchique p**, et l'algorithme utilisant **palm4MSA hiérarchique** par **palm4SMSA hiérarchique h**. Pour une écriture réduit de ces algorithmes, on introduit d'abord la fonction X :

$$X : \{\lambda_l\}_{l=1}^k, \{\mathbf{T}_i\}_{i=1}^k \mapsto \sum_{k=1}^k \lambda_k \prod_{j=1}^J \mathbf{s}_{kj} \qquad (15)$$

Les algorithmes proposés dans cette partie feront l'objet de plusieurs expériences dans la partie suivante.

<footnote n="6">6. Quand $J=2$, les 2 algorithmes sont équivalent</footnote>
## Algorithme 6 : palm4SMSA hiérarchique p

**Entrée :** A : Opérateur.

K : Nombre de termes.

J : Nombre de facteurs.

$\{\mathcal{E}_k\}_{k=1}^K : \text{L'ensemble des contraintes, avec } \mathcal{E}_k = \{\mathcal{E}_{kj}\}_{j=1}^J$.

$\{\mathbf{T}_k^0\}_{k=1}^K : \text{Facteurs initiaux, avec } \mathbf{T}_k^0 = \{\mathbf{S}_{kj}^0\}_{j=1}^J$.

$\{\lambda_k^0\}_{k=1}^K : \text{valeurs initiales.}$

$n_{\max}$ : Nombres d'itérations.

1 pour $k=1$ à $K$ faire

2 $\lambda_k, \mathbf{T}_k \leftarrow \text{Palm4MSA}(\mathbf{A} - X(\{\lambda_l\}_{l=1}^{k-1}, \{\mathbf{T}_l\}_{l=1}^{k-1}), J, \mathcal{E}_k, \mathbf{T}_k^0, \lambda_k^0, n_{\max})$

3 si $k > 1$ alors

4 $\{\lambda_l\}_{l=1}^k, \{\mathbf{T}_l\}_{l=1}^k \leftarrow \text{PALM4SMSA}(\mathbf{A}, k, J, \{\mathcal{E}_l\}_{l=1}^k, \{\mathbf{T}_l\}_{l=1}^k, \{\lambda_l\}_{l=1}^k, n_{\max})$

5 fin

6 fin

Sortie : $\{\lambda_k\}_{k=1}^K, \{\mathbf{T}_k\}_{k=1}^K$

## Algorithme 7 : palm4SMSA hiérarchique h

**Entrée :** A : Opérateur.

K : Nombre de termes.

J : Nombre de facteurs.

$\{\tilde{\mathcal{E}}_k\}_{k=1}^K : \text{L'ensemble des contraintes sur les facteurs, avec}$

$\tilde{\mathcal{E}}_k = \{\tilde{\mathcal{E}}_{kj}\}_{j=1}^{J-1}$.

$\{\mathcal{E}_k\}_{k=1}^K : \text{L'ensemble des contraintes sur les résidus, avec}$

$\mathcal{E}_k = \{\mathcal{E}_{kj}\}_{j=1}^{J-1}$.

$n_{\max}$ : Nombres d'itérations.

1 pour $k=1$ à $K$ faire

2 $\lambda_k, \mathbf{T}_k \leftarrow \text{Palm4MSA hiérarchique}(\mathbf{A} - X(\{\lambda_l\}_{l=1}^{k-1}, \{\mathbf{T}_l\}_{l=1}^{k-1}), J, \tilde{\mathcal{E}}_k, \tilde{\mathcal{E}}_k, n_{\max})$

3 $\tilde{\mathcal{E}}_k = \tilde{\mathcal{E}}_k \cup \{\mathcal{E}_{kJ-1}\}$

4 si $k > 1$ alors

5 $\{\lambda_l\}_{l=1}^k, \{\mathbf{T}_l\}_{l=1}^k \leftarrow \text{PALM4SMSA}(\mathbf{A}, k, J, \{\tilde{\mathcal{E}}_k\}_{k=1}^k, \{\mathbf{T}_k\}_{k=1}^k, \{\lambda_l\}_{l=1}^k, n_{\max})$

6 fin

7 fin

Sortie : $\{\lambda_k\}_{k=1}^K, \{\mathbf{T}_k\}_{k=1}^K$

# 6 Expériences

On présente dans cette partie plusieurs expériences visant à évaluer le modèle **Σ-FAμST** par le biais des algorithmes d'apprentissage proposés dans la section précédente. On propose en tout quatre expériences chacune éclairant un aspect différent de ce modèle. L'expérience 6.1 répond à la question posée dans la section 4 sur une possible meilleure distribution des facteurs. L'expérience 6.2 montre la difficulté d'apprendre des transformées rapides. L'expérience 6.3 propose une méthode d'approximation des matrices de grande taille, et l'expérience 6.4 propose une application sur des
données réelles.

La plupart des expériences sont faites sur des matrices aléatoires, un tel choix est motivé par une
envie de généralisation des conclusions.

## 6.1 Vers une meilleure distribution des facteurs

Cette expérience est conçue pour comparer les performances des deux modèles FAμST et Σ-FAμST sachant qu'on a la même complexité relative : on dispose d'un budget de parcimonie à distribuer, soit en J facteurs, soit en K termes chacun contenant J/K facteurs.

Soit $\mathbf{A}$ une matrice, on propose l'approximation de $\mathbf{A}$ selon les 2 modèles suivants :

$$ \mathbf{A} \approx \prod_{i=1}^{J} \mathbf{S}_i, \quad \mathbf{A} \approx \sum_{k=1}^{K} \prod_{j=1}^{J/K} \mathbf{S}_{kj}. $$

Les 2 modèles ont le même nombre de facteurs, si on pose la même contrainte de parcimonie sur
tous les facteurs lors de l'apprentissage, les approximations obtenues auront la même complexité
relative (CR).

L'expérience est faite sur une matrice aléatoire $\mathbf{A}$ de taille $N \times N$. En faisant varier $K$, on compare
le modèle FAμST contenant $K \log(N)$ facteurs et le modèle Σ-FAμST contenant $K$ termes chacun
contenant $\log(N)$ facteurs. On s'arrête quand $\text{CR} > 1$.

On utilise l'algorithmie palm4MSA et sa version hiérarchique pour obtenir l'approximation FAμST,
les algorithmes palm4SMSA hiérarchique p et h seront utilisés pour obtenir l'approximation
Σ-FAμST. Les initialisations sont toutes aléatoires (tous les facteurs sont aléatoires). Les facteurs
sont tous 2-réguliers.

### Paramètres de l'expérience

- $\mathbf{A}$ : tous les éléments suivent une loi normale $N(0, 1)$

- Taille de la matrice : $N \in \{32, 64\}$

- $K \in \{2, \dots, C\}$ tel que $\frac{C \cdot 2N \cdot \log(N)}{N^2} < 1$

- Tous les facteurs sont de taille $N \times N$

- **Paramètres de l'algorithmie palm4MSA :**

- Nombre de facteurs : $J = K \cdot \log(N)$; Initialisation aléatoire

- **Paramètres de l'algorithmie palm4MSA hiérarchique :**

- L'ensemble des contraintes sur les facteurs : $\forall j \in \mathcal{E}_j = \mathcal{R}_2$

- L'ensemble des contraintes sur les résidus : $\forall j < J - 2 \in \tilde{\mathcal{E}}_j = \mathcal{R}_{\frac{N^2}{2+j}}, \tilde{\mathcal{E}}_{J-1} = \mathcal{R}_2$

- **Paramètres de l'algorithmie palm4MSA hiérarchique p :**

- Nombre de termes : $K$; Nombre de facteurs : $\log(N)$; Initialisation aléatoire

- **Paramètres de l'algorithmie palm4SMSA hiérarchique h :**

- L'ensemble des contraintes sur les facteurs : $\forall k \forall j \in \mathcal{E}_{kj} = \mathcal{R}_2$

- L'ensemble des contraintes sur les résidus : $\forall k \forall j \in \tilde{\mathcal{E}}_{kj} = \mathcal{R}_{\frac{N^2}{2j}}$

**Résultat** : La figure 5 montre les résultats. En rajoutant beaucoup de facteurs à l'algorithmie palm4MSA, on ne fait qu'augmenter l'erreur relative, l'algorithmie palm4MSA hiérarchique affiche des résultats meilleurs, mais le compromis complexité/erreur est très mauvais. L'algorithmie
palm4SMSA hiérarchique h propose la meilleure approximation de A pour tout N et pour tout K, ceci montre qu'il vaut mieux utiliser le modèle Σ-FAμST pour l'approximation des opérateurs au lieu d'augmenter la complexité en rajoutant des facteurs au modèle FAμST. Une autre expérience avec $J = K.\lfloor \frac{\log(N)}{2} \rfloor$ a été faite, nos conclusions restent les mêmes.

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Performance des algorithmes dans le contexte d'approximation d'une matrice aléatoire de taille 32 (à gauche), 64 (à droite). L'abscisse concernant les algorithmes du modèle Σ-FAμST signifie : K termes chacun contenant log(N) facteurs. Les résultats sont en moyenne sur 30 réalisations, la variance étant de l'ordre de 10⁻³. Ces algorithmes montrent une insensibilité à l'initialisation.</figcaption></figure>

## 6.2 Sur la difficulté d'apprendre des transformées rapides

Soit A une matrice dense qui peut s'écrire sous la forme d'une somme d'un produit de matrices parcimonieuses :

$$A = \sum_{k=1}^{K} \prod_{j=1}^{J} S_{kj}, \qquad (16)$$

0 représente alors le minimum global de la fonction objectif (10). Cette expérience essaye de répondre à la question suivante : les algorithmes du modèle Σ-FAμST peuvent-ils atteindre ce minimum global ou au moins proposer de bonnes approximations en terme d'erreur relative ?

L'expérience est faite sur la matrice A :

$$A = F + H, \qquad (17)$$

avec F la matrice de Fourier de taille N et H la matrice de Hadamard de taille N avec N une puissance de 2. Ces deux matrices peuvent s'écrire sous la forme de log(N) matrices 2-régulières, et donc :

$$A = \sum_{k=1}^{2} \prod_{j=1}^{\log(N)} S_{kj}, \qquad (18)$$
## Paramètres de l'expérience

- Taille de la matrice : $N \in \{32, 64, 128\}$

- Tous les facteurs sont de taille $N \times N$

- Nombre de termes : $K = 2$; Nombre de facteurs : $J = \log(N)$

- **Paramètres des algorithmes palm4SMSA et palm4SMSA hiérarchique p :**

- L'ensemble des contraintes : $\forall k \{\mathcal{E}_{kj}\}_{j=1}^J = \mathcal{R}_2$

- Facteurs initiaux : $\forall k \{\mathbf{T}_k^0\} = \{\mathbf{I}\}_{k=0}^k$

- Valeurs initiales : $\forall k \{\lambda_0^k\}_{k=0}^K = \{1\}_{k=0}^K$

- **Paramètres de l'algorithme palm4SMSA hiérarchique h :**

- L'ensemble des contraintes sur les facteurs : $\forall k \forall j \mathcal{E}_{kj} = \mathcal{R}_2$

- L'ensemble des contraintes sur les facteurs : $\forall k \forall j \tilde{\mathcal{E}}_{kj} = \mathcal{R}_{\frac{N^2}{2J}}$

**Résultat** : Le tableau 1 montre les résultats de l'approximation de **A**, l'erreur relative est très grande et ceci indépendamment de la taille, cela montre que les 3 algorithmes ont des difficultés à apprendre la décomposition (18) de **A**. En effet, ces résultats sont attendus : théoriquement, l'algorithme **palm4SMSA** doit atteindre le minimum absolu ou au moins trouver une bonne approximation de **A**, néanmoins le choix d'une bonne initialisation reste un sujet de recherche. Les algorithmes **palm4SMSA hiérarchique p** et **h** ont une structure hiérarchique : on trouve le 1er terme, et après le 2ème et ainsi de suite. Ces algorithmes voient alors **A** comme une seule entité et n'ont aucune information sur sa forme **F** + **H**.

En conclusion les algorithmes proposés dans ce rapport ont des difficultés à apprendre des transformées rapides à cause des minimas locaux concernant **palm4SMSA**, et par conception concernant **palm4SMSA hiérarchique p et h**.

<table><caption>TABLE 1 - Erreur relative des approximations de la matrice $F+H$ de taille $N$ en utilisant les algorithmes **palm4SMSA**, **palm4SMSA hiérarchique p et h**</caption><thead><tr><th>Taille N/Algorithme</th><th>palm4SMSA</th><th>palm4SMSA_h_p</th><th>palm4SMSA_h_h</th></tr></thead><tbody><tr><td>16</td><td>0.91</td><td>0.92</td><td>0.97</td></tr><tr><td>32</td><td>0.96</td><td>0.96</td><td>0.94</td></tr><tr><td>64</td><td>0.98</td><td>0.99</td><td>0.99</td></tr><tr><td>128</td><td>0.99</td><td>0.99</td><td>0.99</td></tr></tbody></table>

## 6.3 Sur l'approximation des opérateurs de grande taille

On propose dans cette expérience une méthode d'approximation d'un opérateur **A** quelconque de grande taille ($M, N > 150$). Cette méthode garantie systématiquement 30% à 10% de réduction de temps de calcul lors d'une multiplication matrice/vecteur entraînant seulement une faible erreur relative. Soit **A** un opérateur, on effectue l'approximation de **A** suivant le $\Sigma$-$FA\mu$ST suivant :

$$ \mathbf{A} \approx \sum_{k=1}^{6} \prod_{j=1}^{2} \mathbf{s}_{kj}. \qquad (19) $$

On utilise l'algorithme **palm4SMSA hiérarchique p** pour l'apprentissage des matrices parcimonieuses. Le choix de (19) est lié à des motivations expérimentales : vue la nature hiérarchique de
**palm4SMSA hiérarchique p**, il est naturel d'obtenir des meilleures approximations en augmentant le nombre de termes *K*. Cependant, en augmentant *K*, on diminue le nombre des non-zéros sur chaque chaque terme, l'approximation $\mathbf{A} \approx \mathbf{S}_1\mathbf{S}_2$ devient alors très pauvre. Il faut alors trouver le meilleur compromis termes/facteurs, et donc le choix de (19) ($K = 6/J = 2$) vient dans ce sens, il est probable de retrouver un compromis meilleur en poussant les expériences.

## Paramètres de l'expérience

- Les approximations sont faites sur les matrices aléatoires suivantes :

- **A**: tous les éléments suivent une loi normale $N(0, 1)$

- **B**: tous les éléments suivent une loi normale $N(10, 25)$

- **C**: tous les éléments suivent une loi uniforme $U([0, 1])$

- **D**: tous les éléments suivent une loi uniforme $U([0, 100])$

- $M = N$ et $N \in \{256, 512, 1024\}$.

- $\mathbf{CR} = \{0.7, 0.8, 0.9\}$

- **paramètres de l'algorithme palm4SMSA hiérarchique p** :

- L'ensemble des contraintes : $\forall k \forall j \quad \mathcal{E}_{kj} = \mathcal{E}_s$ avec $s = \frac{\mathbf{CR} \cdot \|\mathbf{A}\|_0}{12}$. En effet, on distribue le budget total du $\mathbf{\Sigma}-\mathbf{FA}_\mu\mathbf{ST}$ équitablement sur tous les termes et puis équitablement sur les facteurs de chaque terme.

- Initialisation par défaut

**Résultat** : Les tableaux 2,3 et 4 représentent les résultats, il semble que l'erreur relative lors d'une approximation suivant le protocole ci-dessus dépend **uniquement** de la complexité relative. Plusieurs expériences en était faites sur différentes matrices de différente taille ($M, N > 150$) pour renforcer cette conclusion. L'approximation des matrices aléatoires nous laisse penser que ce protocole peut être utilisée sur un grand nombre d'opérateurs pour réduire le coût du stockage et de calcul.

<table><caption>TABLE 2 – $N = M = 256$, l'erreur relative des approximations des matrices aléatoires **A**, **B**, **C** et **D** suivant le $\mathbf{\Sigma}-\mathbf{FA}_\mu\mathbf{ST}$ (19) en utilisant l'algorithme **palm4SMSA hiérarchique p**</caption><thead><tr><th>Matrice/CR</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><td>A</td><td>0.07</td><td>0.04</td><td>0.03</td></tr><tr><td>B</td><td>0.06</td><td>0.04</td><td>0.03</td></tr><tr><td>C</td><td>0.07</td><td>0.05</td><td>0.03</td></tr><tr><td>D</td><td>0.07</td><td>0.05</td><td>0.03</td></tr></tbody></table>

<table><caption>TABLE 3 – $N = M = 512$, l'erreur relative des approximations des matrices aléatoires **A**, **B**, **C** et **D** suivant le $\mathbf{\Sigma}-\mathbf{FA}_\mu\mathbf{ST}$ (19) en utilisant l'algorithme **palm4SMSA hiérarchique p**</caption><thead><tr><th>Matrice/CR</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><td>A</td><td>0.06</td><td>0.04</td><td>0.02</td></tr><tr><td>B</td><td>0.06</td><td>0.04</td><td>0.03</td></tr><tr><td>C</td><td>0.07</td><td>0.05</td><td>0.03</td></tr><tr><td>D</td><td>0.07</td><td>0.05</td><td>0.03</td></tr></tbody></table>
<table><caption>TABLE 4 – N = M = 1024, l'erreur relative des approximations des matrices aléatoires A, B, C et D suivant le Σ-FAμST (19) en utilisant l'algorithme palm4SMSA hiérarchique p</caption><thead><tr><th>Matrice/CR</th><th>0.7</th><th>0.8</th><th>0.9</th></tr></thead><tbody><tr><th>A</th><td>0.06</td><td>0.04</td><td>0.02</td></tr><tr><th>B</th><td>0.06</td><td>0.04</td><td>0.03</td></tr><tr><th>C</th><td>0.07</td><td>0.05</td><td>0.04</td></tr><tr><th>D</th><td>0.07</td><td>0.05</td><td>0.04</td></tr></tbody></table>

## 6.4 Sur la Réduction de coût du stockage d'une base de donnée

Dans cette partie, on utilise le modèle Σ-FAμST pour réduire le coût de stockage de la base de donnée nommée "Absentéisme au travail"[10]. Cette base de donnée est constituée par 740 individus et 21 variables, les variables sont réelles ou entières. On est alors dans le cadre d'une matrice A de taille 740 × 21. On propose l'approximation de A par le Σ-FAμST suivant :

$$A \approx \sum_{k=1}^{K} \prod_{j=1}^{2} S_{kj} \qquad (20)$$

L'approximation est faite par l'algorithme **palm4SMSA hiérarchique p**

### Paramètres de l'expérience

- Nombre de termes : $k \in \{1, ..., 6\}$

- **paramètres de l'algorithme palm4SMSA hiérarchique p :**

- L'ensemble des contraintes : $\forall k \quad \mathcal{E}_{k1} = \mathcal{R}_{2}^{21 \times 21}$, $\mathcal{E}_{k2} = \mathcal{R}_{2}^{740 \times 21}$

- Initialisation par défaut.

**Résultat :** Le Tableau 5 représente les résultats, l'algorithme **palm4SMSA hiérarchique p** arrive à accomplir des bonnes approximations. On peut par exemple économiser 50% du coût de stockage de cette base de donnée en engendrant une faible erreur relative de.

<table><caption>TABLE 5 – Approximation de A suivant le Σ-FAμST (20) en utilisant l'algorithme **palm4SMSA hiérarchique p**</caption><thead><tr><th>Nombre de termes K</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><th>Érreur relative</th><td>0.15</td><td>0.05</td><td>0.03</td><td>0.02</td><td>0.01</td><td>0.007</td></tr><tr><th>Complexité relative</th><td>0.12</td><td>0.24</td><td>0.36</td><td>0.48</td><td>0.59</td><td>0.71</td></tr></tbody></table>

## 7 Discussion

On a proposé dans ce travail le modèle Σ-FAμST qui contribue principalement à l'approximation des matrices par une somme de produits de matrices parcimonieuses. On a motivé le choix de proposer une extension du modèle FAμST dans la section 4 notamment en présentant la proposition 1 et l'expérience 4.2.

Les algorithmes d'apprentissage proposés généralisent les algorithmes proposés dans [3] et [2]. L'algorithme **palm4SMSA** 4 s'avère très décevant par la qualité des minimas locaux obtenus.
On a proposé dans 5.4 une discussion sur l'initialisation et son impact sur les résultats obtenus. La hiérarchisation du problème, nous a permis d'introduire les algorithmes **palm4SMSA hiérarchique p et h** 6 et 7, ces algorithmes montrent leurs efficacités sur plusieurs opérateurs et offrent un bon compromis complexité/précision (Expériences 6.1, 6.3 et 6.4).

On finit ce rapport par une suggestion des prochaines étapes de recherche pour une meilleure exploitation du modèle Σ-FAμST :

1. Les théorèmes 1, 2 et 3 présentés dans la section 4 et les expériences effectuées dans la section 6 montrent que le modèle Σ-FAμST est très expressif, la proposition d'autres expériences, et l'application sur des données réelles, sera un grand pas vers une meilleure exploitation du modèle. Le chercheur peut commencer par raffiner les expériences présentées dans ce rapport notamment la recherche d'un protocole expérimental plus fin en terme d'approximation que (19) et l'appliquer sur un grand nombre de matrices pour consolider les conclusions des expériences 6.3 et 6.4.

2. L'obtention de minima locaux satisfaisants passe par un bon réglage des paramètres des algorithmes proposés (4,6,7), notamment, l'initialisation. Un premier pas vers une meilleure initialisation passera peut-être par comprendre la raison pour laquelle l'initialisation par défaut marche systématiquement mieux que les autres initialisations aléatoires.

3. La proposition de nouveaux algorithmes en se basant sur **PALM** avec un ordre de mise à jour différent peut être bénéfique lors des approximations. On a essayé l'ordre cité dans ce rapport, et l'ordre inverse (de gauche à droite) sans un réel bénéfice, d'autres ordres de mise à jour restent à explorer.

4. Une meilleure exploitation du modèle Σ-FAμST peut passer en dernier recours par d'autres approches d'optimisation que **PALM**. F. Malgouyres propose pour résoudre le problème (10) d'utiliser une variante de **BinaryConnect**[11].

Références

[1] Y. Li, H. Yang, E. Martin, K. L. Ho, and L. Ying. Butterfly factorization. *Multiscale Model. Simul.*, 13 :714–732, 2015.

[2] L. Le Magoarou and R. Gribonval. Flexible multilayer sparse approximations of matrices and applications. *IEEE J. Sel. Topics Signal Process.*, 10(4):688–700, 2016.

[3] J. Bolte, S. Sabach, and M. Teboulle. Proximal alternating linearized minimization for non-convex and nonsmooth problems. *Mathematical Programming*, 146 :459–494, 2014.

[4] F. Malgouyres, O. Chabiron, J. Y. Tourneret, and N. Dobigeon. Toward fast transform learning. *International Journal of Computer Vision*, vol. 114(n° 2) :pp. 195–216, September 2015.

[5] L. Siwei and W. Xin. On algorithms for sparse multi-factor nmf. In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, *Advances in Neural Information Processing Systems*, volume 26. Curran Associates, Inc., 2013.

[6] J. C. Ngouonou. Apprentissage de transformées rapides par décomposition en produit de facteurs 2-uniformes. Master's thesis, Aix-Marseille Univ, 2020.

[7] Q. Le and R. Gribonval. Structured Support Exploration For Multilayer Sparse Matrix Factorization. In *ICASSP 2021 - IEEE International Conference on Acoustics, Speech and Signal*
*Processing*, pages 1–5, Toronto, Ontario, Canada, June 2021. IEEE. Published in ICASSP 2021.

[8] Y. Ito, S. Hattori, and H. Maeda. On the decomposition of a matrix into the sum of stable matrices. *Linear Algebra and its Applications*, 297(1):177–182, 1999.

[9] D. I. Merino. The sum of orthogonal matrices. *Linear Algebra and its Applications*, 436(7):1960–1968, 2012.

[10] A. Martiniano, R. P. Ferreira, R. J. Sassi, and C. Affonso. Application of a neuro fuzzy network in prediction of absenteeism at work. In *7th Iberian Conference on Information Systems and Technologies (CISTI 2012)*, pages 1–4, 2012.

[11] M. Courbariaux, Y. Bengio, and J. P. David. Binaryconnect: Training deep neural networks with binary weights during propagations, 2016.
# Appendices

## A Preuve de la proposition 1

Montrons que :

$$A = BC \iff \forall (i,j) \in [1,4]^2 \exists! k \in [1,4] a_{ij} = b_{ik}c_{kj}. \quad (21)$$

L'existence de $k$ est triviale, montrons l'unicité :

Supposons que $k$ n'est pas unique, il existe alors $i, j$ et $k'$ tel que $a_{ij} = b_{ik}c_{kj} + b_{ik'}c_{k'j}$ (car **B** et **C** ont exactement deux éléments non nuls par ligne et par colonne), et donc, pour tout $j' \neq j$ on a $a_{ij'} = b_{ik}c_{kj'} + b_{ik'}c_{k'j'}$ ou $a_{ij'} = b_{ik}c_{kj'}$ ou $a_{ij'} = b_{ik'}c_{k'j'}$ (car les éléments de **A** sont tous non nuls).

4 cas apparaissent :

- Il existe $j' \neq j$ tels que $a_{ij'} = b_{ik}c_{kj'} + b_{ik'}c_{k'j'}$ et donc, il existe $j'' \neq j, j'$ tels que $a_{ij''} = 0$.

- Il existe $j', j'' \neq j$, tels que $a_{ij'} = b_{ik}c_{kj'}$ et $a_{ij''} = b_{ik}c_{kj''}$, et donc, **C** contient 3 éléments non nuls dans la ligne $k$.

- Il existe $j', j'' \neq j$ tels que $a_{ij'} = b_{ik'}c_{k'j'}$ et $a_{ij''} = b_{ik'}c_{k'j''}$, et donc, **C** contient 3 éléments non nuls dans la ligne $k'$.

- Il existe $j', j'' \neq j$ tel que $a_{ij'} = b_{ik}c_{kj'}$ et $a_{ij''} = b_{ik'}c_{k'j''}$, et donc, il existe $j''' \neq j, j', j''$ tel que $a_{ij'''} = 0$.

En conclusion : $k$ est unique. Montrons maintenant que :

$$\forall (i, j) \in [1, 4]^2 \ k \in [1, 4] \ a_{ij} = b_{ik}c_{kj} \Rightarrow \exists! i' \neq i \ a_{i'j} = b_{i'k}c_{kj}. \quad (22)$$

En effet pour tout $i' \neq i$, $a_{i'j} = b_{i'k}c_{kj}$ ou $a_{i'j} = b_{i'k'}c_{k'j}$ avec $k' \neq k$ (sinon la colonne $j$ de la matrice **B** contiendra plus que 2 éléments). Si pour tout $i' \neq i$ $a_{i'j} = b_{i'k}c_{kj}$, la matrice **B** contient 4 éléments non nuls dans la colonne $j$, et donc, il existe $i' \neq i$ tel que $a_{i'j} = b_{i'k}c_{kj}$.

On déduit de (22) que pour tout $i, j$, il existe $i' \neq i$ tel que $\frac{a_{ij}}{a_{i'j}} = \frac{b_{ik}}{b_{i'k}}$. De même, il existe $j' \neq j$ tel que $\frac{a_{ij'}}{a_{i'j'}} = \frac{b_{ik}}{b_{i'k}}$. Donc $\frac{a_{ij}}{a_{i'j}} = \frac{a_{ij'}}{a_{i'j'}}$.

## B Sur la sensibilité de palm4MSA hiérarchique

L'algorithme **palm4MSA hiérarchique** permet d'obtenir une factorisation exacte de la matrice de Hadamard **H**, et ceci indépendamment de sa taille $N \times N$ avec les paramètres de réglage suivant :

Paramètres de l'algorithme **palm4MSA hiérarchique** :

- L'ensemble des contraintes sur les facteurs : $\forall j \in \mathcal{E}_j = \mathcal{R}_2$

- L'ensemble des contraintes sur les facteurs : $\forall j \in \tilde{\mathcal{E}}_j = \mathcal{R}_{\frac{N^2}{2j}}$

En gardant ces mêmes paramètres de réglage, on explore dans cette expérience le comportement de cet algorithme quand on rajoute un petit bruit à la matrice de Hadamard **H**. Soit la matrice **A** suivante :

$$A = H + \alpha.M$$

avec **M** une matrice aléatoire, et $\alpha$ un scalaire qui contrôle le bruit.
Paramètres de l'expérience

— $N \in \{16, 32\}$

— $\mathbf{M}$ : tous les éléments suivent une loi normale $\mathcal{N}(0, 1)$

— $\alpha \in \{10^{-5}, 10^{-6}, \dots, 10^{-9}, 10^{-10}\}$

<table><caption>TABLE 6 – N = 16. Approximation de $\mathbf{A}$ en utilisant l'algorithme **palm4MSA hiérarchique**</caption><tr><td>α</td><td>10<sup>-5</sup></td><td>10<sup>-6</sup></td><td>10<sup>-7</sup></td><td>10<sup>-8</sup></td><td>10<sup>-9</sup></td><td>10<sup>-10</sup></td></tr><tr><td>Érreur relative</td><td>0.72</td><td>0.72</td><td>0.71</td><td>0.61</td><td>0</td><td>0</td></tr></table>

<table><caption>TABLE 7 – N = 32. Approximation de $\mathbf{A}$ en utilisant l'algorithme **palm4MSA hiérarchique**</caption><tr><td>α</td><td>10<sup>-5</sup></td><td>10<sup>-6</sup></td><td>10<sup>-7</sup></td><td>10<sup>-8</sup></td><td>10<sup>-9</sup></td><td>10<sup>-10</sup></td></tr><tr><td>Érreur relative</td><td>0.82</td><td>0.81</td><td>0.82</td><td>0.68</td><td>0</td><td>0</td></tr></table>

**Résultat** : Les tableaux 6 et 7 montrent les résultats. L'algorithme est très sensible au bruit : quand $\alpha > 10^{-9}$, il est incapable de trouver une bonne approximation, l'erreur relative est très grande.