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

HAL
open science

# Factorisation butterfly par identification algorithmique de blocs de rang un

Léon Zheng, Gilles Puy, Elisa Riccietti, Patrick Pérez, Rémi Gribonval

► To cite this version:

Léon Zheng, Gilles Puy, Elisa Riccietti, Patrick Pérez, Rémi Gribonval. Factorisation butterfly par identification algorithmique de blocs de rang un. XXIXème Colloque Francophone de Traitement du Signal et des Images, Aug 2023, Grenoble, France. hal-04145743

HAL Id: hal-04145743

https://inria.hal.science/hal-04145743

Submitted on 29 Jun 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 Attribution 4.0 International License
# Factorisation butterfly par identification algorithmique de blocs de rang un

Léon ZHENG¹,² Gilles PUY² Elisa RICCIETTI¹ Patrick PÉREZ² Rémi GRIBONVAL¹

¹Univ Lyon, EnsL, UCBL, CNRS, Inria, LIP, F-69342, LYON Cedex 07, France.

²valeo.ai, Paris, France.

**Résumé** – Plusieurs matrices associées à des transformées rapides possèdent une certaine propriété de rang faible qui se caractérise par l'existence de plusieurs partitions par blocs de la matrice, où chaque bloc est de rang faible. À condition de connaître ces partitions, il existe alors des algorithmes, dits de factorisation butterfly, qui approchent la matrice en un produit de facteurs creux, permettant ainsi une évaluation rapide de l'opérateur linéaire associé. Cet article propose une nouvelle méthode pour identifier algorithmiquement les partitions en blocs de rang faible d'une matrice admettant une factorisation butterfly, sans hypothèse analytique sur ses coefficients.

**Abstract** – Many matrices associated with fast transforms possess a certain low-rank property characterized by the existence of several block partitionings of the matrix, where each block is of low rank. Provided that these partitionings are known, there exist algorithms, called butterfly factorization algorithms, that approximate the matrix into a product of sparse factors, thus enabling a rapid evaluation of the associated linear operator. This paper proposes a new method to identify algebraically these block partitionings for a matrix admitting a butterfly factorization, without any analytical assumption on its entries.

## 1 Introduction

L'évaluation rapide d'un opérateur linéaire est un enjeu clé dans de nombreux domaines comme le calcul scientifique, le traitement du signal ou l'apprentissage automatique. Dans des applications mettant en jeu un très grand nombre de paramètres, le calcul direct de la multiplication matrice-vecteur passe difficilement à l'échelle pour cause de complexité quadratique en la taille de la matrice. De nombreux travaux se sont ainsi intéressés à la construction d'algorithmes rapides pour la multiplication matricielle, en s'appuyant typiquement sur des propriétés analytiques ou algébriques des matrices apparaissant dans les problèmes étudiés.

Ces algorithmes rapides sont souvent associés à une factorisation creuse de la matrice correspondante, comme c'est le cas pour la matrice de Hadamard ou de la transformée de Fourier discrète. En effet, à permutation près des lignes et des colonnes, ces matrices de taille $N$ possèdent une factorisation butterfly, dans le sens où elle s'écrivent comme le produit de $O(\log N)$ facteurs ayant chacun $O(N)$ coefficients non nuls, et dont les supports ont une structure particulière illustrée dans la figure 1. Il a été montré dans [4] que la classe des matrices admettant une telle factorisation est expressive, au sens où elle contient plusieurs matrices structurées utilisées en apprentissage ou en traitement du signal. Ce modèle serait donc pertinent pour chercher des factorisations creuses d'opérateurs pour lesquels un algorithme d'évaluation rapide n'est pas connu.

<figure><img src="image_3.png" /><figcaption>FIGURE 1 : Supports des facteurs butterfly, $N = 16$ (cf. section 2).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2 : Histogramme des erreurs relatives pour la résolution de (1) dans le cas sans bruit et $N = 8$, en énumérant et en fixant chaque permutation de lignes et de colonnes, comme expliqué en section 3. Parmi toutes les paires de permutation de lignes et colonnes, seule une paire donne une erreur nulle, comme indiqué par le cercle rouge.</figcaption></figure>

Trouver un algorithme rapide associé à la factorisation butterfly se formalise alors comme un problème d'optimisation, où l'on minimise en norme de Frobenius $\| \cdot \|_F$ l'erreur d'approximation d'une matrice $A \in \mathbb{C}^{N \times N}$ par un produit de facteurs butterfly $X^{(1)}, \dots, X^{(L)}$, à permutations près des lignes et colonnes, encodées par les matrices de permutation $P, Q$ :

$$ \min_{(X^{(\ell)})_{\ell=1}^L, P, Q} \|A - Q^\top X^{(1)} \dots X^{(L)} P\|_F \quad (1) $$

