An Improved Incremental Singular Value Decomposition and New Error Bounds

Résumé technique : Une décomposition en valeurs singulières incrémentale améliorée et de nouvelles bornes d'erreur

Énoncé du problème
La décomposition en valeurs singulières (SVD) incrémentale, initialement proposée par Brand [2, 3] et étendue aux produits scalaires pondérés par Fareed et al. [10], est une méthode standard pour l'approximation de rang faible dans les scénarios de données en flux où la matrice complète ne peut pas être stockée. L'algorithme met à jour une SVD tronquée ($U \approx Q\Sigma R^\top$) à l'arrivée de nouvelles colonnes. Cependant, en arithmétique à virgule flottante, chaque mise à jour multiplie la base singulière courante $Q$ par un petit facteur orthogonal. Bien que les multiplications individuelles préservent l'orthogonalité jusqu'à l'erreur d'arrondi, ces erreurs s'accumulent sur la longueur du flux $n$. Finalement, $Q$ perd son orthogonalité numérique, rendant les valeurs singulières inexactes et les critères de troncature peu fiables.

Le remède standard est une réorthogonalisation périodique (par exemple via le procédé de Gram–Schmidt modifié) déclenchée lorsque l'orthogonalité se dégrade. Cependant, la fréquence requise pour maintenir un niveau de précision prescrit est restée une question ouverte. Comme l'a noté Brand [3, p. 28], il était unclear à quelle fréquence la réorthogonalisation est nécessaire pour garantir un certain niveau global de précision numérique, en particulier dans des contextes pondérés (par exemple, la décomposition orthogonale propre des éléments finis) où la réorthogonalisation est coûteuse en calcul.

Méthodologie
Les auteurs proposent une restructuration de l'algorithme de SVD incrémentale qui diffère les mises à jour des facteurs singuliers courants. La méthodologie centrale repose sur deux observations :

1. Dominance non croissante du rang : Dans les applications de rang faible, la plupart des colonnes entrantes se situent dans le sous-espace courant (résidu < tol). Ces colonnes « non croissantes du rang » n'élargissent pas le sous-espace.
2. Regroupement (Batching) : L'effet cumulatif de $s$ colonnes non croissantes du rang successives peut être reproduit par une seule SVD fine d'une matrice tampon courte et large contenant les coefficients de projection, plutôt que d'effectuer $s$ SVD de matrices bordées distinctes.

L'algorithme proposé (Algorithme 5) maintient un tampon pour les coefficients de projection des colonnes non croissantes du rang. Le facteur courant $Q$ reste inchangé durant ces étapes. Lorsqu'une colonne « croissante du rang » arrive (résidu $\ge$ tol), l'algorithme :

• Vide le tampon en calculant une seule SVD fine des coefficients accumulés et en mettant à jour $Q$ et $R$ par lots.
• Traite la nouvelle colonne, en appliquant potentiellement un second test de troncature sur la plus petite valeur singulière de la matrice bordée locale pour rejeter le rang négligeable.
• Applique un procédé de Gram–Schmidt classique en deux passes (CGS-2) uniquement lors de ces événements de croissance du rang pour restaurer l'orthogonalité $W$ de la colonne ajoutée.

Cette restructuration réduit le nombre de grandes multiplications orthogonales appliquées à $Q$ de $O(n)$ (la longueur du flux) à $O(r + J_{E2})$, où $r$ est le rang numérique et $J_{E2}$ est le nombre d'événements de croissance du rang tronquants.

Contributions clés
L'article fournit deux avancées théoriques primaires et une amélioration algorithmique pratique :

1. Borne d'orthogonalité indépendante de la longueur du flux :
   Les auteurs prouvent (Théorème 4.1) que, sous un modèle probabiliste standard d'erreur d'arrondi, la perte d'orthogonalité $W$ du facteur gauche calculé $Q$ est bornée par :
   $$\|I - Q^\top W Q\|_2 \leq c_1 (r + J_{E2}) r \sqrt{m} \log(2/\delta) \varepsilon_{\text{mach}} + O(\varepsilon_{\text{mach}}^2)$$
   Crucialement, cette borne est indépendante de la longueur du flux $n$. Elle dépend de $m$ (dimension des colonnes) uniquement par un facteur $\sqrt{m}$. Cela résout la question ouverte de Brand : la réorthogonalisation n'est pas nécessaire à chaque étape ; la stratégie de regroupement proposée garantit que l'orthogonalité est maintenue à la précision machine uniformément en $n$.

