# Compression tensorielle d'opérateurs intégraux en électromagnétisme basse fréquence : Premiers essais numériques

A. Torosyan¹, O. Chadebec¹, J.-M. Guichon¹, J.-R. Poirier², G. Meunier¹ et B. Bannwarth¹

Univ. Grenoble Alpes, CNRS, Grenoble INP, G2Elab, 38000 Grenoble, France¹

LAPLACE, Université de Toulouse, CNRS, INPT, UPS, Toulouse, France²

Email : anthony.torosyan@grenoble-inp.fr

**Résumé**—Cet article introduit de manière sommaire les décompositions tensorielles, plus particulièrement le format Tenseur-Train. Ce dernier est ensuite appliqué au cas d'une équation intégrale en électrostatique pour réduire le stockage de la matrice du système obtenu après discrétisation.

## I. INTRODUCTION

Pour la modélisation numérique en électromagnétisme, diverses méthodes de calculs existent pour résoudre les équations mises en jeux, en particulier les méthodes intégrales. Celles-ci présentent l'avantage de ne pas nécessiter une discrétisation de l'air mais l'inconvénient majeur de mener à un système linéaire où la matrice est pleine. Ceci implique que la complexité algorithmique de l'assemblage est en $O(n^2)$ où $n$ est le nombre d'inconnues.

Plusieurs méthodes de compression matricielle ont été mises en place pour réduire la complexité de l'intégration et aussi réduire le stockage en mémoire de l'opérateur intégral. Nous pouvons citer notamment les fast multipole method (FMM) [1] et les matrices hiérarchiques (H-matrices) [4].

Dans ce travail, nous n'adoptons pas le paradigme de la compression matricielle mais plutôt celui de la compression tensorielle. L'idée globale est de réaliser une approximation de bas rang d'un tenseur pour réduire le coût de construction (par exemple d'intégration dans notre cas) et du stockage en mémoire de celui-ci. Divers formats de tenseurs existent, en particulier les décompositions polyadiques canonique (CP) [2] et Tucker [3]. L'arrivée du format Tenseur-Train (TT), popularisé par Oseledets aux alentours des années 2010 [7] permet d'outrepasser certaines difficultés des deux décompositions précédentes, ouvrant la voie à la compression d'équations intégrales par des méthodes tensorielles.

En électromagnétisme haute fréquence, quelques premières tentatives ont été réalisées [6] mais les noyaux de Green intervenant dans ces équations induisent des limitations.

Dans cet article, nous nous plaçons en basse fréquence, où les noyaux sont plus favorables (meilleure régularité notamment) et nous présentons une compression sur un premier exemple simple.

## II. LE FORMAT TENSEUR-TRAIN (TT)

### A. Généralités

Soit $\mathcal{A} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ un tenseur à $d$ dimensions de tailles respectives $n_1, \dots, n_d$ et dont le coefficient relatif au

multi-indice $(i_1, \dots, i_d)$ sera noté $\mathcal{A}(i_1, \dots, i_d)$. L'écriture au format TT de $\mathcal{A}$ est donnée par :

$$ \mathcal{A}(i_1, \dots, i_d) = \sum_{\alpha_1, \dots, \alpha_{d-1}} C_1(\alpha_0, i_1, \alpha_1) \dots C_d(\alpha_{d-1}, i_d, \alpha_d) \quad (1) $$

où les $C_k \in \mathbb{R}^{r_{k-1} \times n_k \times r_k}$ et les $r_k, k \in \{1, \dots, d\}$ (avec $r_0 = r_d = 1$ pour uniformiser les notations des coeurs) sont appelés respectivement les coeurs et les TT-rangs de $\mathcal{A}$. Il existe aussi une écriture équivalente faisant intervenir le produit matriciel :

$$ \mathcal{A}(i_1, \dots, i_d) = C_1(i_1) \dots C_d(i_d) \quad (2) $$

avec $C_k(i_k) \in \mathbb{R}^{r_{k-1} \times r_k}, k \in \{1, \dots, d\}$. Cette dernière écriture montre que le format TT s'interprète comme une généralisation de la séparation des variables : à chaque coeur $C_k$ est associé un des indices $i_k$ du tenseur $\mathcal{A}$ pour $k \in \{1, \dots, d\}$.

Le nombre d'éléments ($N_{TT}$) d'un tenseur au format TT est obtenu en additionnant celui de chaque coeur :

$$ N_{TT} = \sum_{k=1}^{d} r_{k-1} n_k r_k \quad (3) $$