Lorsque les permutations optimales $P, Q$ sont connues, il existe un algorithme hiérarchique efficace en complexité $O(N^2)$ pour trouver des facteurs butterfly $(X^{(\ell)})_{\ell=1}^L$ donnant une faible erreur d'approximation, avec garanties de reconstruction dans le cas du problème sans bruit [6, 13]. Mais lorsque ces permutations ne sont pas connues, le problème (1) est conjecturé comme étant difficile : d'une part, si l'on énumère toutes les permutations possibles pour résoudre (1), on constate numériquement que seule une petite proportion de permutations donne une faible erreur d'approximation, comme illustré en figure 2, ce qui montre la nécessité d'identifier les bonnes permutations afin de résoudre (1); mais d'autre part, une recherche exhaustive de toutes les permutations n'est pas tractable, même en tenant compte de certaines équivalences de permutations vis-à-vis de (1), comme discuté en section 2.
Afin d'identifier les permutations optimales, nous nous ap-
puyons sur le fait qu'une matrice admettant une factorisation
butterfly possède une certaine propriété dite de rang faible
complémentaire [7], dans le sens où il existe des partitions par
blocs de la matrice, où chaque bloc est de rang faible. Il suffit
alors d'identifier ces partitions pour résoudre le problème (1).

Ceci peut se faire *analytiquement* lorsque les coefficients de la matrice s'expriment par un noyau régulier $(\mathbf{x}, \omega) \mapsto K(\mathbf{x}, \omega)$ évalué sur des paramètres $\{\mathbf{x}_i\}_{i=1}^N$, $\{\omega_j\}_{j=1}^N$, par exemple pour des matrices associées à certains opérateurs intégraux [2] ou transformées spéciales de fonctions [8].

En revanche, si la matrice étudiée n'a pas de forme ana-
lytique, ou si celle-ci n'est pas accessible, la littérature ne
propose pas, à notre connaissance, de méthode pour identifier
ces partitions. Aussi, nous proposons une heuristique à base
de partitionnement spectral alterné des lignes et des colonnes
pour identifier les partitions en blocs de rang faible, sans hypo-
thèse analytique. En termes d'applications, cette heuristique
permet de vérifier algorithmiquement qu'un opérateur linéaire,
pour lequel on ignore l'existence d'un algorithme d'évaluation
rapide, possède une propriété de rang faible complémentaire
approchée donnant lieu à une bonne approximation de la ma-
trice associée par un produit de facteurs butterfly.

La section 2 détaille cette propriété de rang faible utilisée
pour identifier les permutations optimales. La section 3 montre
numériquement la nécessité d'identifier ces permutations, ce
qui motive notre méthode expliquée en section 4, et validée
empiriquement en section 5.

2 Formulation du problème

Dans le reste du papier, nous nous plaçons dans le cas des
matrices carrées de taille $N := 2^L$ pour un certain en-
tier $L \ge 2$. Pour une matrice quelconque $\mathbf{M}$, son support
$\text{supp}(\mathbf{M})$ est l'ensemble des indices correspondants aux co-
efficients non nuls de $\mathbf{M}$. Notons $\mathbf{I}_N$ la matrice identité de
taille $N$, et $\otimes$ le produit de Kronecker. Selon [4, 6, 13],
$\mathbf{A} \in \mathbb{C}^{N \times N}$ est une *matrice butterfly* si elle admet une factori-
sation $\mathbf{A} = \mathbf{X}^{(1)} \dots \mathbf{X}^{(L)}$, où chaque facteur $\mathbf{X}^{(\ell)} \in \mathbb{C}^{N \times N}$,
nommé *facteur butterfly*, satisfait la contrainte de support fixe
$\text{supp}(\mathbf{X}^{(\ell)}) \subseteq \text{supp}(\mathbf{S}_{\text{bf}}^{(\ell)})$ pour $\ell \in [L] := \{1, \dots, L\}$, avec
$\mathbf{S}_{\text{bf}}^{(\ell)} := \mathbf{I}_{2^{\ell}-1} \otimes [\frac{1}{1} \frac{1}{1}] \otimes \mathbf{I}_{N/2^{\ell}}$. L'ensemble des matrices admet-
tant une telle factorisation est noté $\mathcal{B}$. Les supports $\mathbf{S}_{\text{bf}}^{(\ell)}$ sont
illustrés par la figure 1 et, par abus de notation, la contrainte
de support fixe sera notée $\text{supp}(\mathbf{X}^{(\ell)}) \subseteq \mathbf{S}_{\text{bf}}^{(\ell)}$ par la suite.