2. Borne d'erreur directe précise :
   L'article établit une nouvelle borne d'erreur directe pour la factorisation dans la norme d'opérateur pondérée $W$ (Théorème 4.5) :
   $$\|U - Q\Sigma R^\top\|_W \leq \sqrt{n} \, \text{tol} + (\text{termes d'arrondi uniformes en } n)$$
   Cela affine la borne existante de Fareed et Singler [9], qui était $n \, \text{tol}$. Les auteurs démontrent via un exemple constructif que le taux $\sqrt{n}$ est atteignable, montrant que la borne précédente était lâche d'un facteur $\sqrt{n}$.

3. Efficacité algorithmique :
   En réduisant le nombre d'étapes coûteuses de réorthogonalisation et en regroupant les mises à jour, l'algorithme proposé réduit considérablement la surcharge de calcul.

Résultats
Des expériences numériques ont été menées sur un problème-test parabolique 2D (contextes d'éléments finis pondérés et non pondérés) et sur des données synthétiques.

• Orthogonalité : Alors que les algorithmes existants (mise à jour directe de Brand et mise à jour réorthogonalisée de Fareed) montraient une dégradation de l'orthogonalité croissante avec $n$ (dérivant vers $10^{-3}$ dans les cas pondérés), l'algorithme proposé a maintenu l'orthogonalité à la précision machine ($\sim 10^{-14}$) sur des longueurs de flux allant jusqu'à $n=10\,000$ et des dimensions de colonnes jusqu'à $m=100\,000$.
• Bornes d'erreur : Les expériences ont confirmé le taux de troncature $\sqrt{n} \, \text{tol}$. Sur un exemple constructif, l'erreur mesurée correspondait à la prédiction $\sqrt{n} \, \text{tol}$, tandis que la borne Fareed–Singler ($n \, \text{tol}$) était lâche d'un facteur croissant comme $\sqrt{n}$.
• Performance : L'algorithme proposé s'est exécuté 4,5 à 34 fois plus vite que ses concurrents les plus proches. Dans le cas pondéré, l'algorithme existant passait 94 % de son temps à la réorthogonalisation, tandis que l'algorithme proposé a réduit ce coût de trois ordres de grandeur.

Signification
L'article prétend résoudre la question ouverte concernant la fréquence de la réorthogonalisation dans la SVD incrémentale. En restructurant l'algorithme pour regrouper les mises à jour non croissantes du rang, il atteint :

1. Orthogonalité uniforme : La base courante reste numériquement orthogonale quelle que soit la longueur du flux, éliminant le besoin de vérifications heuristiques et périodiques de réorthogonalisation qui manquent de garanties a priori.
2. Analyse d'erreur plus serrée : Il fournit une borne d'erreur mathématiquement rigoureuse et précise ($\sqrt{n} \, \text{tol}$) qui améliore l'état de l'art d'un facteur $\sqrt{n}$.
3. Accélération pratique : La méthode offre des accélérations substantielles dans les applications dominées par les coûts de réorthogonalisation, telles que la POD d'éléments finis, sans sacrifier la précision.

Les auteurs notent que l'algorithme a déjà été appliqué dans des versions préliminaires aux EDP à dérivée fractionnaire temporelle, aux équations intégro-différentielles et à l'optimisation sous contraintes d'EDP, s'appuyant sur les garanties d'orthogonalité et d'erreur désormais établies. Un travail futur est suggéré pour affiner l'analyse de la dérive d'orthogonalité (actuellement conservatrice d'un facteur $\sim 10^3$) et pour optimiser l'étape de projection, qui domine actuellement le temps d'exécution.