Fast Evaluation of Truncated Neumann Series by Low-Product Radix Kernels

1. Problématique

L'article aborde le problème de l'évaluation efficace de la série de Neumann tronquée pour une matrice $A$ dont le rayon spectral $\rho(A) < 1$ :
$$S_k(A) = \sum_{j=0}^{k-1} A^j = I + A + A^2 + \dots + A^{k-1}$$
Cette série est fondamentale en algèbre linéaire numérique pour l'inversion de matrices approchée, le préconditionnement polynomial, l'estimation du déterminant logarithmique et la détection de signaux MIMO.

Dans les environnements denses, le coût dominant est la multiplication matrice-matrice (GEMM). L'évaluation naïve nécessite $k-1$ produits. Les méthodes de décomposition classiques (comme la décomposition binaire ou « repeated squaring ») réduisent ce coût à $O(\log k)$, mais avec un coefficient directeur élevé :

• Décomposition binaire (radix-2) : $2 \log_2 k$ produits.
• Décomposition ternaire (radix-3) : $\approx 1.89 \log_2 k$ produits.
• Décomposition quinaire (radix-5) : $\approx 1.72 \log_2 k$ produits (méthode de Gustafsson et al.).

L'objectif est de réduire davantage ce coefficient directeur en explorant des radices supérieurs (radix-9, radix-15, etc.). Cependant, deux obstacles majeurs existent :

1. L'absence de constructions explicites de « noyaux » (kernels) exacts pour des radices supérieurs à 5.
2. L'incertitude sur la possibilité de construire des noyaux exacts pour des radices très élevés, et la gestion des erreurs si l'on utilise des approximations.

2. Méthodologie

L'auteur propose une approche basée sur la construction de noyaux radiciels (radix kernels) et un cadre général pour gérer les approximations.

A. Définition du Noyau

Le calcul de $S_{mn}(A)$ peut être décomposé en $S_n(A) \cdot T_m(A^n)$, où $T_m(B) = I + B + \dots + B^{m-1}$. Le coût d'une mise à jour radice-$m$ est donné par $C(m) = \mu_m + 2$, où $\mu_m$ est le nombre minimal de produits matriciels nécessaires pour évaluer le noyau $T_m$. L'efficacité asymptotique est mesurée par le coefficient $\frac{\mu_m + 2}{\log_2 m}$.

B. Construction du Noyau Radix-9 (Exact)

Pour le radix-9, l'auteur dérive un noyau exact utilisant uniquement 3 produits matriciels ($\mu_9 = 3$), atteignant la borne inférieure théorique $\lceil \log_2(9-1) \rceil = 3$.

• Technique : Factorisation de $T_9(B)$ sous la forme d'un produit $P \cdot Q$ où $P$ et $Q$ sont des combinaisons linéaires de termes de base ($B, B^2, B^3+2B^4$).
• Résultat : Les coefficients sont rationnels, permettant un calcul exact (arithmétique rationnelle ou fixe). Cela permet d'obtenir un coefficient global de $5 \log_9 k \approx 1.58 \log_2 k$.

C. Construction du Noyau Radix-15 (Approximatif)

Pour le radix-15, la recherche d'un noyau exact à 4 produits ($\mu_{15}=4$) s'avère impossible par des méthodes symboliques. L'auteur utilise une optimisation numérique (L-BFGS-B) pour trouver un noyau qui correspond aux termes de degré 0 à 14, mais qui génère des termes d'ordre supérieur non nuls (phénomène de « débordement » ou spillover).

• Problème : Le débordement brise l'identité télescopique standard nécessaire pour extraire la puissance suivante $A^{mn}$.
• Solution (Cadre Général) : L'auteur introduit un cadre basé sur le résidu. Au lieu d'appliquer le noyau directement aux puissances de $A$, on l'applique au résidu de l'approximation inverse.
  • Soit $f(z)$ le noyau approximatif et $E(z) = 1 - (1-z)f(z)$ l'application d'erreur.
  • L'itération est définie par : $Y_{n+1} = Y_n \cdot f(R_n)$ et $R_{n+1} = I - (I-A)Y_{n+1}$.
  • Cette méthode permet de tolérer le débordement tout en garantissant que la précision de la série tronquée croît comme $m^n$.

3. Contributions Clés

1. Premier noyau exact au-delà du radix-5 : Construction d'un noyau radix-9 exact avec $\mu_9=3$ et coefficients rationnels. C'est la première construction explicite de ce type.
2. Noyau radix-15 optimisé : Découverte d'un noyau à 4 produits pour le radix-15 via optimisation numérique, atteignant le coefficient théorique optimal de $6/\log_2 15 \approx 1.54$, bien qu'il soit approximatif.
3. Cadre de noyau radiciel général : Développement d'un algorithme itératif basé sur le résidu qui intègre les noyaux approximatifs (avec débordement) sans perdre l'avantage asymptotique. Ce cadre unifie les méthodes exactes (télescopage) et approximatives.

4. Résultats

Les résultats théoriques et numériques démontrent une réduction significative du nombre de multiplications matricielles :

• Coefficients asymptotiques :
  • Décomposition binaire : $2.00$
  • Quinaire (radix-5) : $1.72$
  • Radix-9 (exact) : $1.58$ (réduction de 21 % par rapport au binaire).
  • Radix-15 (approximatif) : $1.54$ (le meilleur taux connu).
• Expériences numériques :
  • Pour des matrices aléatoires de taille $d=500$, les méthodes radix-9 et radix-15 atteignent la même précision que la méthode binaire avec environ 25 % de produits en moins.
  • Les méthodes exactes (binaire, radix-9) atteignent la précision machine ($\approx 10^{-14}$).
  • La méthode radix-15 converge vers un « plancher d'erreur » d'environ $10^{-6}$, déterminé par la précision du noyau approximatif, ce qui reste acceptable pour de nombreuses applications d'inversion approchée.
  • Les gains de temps d'exécution (wall-clock time) suivent directement la réduction du nombre de produits GEMM.

5. Signification et Impact

Ce travail repousse les limites de l'évaluation des séries de Neumann en démontrant que :

• L'utilisation de radices supérieurs est non seulement possible mais supérieure aux méthodes hybrides {2,3} précédemment considérées comme optimales.
• La transition des noyaux exacts (radix-9) aux noyaux approximatifs (radix-15) ouvre une nouvelle voie où l'on accepte un léger débordement d'erreur pour gagner en efficacité algorithmique.
• Le cadre proposé permet d'appliquer ces noyaux à des problèmes réels sans nécessiter d'extraction exacte des puissances, rendant les algorithmes plus flexibles.

En conclusion, l'article établit un nouveau standard de complexité pour l'inversion de matrices approchée et le préconditionnement polynomial, offrant des gains de performance substantiels pour les calculs haute performance (HPC) sur des matrices denses.