Le problème que nous souhaitons résoudre est donc (1),
sous la contrainte $\text{supp}(\mathbf{X}^{(\ell)}) \subseteq \mathbf{S}_{\text{bf}}^{(\ell)}$ pour tout $\ell \in [L]$.
Pour trouver une bonne solution de l'instance du problème
où les matrices de permutations $\mathbf{P}$ et $\mathbf{Q}$ sont fixées, on ap-
plique l'algorithmhe hiérarchique de [6, 13]. Par la suite, on
note $E_{\text{bf}}(\mathbf{A}, \mathbf{P}, \mathbf{Q}) := \|\mathbf{A} - \mathbf{Q}^{\top}\tilde{\mathbf{X}}^{(1)} \dots \tilde{\mathbf{X}}^{(L)}\mathbf{P}\|_{F}$ l'erreur
d'approximation donnée par la sortie $(\tilde{\mathbf{X}}^{(1)}, \dots, \tilde{\mathbf{X}}^{(L)})$ de l'al-
gorithme 1, quand $\mathbf{P}$ et $\mathbf{Q}$ sont fixées. En revanche, le cas
difficile est celui où les permutations ne sont pas fixées : le
reste de la section explique notre approche qui s'appuie sur la
propriété de rang faible complémentaire.

**Propriété de rang faible complémentaire** [7] Celle-ci se définit à l'aide de deux « cluster tree » $T^X$ et $T^\Omega$, qui sont des arbres binaires avec $L = \log_2(N)$ niveaux (sans compter la

**Algorithmhe 1 :** Algorithme hiérarchique pour ré-soudre l'instance de (1) où **P** et **Q** sont fixées.

Entrées : $\mathbf{A} \in \mathbb{C}^{N \times N}$, permutations $\mathbf{P}, \mathbf{Q}$.

1 Résoudre $\min \|\mathbf{Q}\mathbf{A}\mathbf{P}^\top - \mathbf{X}^{(1)} \dots \mathbf{X}^{(L)}\|_F$ sous la contrainte $\text{supp}(\mathbf{X}^{(\ell)}) \subseteq \mathbf{S}_{\text{bf}}^{(\ell)}$ pour tout $\ell \in [L]$ via [13, Algorithme 3.2]

2 retourner les facteurs butterfly $(\tilde{\mathbf{X}}^{(1)}, \dots, \tilde{\mathbf{X}}^{(L)})$

racine) dont chaque nœud est un sous-ensemble non vide de
$[[N]]$, avec comme racine $[[N]]$ au niveau 0, et où les enfants
constituent une partition de leur parent en deux sous-ensembles
de même cardinal [5]. Une matrice **A** possède la propriété de
rang faible complémentaire pour $T^X$ et $T^\Omega$ si, pour chaque $\ell \in$
$[[L-1]]$, pour chaque nœud $R$ au niveau $L-\ell$ de $T^X$ et pour
chaque nœud $C$ au niveau de $\ell$ de $T^\Omega$, la restriction $\mathbf{A}_{R,C}$ de **A**
sur les lignes et colonnes indexées par $R$ et $C$ est de rang faible.
Sachant que les nœuds d'un même niveau de $T^X$ et ceux de
$T^\Omega$ forment respectivement une partition des indices de lignes
et de colonnes, la propriété de rang faible complémentaire
impose que les blocs $\{\mathbf{A}_{R,C}\}_{R,C}$ des partitions décrites par
les niveaux de $T^X$ et $T^\Omega$ soient de rang faible.

Il a été montré dans [13] que toute matrice butterfly $\mathbf{A} \in \mathcal{B}$
satisfait la propriété de rang faible complémentaire (de rang 1)
pour une paire d'arbres $(T_{\text{bf}}^X, T_{\text{bf}}^{\Omega})$ spécifique, comme expliqué
ci-dessous pour rappeler les concepts qui nous seront utiles
par la suite. Notons $\mathbf{S}_{\text{bf}}^{(p:q)} := \mathbf{S}_{\text{bf}}^{(p)} \dots \mathbf{S}_{\text{bf}}^{(q)}$ ($1 \le p < q \le L$),
et définissons la $\ell$-ème classe Monarch [3] pour $\ell \in [[L-1]]$
comme étant l'ensemble

$$
\mathcal{M}^{(\ell)} := \{ \mathbf{X}\mathbf{Y}, \text{supp}(\mathbf{X}) \subseteq \mathbf{S}_{\text{bf}}^{(1:\ell)}, \text{supp}(\mathbf{Y}) \subseteq \mathbf{S}_{\text{bf}}^{(\ell+1:L)} \}.
$$

