# Approximation de faible rang de l'inverse d'un opérateur pour le pré-conditionnement de systèmes linéaires sur un espace produit tensoriel

Loïc Giraldi, Anthony Nouy, Grégory Legrain

► To cite this version:

Loïc Giraldi, Anthony Nouy, Grégory Legrain. Approximation de faible rang de l'inverse d'un opérateur pour le pré-conditionnement de systèmes linéaires sur un espace produit tensoriel. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717101

HAL Id: hal-01717101

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

Submitted on 25 Feb 2018

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_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Approximation de faible rang de l'inverse d'un opérateur pour le pré-conditionnement de systèmes linéaires sur un espace produit tensoriel

Loïc GIRALDI, Anthony NOUY, Grégory LEGRAIN

LUNAM Université, GeM, UMR CNRS 6183, École Centrale de Nantes, Université de Nantes,
{Loic.Giraldi, Anthony.Nouy, Gregory.Legrain}@ec-nantes.fr

**Résumé** — Des préconditionneurs sont proposés pour les systèmes linéaires sur un espace produit tensoriel. Ils sont définis comme des approximations de faible rang de l'inverse de l'opérateur et obtenus par la formulation d'un problème de meilleure approximation au sens d'une norme adaptée. Les pré-conditionneurs sont appliqués aux solveurs itératifs couplés aux méthodes d'approximation de tenseurs, ainsi qu'aux méthodes d'approximation directe de la solution dans des ensembles de tenseurs de faible rang. Cette stratégie est illustrée sur des exemples numériques tels que la résolution d'une équation de diffusion en grande dimension ou d'une équation paramétrique pour la quantification d'incertitudes.

**Mots clés** — Approximation de tenseurs, Solveurs itératifs, Préconditionneurs, Réduction de modèle, Proper Generalized Decomposition

## 1 Introduction

On s'intéresse ici aux systèmes linéaires issus de problèmes variationnels formulés sur des espaces produit tensoriel de grande dimension, tels que les problèmes régis par des équations aux dérivées partielles paramétriques ou stochastiques [14], de la forme

$$Au = b, \tag{1}$$

où $u$ et $b$ sont dans l'espace de dimension finie $V = \bigotimes_{\mu=1}^{d} V^{\mu} \simeq \bigotimes_{\mu=1}^{d} \mathbb{R}^{n_{\mu}}$ et $A \in \mathcal{L}(V) = \bigotimes_{\mu=1}^{d} \mathcal{L}(V^{\mu}) \simeq \bigotimes_{\mu=1}^{d} \mathbb{R}^{n_{\mu} \times n_{\mu}}$. Des approches non standard doivent être employées pour contourner la malédiction de la dimensionnalité associée à cette classe de problèmes. Des méthodes basées sur l'approximation de faible rang de la solution ont été récemment développées. On peut distinguer deux manières de construire une telle approximation. La première consiste à introduire des approximations de tenseurs à l'intérieur des solveurs itératifs [1, 12]. La méthode est simple à mettre en oeuvre, cependant les propriétés de convergence et de robustesse des solveurs sont souvent perdues. La seconde approche, parfois appelée Proper Generalized Decomposition (PGD) [2, 4], recherche directement l'approximation de la solution dans un ensemble de tenseurs de faible rang.

Afin d'utiliser ces solveurs, le système linéaire doit être préconditionné pour plusieurs raisons. Tout d'abord, les convergences des solveurs itératifs dépendent fortement du conditionnement de l'opérateur. De plus, étant donné que l'on introduit des approximations dans les solveurs, le préconditionnement a pour effet de diminuer la sensibilité du système vis-à-vis des perturbations. Dans le cadre des méthodes d'approximation directe de la solution basée sur une minimisation de résidu, le préconditionnement a pour effet d'améliorer la qualité de l'approximation obtenue. L'opérateur A s'écrivant usuellement

$$A = \sum_{i=1}^{r_A} A_i^1 \otimes \dots \otimes A_i^d,$$

