Ce cours est une adaptation des «slides» que j'avais présentés à mon cours sur le calcul analytique aux Journées Nationales de Calcul Formel, en novembre 2011 à Luminy. J'ai rajouté un peu de matériel supplémentaire et quelques explications plus détaillées par rapport à la présentation originale, sans toutefois avoir cherché la perfection d'un cours sur papier habituel.

Pourquoi un cours sur le calcul analytique ? Traditionnellement, on choisit le calcul formel pour la beauté du calcul symbolique! Or, les techniques du calcul formel admettent de nombreuses applications en analyse, et le calcul numérique peut parfois permettre une résolution plus efficace ou directe de problèmes en calcul formel. Toutefois, quand j'ai commencé à m'intéresser à des problèmes plus analytiques, il me manquait la culture nécessaire en analyse calculable et en arithmétique d'intervalles. Ce cours cherche à exposer les techniques les plus utiles de ces domaines, tout en gardant le lien avec le calcul formel en filigrane.

Le document inclut un certain nombre de sessions du système . Il s'agit de la version interprétée mmx-light du langage de novembre 2011. D'ici quelques années, quand l'interpréteur sera remplacé par un compilateur, il pourra être nécessaire de typer plus précisément les déclarations afin de faire fonctionner les exemples.

Le projet vise la conception et le développement d'un nouveau système de calcul formelet analytique. Àl'heure actuelle, deux types de systèmes de calcul mathématique ont connu un grand succès. De manière historique, on trouve les systèmes de calcul numérique, comme , ou . Ces logiciels permettent la résolution approchée de problèmes analytiques, comme l'intégration d'équations différentielles. D'autre part, on dispose des systèmes de calcul formel, comme , , , , et plus récemment . Ces logiciels visent la manipulation exacte d'objets mathématiques, en calculant directement sur des formules générales plutôt que sur des instances numériques particulières.

L'approche du calcul formel est particulièrement bien adaptée pour des objets de nature algébrique ou combinatoire algèbre linéaire, élimination polynomiale, énumération combinatoire, théorie de Galois, etc. En revanche, les systèmes numériques classiques sont beaucoup plus efficaces pour des problèmes plus analytiques, comme l'intégration d'équations différentielles. Malheureusement, les systèmes numériques actuels opèrent principalement en précision fixe et incorporent peu d'outils pour certifier les résultats. La session suivante en illustre bien ce problème :

Il est donc tentant de transposer l'esprit de calcul exact et garanti du calcul formel à l'analyse. Or le calcul formel et le calcul numérique sont encore dans une grande mesure des mondes distincts. D'une part, il y a peu d'intersections quant aux algorithmes de base utilisés (algèbre linéaire, polynômes, séries, etc.). D'autre part, il n'existe pas d'environnement de programmation convivial qui soit adapté aux deux mondes. En effet, les interactions actuelles se limitent surtout àdes interfaces de type « boîte noire ».

Le calcul analytique se situe à la confluence de plusieurs domaines. L'analyse calculable, qui sera discutée dans le chapitre 2, propose une fondation réaliste au regard des architectures existantes d'ordinateurs. Un deuxième point théorique important concerne le développement d'algorithmes numériques certifiés, poursuivant des travaux classiques sur l'arithmétique des intervalles. Dans le chapitre 3, nous exposerons une variante de cette théorie: l'arithmétique de boules.

Évidemment, dans un système de calcul analytique, il est impératif de pouvoir approcher des nombres réels aussi précisément que nécessaire, donc de travailler en précision multiple. Ceci conduit à une pénalité énorme en efficacité par rapport à l'analyse numérique classique, et il convient de réduire cette pénalité autant que possible. Dans le chapitre 4, nous rappelons quelques résultats classiques en théorie de complexité. Dans les derniers chapitres5,6 et7, nous combinons les différentes techniques en les appliquant à l'intégration de systèmes dynamiques et à la résolution polynomiale.

Pendule

Pour , , , calculer (ou ) avec

Systèmes dynamiques plus généraux

Soit et considérons

Pour , et , calculer (si possible) avec

Fonctions spéciales

Étant donné un système zéro-dimensionnel

avec , trouver les zéros “dans ”. Plus précisément et ici encore, il faudrait pouvoir approcher les zéros avec autant de précision que souhaitée.

Variante : prendre , où désigne l'ensemble des “constantes exp-log”, construits à partir de en utilisant les opérations et .

Variante : chercher les racines dans .

avec zéro-dimensionnel en pour et fini.

Étant donné avec et un chemin qui évite , calculer le chemin avec pour tout .

Problème 1 : calculer avec une erreur relative de .

Problème 2 : trouver automatiquement le développement asymptotique

Problème 3 : pour , trouver une constante explicite pour le .

La technique incontestablement la plus importante en calcul analytique est la déformation. Une application astucieuse est la résolution de systèmes polynomiaux par homotopie. Supposons par exemple que l'on veuille résoudre un système

Ce système est suffisamment «générique» pour que toutes ses solutions soient simples et pour qu'il y ait exactement solutions (comme prédit par la borne de Bézout). En particulier, le système «ressemble» (en degré, nombre et nature des solutions) au système beaucoup plus simple

On considère maintenant l'homotopie

qui donne le système en et le système du départ en . Partant des solutions en , on obtient donc les solutions de par déformation, en suivant les solutions quand on fait bouger de vers (voir la section1.2.3 et le chapitre 7).

La technique de suivi de chemin est aussi utile pour le calcul du groupe de Galois d'un polynôme , vu comme polynôme en sur . On procède comme suit :

Pour tout :

On montre que les génèrent .

L'algorithme pour calculer le groupe de Galois d'une fonction algébrique s'adapte au cas des opérateurs différentiels . Considérons par exemple

Base de solutions pour en :

Prolongement analytique de la solution autour de :

La matrice de monodromie joue un rôle analogue à la permutation de l'exemple1.2.

Soit un opérateur différentiel linéaire qui n'admet que des singularités régulières. Alors son groupe de Galois différentiel est généré par ces matrices de monodromie en tant que groupe de algébrique fermé ; c'est le théorème de densité de Schlesinger . Dans le cas général, il faut rajouter les matrices de Stokes et les matrices exponentielles ; c'est le théorème de Martinet-Ramis .

Toutes ces matrices de monodromie peuvent se calculer de façon certifiée et en temps presque linéaire . Ceci reste vrai pour les matrices exponentielles et les matrices de Stokes, grace à la théorie d'accéléro-sommation d'Écalle. Pour une précision suffisamment grande de calcul, ceci conduit a un algorithme pour calculer le groupe de Galois différentiel. En particulier, on obtient un algorithme de factorisation de .

Prétraitement symbolique meilleure efficacité

Plus généralement: asymptotique formelle évaluation rapide de fonctions

De façon similaire, la résolution de systèmes algébriques surdéterminés peut bénéficier d'un prétraitement pour calculer une base de Gröbner ; ensuite on peut par exemple chercher les solutions numériques par méthode d'homotopie.

Calcul de schémas numériques rapides

Arithmétique rapide

Besoin d'un bon langage

Comme l'ensemble des nombres réels est non dénombrable, il est impossible de construire un type de données susceptible de pouvoir représenter n'importe quel nombre réel. Toutefois, au moins d'un point de vue théorique, il serait utile de pouvoir calculer directement avec des nombres réels, quitte à se restreindre à une sous-classe de . Le choix le plus naturel consiste à prendre la sous-classe des nombres que l'on peut approcher aussi précisément que nécessaire par des nombres rationnels:

Ici, il est important d'intégrer le fait qu'un nombre est en réalité un programme ou encore la promesse de pouvoir trouver une approximation aussi fine que souhaitée. Bien sûr, la définition admet de nombreuses variantes, en jouant sur la façon d'approcher . Par exemple :