On vérifie alors que $\mathcal{B} \subseteq \bigcap_{\ell=1}^{L-1} \mathcal{M}^{(\ell)}$, et qu'une matrice $\mathbf{A}$ appartient à $\mathcal{M}^{(\ell)}$ si, et seulement si, $\mathbf{A}_{R,C}$ est de rang au plus 1 pour tout $(R, C) \in \mathcal{P}^{(\ell)} := \{R_i^{(\ell)}\}_{i=1}^{N/2^{\ell}} \times \{C_j^{(\ell)}\}_{j=1}^{2^{\ell}}$, où

$$
R_i^{(\ell)} := \{i + (k - 1)N/2^{\ell}, k \in [2^\ell]\},
$$

$$
C_j^{(\ell)} := \{(j-1)N/2^{\ell} + k, k \in [N/2^{\ell}]\}.
$$

On observe que $\{R_i^{(\ell)}\}_{i=1}^{N/2^\ell}$ et $\{C_j^{(\ell)}\}_{j=1}^{2^\ell}$ constituent chacun
une partition de $[[N]]$ : ainsi, lorsque $\mathbf{A} \in \mathcal{M}^{(\ell)}$, ces partitions
de lignes et de colonnes décrivent une partition de $\mathbf{A}$ en blocs
de rang 1. On définit alors $T_{\text{bf}}^X$ et $T_{\text{bf}}^\Omega$ comme étant les deux
arbres pour lesquels les nœuds de $T_{\text{bf}}^X$ au niveau $L-\ell$ et
les nœuds de $T_{\text{bf}}^\Omega$ au niveau $\ell$ sont précisément $\{R_i^{(\ell)}\}_{i=1}^{N/2^\ell}$ et
$\{C_j^{(\ell)}\}_{j=1}^{2^\ell}$, pour chaque $\ell \in [L-1]$. Ainsi, $\mathbf{A} \in \bigcap_{\ell=1}^{L-1} \mathcal{M}^{(\ell)}$
est précisément une reformulation de la propriété de rang faible
complémentaire (de rang 1) pour les arbres $T_{\text{bf}}^X$ et $T_{\text{bf}}^\Omega$.

**Permutation d'indices** Pour tout « cluster tree » $T$ dont la racine est $[[N]]$, et pour toute permutation $\sigma : [[N]] \to [[N]]$, on définit le « cluster tree » $\sigma(T)$ obtenu en permutant selon $\sigma$ les indices dans $T$. En particulier, on remarque que pour tout arbre $T^X$ et $T^\Omega$, il existe plusieurs matrices de permutation $\mathbf{P}$ et $\mathbf{Q}$ pour lesquelles $\sigma_Q(T_{\text{bf}}^X) = T^X$ et $\sigma_P(T_{\text{bf}}^\Omega) = T^\Omega$, où $\sigma_P, \sigma_Q$ sont les permutations associées aux matrices $\mathbf{P}, \mathbf{Q}$. Ceci définit alors des classes d'équivalence de permutations de lignes $[\mathbf{Q}_{T \times}]$ et de colonnes $[\mathbf{P}_{T^\Omega}]$.

Approche pour résoudre (1) Supposons que la matrice
cible dans (1) soit de la forme $\mathbf{A} := \tilde{\mathbf{Q}}^{\top}\tilde{\mathbf{A}}\tilde{\mathbf{P}}$, où $\tilde{\mathbf{A}} \in \mathcal{B}$
est une matrice butterfly, et $\tilde{P}, \tilde{Q}$ sont deux matrices de permutations arbitraires inconnues. Puisque $\tilde{A}$ satisfait la propriété de rang faible complémentaire pour les arbres $T_{\text{bf}}^X$ et $T_{\text{bf}}^{\Omega}$, la matrice $\tilde{Q}^{\top}\tilde{A}\tilde{P}$ la satisfait également mais pour les arbres $T^X := \sigma_{\tilde{Q}}(T_{\text{bf}}^X)$ et $T^{\Omega} := \sigma_{\tilde{P}}(T_{\text{bf}}^{\Omega})$. Si l'on parvient à reconstruire $T^X$ et $T^{\Omega}$ à partir de l'observation de $\mathbf{A}$, alors (1) peut se résoudre en choisissant une paire quelconque $(\mathbf{P}, \mathbf{Q}) \in [\mathbf{P}_{T^{\Omega}}] \times [\mathbf{Q}_{T^X}]$, et en appliquant l'algorithme 1 avec ces permutations fixées. En effet, un tel choix suffit pour garantir que $\mathbf{QAP}^{\top}$ satisfasse la propriété de rang faible complémentaire pour $T_{\text{bf}}^X$ et $T_{\text{bf}}^{\Omega}$, i.e., $\mathbf{QAP}^{\top} \in \bigcap_{\ell=1}^{L-1} \mathcal{M}^{(\ell)}$, puisque les matrices $\mathbf{P}^{\top}$ et $\mathbf{Q}^{\top}$ sont associées aux permutations inverses $\sigma_{\tilde{P}}^{-1}$ et $\sigma_{\tilde{Q}}^{-1}$, et par définition des classes d'équivalence, on a bien $\sigma_{\tilde{P}}^{-1}(T^{\Omega}) = T_{\text{bf}}^{\Omega}$ et $\sigma_{\tilde{Q}}^{-1}(T^X) = T_{\text{bf}}^X$. En conclusion, afin de résoudre (1), il est suffisant d'identifier les arbres $T^X$ et $T^{\Omega}$ pour lesquelles la matrice cible satisfait la propriété de rang faible complémentaire, ce qui se ramène, par définition, à identifier des partitions de $\mathbf{A}$ en blocs de rang 1.