les méthodes classiques de préconditionnement (e.g. factorisations incomplètes) conduisent à des opérateurs préconditionnés ne possédant plus une structure de faible rang. Pour palier ceci, des approches ont récemment été proposées où le préconditionneur est donné sous une forme tensorielle. On peut distinguer alors les préconditionneurs de rang 1 et de rang r. Dans la première catégorie, trois approches ont été proposées. Si le système linéaire provient d'une équation paramétrique stochastique, un préconditionneur
classiquement utilisé [8] est défini par l'inverse de l'opérateur moyen. On peut aussi utiliser l'inverse de la meilleure approximation de rang 1 de A [13, 18], ou directement la meilleure approximation de rang 1 de $A^{-1}$ [19] en minimisant la norme de $I - PA$. En ce qui concerne les préconditionneurs de rang supérieur, on retrouve deux approches dans la littérature consistant à s'appuyer sur une expression analytique de l'inverse [3, 11] (applicable à des problèmes spécifiques) ou à supposer une structure particulière du préconditionneur [17], typiquement $P_1 \otimes I_2 \otimes \dots \otimes I_d + \dots + I_1 \otimes \dots \otimes I_{d-1} \otimes P_d$.

On propose ici de construire un préconditionneur de manière automatique, avec un faible coût de calcul, sans hypothèse sur le rang ou la structure du préconditionneur, et avec des propriétés imposées, comme un caractère creux ou des symétries. La démarche introduite peut être utilisée pour corriger un autre préconditionneur. La méthode est appliquée à un problème de diffusion en grande dimension ainsi qu'à une équation paramétrique pour la quantification d'incertitudes.

# 2 Approximation de faible rang de l'inverse

## 2.1 Choix de la norme et imposition de propriétés

Afin de construire un préconditionneur avec les propriétés désirées, on s'appuie sur une approximation $P \in \mathcal{L}(\mathcal{V})$ de $A^{-1}$ telle que

$$ \|A^{-1} - P\| = \min_{B \in \mathcal{M}} \|A^{-1} - B\|, $$

où $\mathcal{M} \subset \mathcal{L}(\mathcal{V})$ est un sous-ensemble de tenseurs et $\| \cdot \|$ est une norme sur $\mathcal{L}(\mathcal{V})$. Afin d'approcher $A^{-1}$ a priori, on introduit une norme appropriée, pouvant différer suivant les propriétés de A. Si A est symétrique défini positif, on considère la norme induite par l'opérateur, c'est-à-dire $\|X\|_A^2 = \|XA^{1/2}\|_A^2 = \langle X, XA \rangle$. Si A est défini, on considère $\|X\|_{AA^*}^2 = \|XA\|_A^2 = \langle X, XAA^* \rangle$. On notera par la suite $\| \cdot \|_{AB}$ avec $B \in \{I, A\}$ suivant la norme choisie. Les normes proposées sont similaires aux normes vectorielles introduites dans différentes définitions de la PGD [15], ou encore aux normes utilisées pour les solveurs itératifs [16, section 6].

On montre que l'on peut aussi imposer différentes propriétés à P, comme des symétries ou un caractère creux, si ces dernières peuvent être caractérisées par l'appartenance de P à un sous-espace $U = \bigotimes_{\mu=1}^d U^\mu$, avec $U^\mu \subset \mathcal{L}(\mathcal{V}^\mu) = \mathbb{R}^{n_\mu \times n_\mu}$ [6]. On peut considérer par exemple les espaces de matrices symétriques pour $U^\mu$, ce qui donnera une approximation P symétrique.

## 2.2 Algorithmes de construction

### 2.2.1 Approche gloutonne

Afin que le coût de calcul de l'approximation reste raisonnable tout en controlant la qualité de l'approximation, on utilise une démarche gloutonne pour approximer $A^{-1}$. Cela consiste à construire une séquence d'approximations obtenues par corrections successives dans l'ensemble des tenseurs canoniques de rang un

$$ C_1 = \left\{ Q = \bigotimes_{\mu=1}^d Q^\mu; Q^\mu \in \mathcal{L}(\mathcal{V}^\mu) \right\}. $$

A l'itération r de cette construction, on trouve une approximation $P_r$ de rang-r telle que

$$ P_r = \sum_{i=1}^{r} Q_i, \quad Q_i = \bigotimes_{\mu=1}^{d} Q_i^\mu \in C_1. $$

La correction $Q_r$ telle que $P_r = P_{r-1} + Q_r$ est trouvée en résolvant le problème d'optimisation

$$ \min_{Q \in C_1} \|A^{-1} - P_{r-1} - Q\|_{AB}. \qquad (2) $$

