Accelerating operator Sinkhorn iteration with overrelaxation

Résumé technique : Accélération de l'itération de Sinkhorn des opérateurs par sur-relaxation

Énoncé du problème
L'article aborde le problème de mise à l'échelle des opérateurs, qui cherche à trouver des matrices carrées inversibles $L \in GL_m(\mathbb{R})$ et $R \in GL_n(\mathbb{R})$ transformant conjointement une collection donnée de matrices $A_1, \dots, A_k \in \mathbb{R}^{m \times n}$ en une forme « mise à l'échelle » $\bar{A}_i = L A_i R^\top$. L'objectif est de satisfaire deux conditions simultanées :

1. $\sum_{i=1}^k \bar{A}_i \bar{A}_i^\top = \frac{1}{m} I_m$
2. $\sum_{i=1}^k \bar{A}_i^\top \bar{A}_i = \frac{1}{n} I_n$

Ce problème généralise la mise à l'échelle classique des matrices et apparaît dans des applications telles que le test d'identité de polynômes non commutatifs, la théorie des invariants et le traitement du signal. La méthode de résolution standard est l'itération de Sinkhorn des opérateurs (OSI), un algorithme itératif à point fixe alterné. Bien que l'OSI soit globalement convergente sous certaines conditions, sa convergence peut être lente. Les auteurs visent à accélérer ce processus en utilisant la sur-relaxation successive (SOR), une technique bien établie pour les systèmes linéaires mais moins explorée dans ce contexte spécifique de mise à l'échelle non linéaire d'opérateurs.

Méthodologie
Les auteurs proposent trois variantes distinctes de sur-relaxation pour l'itération de Sinkhorn des opérateurs, toutes basées sur l'interprétation du problème comme une minimisation alternée sur des cônes de matrices définies positives (PD).

1. Sur-relaxation euclidienne (PD-OR) : Cette approche applique des combinaisons affines standard aux itérés $(X_t, Y_t)$, où $X_t = L_t^\top L_t$ et $Y_t = R_t^\top R_t$. La règle de mise à jour est $(1-\omega)X_t + \omega S_1(Y_t)$, où $S_1$ est l'application de mise à l'échelle. Bien que computationnellement simple, cette méthode manque de garanties théoriques assurant que les matrices résultantes restent définies positives lorsque le paramètre de relaxation $\omega > 1$.
2. Sur-relaxation par transformation de coordonnées (Cholesky-OR) : Pour garantir que les itérés restent dans le cône PD, les auteurs suggèrent d'appliquer la sur-relaxation dans un espace de coordonnées transformé. Ils utilisent spécifiquement les facteurs de Cholesky ($L$ et $R$) des matrices PD. Puisque les facteurs de Cholesky appartiennent à des espaces linéaires (matrices triangulaires), les combinaisons affines sont bien définies et préservent la structure requise pour l'algorithme. Cette variante est présentée comme la plus efficace numériquement.
3. Sur-relaxation géodésique (Geodesic-OR) : Cette méthode respecte la géométrie hyperbolique des cônes PD induite par la métrique de Hilbert. Au lieu de combinaisons affines, les auteurs définissent des mises à jour le long des géodésiques en utilisant la moyenne géométrique des opérateurs $X \#_\omega \tilde{X}$. Cette approche est théoriquement robuste, préservant l'équivariance orbitale du problème et assurant que les itérés restent dans le cône PD.

Contributions clés et résultats théoriques

• Analyse de convergence locale : Les auteurs démontrent qu'à proximité d'un point fixe, les trois variantes de sur-relaxation sont équivalentes au premier ordre. En linéarisant les itérations et en utilisant l'équivariance orbitale (invariance sous le changement d'échelle de l'orbite de la solution), ils dérivent une formule unifiée pour le taux de convergence locale asymptotique.

  • Ils établissent que le rayon spectral de l'itération SOR linéarisée suit le motif classique du théorème de Young pour la SOR.
  • Ils fournissent une formule pour le paramètre de relaxation asymptotiquement optimal $\omega_{opt}$, qui dépend du rayon spectral $\beta$ de l'itération standard (non relaxée) de Sinkhorn des opérateurs : $\omega_{opt} = \frac{2}{1 + \sqrt{1-\beta^2}}$.
  • Cette analyse repose sur le fait que le Hessien de la fonction objectif géodésiquement convexe sous-jacente possède un espace nul unidimensionnel correspondant à l'orbite de mise à l'échelle, une propriété qu'ils prouvent rigoureusement.