### 3 Nécessité de retrouver les partitions

Nous montrons à présent empiriquement que l'identification des arbres $T^X$ et $T^\Omega$ est en fait nécessaire pour résoudre (1) avec $\mathbf{A} := \tilde{Q}^\top \tilde{A} \tilde{P}$. Considérons $\tilde{A} \in \mathcal{B}$ dont les facteurs ont des coefficients non nuls tirés selon une gaussienne centrée réduite. Puis, nous énumérons tous les arbres $T^X, T^\Omega$ possibles, calculons une solution via l'algorithme 1 en fixant une paire arbitraire $(\mathbf{P}, \mathbf{Q}) \in [\mathbf{P}_{T^\Omega}] \times [\mathbf{Q}_{T^X}]$, et vérifions que les seuls arbres donnant une erreur $E_{\text{bf}}(\mathbf{A}, \mathbf{P}, \mathbf{Q})$ faible sont $T^X$ et $T^\Omega$. Un dénombrement de tous les arbres montre en revanche que cette expérience n'est pas tractable pour une grande taille $N$. En effet, le nombre $u_N$ de « cluster tree » pour une même racine de cardinal $N$ satisfait la relation de récurrence $u_N = \frac{1}{2}(\frac{N}{2})(u_N)^2$ avec $u_2 = 1$, car il y a $\frac{1}{2}(\frac{N}{2})$ paires d'enfants possibles pour la racine, et chaque enfant est un « cluster tree » dont la racine est de cardinal $N/2$. Ainsi, nous considérons $N = 8$ à titre illustratif, ce qui donne $u_8 = 315$.

Dans la figure 2, la recherche exhaustive sur toutes les paires d'arbres montrent que l'erreur est nulle seulement pour une seule paire, et que les autres paires échouent à la résolution de (1). Ceci illustre donc empiriquement la nécessité d'identifier les bons arbres $T^X$ et $T^\Omega$ pour résoudre (1). Démontrer formellement une telle nécessité pour n'importe quelle taille de matrices pourra faire l'objet de futurs travaux.

### 4 Partitionnement spectral alterné

Nous proposons l'algorithme 3 à base de partitionnement spectral alterné pour identifier les arbres $T^X$ et $T^\Omega$ pour lesquels la matrice $\mathbf{A}$ satisfait la propriété de rang faible complémentaire. Pour décrire notre approche, commençons par expliquer la résolution du problème $\min_{\mathbf{M} \in \mathcal{M}^{(\ell)}, \mathbf{P}, \mathbf{Q}} \|\mathbf{A} - \mathbf{Q}^{\top}\mathbf{MP}\|_F^2$ pour chaque $\ell \in \llbracket L-1 \rrbracket$. Étant donné qu'une matrice $\mathbf{M}$ appartient à $\mathcal{M}^{(\ell)}$ si et seulement si le bloc $\mathbf{M}_{R,C}$ est de rang au plus 1 pour tout jeu de lignes et colonnes $(R, C) \in \mathcal{P}^{(\ell)}$, ce problème équivaut à

$$ \min_{\{R_i\}_{i=1}^{N/2^\ell}, \{C_j\}_{j=1}^{2^\ell}} \sum_{i=1}^{N/2^\ell} \sum_{j=1}^{2^\ell} \min_{\mathbf{x}, \mathbf{y}} \|\mathbf{A}_{R_i, C_j} - \mathbf{x}\mathbf{y}^*\|_F^2, \quad (2) $$

**Algorithme 2 :** Partitionnement spectral alterné pour résoudre (2) durant $n$ itérations, à $\ell \in \llbracket L-1 \rrbracket$ fixé.

**Entrées :** $\mathbf{A}, \ell, \alpha > 0$, graine aléatoire seed.

1. $\{C_j\}_{j=1}^{2^\ell}$ ← partition aléatoire de $\llbracket N \rrbracket$ selon seed

2. **pour** $i = 1, \dots, n$ **faire**

