Subspace gradient descent method for linear tensor equations

1. Problématique

L'article s'intéresse à la résolution numérique d'équations tensorielles linéaires de la forme :
$$\mathcal{A}(x) = c$$
où $\mathcal{A}$ est un opérateur linéaire agissant sur un espace de tenseurs $\mathbb{R}^{n_1 \times \dots \times n_d}$, supposé symétrique et défini positif. Les tenseurs $x$ (solution) et $c$ (second membre) sont d'ordre $d$.

Ces équations apparaissent dans divers domaines tels que la chimie quantique, les mathématiques financières, l'imagerie et les réseaux de neurones profonds, ainsi que dans la discrétisation d'équations aux dérivées partielles (EDP) de dimension trois ou supérieure.

Le défi principal réside dans la « malédiction de la dimensionnalité » : la taille des tenseurs croît exponentiellement avec le nombre de modes $d$. Pour contourner cela, les méthodes existantes utilisent des formats de faible rang (Tucker, CP, TT, HT). Cependant, les solveurs actuels, comme l'algorithme AMEn (Alternating Minimal Energy) basé sur le format Tensor-Train (TT), peuvent souffrir de problèmes de convergence (minima locaux) ou de gestion complexe de la mémoire et de la précision.

2. Méthodologie

Les auteurs proposent deux nouvelles méthodes de type descente de gradient adaptées aux équations tensorielles, en généralisant l'approche « Subspace Conjugate Gradient » (Ss-cg) développée pour les matrices.

A. Format de représentation : Tucker

Contrairement à l'approche AMEn qui utilise le format Tensor-Train (TT), cette méthode adopte le format de décomposition de Tucker. Ce choix est motivé par :

• L'efficacité du format Tucker pour un nombre modéré de modes (ici $d=3$).
• La facilité d'identification des sous-espaces de projection nécessaires à la mise à jour des itérés.
  Un tenseur $x$ est représenté par un tenseur cœur $C$ et des matrices factorielles $U_k$ : $x = C \times_1 U_1 \times_2 \dots \times_d U_d$.

B. Algorithmes proposés

Les auteurs dérivent deux algorithmes itératifs basés sur la minimisation d'une fonctionnelle convexe $\Phi(x) = \frac{1}{2}\langle x, \mathcal{A}(x) \rangle - \langle x, c \rangle$ :

1. Tk–ss–cg (Subspace Conjugate Gradient Tensor) :

   • Généralisation tensorielle de la méthode du gradient conjugué.
   • À chaque itération $k$, une nouvelle direction de recherche est construite en combinant le résidu courant et la direction précédente, avec une condition d'orthogonalité par rapport à l'opérateur $\mathcal{A}$ (A-orthogonalité).
   • Cela nécessite la résolution d'équations tensorielles réduites pour trouver les coefficients de mise à jour ($\alpha_k$ et $\beta_k$).

2. Tk–ss–sd (Subspace Steepest Descent Tensor) :

   • Généralisation tensorielle de la méthode du gradient à pente la plus raide.
   • La direction de recherche est simplement le résidu courant (ou son approximation préconditionnée).
   • Cela évite le calcul du coefficient $\beta_k$, réduisant ainsi le coût par itération, au prix d'une convergence asymptotiquement plus lente.

C. Stratégies d'implémentation et d'optimisation

Pour rendre ces méthodes viables et efficaces, plusieurs techniques avancées sont employées :

• Troncature de rang : Les opérations (sommes, produits) augmentent le rang des tenseurs. Des stratégies de troncature (SVD tronquée, ST-HOSVD) sont appliquées à chaque étape pour maintenir le rang multilinéaire sous une limite fixée (maxrank).
• Calcul en précision mixte : La résolution des systèmes réduits pour les coefficients ($\alpha_k, \beta_k$) et l'application des préconditionneurs sont effectuées en précision simple (single precision), tandis que le reste du calcul reste en double précision. Cela réduit l'utilisation de la mémoire et accélère les calculs sans perte significative de précision finale.
• Préconditionnement : Pour accélérer la convergence, plusieurs préconditionneurs sont testés :
  • P–Eig : Basé sur la factorisation spectrale de l'opérateur Laplacien discret.
  • P–FFT : Une version optimisée utilisant la Transformée en Sinus Discrète (DST) pour exploiter la structure de Kronecker de l'opérateur Laplacien.
  • P–InnOut : Une stratégie « intérieur-extérieur » où quelques itérations d'un solveur rapide (Tk–ss–sd) sont utilisées pour approximer l'action de l'inverse de l'opérateur.

3. Contributions Clés

1. Extension des méthodes de sous-espace : Transfert réussi de la méthode Ss-cg (initialement pour les matrices) au cadre tensoriel avec le format Tucker.
2. Nouvelle variante Steepest Descent : Proposition d'une version tensorielle du gradient à pente la plus raide, offrant un compromis coût/convergence intéressant.
3. Intégration de stratégies hybrides : Combinaison efficace du format Tucker, de la précision mixte et de préconditionneurs adaptés aux EDP elliptiques.
4. Comparaison rigoureuse : Évaluation systématique contre l'état de l'art (AMEn) sur des problèmes de dimension 3.

4. Résultats Expérimentaux

Les expériences ont été menées sur des discrétisations d'EDP elliptiques en 3D (Poisson, diffusion avec coefficients variables, diffusion-réaction, coefficients discontinus) avec des grilles allant jusqu'à $n=1000$.

• Performance globale : Les méthodes proposées (notamment Tk–ss–sd préconditionné) sont compétitives, voire supérieures, à l'algorithme AMEn en termes de temps CPU total.
• Efficacité du préconditionnement P–FFT : Le préconditionneur basé sur les FFT s'avère le plus performant, réduisant le nombre d'itérations à 2 ou 3, ce qui rend la méthode extrêmement rapide (ex: 0.05s pour $n=1000$ contre 0.75s pour AMEn sur le problème de Poisson).
• Comparaison CG vs SD :
  • Sans préconditionnement, Tk–ss–cg converge en moins d'itérations mais Tk–ss–sd peut être plus rapide en temps CPU si le nombre d'itérations reste faible.
  • Avec préconditionnement, Tk–ss–sd devient la méthode de choix car la réduction drastique du nombre d'itérations compense sa convergence théorique plus lente.
• Robustesse : Les méthodes montrent une bonne robustesse face à des coefficients non constants et discontinus, bien que le préconditionneur P–Eig montre une certaine sensibilité à la troncature du rang pour les résidus très faibles.

5. Signification et Conclusion

Ce travail démontre que les méthodes de descente de gradient dans des sous-espaces, couplées au format Tucker, constituent une alternative puissante aux méthodes d'optimisation alternée (comme AMEn) pour les équations tensorielles symétriques.

L'approche proposée offre :

• Une gestion transparente de la mémoire et de la précision.
• Des garanties de convergence plus claires que les méthodes basées sur l'optimisation non convexe (risque de minima locaux).
• Une compétitivité en temps de calcul, en particulier lorsque des préconditionneurs adaptés à la structure du problème (comme P–FFT) sont utilisés.

Les auteurs concluent que cette classe de méthodes est prête pour des applications pratiques sur des problèmes de grande dimension issus de la physique et de l'ingénierie, offrant un équilibre optimal entre coût computationnel et précision.