En optant pour $r_k = r$ et $n_k = n$, $k \in \{1, \dots, d\}$, le nombre d'éléments au format TT est en $O(dnr^2)$ là où celui d'un tenseur classique (noté $N_{tens}$) est en $O(n^d)$.

### B. Algorithme de construction

La construction d'une approximation de $\mathcal{A}$ à l'aide de son écriture en format TT est donnée dans [7] :

- La première étape est de représenter $\mathcal{A}$ par $\mathbf{A}^{(1)} \in \mathbb{R}^{n_1 \times n_2 \times \dots \times n_d}$ et de chercher une approximation de rang $r_1$ de $\mathbf{A}^{(1)}$. En général, $\mathbf{A}^{(1)} \approx \mathbf{X}_1 \mathbf{Y}_1^T$ avec $\mathbf{X}_1 \in \mathbb{R}^{n_1 \times r_1}$ et $\mathbf{Y}_1 \in \mathbb{R}^{n_2 \times n_d \times r_1}$. Le coeur $C_1$ est obtenu en représentant $\mathbf{X}_1$ par $\mathcal{X}_1 \in \mathbb{R}^{r_0 \times n_1 \times r_1}$ ($r_0 = 1$) et en posant $C_1 = \mathcal{X}_1$. L'algorithme se poursuit avec le tenseur à $d-1$ dimensions $W_2 \in \mathbb{R}^{r_1 n_2 \times n_3 \dots n_d}$ constitué des termes de $\mathbf{Y}_1$.

- A l'étape $k$, le tenseur $W_k \in \mathbb{R}^{r_{k-1} n_k \times n_{k+1} \times \dots \times n_d}$ est représenté par $\mathbf{W}_k^{(1)} \in \mathbb{R}^{r_{k-1} n_k \times n_{k+1} \times \dots \times n_d}$. L'approximation de rang $r_k$ donne $\mathbf{W}_k^{(1)} \approx \mathbf{X}_k \mathbf{Y}_k^T$ avec $\mathbf{X}_k \in \mathbb{R}^{r_{k-1} n_k \times r_{k+1} \times \dots \times r_d}$ et $\mathbf{Y}_k \in \mathbb{R}^{n_{k+1} \times n_d \times r_{k+1} \times \dots \times r_d}$. Le coeur $C_k$ est obtenu en représentant $\mathbf{X}_k$ par $\mathcal{X}_k \in \mathbb{R}^{r_{k-1} n_k \times r_{k+1} \times \dots \times r_d}$ et en posant $C_k = \mathcal{X}_k$. L'algorithme se poursuit avec le tenseur à $d-k$ dimensions $W_{k+1} \in \mathbb{R}^{r_k n_{k+1} \times r_{k+2} \times \dots \times n_d}$ constitué des termes de $\mathbf{Y}_k$.
Un choix de décomposition classique serait de déterminer l'approximation de bas rang en tronquant la Décomposition en Valeurs Singulières (SVD) de chacune des matrices successives. Ceci n'est pas judicieux car cela nécessite de calculer entièrement A. Une solution à ce problème consiste à utiliser une décomposition CUR où C, R contiennent respectivement r colonnes et r lignes linéairement indépendantes d'une matrice et U est l'intersection entre C et R. Cet algorithme se nomme le TT-cross [8].

### C. Tensorisation

Après la discrétisation d'une équation intégrale, intervient la résolution d'un système linéaire $Ax = y$ avec $A \in \mathbb{R}^{N \times N}$, $x \in \mathbb{R}^N$ et $y \in \mathbb{R}^N$ désignant respectivement la matrice de l'opérateur intégral, les inconnues et le vecteur contenant les données (la source). Nous constatons que $A$ ne permet pas de tirer profit du format TT. Pour cela, il faut tensoriser le système linéaire, c'est-à-dire transformer les matrices et les vecteurs en tenseurs de dimensions supérieures à deux. Ceci est possible en supposant que $N = \prod_{k=1}^{d} n_k$. Ainsi, un isomorphisme [5] est construit entre les espaces $\mathbb{R} \prod_{k=1}^{d} n_k$ et $\mathbb{R}^{n_1 \times \cdots \times n_d}$ pour les vecteurs ainsi qu'entre $\mathbb{R} \prod_{k=1}^{d} n_k \times \prod_{k=1}^{d} n_k$ et $\mathbb{R}^{n_1^2 \times \cdots \times n_d^2}$ pour les matrices.

### III. COMPRESSION D'UNE ÉQUATION INTÉGRALE ISSUE DE L'ÉLECTROSTATIQUE

Soit le cercle unité $\mathbb{S} \subset \mathbb{R}^2$ et l'équation intégrale suivante dans le cas électrostatique [10] :