L'approximation $P_r$ appartient alors à l'ensemble $C_r = \{\sum_{i=1}^r X_i; X_i \in C_1\}$ des tenseurs canoniques de rang $r$. Cet ensemble est fermé, et donc l'existence d'une meilleure approximation dans $C_r$ est garantie, si et seulement si $d=2$ ou $r=1$. C'est pour cette raison que l'on considère par la suite l'ensemble $C_1$.
### 2.2.2 Méthode de projection

La construction gloutonne pouvant être très sous-optimale, on propose alors une série de mises à jour de $P_r$ afin de diminuer l'erreur $\|A^{-1} - P_r\|_{AB}$ en tirant partie de formats alternatifs de tenseurs [7, 10]. Une méthode naturelle consiste à utiliser les tenseurs de Tucker. On projette alors $A^{-1}$ sur l'espace $\mathcal{U}_r = \bigotimes_{\mu=1}^{d} \mathcal{U}_r^\mu$ avec $\mathcal{U}_r^\mu = \text{span}\{Q_i^\mu; 1 \le i \le r\}$. On obtient alors une approximation de la forme

$$P_r = \sum_{i_1=1}^{r_1} \dots \sum_{i_d=1}^{r_d} \alpha_{i_1 \dots i_d} \bigotimes_{\mu=1}^{d} Q_{i_\mu}^{\mu}, \quad \alpha \in \mathbb{R}^{r_1 \times \dots \times r_d}.$$

La construction de $P_r$ est alors réalisée à partir de $P_{r-1}$ en deux étapes : on trouve le tenseur $Q_r$ solution de (2) puis on projette $A^{-1}$ sur $\mathcal{U}_r$. On remarque alors que la phase de projection souffre encore de la malédiction de la dimensionnalité. Celle-ci peut être contournée en approchant la projection grâce aux tenseurs hiérarchiques de Tucker.

### 2.2.3 Algorithme de projection approchée

La difficulté liée à la projection est contournée grâce aux tenseurs hiérarchiques de Tucker [9]. Celle-ci est approchée en exploitant une structure hiérarchique de l'espace $\mathcal{V} = \bigotimes_{\mu=1}^{d} \mathcal{V}^\mu$ comme illustrée en figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Représentation hiérarchique de l'espace produit tensoriel $\mathcal{V} = \mathcal{V}^1 \otimes \mathcal{V}^2 \otimes \mathcal{V}^3 \otimes \mathcal{V}^4$ et arbre des dimensions</figcaption></figure>

En notant $T$ l'arbre des dimensions, c'est-à-dire l'arbre supportant la décomposition de l'espace produit tensoriel choisie, on note $\mathcal{H}_r^T$ l'ensemble des tenseurs hiérarchiques de Tucker $v$ de rang $r = (r_t)_{t \in T}$ qui s'écrivent

$$v = \sum_{i=1}^{r_{D_1}} \sum_{j=1}^{r_{D_2}} \alpha_{ij}^{D_1} v_i^{D_1} \otimes v_j^{D_2}, \quad D = D_1 \cup D_2, D_1 \cap D_2 = \emptyset,$$

$$v_k^t = \sum_{i=1}^{r_1} \sum_{j=1}^{r_2} \alpha_{ijk}^t v_i^{t_1} \otimes v_j^{t_2}, \quad t = t_1 \cup t_2, t_1 \cap t_2 = \emptyset, \quad \forall t \in I(T),$$

$$v_k^\mu \in \mathcal{V}^\mu, \quad \forall \mu \in L(T),$$

où $D = \{1, \dots, d\}$, $I(T)$ représente les noeuds intérieurs de l'arbre des dimensions et $L(T)$ ses feuilles.

L'ensemble $\mathcal{H}_r^T$ étant un ensemble fermé [5], la phase de projection présentée en section 2.2.2 peut être remplacée par une optimisation sur la famille de tenseurs $(\alpha')_{t \in T}$ (d'ordre 2 ou 3). Pour cela on applique un algorithme de minimisations alternées sur $(\alpha')_{t \in T}$. On arrive alors à chaque itération à la résolution d'un système

$$A^t \alpha^t = b^t,$$

qui peut s'avérer être singulier. Dans ce cas, on prend une solution non nulle du problème de minimisation $\min_{\alpha'} \|A^t \alpha^t - b^t\|$.
# 3 Illustration