3. $\{R_i\}_{i=1}^{N/2^\ell}$ ← partitionnement en fixant $\{C_j\}_{j=1}^\ell$

4. $\{C_j\}_{j=1}^{2^\ell}$ ← partitionnement en fixant $\{R_i\}_{i=1}^{N/2^\ell}$

5. **fin**

6. $E_\ell \leftarrow \sum_{i,j} \min_{\mathbf{x}, \mathbf{y}} \|\mathbf{A}_{R_i, C_j} - \mathbf{x}\mathbf{y}^*\|_F^2$

**Sorties :** $E_\ell, \{R_i\}_{i=1}^{N/2^\ell}, \{C_j\}_{j=1}^{2^\ell}$.

où $\{R_i\}_{i=1}^{N/2^\ell}, \{C_j\}_{j=1}^{2^\ell}$ sont des partitions de lignes et de colonnes en sous-ensembles de même cardinal, et $\min_{\mathbf{x}, \mathbf{y}} \|\mathbf{A}_{R,C} - \mathbf{x}\mathbf{y}^*\|_F^2$ calcule la meilleure approximation de rang 1 de $\mathbf{A}_{R,C}$. Le symbole $*$ désigne la matrice adjointe.

**Partitionnement des lignes** Pour qu'une optimisation alternée fonctionne, il est nécessaire de pouvoir résoudre le problème 2 lorsqu'une des deux partitions est connue. Fixons ainsi sans perte de généralité une partition de colonnes $\{C_j\}_{j=1}^{2^\ell}$, et cherchons une partition de lignes $\{R_i\}_{i=1}^{N/2^\ell}$ qui minimise (2).

Pour cela, nous nous inspirons des méthodes existantes pour le problème de « subspace clustering » [11]. Définissons $2^\ell$ graphes $\{\mathcal{G}_j\}_{j=1}^{2^\ell}$, où les $N$ nœuds du graphe $\mathcal{G}_j$ sont les $N$ lignes de $\mathbf{A}$ restreintes aux colonnes $C_j$ –notées $\mathbf{A}_{k,C_j}$ pour $k \in \llbracket N \rrbracket$– et les poids des arêtes de $\mathcal{G}_j$ sont donnés par la matrice de similarité $\mathbf{W}^{(j)} \in \mathbb{R}^{N \times N}$ définie par

$$ \mathbf{W}_{k,l}^{(j)} := \left( \frac{\| \mathbf{A}_{k,C_j} * \mathbf{A}_{l,C_j} \|}{\| \mathbf{A}_{k,C_j} \|_2 \| \mathbf{A}_{l,C_j} \|_2} \right)^{\alpha} \quad \forall k, l \in \llbracket N \rrbracket, $$

avec $\alpha > 0$ un paramètre contrôlant le contraste entre les poids des arêtes. Intuitivement, un groupe de lignes restreintes aux colonnes $C_j$ est inter-connecté par des poids de fortes valeurs dans le graphe $\mathcal{G}_j$ lorsque les lignes correspondantes sont corrélées, ce qui est le cas lorsqu'elles forment un bloc de rang 1. Inversement, deux lignes non corrélées ont un poids de faible valeur. Ainsi, en résolvant un problème de coupe minimale sur le graphe $\mathcal{G}$ dont la matrice de similarité est $\mathbf{W} := \sum_{j=1}^{2^\ell} \mathbf{W}^{(j)}$, les groupes de nœuds obtenus doivent correspondre à un partitionnement des lignes qui minimise (2).

Concrètement, à partir de la matrice de similarité $\mathbf{W}$, un partitionnement spectral [12] du graphe $\mathcal{G}$ est effectué en calculant la décomposition en vecteurs propres du Laplacien non normalisé $\mathbf{L} := \mathbf{D} - \mathbf{W}$, où $\mathbf{D}$ est la matrice de degrés dont les coefficients diagonaux sont $\mathbf{W}(1 \dots 1)^{\top}$. Afin de garantir un partitionnement des lignes en groupe de même taille, l'étape de partitionnement $k$-moyenne sur les représentations spectrales est implémentée selon la méthode de [1].

**Optimisation alternée** Quand la partition de colonnes n'est plus fixée, la résolution de (2) suit l'algorithme 2 de partitionnement alterné, où l'on initialise aléatoirement une partition des colonnes et, à chaque itération, un partitionnement spectral des lignes est effectuée en fixant la partition de colonnes de l'itération précédente, et vice-versa en échangeant le rôle des lignes et des colonnes. Sans garanties de réussite, l'algorithme 2 peut nécessiter plusieurs réinitialisations pour trouver une solution.
**Algorithme 3 :** Identifier $T^X$, $T^\Omega$ pour résoudre (1).