L'analyse calculable est le sujet qui se propose de «réécrire» l'analyse sous cet angle de la calculabilité. La encore, il y a des sujets voisins, comme l'analyse constructive.

Une autre approche est de faire «comme si» l'on pouvait représenter n'importe quel nombre réel et effectuer des opérations habituelles en temps constant. Évidemment, ce modèle ne correspond encore moins à des ordinateurs concrets que les machines de Turing. Toutefois, ce modèle simplificateur n'est pas dénué de sens si on calcule en précision fixe et il n'est pas inutile d'adopter ce point de vue par moments, en particulier pour la résolution de systèmes polynomiaux (voir les chapitres 6 et 7).

Le premier théorème de non calculabilité parait naturel aujourd'hui. Or c'est ce théorème qui est à l'origine des machines de Turing et des premiers résultats de non calculabilité !

Asymmétrie. Soit .

Restriction à des sous-classes.

Le deuxième théorème de non calculabilité est un peu plus surprenant, du moins pour moi-même: ce n'est pas la peine de chercher à calculer des fonctions discontinues, puisque c'est impossible !

Le troisième théorème de non calculabilité est surprenant dans ce sens que les données ne sont pas des nombres ou fonctions calculables, mais plutôt des simples vecteurs et polynômes à coefficients dans . Le théorème pose clairement des limites à ce que l'on pourra calculer dans la théorie des systèmes dynamiques.

Comme exercice, le lecteur pourra essayer de déterminer lesquels des problèmes suivants sont calculables :

Les théorèmes de non calculabilité de la section 2.2 sont inutilement pessimistes. Par exemple, le théorème de Turing ne rend pas compte du fait suivant: étant donné un nombre , on peut démontrer en temps fini que , si tel est en effet le cas. Pour mieux cerner ce qui reste calculable, il est utile d'introduire les ensembles de nombres réels calculables à gauche et à droite, et qui reflètent mieux l'asymmétrie profonde de l'analyse calculable entre l'égalité et l'inégalité.

Calculabilité à gauche

s'il existe calculable avec et

Calculabilité à droite

si . On a

On peut introduire une ribambelle de notions de calculabilité en poursuivant sur cette voie. Par exemple :

Inférieurement calculable

s'il existe calculable avec , où chaque est une réunion finie de « blocks » fermés à extrémités dans

Supérieurement calculable

si

Il est utile de considérer des ensembles comme , , , , etc. comme des types de données concrets, pouvant intervenir dans des implantations. On peut d'ailleurs construire plein d'autres types du même acabit, comme les booléens à droite ou les «booléens modulo la conjecture de Schanuel» .

Ces types de données plus fins permettent de transformer les énoncés quelque peu négatives de la section2.2 en assertions plus positives :

D'un point de vue logique, les types de données peuvent être vu comme des «structures effectives» sur un ensemble . Chaque élément de admet alors une représentation dans un ensemble de représentations. Chaque représentation peut s'encoder à son tour par un entier dans , et donc être manipulée par des machines de Turing.

Par exemple, un nombre dans pourrait être représenté par une suite calculable croissante :

Il est à noter que, de même que l'on peut définir plusieurs structures de groupe sur un ensemble à quatre éléments, on peut définir plusieurs structures effectives sur : on a bien sûr , mais aussi et . La représentation d'un élément de pourrait être une fonction qui retourne un élément de si la conjecture de Schanuel est vraie.

En analyse, les objets que l'on manipule correspondent généralement à des résultats de processus d'approximation. Il s'agit donc de structures effectives d'une nature particulière :

Le plus souvent, est constitué de parties de (comme des intervalles ou des boules dans le cas où ) et la notion de limite correspond à prendre une intersection. Par exemple :

Toutefois, il faut parfois tordre un peu la notion de limite :

L'intérêt de cette formalisation est une certaine fonctorialité. Par exemple, si on peut approcher les éléments de et de , il en est de même pour les éléments de et les fonctions dans de vers :

si

Dans la définition2.1, nous avions choisi d'approcher les nombres réels par des rationnels. Malheureusement, de telles approximations ne sont pas des plus judicieuses si nous voulons borner de façon automatique les erreurs intervenant dans un calcul complexe.

Pour cette raison (voir aussi la section2.3.4), il est plus commode d'approcher un nombre réel donné non pas par d'autres nombres avec «petit», mais plutôt par des boules fermées de la forme

avec et . Ainsi, lorsque est le résultat d'un calcul par arithmétique de boules, on est sûr que le «vrai résultat» est dans . On peut donc interpréter une boule comme un «nombre que l'on ne connait pas, mais pour lequel on est sûr que».

Bien sûr, on peut prendre des boules dans des espaces métriques plus généraux que , en commençant par . En fait, il n'est même pas impératif que les rayons des boules vivent dans un espace métrique au sens classique. Plus tard, on verra par exemple l'utilité de considerer des boules dont les centres sont des matrices dans et dont les rayons sont également des matrices àcoefficients positifs dans .

Dans la suite, étant donnés un ensemble de centres et un ensemble de rayons, on notera par l'ensemble des boules avec des centres dans et des rayons dans .

Soit une fonction. Une fonction est appelée une extension de si

pour tous les . Ainsi, si , alors on peut être sûr que , conformément à la sémantique décrite dans la section précédente.

On peut utiliser les formules suivantes pour les opérations élémentaires :

Il est à noter que la dernière formule est simple, mais pas nécessairement optimale : par exemple, on trouve , alors que la boule vérifie bien .

Pour des calculs concrets sur machine, on ne peut prendre les centres et les rayons dans. Dans ce cas, on considère plutôt des centres et des rayons dans , voire dans l'ensemble

des nombres flottants avec un mantisse de bits et un exposant de bits. Dans ce dernier cas, il est à noter qu'il est impératif de rajouter la possibilité de prendre des rayons infinis, afin de représenter le résultat d'un overflow (dépassement de précision).

Quand on calcule avec en précisions et fixes, le résultat exact d'une opération sur des nombres dans n'est pas nécessairement dans et son approximation par un élément dans induit donc une erreur dont il faut tenir compte dans l'implantation d'une arithmétique de boules.

La norme IEEE pour le calcul avec des nombres flottants spécifie que l'approximation est obtenue en arrondissant le résultat exact suivant un mode d'arrondi à spécifier. Nous utiliserons les notations pour les modes d'arrondi vers le haut, vers le bas et au plusprès. En notant , nous notons que est une borne supérieure pour l'erreur quelque soit le mode d'arrondi choisi. Les formules(3.2),(3.3) et(3.4) peuvent donc adaptées au cas du calcul en précision limitée, en les remplaçant par

Il est à noter que la norme IEEE est désormais suivie par la plupart des constructeurs de micro-processeurs. fournit aussi une bonne implantation en précision arbitraire.

On a déjà mentionné le fait que l'arithmétique de boules est classiquement considérée comme une variante de l'arithmétique d'intervalles. Pourquoi choisir l'une ou l'autre? Voici quelques avantages et inconvénients des deux approches.

Une question intéressante concerne la sémantique des prédicats comme . Si on interprète les boules comme un ensemble de «valeurs possibles», alors le résultat d'un test devrait être un intervalle de avec

En revanche, si on interprète et comme des approximation d'un vrai nombre réel à une certaine précision, il est plus naturel de prendre

En effet, si sont des nombres calculables représentés par des suites calculables décroissantes de boules avec et , alors la suite à valeurs dans représente le résultat du test d'égalité dans . La deuxième définition reflète donc l'asymmétrie profonde de l'égalité dans l'analyse calculable.

L'arithmétique des boules est basée sur une analyse a posteriori des erreurs : on effectue une opération en arithmétique de boules en utilisant une certaine précision que l'on augmentera jusqu'au moment où le rayon du résultat est suffisamment petit. Parfois, bien que rarement, on peut préférer une estimation a priori de l'erreur.

Exemple : calcul d'une -approximation de