• Convergence globale pour la SOR géodésique : Bien que la convergence globale pour la SOR non linéaire soit généralement difficile, les auteurs prouvent la convergence globale pour la variante Geodesic-OR dans une plage spécifique de $\omega$.

  • Ils utilisent la métrique de Hilbert et la contractivité de l'itération standard de Sinkhorn.
  • Une contribution novatrice est une inégalité métrique pour la métrique de Hilbert sur les matrices PD (Lemme 3.4), qui borne la distance d'un point géodésique par rapport à un troisième point. Cette inégalité leur permet d'étendre les résultats de contractivité au cas sur-relaxé, prouvant la convergence pour $\omega \in (0, \frac{2}{1+\sqrt{\Lambda_1 \Lambda_2}})$, où $\Lambda_1, \Lambda_2$ sont les constantes de Lipschitz des applications de mise à l'échelle.

• Sélection adaptative du paramètre : Reconnaissant que le $\omega$ optimal dépend du taux de convergence inconnu de l'instance spécifique, les auteurs adoptent une stratégie adaptative. Ils estiment le taux de convergence $\beta$ à partir des itérations initiales de la méthode standard et mettent à jour dynamiquement $\omega$ vers $\omega_{opt}$, une technique adaptée d'un travail récent sur la mise à l'échelle des matrices.

Résultats numériques
Les auteurs ont mené des expériences sur deux types d'instances :

1. Mise à l'échelle de cadres (Frame Scaling) : Dans cette application, les méthodes SOR (en particulier Cholesky-OR et Geodesic-OR) ont démontré une convergence nettement plus rapide que l'OSI standard. Les méthodes ont atteint des erreurs de l'ordre de $10^{-14}$ en environ la moitié des itérations requises par l'OSI pour atteindre $10^{-8}$. Les taux de convergence observés correspondaient étroitement aux prédictions théoriques.
2. Opérateurs mal conditionnés : En utilisant des matrices de Hilbert comme entrées, les auteurs ont observé que, bien que la SOR accélère la convergence initialement, elle stagne à une précision inférieure ($10^{-6}$) par rapport à l'OSI standard ($10^{-11}$). Ils attribuent cela au fait que les variantes SOR basées sur des matrices intermédiaires « demi-mises à l'échelle » ne compensent pas complètement le mauvais conditionnement initial aussi efficacement que les mises à jour alternées complètes de l'OSI standard.

Importance et affirmations
L'article prétend généraliser les résultats existants sur la mise à l'échelle des matrices (spécifiquement les travaux de Thibault et al. et Lehmann et al.) au cadre plus complexe de la mise à l'échelle des opérateurs.

• Insight théorique : Les auteurs soulignent que leur analyse rigoureuse de la convergence locale, qui gère la nature non définie positive du Hessien due à l'équivariance de groupe, apporte un nouvel éclairage sur la structure des problèmes de mise à l'échelle des opérateurs.
• Utilité pratique : Les expériences numériques suggèrent que la sur-relaxation peut accélérer la mise à l'échelle des opérateurs de plusieurs ordres de grandeur dans des scénarios bien conditionnés (comme la mise à l'échelle de cadres) avec une surcharge computationnelle négligeable, en particulier lors de l'utilisation de l'approche basée sur Cholesky.
• Limites : Les auteurs reconnaissent modestement que leur preuve de convergence globale pour la version géodésique couvre une plage de $\omega$ probablement trop pessimiste (souvent proche de 1), et que les implémentations SOR actuelles peuvent souffrir d'une précision numérique réduite sur des entrées mal conditionnées. Ils identifient l'établissement d'une convergence globale pour une plage plus large de $\omega$ et l'amélioration de la précision sur les problèmes mal conditionnés comme des pistes de travail futur.