Entrées : $\mathbf{A}$, $\{\alpha_k\}_{k=1}^K$, $\{\text{seed}_m\}_{m=1}^M$.

1. pour $\ell = 1, \dots, L-1$ faire

2. pour $k = 1, \dots, K$, $m = 1, \dots, M$ faire

3. $E_\ell$, $\{\tilde{R}_i^{(\ell)}\}_i$, $\{\tilde{C}_j^{(\ell)}\}_j \leftarrow$ **algorithme 2**
appliqué à $(\mathbf{A}, \ell, \alpha_k, \text{seed}_m)$

4. fin

5. Garder $\{\tilde{R}_i^{(\ell)}\}_i$, $\{\tilde{C}_j^{(\ell)}\}_j$ avec le plus petit $E_\ell$

6. fin

7. Vérifier que $\{\tilde{R}_i^{(\ell)}\}_i$, $\{\tilde{C}_j^{(\ell)}\}_j$ pour $\ell \in [L-1]$ forment chacun des « cluster tree » $T^X$, $T^\Omega$ valides

8. **si un des deux arbres n'est pas valide alors**

9. retourner échec

10. sinon

11. $(\mathbf{P}, \mathbf{Q}) \in [\mathbf{P}_{T^\Omega}] \times [\mathbf{Q}_{T^X}]$

12. retourner **succès**, $\mathbf{P}, \mathbf{Q}, E_{bf}(\mathbf{A}, \mathbf{P}, \mathbf{Q})$

13. fin

**Résolution finale de (1)** Étant donné la matrice cible $\mathbf{A} := \tilde{\mathbf{Q}}^\top \tilde{\mathbf{A}} \tilde{\mathbf{P}}$, l'algorithme 3 résout indépendamment chaque problème (2) pour $\ell \in [L-1]$ via l'algorithmhe 2. Si chaque problème (2) est bien résolu, et que les partitions trouvées forment des arbres $T^X$ et $T^\Omega$ valides dans le sens où ils satisfont les axiomes d'un « cluster tree », alors $\mathbf{A}$ possède la propriété de rang faible complémentaire pour $T^X$ et $T^\Omega$. On résout alors (1) via l'algorithmhe (1), en fixant $(\mathbf{P}, \mathbf{Q}) \in [\mathbf{P}_{T^\Omega}] \times [\mathbf{Q}_{T^X}]$.

# 5 Expériences

Nous évaluons les performances empiriques de notre méthode pour factoriser $\mathbf{A} := \tilde{\mathbf{Q}}^\top \tilde{\mathbf{A}} \tilde{\mathbf{P}} + \epsilon (\|\tilde{\mathbf{A}}\|_F / \|\mathbf{N}\|_F)\mathbf{N}$, où $\tilde{\mathbf{P}}, \tilde{\mathbf{Q}}$ sont des permutations aléatoires, $\mathbf{N}$ est une matrice avec des coefficients suivant une loi gaussienne centrée réduite, $\epsilon \ge 0$ contrôle le niveau relatif de bruit, et $\tilde{\mathbf{A}}$ correspond soit à une matrice butterfly orthogonale aléatoire définie par [9], soit à la matrice de la transformée de Fourier discrète (TFD). Nous appliquons l'algorithmhe 3 avec $\{\alpha_k\}_{k=1}^K := \{10^p\}_{p \in \{-2, -1, 0, 1, 2\}}$ et $M = 5$ sur 20 instances du problème pour $\epsilon \in \{0, 0.01, 0.03, 0.1\}$, et pour $N \in \{2^L\}_{L \in \{2, \dots, 7\}}$. La taille $N > 128$ n'est pas considérée car l'algorithmhe 3 a une complexité cubique en $N$ : une exécution prend quelques minutes pour $N = 64$, et une heure pour $N = 128$.

Lorsque $\tilde{\mathbf{A}}$ est une matrice butterfly orthogonale aléatoire, l'algorithmhe 3 a 100 % de succès sur les 20 instances du problème, pour tous les niveaux de bruit et toutes les tailles considérés, ce qui veut dire que l'algorithmhe 2 répété avec suffisamment de $\alpha$ et de graines aléatoires permet de résoudre indépendamment chaque problème (2) pour $\ell \in [L-1]$, et que les partitions trouvées forment des arbres $T^X$ et $T^\Omega$ valides. La figure 3 illustre les résultats de l'algorithmhe 2. En cas de succès, l'algorithmhe 3 retourne la même erreur d'approximation que $E_{bf}(\mathbf{A}, \tilde{\mathbf{P}}, \tilde{\mathbf{Q}})$ obtenue en connaissant $\tilde{\mathbf{P}}, \tilde{\mathbf{Q}}$. Dans le cas sans bruit, l'erreur relative atteint la précision machine. Dans le cas bruité, la figure 4 montre qu'elle est de l'ordre de $\epsilon$.