Calculer des -approximations et de et

Retourner

Avantage : précision de calcul adaptatif dans le case

En effet, dans l'arithmétique de boules classique, on calcule et en utilisant la même précision, alors que l'on peut simplement utiliser comme approximation de tant que . Supposant que est un nombre dont l'approximation est très couteuse en temps, une estimation a priori de l'erreur est donc plus efficace.

Problème : évaluation de par Horner en

Calcul d'une -approximation de pour

Problème : nécessite des « dags » (directed acyclic graphs) pour tous les résultats intermédiaires

Le problème principal de l'arithmétique de boules concerne la surestimation des erreurs. Par exemple, l'évaluation de la fonction en donne

Plus généralement, ce problème intervient chaque fois qu'il y a des quantités qui s'annulent lors d'un calcul, sans que la méthode de calcul s'en rend compte.

Dans certain cas, la surestimation est due à notre façon de représenter les encadrements. Par exemple, en arithmétique d'intervalles standard, on encadre les nombres complexes par des rectangles de la forme . Ces rectangles ont tendance à se tourner lors de la multiplication par un nombre complexe non réel, comme dans l'exemple

Ici l'inclusion du résultat dans un autre rectangle de la forme requise implique une perte automatique de précision. On appelle ceci l'effet d'enveloppement. Sur cet exemple précis, l'inconvénient disparaît lorsque l'on utilise l'arithmétique de boules, puisqu'une boule reste une boule si on la tourne. Toutefois, on verra un exemple plus complexe dans la section5.5.2, où il faudra travailler plus pour résoudre le problème.

Illustrons la perte de précision lorsque l'on calcule pour de façon naïve en utilisant l'arithmétique d'intervalles. Dans la sortie, on utilise la notation scientifique. Par exemple désigne un nombre compris entre et .

Comme on l'avait prédit, la perte de précision n'apparait pas lorsque l'on remplace l'arithmétique d'intervalles par l'arithmétique de boules.

Une autre technique pour limiter l'effet d'enveloppement consiste à utiliser des algorithmes qui minimisent la profondeur du calcul. En effet, dans la section 3.2.2, l'effet est amplifié par le fait que l'on réinjecte fois le résultat du calcul précédent dans l'étape d'après, perdant ainsi bits de précision. Si on utilise un algorithme diviser pour régner pour calculer des puissances, la perte de précision se limite à bits. À noter que ce genre d'algorithmes sont intéressants de toute façon, car ils se parallèlisent généralement mieux et ils se comportent mieux vis à vis de la mémoire cache.

Nous avons déjà souligné l'importance des méthodes perturbatives en calcul analytique. Elles interviennent notamment lorsque l'on cherche à résoudre des équations et que l'on a déjà un moyen pour résoudre l'équation de façon approchée.

Considérons par exemple l'inversion d'une matrice , qui correspond à la résolution de l'équation . Si est une matrice de boules, réécrite comme une boule de centre et de rayon , le calcul de par pivot de Gauss naïf produit une surestimation énorme. Or il existe de bons algorithmes numériques pour inverser le centre de façon approchée, produisant une matrice avec . Dès lors, il suffit d'étudier de combien l'inverse de peut bouger lorsque l'on fait varier . Ceci conduit à l'algorithme suivant :

Bien sûr, la méthode perturbative s'applique à d'autres problèmes similaires, comme le calcul de la forme échelon. Montrons d'abord le calcul quand on utilise l'algorithme naïf qui consiste à directement appliquer le pivot de Gauss sur une matrice dont les coefficients sont des intervalles :

On observe bien une déperdition de la précision au cours de l'algorithme. Lorsque les coefficients de la matrice sont remplacés par des boules, utilise la méthode perturbative :

En calcul numérique, la précision relative du résultat est généralement moindre que la précision relative des données. On pourrait appeler la différence entre ces deux précisions la «déperdition de précision». Cette déperdition admet deux sources tout à fait distinctes:

Fixons une norme sur et soit une fonction. On définit le nombre de conditionnement de en par

Le logarithme en base deux mesure donc la déperdition intrinsèque de précision pour l'évaluation de en . En algèbre linéaire, on définit aussi le nombre de conditionnement d'une matrice par

Ceci correspond au maximum des pour , où .

Supposons maintenant que l'on cherche à approcher par , en calculant avec une précision de bits. On définit le facteur d'éloignement de en par

Quand on calcule en on perd donc environ bits de précision de façon intrinsèque et environ bits additionnels à cause de l'algorithme employé.

Essayons maintenant de transposer l'esprit de la section précédente à l'arithmétique de boules. Le nombre de conditionnement n'a pas de contre-partie directe. En revanche, on la notion d'une «extension optimale» :

Soit une fonction. Nous avons déjà signalé l'existence a priori d'une multitude d'extensions de vérifiant la relation

pour tous les . Néanmoins, si on impose la contrainte supplémentaire que est de la forme pour un certain , alors il existe une unique extension optimale , définie par

Ici, on utilise la notation vectorielle. Par exemple .

Maintenant, nous pouvons définir la surestimation d'une extension arbitraire de façon naturelle par rapport à l'extension optimale . Plus précisément, étant donnée une boule , on définit la surestimation de en par

On définit également la surestimation ponctuelle de en par

Nous allons voir maintenant que cette surestimation ponctuelle est souvent simple à calculer. De la même manière qu'il est bon pour un algorithme numérique de surveiller le nombre de conditionnement et le facteur d'éloignement, c'est généralement une bonne idée pour un algorithme certifié de surveiller la surestimation. En effet, si la surestimation devient trop importante, il est souvent possible de déclencher un autre algorithme a priori plus couteux, mais plus précis.

Supposons que l'on ait une expression construite à partir de constantes dans , d'un nombre fini d'indéterminées et les opérations . D'une part, on peut interpréter comme une fonction . D'autre part, en utilisant l'arithmétique de boules standard(3.23.4), l'expression donne naturellement lieu à une fonction . Il n'est pas difficile de montrer que la surestimation ponctuelle de se calcule explicitement par

où l'opérateur de «gradient majoré» est défini par

Lorsque , la formule(3.6) se simplifie en

Toujours pour l'arithmétique de boules standard, une problème intéressant est de pouvoir calculer ou au moins estimer la surestimation sur une boule générale par rapport à la surestimation ponctuelle. Voici une question plus précise qui semble abordable:

Ayant défini la surestimation d'un algorithme de façon précise, la question est maintenant comment construire des algorithmes à la fois efficaces et de bonne qualité (c'est à dire, avec une surestimation proche de ). Bien évidemment, il s'agit de deux objectifs contraires, donc il faudra chercher un compromis.

Par exemple, en rendant systématiquement la boule comme résultat, on obtient un algorithme très efficace, mais sans intérêt. Réciproquement, dans de nombreux cas, l'obtention de bornes optimales est NP complet ou pire . Généralement, les algorithmes qui nous intéressent introduisent un facteur constant dans la complexité pour calculer des bornes de qualité acceptable, ou des petits facteurs logarithmiques ou polynomiaux pour avoir un algorithme avec une surestimation proche de .

On rappelle aussi qu'il est souvent possible d'appliquer la stratégie de la «terminaison précoce» : on commence avec un algorithme rapide et a priori de faible qualité. Seulement si la qualité du résultat est jugée insuffisante, on lance des algorithmes plus lents, mais de meilleure qualité.

La recherche d'un compromis entre efficacité et qualité s'illustre bien sur le problème de la multiplication de deux matrices .

Stratégie naïve

On calcule le produit par l'algorithme classique en utilisant opérations en arithmétique de boules. C'est un algorithme de bonne qualité, mais assez lent.

Emploi de boules matricielles

On réinterprète comme des boules avec des centres et des rayons dans . Dès lors, on peut utiliser la formule(3.4) pour la multiplication, quitte à rajouter le nécessaire pour les erreurs d'arrondi.