$$ \int_{\mathbb{S}} G(\mathbf{x}, \mathbf{y}) q(\mathbf{y}) \, d\gamma(\mathbf{y}) = \phi(\mathbf{x}) \quad (4) $$

avec $\phi$ le potentiel électrique sur le cercle (imposé constant dans notre cas), $q$ la densité linéique de charges solution du problème et $G(\mathbf{x}, \mathbf{y}) = -\ln(||\mathbf{x} - \mathbf{y}||_2)/(2\pi)$ la fonction de Green du Laplacien en dimension deux ($|| \cdot ||_2$ étant la norme euclidienne).

La discrétisation est obtenue en considérant que $\mathbb{S} = \bigcup_{k=1}^{N} \Gamma_k$ où les $\Gamma_k$ sont des segments de longueurs identiques. La solution $q$ est recherchée dans une base de fonctions $u_k$, $k \in \{1, \dots, N\}$ d'ordre 0 (c'est-à-dire que $q = \sum_{k=1}^{N} \lambda_k u_k$ où les $\lambda_k \in \mathbb{R}$).

Par une méthode de collocation réalisée au barycentre de chaque élément, le système discrétisé devient :

$$ \sum_{k=1}^{N} \int_{\Gamma_k} G(\mathbf{x}_k, \mathbf{y}) u_k(\mathbf{y}) \, d\gamma(\mathbf{y}) \lambda_k = \phi(\mathbf{x}_k), \quad \forall \mathbf{l} \in \{1, \dots, N\} \quad (5) $$

où les $\mathbf{x}_k$ sont les points milieux des segments $\Gamma_k$.

Pour l'expérience numérique, la librairie TTpy [9] est utilisée pour compresser la matrice de l'opérateur intégral (l'algorithme du TT-cross présenté précédemment) et le nombre d'connues dans le maillage $N$ est choisi comme une puissance de 2 [6] (c'est-à-dire que $N = 2^d$ où $d \in \mathbb{N}$).

En faisant le rapport $N_{TT}/N_{tens}$ à N fixé (avec $N_{tens} = N^2$), nous obtenons des taux de compression de l'ordre de $10^{-6}$ pour le plus grand N (Voir Fig.1).

<figure><img src="image_1.png" /><figcaption>FIGURE 1. Comparaison du nombre de coefficients entre la matrice issue de l'opérateur intégral et son approximation au format TT.</figcaption></figure>

### IV. CONCLUSION

Cette première expérience numérique sur une géométrie simple est prometteuse, les taux de compression sont élevés pour les plus grands $N$. Cependant, un point reste limitant : la contrainte du choix du nombre d'éléments.

L'étape suivante sera de résoudre le système linéaire. Le format TT donne accès à un produit matrice-vecteur, rendant possible l'utilisation de GMRES, mais il existe aussi l'algorithme itératif AMEN [11] conçu pour le format TT.

### REMERCIEMENT

Ce travail est supporté par le projet ANR TensorVim 22-CE50-0028.

### RÉFÉRENCES

[1] Rokhlin, V. (1985). Rapid solution of integral equations of classical potential theory. Journal of computational physics, 60(2), 187-207.

[2] Hitchcock, F. L. (1927). The expression of a tensor or a polyadic as a sum of products. Journal of Mathematics and Physics, 6(1-4), 164-189.

[3] Tucker, L. R. (1966). Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3), 279-311.

[4] Hackbusch, W. (2015). Hierarchical matrices : algorithms and analysis (Vol. 49). Heidelberg : Springer.

[5] Hackbusch, W. (2012). Tensor spaces and numerical tensor calculus (Vol. 42, pp. p-xxiv). Berlin : Springer.

[6] Poirier, J. R., Coulaud, O., Kaya, O. (2020). Fast BEM solution for 2-D scattering problems using quantized tensor-train format. IEEE Transactions on Magnetics, 56(3), 1-4.

[7] Oseledets, I. V. (2011). Tensor-train decomposition. SIAM Journal on Scientific Computing, 33(5), 2295-2317.

[8] Oseledets, I., Tyryshnikov, E. (2010). TT-cross approximation for multi-dimensional arrays. Linear Algebra and its Applications, 432(1), 70-88.

[9] I. V. Oseledets. « tpy », 2012. https://github.com/oseledets/tpy.

[10] Nédélec, Jean-Claude. Acoustic and electromagnetic equations : integral representations for harmonic problems. Vol. 144. New York : Springer, 2001.

[11] Dolgov, S. V., Savostyanov, D. V. (2014). Alternating minimal energy methods for linear systems in higher dimensions. SIAM Journal on Scientific Computing, 36(5), A2248-A2271.