Lorsque $\tilde{\mathbf{A}}$ est la matrice TFD, l'algorithmhe 3 a également 100 % de succès pour tous les niveaux de bruit et pour $N \le 64$. Pour $N = 128$, la taux de succès est à 100 % pour le cas sans bruit, mais se dégrade dans le cas bruité comme illustré dans la

<table><caption>TABLE 1 : Taux de succès de l'algorithmhe 3 sur 20 instances.</caption><thead><tr><th>ε</th><th>0</th><th>0.01</th><th>0.03</th><th>0.1</th></tr></thead><tbody><tr><td>Butterfly orthogonale aléatoire (N = 128)</td><td>100 %</td><td>100 %</td><td>100 %</td><td>100 %</td></tr><tr><td>Transformée de Fourier discrète (N = 128)</td><td>100 %</td><td>95 %</td><td>90 %</td><td>50 %</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>FIGURE 3 : Résolution de (2) par l'algorithmhe 2 (avec $n = 50$) pour une matrice butterfly orthogonale aléatoire de taille 128 (orange, plein). Les barres d'erreur (négligeables) affichent les extrema, la croix indique la médiane. Bleu, tiret : erreur avec $T^X$, $T^\Omega$ connus. Vert, pointillé : erreur minimale sur 1000 tirages aléatoires de partitions, où l'on résout (2) en fixant les partitions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4 : Erreur relative divisée par $\epsilon$ de l'algorithmhe (3) pour approcher une matrice butterfly orthogonale aléatoire bruitée (moyennes, écart-type).</figcaption></figure>

table 1. La robustesse pourrait être améliorée par exemple en résolvant les problèmes (2) pour $\ell \in [L-1]$ conjointement, et non pas indépendamment.

**Conclusion** Nous avons proposé une heuristique pour identifier sans hypothèse analytique les partitions d'une matrice en blocs de rang faible permettant une factorisation butterfly. Lever les verrous de robustesse et de passage à l'échelle de l'heuristique permettrait à terme de chercher une factorisation butterfly pour des opérateurs utilisés en traitement de signal ou en apprentissage, comme les transformées de Fourier sur graphes [10] ou les couches de réseaux de neurones [3].

**Remerciements** Ce travail a été soutenu par le projet ANR AllegroAssai ANR-19-CHIA-0009.

# Références

[1] P. S. BRADLEY, K. P. BENNETT et A. DEMIRIZ : Constrained k-means clustering. *Microsoft Research, Redmond*, 2000.

[2] E. CANDES, L. DEMANET et L. YING : A fast butterfly algorithm for the computation of Fourier integral operators. *Multiscale Modeling & Simulation*, 2009.

[3] T. DAO, B. CHEN, N. S. SOHONI, A. DESAI, M. POLI, J. GROGAN, A. LIU, A. RAO, A. RUDRA et C. RÉ : Monarch : Expressive structured matrices for efficient and accurate training. *In ICML*, 2022.

[4] T. DAO, A. GU, M. EICHHORN, A. RUDRA et C. RÉ : Learning fast algorithms for linear transforms using butterfly factorizations. *In ICML*, 2019.

[5] W. HACKBUSCH : *Hierarchical matrices : algorithms and analysis*. Springer, 2015.

[6] Q.-T. LE, L. ZHENG, E. RICCIETTI et R. GRIBONVAL : Fast learning of fast transforms, with guarantees. *In ICASSP*, 2022.

[7] Y. LI, H. YANG, E. R. MARTIN, K. L. HO et L. YING : Butterfly factorization. *Multiscale Modeling & Simulation*, 2015.

[8] M. O'NEIL, F. WOOLFE et V. ROKHLIN : An algorithm for the rapid evaluation of special function transforms. *Applied and Computational Harmonic Analysis*, 2010.

[9] D. S. PARKER : *Random butterfly transformations with applications in computational linear algebra*. UCLA Comput. Sci. Department, 1995.

[10] D. I. SHUMAN, S. K. NARANG, P. FROSSARD, A. ORTEGA et P. VANDERGHEYNST : The emerging field of signal processing on graphs : Extending high-dimensional data analysis to networks and other irregular domains. *IEEE Signal Processing Magazine*, 2013.

[11] R. VIDAL, Y. MA et S.S. SASTRY : *Generalized Principal Component Analysis*. Springer New York, 2016.

[12] U. VON LUXBURG : A tutorial on spectral clustering. *Statistics and computing*, 2007.

[13] L. ZHENG, E. RICCIETTI et R. GRIBONVAL : Efficient identification of butterfly sparse matrix factorizations. *SIAM Journal on Mathematics of Data Science*, 2023.