L'avantage principal de cette méthode est que l'on s'appuie directement sur la multiplication matricielle dans au lieu de faire toutes les opérations sur les coefficients un par un dans l'arithmétique des boules. En effet, la multiplication dans est généralement hautement optimisée (en précision machine, on peut par exemple utiliser des ), donc on gagne au moins une grosse constante par rapport à la stratégie naïve.

Majorations brutales

Dans la méthode précédente, on majore le rayon du produit par la somme de deux produits de matrices dans . Étant données deux matrices , il est souvent possible d'obtenir une majoration raisonnable du produit en faisant seulement opérations.

En effet, supposons qu'après préconditionnement, on peut s'arranger pour que les entrées de chaque ligne de et de chaque colonne de soient du même ordre de grandeur. Dans ce cas, on peut utiliser la formule

Autres compromis

Dans la pratique, il arrive souvent que les entrées d'une matrice soient localement du même ordre de grandeur, mais pas globalement, même après préconditionnement. Dans ce cas, on peut découper la matrice en petits blocs de matrices , utiliser des majorations brutales sur ces petits blocs et un algorithme plus naïf pour la matrice toute entière. Ceci conduit à un algorithme de coût supérieur à la multiplication dans , avec .

De manière générale, on observe aussi que le coût de certification des calculs tend souvent à devenir négligeable pour des gros problèmes. Par exemple, pour la multiplication dans , l'estimation des erreurs se fait généralement en simple précision. En grande précision, le coût de deux multiplications en simple précision devient négligeable devant le coût d'une grosse multiplication dans . De même, dans des cas bien conditionnés où on peut utiliser la méthode des majorations brutales pour la multiplication des matrices, le coût des opérations pour estimer les erreurs est négligeable devant le coût pour la multiplication des centres.

À présent, nous avons vu les types de base pour le calcul analytique. D'un point de vue haut niveau, nous avons les nombres calculables de, , , etc. Juste en-dessous suit l'arithmétique de boules qui permet le calcul certifié avec des approximations. Tout en bas, nous avons le calcul numérique classique avec des nombres en virgule flottante et l'arithmétique rapide pour les mantisses en précision arbitraire.

Cette division en quatre niveaux de la «hiérarchie numérique» est pertinente pour la plupart des algorithmes en calcul analytique. Dans le cas d'une multiplication de matrices par exemple, on procède comme suit :

L'algorithmique rapide est un sujet classique en calcul formel. Nous ferons quelques rappels, sans chercher à être exhaustif. Voir par exemple pour deux références classiques.

L'opération clef à comprendre pour des analyses en complexité est la multiplication, que ça soit sur les entiers, sur les nombres flottants, les polynômes, les matrices, ou encore les séries. Les complexités d'autres opérations (division, racine carrée, etc.) s'expriment généralement en fonction du coût de la multiplication. Voici quelques complexités classiques concernant la multiplication:

Il existe deux variantes pour multiplier deux polynômes de degrés . On peut couper les polynômes en deux à l'exposant , ou par exposants pairs-impairs. L'algorithme pair-impair est généralement plus stable d'un point de vue numérique:

La multiplication de Karatsuba est un algorithme «multi-modulaire» classique, procédant par évaluation-interpolation dans les points (en projectif). Si admet une racine primitive de l'unité d'ordre (donc ), une stratégie multi-modulaire plus efficace est basée sur la transformation de Fourier discrète rapide (FFT). Pour un polynôme avec , on définit par

et il est classique que l'on peut calculer et son inverse en temps . Étant donné deux polynômes avec on peut donc calculer leur produit en temps par la formule

Un point de vue naturel pour calculer avec des séries dans est de calculer systématiquement avec des séries à l'ordre . Pour multiplier deux telles séries, il suffit de multiplier les séries tronquées en tant que polynômes et de tronquer le résultat. Toutefois, si on utilise la multiplication naïve, ceci conduit à faire environ deux fois trop de travail. C'est un problème ouvert de savoir si on peut faire mieux en général:

On a déjà signalé le fait que la complexité de la plupart des opérations plus complexes, comme la division, la racine carrée, l'exponentielle, etc. s'expriment en fonction de la complexité de la multiplication. Une technique puissante pour démontrer ceci est la méthode de Newton.

Supposons par exemple que l'on veuille inverser une série avec à l'ordre. En d'autres mots, étant donnés les coefficients , on voudrait calculer les coefficients de . La méthode de Newton appliquée à l'équation donne l'itération

Comme on ne connait pas encore dans cette itération, on utilisera plutôt l'itération

Si est bon jusqu'à l'ordre :

la nouvelle valeur sera bonne jusqu'à l'ordre , puisque

Si désigne le coût de la division à l'ordre , ceci donne

En supposant que est croissante, on laisse au soin du lecteur de vérifier que ceci implique . On déduit aussi aisément que le calcul du quotient et du reste de la division euclidienne d'un polynôme de degré par un polynôme de degré se calcule en temps .

L'utilisation de la méthode de Newton dans ce cadre remonte à Hensel et se généralise pour des équations polynomiales plus générale (à la fois dans les séries et les nombres -adiques d'ailleurs). Brent et Kung ont remarqué que la même technique peut être utilisée pour composer des séries formelles, calculer leurs inverses fonctionnelles, ou résoudre des équations différentielles. Par exemple, l'exponentielle de avec se calcule par l'itération

Une méthode efficace pour la résolution d'équations différentielles algébriques plus générales a été proposée par Sedoglavic.

Il y a de nombreux autres opérations sur les polynômes (racine carrée, pgcd, ppcm, etc.) qui peuvent se calculer en temps presque linéaire. Nous donnons un dernier exemple qui est utile pour le calcul des zéros d'un polynôme. En effet, si on dispose déjà de bonnes approximations pour les racines, l'évaluation multi-point rapide permet d'appliquer la méthode de Newton pour améliorer toutes ces approximations de façon simultanée.

On remarque que la complexité du précalcul et de la dichotomie proprement dite vérifient une estimation du type

En supposant que est croissante, on laisse au soin du lecteur de vérifier que ceci implique . On peut encore gagner un facteur constant sur l'algorithme présenté ici.

Nous avons vu comment calculer avec des séries en les tronquant à un certain ordre . Une autre approche consiste à considérer les séries comme des «flots de coefficients» qui sont calculés un par un. Ce point de vue détendu impose une contrainte sur notre façon de concevoir les opérations sur les séries. Par exemple, la multiplication détendue de deux séries impose la sortie de dès que sont connus. Par conséquent, on ne peut pas directement appliquer certains algorithmes rapides, comme la multiplication FFT. L'avantage du calcul détendu est qu'il permet naturellement de résoudre des équations «récursives»

où l'extraction du coefficient en de ne fait intervenir que les coefficients de . Voici deux exemples:

La stratégie la plus naïve pour calculer le coefficient d'un produit est d'utiliser la formule classique de convolution

Pour un calcul jusqu'à l'ordre , cela donne une complexité . Dans la figure4.1, on montre les étapes successives du calcul détendu de l'exponentielle .

Lorsque l'on applique l'algorithme de Karatsuba sur des polynômes et avec , où on considère comme des paramètres formels, on observe que la formule pour ne dépend que de pour chaque . En d'autres termes, l'algorithme de Karatsuba est naturellement détendu, si on fait les calculs dans le bon ordre. Ceci montre qu'il existe un algorithme de multiplication détendue de complexité .

On peut encore faire mieux en anticipant des calculs à venir. Par exemple lorsque et sont connus, on peut calculer la contribution de au produit par n'importe quel algorithme de multiplication rapide pour les polynômes. De même, lorsque l'on connait et , on peut calculer la contribution de au produit de façon rapide. En exploitant cette idée, les contributions de tous les grands carrés dans la figure4.2 peuvent se calculer par un algorithme rapide. Il s'ensuit que