## 3.1 Problème de diffusion

On applique cette méthode au problème de Poisson en dimension $d$ sur un domaine $\Omega = \omega \times \dots \times \omega$ avec $\omega = [0, 1]$ dont la formulation variationnelle s'écrit

$$u \in H_0^1(\Omega), \quad \int_{\Omega} \nabla \delta_v \cdot \nabla u \, d\Omega = \int_{\Omega} \delta_v \, d\Omega, \quad \forall \delta_v \in H_0^1(\Omega).$$

On introduit des espaces éléments finis linéaires $V^\mu \subset H_0^1(\omega)$ avec $V^\mu = \operatorname{span}(\varphi_i)_{1 \le i \le n+1}$. Le problème discrétisé s'écrit alors

$$Au = b,$$

avec

$$A = \sum_{i=1}^{d} \bigotimes_{\mu=1}^{d} (\delta_{i\mu}K + (1 - \delta_{i\mu})M), \quad b = \bigotimes_{\mu=1}^{d} b_n,$$

$$K_{ij} = \int_{\omega} \nabla \varphi_i \cdot \nabla \varphi_j \, d\omega, \quad M_{ij} = \int_{\omega} \varphi_i \varphi_j \, d\omega, \quad (b_n)_i = \int_{\omega} \varphi_i d\omega,$$

où $\delta_{i\mu}$ est le symbole de Kronecker.

Le solveur utilisé est le gradient conjugué préconditionné tronqué dans l'ensemble des tenseurs hiérarchiques de rang 15 [12]. On considère ici un problème de grande dimension avec $d = 20$ et $n = 100$. Les symétries sont imposées au préconditionneur pour qu'il puisse être utilisé avec le gradient préconditionné. Les résultats sont présentés en figure 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Convergence du gradient conjugué préconditionné par une approximation symétrique de l'inverse avec un rang allant de 1 à 10, et un préconditionneur construit de manière gloutonne (+) ou par la méthode de projection approchée (○).</figcaption></figure>

On observe en figure 2 que tout d'abord, le fait de considérer des préconditionneurs de rang supérieur à 1 améliore grandement la qualité de l'approximation atteignable par l'algorithme du gradient conjugé préconditionné perturbé. De plus la phase de projection approchée sur l'espace réduit $U_r$ permet d'améliorer la qualité de l'approximation de l'inverse, ce qui se traduit par un préconditionnement de bonne qualité dans la figure 2.
## 3.2 Problème de transport stochastique

On considère un domaine spatial $\Omega = [0, 1] \times [0, 1]$, un domaine temporel $\Omega_t = [0, 3]$ et on considère le problème suivant :

$$ \frac{\partial u}{\partial t} - \kappa \Delta u + c(D \cdot \nabla u) + \sigma u = 0 \quad \text{sur } \Omega \times \Omega_t, $$

$$ u = 1 \quad \text{sur } \Gamma_1 \times \Omega_t, $$

$$ u = -1 \quad \text{sur } \Gamma_2 \times \Omega_t, $$

$$ u = 0 \quad \text{sur } \partial\Omega \setminus (\Gamma_1 \cup \Gamma_2) \times \Omega_t $$

$$ u = 0 \quad \text{sur } \Omega \times \{0\}, $$

où $\Gamma_1 = \{0\} \times [0, 1]$ et $\Gamma_2 = \{1\} \times [0, 1]$, et où $\kappa$ (resp. $c$, $\sigma$) est une variable aléatoire uniforme sur $[0.01, 1]$ (resp. $[0, 1]$, $[0.5, 1.5]$). D est un champ vectoriel tel que

$$ D(x_1, x_2) = \frac{1}{\|(x_1, x_2) - (0.5, 0.5)\|} \begin{pmatrix} (x_2 - 0.5) \\ -(x_1 - 0.5) \end{pmatrix} \quad \text{sur } \Omega \setminus \{(0.5, 0.5)\} \quad \text{et} \quad D(0.5, 0.5) = 0. $$