Est-ce que l'on peut faire encore mieux ? Si admet suffisamment de racines -ièmes de l'unité, on montre que:

Dans le cas général, la question reste ouverte:

La multiplication FFT de la section4.2.1 est un exemple classique d'algorithme multi-modulaire. En effet, elle repose en réalité sur l'isomorphisme

entre les polynômes modulo et les réductions de modulo les avec . De la même manière, étant donnés nombres premiers mutuellement distincts, le théorème des restes chinois nous fournit un isomorphisme

Dans chaque cas, le calcul multi-modulaire repose sur un changement de représentation, où un objet dans ou où la multiplication coûte cher est remplacé par un vecteur d'objets pour lesquels la multiplication coûte moins cher.

Comme les changements de représentation coûtent généralement cher aussi, il est recommandé d'effectuer un maximum de calculs dans le modèle multi-modulaire. Supposons par exemple que l'on veuille multiplier deux matrices et dans . Au lieu de faire opérations dans , de coût , il vaut mieux utiliser la formule

qui correspond à mettre tous les coefficients de et dans le modèle FFT, de multiplier matrices scalaires dans , et de faire la transformation inverse. En effet, cet algorithme admet une complexité bien meilleure.

Il y a plusieurs variantes pour les méthodes multi-modulaires. Chaque variante admet ses caractéristiques propres quant au nombre de moduli utilisés et quant au coût de la réduction et de la reconstruction. Les méthodes les plus utilisées sont les suivantes:

Restes chinois

Cette méthode permet grosso modo d'encoder un entier d'environ mots machines en utilisant seulement moduli. En revanche, les coûts de la réduction et de la reconstruction sont en pour un entier de taille .

FFT sur avec beaucoup de racines -ièmes de l'unité

Certains nombres premiers comme admettent des racines primitives de l'unité d'un ordre avec grand. Cela permet d'utiliser la multiplication FFT pour des polynômes de degrés importants à coefficients dans . Par ailleurs, tout nombre dans peut se coder par avec . En prenant , on peut aussi ramener la multiplication dans à la multiplication dans pour des nombres pas trop grands. Par rapport aux restes chinois, ceci revient à prendre plus de moduli, mais aussi à gagner sur le coûts de réduction et de reconstruction. On peut d'ailleurs combiner les deux approches.

FFT numérique

Lorsque l'on calcule avec des nombres flottants complexes, on peut directement calculer avec une racine -ième de l'unité approchée. En encadrant les erreurs d'arrondi avec soin, ceci conduit à des algorithmes de multiplication rapide dans et qui ont des avantages et inconvénients semblables par rapport à la méthode précédente.

FFT synthétique de Schönhage-Strassen

En cas d'absence de racines -ièmes de l'unité pour suffisamment grand, on peut les synthétiser. Supposons par exemple, que l'on veuille multiplier des polynômes dans et prenons . Alors est une racine -ième primitive de l'unité dans et montre que la multiplication dans se ramène à une multiplication de polynômes à coefficients dans . En outre, la FFT à coefficients dans est très efficace, car elle ne nécessite que des additions et des soustractions dans . Par ailleurs, la méthode est parfaitement générale, et marche encore pour des grands entiers. Le seul point noir de l'approche est qu'elle nécessite plus de moduli que les méthodes précédentes.

Un problème qui revient fréquemment est l'évaluation rapide d'un vecteur de polynômes en plusieurs variables dans une -algèbre . En effet, ce problème intervient dans l'intégration de système dynamiques ou encore dans des méthodes algébriques par homotopie. En outre, il n'est pas rare que l'on veuille évaluer à la fois et sa matrice jacobienne .

Nous allons nous concentrer sur le cas , en supposant que la multiplication dans s'effectue par un algorithme multi-modulaire avec moduli et avec un coût de réduction-reconstruction de opérations dans .

Considérons d'abord le cas où les polynômes sont donnés par des expressions avec éventuellement des sous-expressions communes. En d'autres termes, est représenté par un «dag» (directed acyclic graph). La figure4.3 montre un dag typique, correspondant au vecteur avec

Un sous-dag de est dit scalaire s'il ne fait pas intervenir les . Dans la figure4.3, les seuls sous-dags scalaires sont et . Un sous-dag multiplicatif de de la forme est dit

Dans notre exemple, les ensembles de sous-dags multiplicatifs homothétiques, terminaux et actifs sont respectivement , et . On a:

Si les polynômes sont donnés comme combinaisons linéaires de monômes dans , alors il faut construire un dag pour avant d'appliquer la proposition4.1. Ceci conduit au problème comment trouver un dag pour lequel le nombre de sous-dags multiplicatifs actifs est minimal. Bien sûr, ce problème de nature plutôt combinatoire est assez délicat, donc on cherche plutôt de bonnes stratégies pour s'approcher du minimum.

Soit l'ensemble des monômes dans apparaissant dans . La question se réduit essentiellement à la nouvelle question comment trouver deux sous-ensembles tels que et tels que le cardinal soit minimal. Deux approches heuristiques semblent raisonnables:

Pour chaque partie de fonction caractéristique , on a une projection naturelle

La projection peut se calculer en temps linéaire. L'idée est maintenant de prendre et pour un qui minimise . Si est grand, ce calcul peut encore devenir trop cher, auquel cas on peut se restreindre aux parties de la forme avec .

Étant donnés des entiers , on peut considérer les projections

avec . L'idée est alors de prendre et pour des entiers bien choisis. On peut trouver ces entiers en commençant avec et en doublant l'un des entiers tant que cela fait chuter la cardinalité de .

Malheureusement, les algorithmes asymptotiquement rapides pour multiplier des polynômes à coefficients flottants sont numériquement instables dans le cas général. Par exemple, considérons le calcul suivant d'un carré par la méthode de Karatsuba, en utilisant une précision de chiffres binaires:

On voit bien que le fait d'avoir ajouté et soustrait des coefficients d'ordres de grandeur différents a pollué le calcul, de manière à rendre le coefficient en du résultat incorrect.

Une manière simple pour régler le problème de la section précédente est de rendre tous les coefficients de de même ordre de grandeur modulo une postcomposition par pour une«échelle» bien choisie:

Malheureusement, cela ne règle pas le problème dans le cas général, car le choix d'un bon peut être ambigu. Par exemple, si

il faudrait prendre pour que et soient du même ordre de grandeur, mais pour que et le soient.

Heureusement, dans le cas d'une série , le comportent des est généralement assez régulière pour , car dicté par la nature de la singularité de qui est la plus proche de l'origine. Dans de nombreux cas (et notamment quand est algébrique avec une unique singularité dominante), on a par exemple

est le rayon de convergence de . Dès lors, l'existence d'une bonne échelle est garantie pour tout polynôme avec et suffisamment grands. D'un point de vue pratique, s'approche bien, en considérant la pente au milieu du polygone numérique de Newton (voir la section6.3).

Afin de vérifier la stabilité numérique de la méthode de préconditionnement par mise à l'échelle, il suffit d'effectuer le même calcul avec des précisions différentes. Bien entendu, cela ne garantit rien, mais ça indique que tout se passe comme prévu.

Pour une séries à coefficients dans , obtenue comme le résultat d'un calcul naïf comme

il est instructif d'étudier le comportement des rayons des par rapport àcelui des centres . La différence principale entre les calculs de et est que «tous les ont été remplacés par des » dans les formules de récurrence. En effet, dans notre exemple, on a

Par conséquent, si on calcule avec une précision de chiffres binaires, on obtient

pour un certain . Comme le rayon de convergence de la série est plus petite que le rayon de convergence de , il s'en suit qu'il existe une constante telle que l'algorithme d'inversion naïf perd toute sa précision pour .

Le remède est pareil que pour l'inversion des matrices dans la section3.2.5: supposons que l'on veuille inverser une série . Alors on calcule un inverse numérique approché de , et on prend , où est calculé utilisant une méthode détendue générique.

Considérons maintenant le cas du calcul certifié d'une exponentielle . Si on fait une analyse similaire à celle de la section précédente, on se rend compte que les séries génératrices et formées par les centres et rayons des coefficients vérifient cette fois-ci

est la série génératrice avec . Or les rayons de convergence de et sont les mêmes et les singularités dominantes généralement de même nature. Par magie, on observe donc rarement de la surestimation dans ce cas! Cette observation s'étend d'ailleurs aux intégrales de systèmes dynamiques plus générales.

Ceci dit, dans les rares cas où on observe quand même de la surestimation, on peut bien sûr adapter la méthode perturbative. En effet, pour le calcul de avec , il suffit de prendre , après quoi la fonction se calcule en intégrant le système dynamique .

Pour l'instant, on a surtout vu des boules scalaires dans ou et des boules qui opèrent coefficient par coefficient, comme des boules dans ou . Il est maintenant temps d'étendre les principes du calcul à des espaces un peu plus sophistiqués comme l'ensemble des fonctions réelles analytiques sur la boule. La norme naturelle sur cet espace est donnée par

Pour et , on peut donc considérer la boule fermée

Étant donné un ordre , on peut ensuite prendre les centres dans l'ensemble

Une boule avec et est appelé un modèle de Taylor d'ordre sur la boule . On note par l'ensemble des modèles de Taylor de ce type. Pour des calculs sur machine, on peut évidemment définir la variante .

La multiplication doit être définie avec un peu plus de soin pour les modèles de Taylor:

Si on calcule avec des flottants de précision , il faut en plus prendre en compte les erreurs d'arrondi, par exemple en rajoutant

au rayon en utilisant systématiquement le mode d'arrondi vers le haut pour les estimations d'erreur. Dans tous les cas, on observe que le calcul des erreurs est négligeable devant le coût de la multiplication polynomiale tronquée pour pas trop petit.

Contrairement au boules traditionnelles, les modèles de Taylor permettent aussi l'implantation d'opérations fonctionnelles comme l'intégration:

Ici, on utilise le fait que pour tout .

Les définitions se généralisent naturellement aux cas des fonctions analytiques complexes et des fonctions analytiques en variables sur un polydisque avec . Les centres d'un modèle de Taylor en plusieurs variables vivent généralement dans un ensemble du type

avec , et .

Une application importante des modèles concerne la réduction de la surestimation quand on évalue une fonction en arithmétique de boules. Revenons à l'exemple3.3:

Considérons maintenant l'évaluation de l'expression en un modèle de Taylor d'ordre (ou plus):

Pour tout on a donc , ce qui nous permet de définir une meilleure extension de par

Il est facile de vérifier que cette extension est ponctuellement optimale dans le sens que

pour tout . Bien sûr, il faudra recourir à des modèles de Taylor d'ordre au moins si présente des racines de multiplicité (ou lorsqu'une petite perturbation de présente de telles racines).

Considérons un système de polynômes réels

Un algorithme classique pour trouver les solutions de ce système dans une «boîte» est de la découper en des boîtes de plus en plus petites jusqu'au moment où l'une des deux situations suivantes se présente:

Si n'admet que des racines simples dans , alors cet algorithme simpliste termine. Toutefois, en présence de racines multiples, ou racines proches (éventuellement dans des voisinages imaginaires), l'évaluation de présentera une surestimation qui rend la méthode inutilisable.

En utilisant les modèles de Taylor, on aimerait donc réduire cette surestimation. Le principal problème est alors de trouver un ordre de troncature adéquat. Ceci peut se faire de façon heuristique. Afin de simplifier l'exposition, considérons le cas de dimension un. Supposons que l'on a calculé un encadrement de sur en utilisant des modèles de Taylor à l'ordre. Pour un coût plus grand on peut aussi calculer un encadrement en utilisant des modèles de Taylor à l'ordre . Même si on ne connait pas les surestimations et de et dans l'absolu, le quotient entre les rayons de et de nous fournit le ratio . On s'attend donc à ce que l'on doive découper la boîte en des boîtes fois plus petites en utilisant des modèle de Taylor d'ordre plutôt que . Lorsque , on a donc intérêt à prendre des modèles de Taylor d'ordre plutôt que des modèles d'ordre . Ce critère nous permet de trouver un ordre à peu près optimal.

Nous notons en outre que les modèles de Taylor en plusieurs variables nous permettent en principe des prendre des boîtes qui ne sont pas forcément parallèles aux axes (voir aussi la section7.2.3) et donc de mieux suivre la géométrie du problème. En fait, les boîtes ont même le droit de devenir courbées, si cela peut arranger les estimations d'erreur.

L'intégration certifiée de systèmes dynamiques utilisant la technique des modèles de Taylor a été proposé par Berz et Makino. Un grand avantage de ces méthodes est qu'ils permettent d'étudier la dépendance du flot par rapport aux conditions initiales pour des conditions initiales dans des domaines relativement grands. Bien sûr, ceci nécessite de recourir à des modèles de Taylor en plusieurs variables, c'est à dire par rapport au temps et par rapport aux conditions initiales. Dans cette section, nous allons plutôt nous intéresser à des techniques qui continuent à marcher en haute précision. Dans ce cadre, les conditions initiales sont forcément spécifiées de façon très précise, donc on peut se contenter de modèles de Taylor en une variable. Dans la section5.5.2, on aura juste besoin d'étudier la dépendance du flot par rapport aux conditions initiales à l'ordre un.

Supposons que l'on veuille intégrer le système dynamique

avec , , et . Plus précisément, supposons que l'on veuille calculer la valeur de en . Alors on peut utiliser la méthode suivante:

Vérifions que cette méthode est correcte. En effet, si on a l'inclusion(5.1), alors l'opérateur de envoie la boule dans lui même. Comme est une boule compacte dans un espace de Banach, il s'en suit que admet un point fixe. Or il existe une unique façon de calculer les coefficients en série d'un tel point fixe, donc ce point fixe est nécessairement unique.

Il nous reste à montrer comment calculer la constante . L'idée est tout simplement de commencer avec et d'itérer quelque fois le processus suivant: on calcule et on remplace par le supremum de et le rayon du résultat. Si cette itération diverge numériquement, alors on arrête et (5.1) ne sera pas vérifier. Sinon, on continue jusqu'au moment où ne bouge plus trop, disons , et on refait une nouvelle fois avant de finir.

Si on applique l'algorithme de la section précédente tel quel dans le cas du système

on observe un effet d'enveloppement similaire au calcul de puissances successives d'un nombre complexe (voir la section3.2.2, ainsi que la figure5.1). Pour y remédier, on utilise les idées suivantes.

Idée 1 : calculer la dependence en fonction de la condition initiale

Idée 2 : encadrements de l'erreur dans des coordonnées déterminées par

Il est naturel de se demander s'il ne faudrait pas remplacer la définition de l'ensemble des modèles de Taylor sur le disque et d'ordre par quelque chose d'un peu plus fin, comme

Bien sûr, la définition classique est plus simple, ce qui rend les opérations de base plus efficaces. Toutefois, en grande précision, le calcul avec des boules dans ne coûte pas sensiblement plus cher que le calcul avec des flottants dans . Par ailleurs, la définition alternative admet aussi quelques avantages:

Examinons plus en détails la qualité des estimations lorsque l'on utilise la variante des modèles de Taylor de la section précédente. On suppose que l'on a calculé des encadrements pour les coefficients de la solution jusqu'à l'ordre et on pose

Notre but est d'obtenir un encadrement pour le reste

Soit la fonctionnelle définie par

En calculant en utilisant l'arithmétique pour les modèles de Taylor, on trouve une boule avec

Rappelons que notre but était de trouver un encadrement pour avec

Or

Nous en sommes donc réduit à trouver un avec

Supposant que