Le problème étant non symétrique, on utilise un solveur GMRES réinitialisé toutes les 5 itérations (GMRES(5)), dans lequel les itérés sont approchés dans l'ensemble $H_{20}^T(\mathcal{V})$ des tenseurs hiérarchiques de Tucker de rang 20 [1]. Le préconditionneur est construit à partir de la norme $\|\cdot\|_{AA^*}$. La convergence du solveur pour un préconditionneur de rang 5 construit de manière gloutonne ou grâce à la méthode de projection approchée est illustrée en figure 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Convergence de GMRES(5) préconditionné par une approximation de l'inverse de rang 5 construite de manière gloutonne ou par la méthode de projection approchée</figcaption></figure>

On observe d'une part que le solveur converge très rapidement vers le voisinage de la solution pour les deux préconditionneurs considérés, et d'autre part le gain substantiel apporté par la phase de projection approchée sur la précision de la solution.

## 4 Conclusion

On a proposé une méthode pour approcher l'inverse d'un opérateur. Cette approximation peut avoir des propriétés imposées telles que des symétries ou un caractère creux. De plus l'algorithme proposé améliore grandement la qualité d'une construction purement gloutonne pour un faible coût supplémentaire. Cette approximation de l'inverse est utilisée pour préconditionner des opérateurs définis possiblement symétriques positifs.
Références

[1] J. Ballani, L. Grasedyck. *A projection method to solve linear systems in tensor format*, Numerical Linear Algebra with Applications 20(1), 27-43, 2013.

[2] F. Chinesta, P. Ladevèze, E. Cueto. *A short review on model order reduction based on proper generalized decomposition*, Archives of Computational Methods in Engineering 18(4), 395-404, 2011.

[3] S. Dolgov, B. N. Khoromskij, I. Oseledets, E. Tyrtyshnikov. *A reciprocal preconditioner for structured matrices arising from elliptic problems with jumping coefficients*, Linear Algebra and its Applications 436(9), 2980-3007, 2012.

[4] A. Falcó, A. Nouy. *Proper generalized decomposition for nonlinear convex problems in tensor Banach spaces*, Numerische Mathematik 121(3), 503-530, 2012.

[5] A. Falcó, W. Hackbusch. *On Minimal Subspaces in Tensor Representations*, Foundations of Computational Mathematics 12(6), 765-803, 2012.

[6] L. Giraldi. *Contributions aux méthodes de calcul basées sur l'approximation de tenseurs et applications en mécanique numérique*, PhD thesis, École Centrale Nantes, 2012.

[7] L. Giraldi, A. Nouy, G. Legrain, P. Cartraud. *Tensor-based methods for numerical homogenization from high-resolution images*, Computer Methods in Applied Mechanics and Engineering 254, 154-169, 2013.

[8] R. G. Ghanem, R. M. Kruger. *Numerical solution of spectral stochastic finite element systems*, Computer Methods in Applied Mechanics and Engineering, 129(3), 289-303, 1996.

[9] W. Hackbusch, S. Kühn. *A new scheme for the tensor representation*, Journal of Fourier Analysis and Applications 15(5), 706-722, 2009.

[10] W. Hackbusch. *Tensor spaces and numerical tensor calculus*, Springer series in computational mathematics 42, 2012.

[11] B. Khoromskij. *Tensor-structured preconditioners and approximate inverse of elliptic operators in $\mathbb{R}^d$, Constructive Approximation* 34(2), 599-620, 2009.

[12] D. Kressner, C. Tobler. *Low-rank tensor krylov subspace methods for parametrized linear systems*, SIAM Journal on Matrix Analysis and Applications 32(4), 1288-1316, 2011.

[13] A. N. Langville, W. J. Stewart. *A kronecker product approximate preconditioner for SANs*. Numerical Linear Algebra with Applications 11 (8-9), 723-752, 2004.

[14] A. Nouy. *Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations*, Archives of Computational Methods in Engineering 16(3), 251-285, 2009.

[15] A. Nouy. *A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations*, Computer Methods in Applied Mechanics and Engineering 199 (23-24), 1603-1626, 2010.

[16] Y. Saad. *Iterative methods for sparse linear systems*, 2nd edition, SIAM, 2003.

[17] A. Touzene. *A tensor sum preconditioner for stochastic automata networks*. INFORMS Journal on Computing 20(2), 234-242, 2008.

[18] E. Ullmann. *A kronecker product preconditioner for stochastic galerkin finite element discretizations*. SIAM Journal on Scientific Computing 32(2), 923-946, 2010.

[19] E. Zander. *Tensor approximation methods for stochastic problems*. PhD thesis, TU Braunschweig, 2012.