il suffit alors de prendre

L'intérêt de majorer le reste par une expression de la forme et non pas par une boule de la forme saute désormais au yeux: dans le deuxième cas, il aurait fallu demander que

et on aurait dû se résoudre à prendre

La taille de nos pas aurait donc été environ fois trop pessimiste.

La formule(5.3) montre aussi qu'il est intéressant de rendre aussi petit que possible. Ceci suggère de ne pas chercher des majorations directes dans les coordonnées d'origine, mais plutôt des majorations de la forme tel que soit petit.

On observe enfin que la condition(5.2) est vérifiée en prenant suffisamment grand. À partir de cette observation, on montre que le rayon pour lequel on peut certifier l'existence d'une borne tend vers le rayon de convergence de lorsque l'on fait tendre et vers l'infini. Dans la terminologie du chapitre2, cela implique

L'équation de Volterra est souvent utilisée comme benchmark pour des systèmes d'intégration certifiés. L'implantation des modèles de Taylor est en cours dans et dans la session qui suit nous trichons encore un peu quant à la «certification». Toutefois, l'exemple montre bien ce que nous cherchons à faire un jour.

D'une part, il existe des algorithmes d'une bonne complexité théorique pour calculer les racines d'un polynôme en une variable. D'autre part, il existe des bons implantations. La messe est-elle dite? Malheureusement, l'état de l'art présente encore plusieurs défauts.

D'une part, les algorithmes asymptotiquement rapides marchent bien quand la précision de calcul excède le degré du polynôme. Or une précision bien moindre suffit généralement pour isoler la plupart des racines. Par ailleurs, il est souvent plus rapide d'appliquer des algorithmes numériques brutaux en précision machine, plutôt que de payer un grand facteur constant pour l'utilisation d'une arithmétique «rapide» en précision multiple. D'autant plus que ces algorithmes brutaux sont généralement plus faciles à parallèliser. Manifestement le passage des algorithmes naïfs mais asymptotiquement mauvais à des algorithmes intelligents mais asymptotiquement bons pose encore un problème dans ce domaine.

En outre, trouver toutes les racines dans le cas général pose de sérieux problèmes techniques: souvent les bonnes idées ne s'appliquent que dans certains cas et il est difficile de les combiner. D'ailleurs, les algorithmes de Pan attendent toujours à être implantés.

Enfin, il est fort probable que les mesures de complexité traditionnels ne soient pas tout à fait opérationnels pour ce problème. En effet, il existe de nombreux types de polynômes(voir la section6.1) et l'efficacité et la qualité des solveurs dépend grandement de la répartition géométrique des racines. Faudrait-il des nombres de conditionnement plus fins pour exprimer la vraie complexité? Est-ce qu'il vaudrait mieux considérer un autre problème, comme la factorisation en deux facteurs de degrés deux fois moindre? Je ne connais pas les réponses pour l'instant, d'où le caractère plus «expérimental» de ce chapitre.

Meilleur conditionnement :

Pire conditionnement :

Pour différentes multiplicités :

On pourrait penser que le mauvais conditionnement est lié à la présence d'une racine de grande multiplicité (modulo petites perturbations du moins). L'exemple qui suit montre qu'il n'en est rien: des polynômes dont toutes les racines sont similaires en norme et dans le même demi plan sont à peu près aussi mal conditionnés qu'un polynôme de la forme.

Un algorithme numérique classique et efficace pour trouver les racines d'un polynôme

est la méthode d'Aberth. On part d'un ansatz pour les racines, comme

L'idée est ensuite d'appliquer pour chaque la méthode de Newton à la fonction

ce qui conduit à l'itération

Une propriété intéressante de cet algorithme est le fait que dès que l'on a trouvé des bonnes approximations de racines de , disons , alors les itérations suivantes appliquent essentiellement la méthode d'Aberth pour le polynôme

et les racines restantes . La méthode est particulièrement efficace quand on l'applique de façon naïve en précision machine (avec une complexité de calcul en ); c'est une des raisons de l'efficacité du logiciel .

La méthode par homotopie est un autre exemple d'une méthode qui est efficace quand on l'applique de façon naïve en précision machine. L'itération est un peu différente que pour la méthode d'Aberth, mais son efficacité est du même ordre de grandeur. On y reviendra dans un contexte plus général dans le chapitre7. Voir aussi.

La plupart des algorithmes asymptotiquement efficaces font intervenir la transformation de Graeffe , dont voici la définition:

On vérifie que

Il s'en suit la propriété fondamentale que et

En d'autres termes, la transformation est un séparateur puissant de racines, au moins en norme: si pour deux racines et de , alors pour les racines et de . Or, lorsque les racines de sont très différents en module, disons , alors on a pour chaque . Si les racines de sont toutes différentes en module, il suffit donc de quelques transformations de Graeffe pour obtenir pour un certain .

Reste le problème comment retrouver à partir de . On peut utiliser une astuce pour cela qui consiste à calculer avec des «nombres tangents» dans au lieu de simples coefficients dans . En effet, lorsque l'on fait le remplacement

et que l'on calcule les puissances -ièmes des racines de en faisant quelques transformations de Graeffe, on obtient pour chaque racine :

Nous venons de voir l'importance de l'étude des normes des racines du polynôme au delà de leur calcul à proprement parler. Le polygone numérique de Newton nous fournit les valeurs approximatives de ces normes, directement à partir des coefficients du polynôme.

Le polygone de Newton intervient aussi lorsque l'on veut développer une arithmétique efficaceet numériquement stable pour les polynômes. Une telle arithmétique reste typiquement précis pour le calcul de avec . Un tel algorithme de multiplication est donné dans , mais la borne(3) dans ce papier est fausse.

Le polygone numérique de Newton de est défini comme l'enveloppe convexe de l'ensemble

Moralité : Polygone de Newton Comportement en évaluation

En particulier : Pentes du polygone de Newton Modules des racines

La question suivante semble ouverte et abordable:

L'espoir

Plus réaliste

Pour un certain facteur dépendant du conditionnement et de :

Exemple d'un théorème précis:

C'est le problème inverse de la recherche des racines. En effet, ce problème est important lorsque l'on veut raffiner une approximation pour une précision en une meilleure approximation pour la précision (par exemple, en utilisant la méthode de Newton).

Arithmétique naïve en précision machine

Coût , en appliquant fois Horner.

Arithmétique rapide en précision multiple

Algorithme intermédiaire

, en faisant des FFT massives sur différents cercles de centre . Cet algorithme est meilleur que les précédents lorsque .

Les algorithmes de Schönhage et de Pan, qui admettent les meilleures complexités actuelles, certifient déjà les résultats. Toutefois, en suivant la même philosophie que dans les autres chapitres, il peut être utile de découpler les étapes de calcul et la certification des résultats.

Certification par la méthode de Krawczyk, présentée pour le cas de fonctions en plusieurs variables:

Ce théorème ne garantit pas l'unicité de la racine. Rump a donné une version plus fine du théorème qui garantit aussi l'unicité:

Pour trouver , faire

Méthode de -inflation de Rump.

Prendre pour et évaluer

Montrer que tout admet racines dans .

Par exemple, il suffit de vérifier que

et d'appliquer le théorème de Rouché. Voir pour des variantes de cette méthode.

Une autre approche consiste à examiner de combien il faut déformer le polynôme pour qu'il admette exactement une racine de multiplicité . Cette approche«par déflation» permet de se ramener à la théorie de la section6.6.1 en cherchant à résoudre simultanément les équations . Toutefois, l'astuce introduit un facteur dans la complexité en temps.

Dans la section1.3.1, nous avons esquissé comment résoudre des systèmes polynomiaux par homotopie. Cette technique ancienne remonte à la démonstration par Gauss que est algébriquement clos. En relation avec des implantations concrètes, les méthodes par homotopie apparaissent dans plusieurs contextes différents mais liés:

Homotopies numériques

Les homotopies numériques non certifiés ont été étudié par de nombreux auteurs. Voir par exemple et les références dans ces documents. Quelques logiciels connus pour des homotopies numériques sont et .

Complexité théorique des méthodes par homotopie

La complexité théorique des méthodes par homotopie a été analysé dans. Ces travaux donnent une justification théorique pourquoi les méthodes par homotopie marchent aussi vite, tout en faisant des hypothèses assez fortes sur le modèle de calcul.

Homotopies algébriques

Au lieu de calculer avec des nombres flottants, on peut aussi définir les trajectoires de manière algébrique. Ce point de vue a donné lieu à la méthode de Kronecker, qui admet également une implantation concrète. Par construction, cette méthode est algébrique donc certifiée, bien qu'elle ne tire pas profit de l'efficacité du calcul numérique pur. Par ailleurs, ce point de vue évite les chemins partant vers l'infini; voir la section7.1.2 plus bas.

Homotopies numériques certifiées

Il est naturel de chercher à certifier les méthodes numériques par homotopie, mais curieusement, peu d'efforts ont été déployés dans ce sens. Quelques exceptions: . Des implantations partielles sont présentes dans et depuis 2009.

Les méthodes par homotopie ne viennent pas sans leur lot de problèmes. D'une part, si le système est surdéterminé, le problème est mal posé d'un point de vue numérique. Dans ce cas, le problème est intrinsèquement irrésoluble par des voies exclusivement numériques, du moins si on souhaite certifier les solutions. Or, même en excluant des systèmes surdéterminés, et en considérant exclusivement des systèmes zéro dimensionnels avec équations et inconnus, les méthodes par homotopie comportent plusieurs inconvénients.

Premièrement, il est crucial que le système initial «facile à résoudre» ressemble autant que possible au système final qui nous intéresse réellement. En particulier, en prenant un système du même multi-degré, on espère que les deux systèmes admettent le même nombre de solutions. Malheureusement, tel n'est pas toujours le cas, comme on le voit sur l'exemple

qui admet seulement une solution au lieu de . Dans ce cas, trois des quatre chemins de l'homotopie partent vers l'infini. On peut remédier au problème des chemins qui partent vers l'infini en homogénéisant le système et en coupant par un hyperplan affine générique:

Pour ce nouveau système, il n'y a plus de chemins qui partent à l'infini, mais il existe toujours trois chemin qui tendent vers une racine multiple qui nous intéresse pas.

Dans ce cas précis, il aurait été possible d'éviter ces trois chemins parasites, en prenant des systèmes initiaux et finaux qui ne se ressemblent pas seulement en degré mais aussi quant à leurs polytopes de Newton. Pour plus de détails, voir par exemple et des travaux plus récents du même auteur. En général, il est difficile de trouver un système initial à résoudre qui ressemble suffisamment au système final pour que chaque solution du système initial se déforme en une unique solution du système final.

Une deuxième difficulté concerne les systèmes dégénérés avec des racines multiples. Reprenons l'exemple de la section précédente dans :

Pour une précision fixée, l'algorithme peine à converger vers la solution triple , et s'arrête bien avant d'avoir atteint une bonne qualité d'approximation. Ceci tient au fait que les méthodes naïves de suivi de chemin sont à convergence quadratique pour des solutions simples, mais seulement à convergence linéaire pour des solutions multiples. Une technique pour remédier à ce problème sera discutée dans la section7.3.

Méthode « prédicteur-correcteur »

Contrôle de pas

Critère numérique d'acceptation

Contrôle de la précision

Surestimation pour un certain seuil .

Toute méthode de certification d'une racine simple d'un système polynomial donne lieu àune méthode de certification pour des algorithmes de suivi de chemin en remplaçant les scalaires par des petite boules qui bornent la dépendance en de façon uniforme. On peut par exemple utiliser la méthode de Krawczyk-Rump de la section6.6.1 ou . Géométriquement parlant, ceci revient à certifier que, pour des dans une certaine boule, le chemin à suivre reste dans une certaine boîte comme dans la figure7.1.

La méthode de la section précédente conduit à des tailles de pas inutilement pessimistes. En effet à cause de la forme géométrique des encadrements, on ne souffre pas seulement de la surestimation de l'erreur due à la variation de , mais aussi de la surestimation en . Afin de réduire la surestimation en , il vaut mieux encadrer le chemin par un ensemble qui suit mieux la courbe. Dans , je montre comment faire cela en évaluant l'homotopies non pas en des boules, mais en des modèles de Taylor d'un type particulier.

Considérons une homotopie générale

avec et . Pour , toute solution est donnée par une série de Laurent algébrique

En particulier, en remplaçant par avec , on obtient un «troupeau» de solutions. Une méthode pour approcher des solutions multiples consiste à considerer tout le troupeau comme un unique chemin dans un nouvel espace. Cette méthode admet à la fois les avantages d'une déflation et des éclatements (des chemins distincts modulo et de même multiplicité ont des limites distinctes dans le nouvel espace). En utilisant la FFT, les troupeaux se suivent également bien d'un point de vue numérique.

Idée : Considérer comme « chemin » avec idéal

Représentation de (en position générale)

Relever l'homotopie

Évaluer en avec

Une question naturelle est de savoir quelle est l'efficacité des méthodes par homotopie par rapport à des algorithmes rapides pour calculer des bases de Gröbner. La session en dessous est éclairante à cet égard. Sur un problème plus ou moins générique (et de ce fait favorisant les méthodes par homotopie) avec solutions, on observe que le résultat du calcul des solutions certifiées prend environ une page, alors que la base de Gröbner en prend plus de (pour des raisons d'espace, nous avons tronqué les résultats ici). Plus qu'il y a de solutions, plus que ce ratio entre les tailles tend à devenir grand.

Bien sûr, la comparaison que nous venons de faire est favorable aux méthodes par homotopie. Toutefois, elle met en lumière un problème de fond concernant les bases de Gröbner: au moins sur certains exemples, la taille du résultat est beaucoup plus grande que nécessaire.

En outre, si on souhaite avoir une solution numérique, le simple calcul d'une base de Gröbner n'est pas suffisant. Même si on dispose n'une base pour l'ordre lexicographique, il faudra encore trouver les racines d'un gros polynôme en une variable. Or ce dernier problème pourrait se révéler plus dur que le problème de départ. En effet, il y a moins d'espace pour les racines dans que dans , ce qui nous impose de calculer avec une plus grande précision , comme on l'avait déjà constaté dans le chapitre 6. Par contraste, la précision machine suffit dans de nombreux cas, lorsque l'on applique des méthodes par homotopie sur le système de départ.

Toutefois, les systèmes surdéterminés et les chemins partant vers l'infini forment deux problèmes importants pour les méthodes par homotopie. Il est donc fort probable qu'il n'existe pas une technique unique idéale pour la résolution de systèmes polynomiaux.

Ce qui soulève les questions de pouvoir combiner plusieurs méthodes, et, plus modestement, de rendre les résultats d'une méthode utilisables pour d'autres méthodes. À l'intérieur des systèmes de calcul de bases de Gröbner, on dispose par exemple de l'algorithme FGLM pour passer d'un ordre sur les monômes à un autre. De la même manière, on pourra vouloir faire des conversions entre bases de Gröbner, représentations de Kronecker, systèmes triangulaires et des représentations numériques certifiées des variétés de solutions.

Par exemple, lorsque l'on dispose des racines approchées d'un polynôme unitaire , il est facile de reconstituer le polynôme en calculant et en retrouvant les coefficients par développement en fractions continues. Cette approche se généralise aux polynômes en plusieurs variables, et permet le calcul d'une représentation de Kronecker à partir de l'ensemble des solutions d'un système zéro-dimensionnel. Ceci permet en particulier l'utilisation de méthodes numériques par homotopie pour la résolution de systèmes surdéterminés.