Télécharger le fichier pdf d’un mémoire de fin d’études
Introduction au calcul scientifique et au GPGPU
L’apparition de l’informatique et sa démocratisation croissante depuis la deuxième moitié du XXe siècle constitue sans doute le changement le plus important dans la démarche scientifique et ses applications industrielles. La faculté des ordinateurs à effectuer un très grand nombre d’opéra-tions selon des instructions prédéfinies est en effet un atout majeur pour résoudre une multitude de problèmes formulés en langage scientifique. Si ses frontières ne sont pas immédiates à définir, le calcul scientifique désigne en général l’utilisation des ordinateurs pour résoudre un problème mathématique qui ne possède le plus souvent pas de solution analytique exacte, mais dont la solution peut s’exprimer au moyen d’un certain nombre d’équations à résoudre. Ces problèmes concernent en général la simulation, la compréhension ou le contrôle d’un système gouverné par les lois de la physique, et peuvent relever de la physique elle-même mais aussi de la biologie, de la chimie ou des sciences de l’ingénieur. La qualité et la précision de la solution obtenue est natu-rellement liée à la puissance de calcul disponible. Le plus souvent, une grandeur physique définie continûment dans la réalité est discrétisée en N morceaux, et le temps de résolution du problème est directement proportionnel à la finesse de cette discrétisation. Dans le cas de problèmes intrin-sèquement discrets comme on peut les trouver en biologie moléculaire ou en mécanique quantique, la simulation du comportement d’un grand nombre de particules est souvent désirée et une fois encore une grande puissance de calcul est importante. Pour tous ces problèmes, un compromis entre qualité de la solution et temps de calcul associé est nécessaire. Par ailleurs, l’augmentation de la fidélité de la modélisation passe aussi par la prise en compte d’un phénomène observé ex-périmentalement et dont on va rendre compte par l’addition d’équations supplémentaires, ce qui contribue à l’augmentation du nombre d’opérations nécessaires à la résolution du problème. Ainsi donc, une puissance de calcul supérieure permet de résoudre plus fidèlement un problème donné, ou peut rendre accessible une méthode considérée comme trop longue auparavant, ce qui justifie l’intérêt de son augmentation continue.
Évolution de la puissance de calcul des processeurs
Depuis le début des années 70, la puissance de calcul des processeurs, ou CPU (Central Proces-sing Unit) n’a cessé de s’accroître. Très rapidement, plusieurs conjectures ont été formulées pour apprécier cet accroissement exponentiel. La plus célèbre d’entre elles, la loi de Moore [1.10(a)], qui a été énoncée en 1965 puis corrigée en 1975, prédit un doublement du nombre de transistors sur une puce de silicium tous les 2 ans. Cette conjecture s’est révélée exacte depuis 1971 jusqu’à ces dernières années (2016), même si l’industrie des semi-conducteurs reconnaît que cette cadence ne peut plus être maintenue à présent. Pour mieux appréhender ce qu’implique réellement cette loi en termes de performances, on peut exprimer le temps d’exécution d’un programme de la façon suivante T = N × CHPI / F où N est le nombre d’instructions du programme, CHPI repré-sente le nombre de Cycles d’Horloge Par Instruction, et F est la fréquence du processeur, qui est inversement proportionnelle au nombre de cycles d’horloge.
L’augmentation du nombre de transistors a notamment permis la diminution des CHPI, en ajoutant de l’intelligence dans la gestion des différentes instructions à exécuter. Les micro-processeurs ont introduit le parallélisme du jeu d’instruction, qui permet d’exécuter certaines instructions indépendantes en simultané, ou d’en modifier d’ordre si les valeurs des quantités à additionner ou multiplier sont immédiatement accessibles. Pour la plupart des processeurs actuels, il est possible d’exécuter 4 Floating-Point OPerations (FLOP) lors d’un même cycle d’horloge sur le cœur de calcul d’un CPU. Les modèles dédiés au calcul comme les Intel Xeon Phi atteignent jusqu’à 32 FLOP par cœur. Des algorithmes de prédiction des instructions à opérer peuvent être utilisés pour augmenter le nombre d’opérations exécutées en parallèle. De même, un grand nombre de transistors du processeur sont associés à la gestion de la mémoire cache et à la prédiction des données à utiliser, car l’accès à la mémoire RAM nécessite plusieurs centaines de cycles d’hor-loge et peut être la source d’un ralentissement considérable du programme. Ainsi, l’utilisation de l’augmentation du nombre de transistors T a essentiellement contribué à diminuer le nombre de cycles par instruction, et le gain en performance associé est estimé à √T selon la loi de Pollack. Par ailleurs, grâce à la miniaturisation des transistors, la fréquence des processeurs a également augmenté de façon exponentielle, mais à un rythme inférieur car ces transistors sont des com-posants gravés√par lithographie sur des plaques à deux dimensions, ce qui implique aussi une croissance en T . Cependant, la fréquence d’horloge ne suit plus cette tendance depuis 2004 : la miniaturisation de ces composants a mis à mal la faculté à évacuer la chaleur que dégage ces transistors, et aucune solution n’a permis une réelle augmentation en fréquence, qui stagne de-puis autour de 3GHz. De plus, un lien très fort unit fréquence et consommation énergétique :
C ∝ U(f)2 × f. Pour des fréquences de l’ordre du GigaHertz, la tension U(f) est une fonction linéaire de la fréquence, ce qui induit une relation cubique entre fréquence d’horloge et consom-mation du processeur. La consommation électrique des processeurs a ainsi également augmenté continuellement pendant plusieurs décennies. Paradoxalement, l’efficacité énergétique a quant à elle aussi augmenté exponentiellement d’un facteur 1.57 tous les deux ans depuis 60 ans, comme le montre la loi de Koomey illustrée sur la figure 1.10(b).
Pour faire face à la limitation en fréquence, les processeurs sont actuellement conçus avec plusieurs cœurs de calcul indépendants. D’un point de vue d’une utilisation grand public, cette approche permet surtout l’exécution simultanée de plusieurs programmes potentiellement exi-geants. Elle permet également une véritable accélération des performances pour des programmes capables de se scinder en plusieurs processus, ce qui demande une réécriture de l’algorithme. Néanmoins, pour l’essentiel des applications, cette approche est sans effet sur ces programmes qui sont écrits de façon séquentielle, ce qui explique la multiplication très tempérée du nombre de cœurs des processeurs sur ces dix dernières années.
Ainsi, pendant près de 40 ans, un programme donné pouvait sans aucune modification néces-saire bénéficier d’une accélération significative en utilisant un processeur de nouvelle génération. Le changement de paradigme imposé par la programmation parallèle requiert quant à lui une réécriture au moins partielle d’un programme, en suivant des règles plus exigeantes, ce qui peut expliquer la relative inertie à la transcription de ces programmes, y compris dans le domaine scientifique. Cet effort de transition a toutefois été adopté par plusieurs solutions commerciales à but scientifique, et en particulier celles impliquant des modèles éléments finis, coûteux en temps de calcul. Le gain en temps d’exécution, restant tributaire du nombre de cœurs disponibles, de-meure modéré. C’est pourquoi l’intérêt des scientifiques désireux d’une puissance de calcul plus importante, et tirant davantage parti de l’augmentation exponentielle du nombre de transistors, s’est porté vers un autre composant informatique, le processeur graphique.
General-Purpose computing on Graphics Processing Units (GPGPU)
Les processeurs graphiques, ou GPUs ( Graphics Processing Unit ) constituent le moteur de calcul des cartes graphiques. L’objectif de ces cartes est d’alléger le CPU des tâches relatives à l’affichage, lequel a beaucoup évolué depuis les années 80. Initialement, il se résumait à l affi-chage de caractères ASCII, avant de pouvoir contrôler l’allumage de chacun des pixels de l’écran. Avec l’apparition de la couleur, codée sur un nombre croissant de niveaux au cours du temps, l’augmentation de la résolution écran, et enfin la gestion autonome de la représentation d’objets à 3D, la tâche de l’affichage s’est rapidement complexifiée. Ceci s’est traduit par une évolution hardware suivant le modèle Single Instruction Multiple Data (SIMD) conduisant en particulier à l’augmentation du nombre d’unités logiques de calcul et de la bande passante mémoire. De plus, la complexification progressive du jeu d’instructions à exécuter a également influé sur l’architecture matérielle, se rapprochant un peu plus de celle d’un CPU avec l’apparition de plusieurs niveaux de mémoire cache. Cette modification hardware a également permis une évolution très impor-tante côté software. Du côté de l’affichage, les fonctionnalités proposées par l’API OpenGL se sont assouplies, permettant au programmeur de définir lui-même le mouvement et la modification structurelle des objets à 3D à représenter. Ces possibilités de programmation ont rapidement été exploitées pour des objectifs autres que l’affichage écran : au début des années 2000, le terme General-Purpose computing on GPU (GPGPU), a été inventé par Mark Harris pour désigner cette pratique. En 2007 puis en 2009, l’apparition des plate-formes CUDA, développée par le fabriquant de GPU Nvidia, et OpenCL ont grandement contribué à faciliter l’écriture d’un code scientifique utilisant le GPU comme outil de calcul. De nombreuses publications faisant état des gains en performance [20, 21] sont alors parues, et le GPU s’est depuis positionné comme un acteur incontournable du calcul haute-performance. Depuis 10 ans, l’écart de performance entre CPU et GPU n’a cessé de s’accroître, comme on peut le voir sur la figure 1.11(a). On peut l’expliquer par la demande du principal acteur de ce domaine, le secteur des jeux vidéos, qui souhaite obtenir un rendu toujours plus réaliste et s’adapter aux normes de haute définition des écrans, ce qui passe inéluctablement par une augmentation du nombre d’opérations à effectuer par seconde. Du point de vue du composant, cette accélération est due à une augmentation drastique du nombre d’unités logiques de calcul présentes sur le GPU : les GPU actuels (2016) en contiennent plusieurs milliers. Une autre caractéristique non moins importante à observer est la bande passante du composant, qui va conditionner la vitesse à laquelle les données vont être acheminées vers l’unité logique de calcul, et est le principal frein à l’exploitation de la puissance de calcul du matériel. Elle affecte en particulier les problèmes qui nécessitent l’accès à de larges tableaux, ce qui inclut la plupart des situations de simulation numérique, dépendantes de plusieurs variables d’entrée. Une fois encore, l’avantage est donné au GPU comme montré figure [1.11(b)], ce qui était prévisible dans l’optique de l’exploitation du très grand nombre d’unités de calcul qui le constitue. Enfin, la performance énergétique entre CPU et GPU, mesurée en nombre d’opérations flottantes par watt est aussi à la faveur des GPU [1.11(c)]. La fréquence d’horloge en moyenne trois fois moins élevée pour les GPU explique en grande partie ce phénomène, du fait de la dépendance cubique entre fréquence et consommation.
Un autre avantage significatif des GPU est qu’en plus d’offrir une puissance de calcul très supérieure, d’environ 2 ordres de grandeur actuellement, le coût de ces deux composants est équi-valent, pour des produits d’entrée, milieu ou haut de gamme. C’est sans doute cette caractéristique a priori surprenante qui est à l’origine du changement progressif de pratique dans le calcul scien-tifique. Et si de plus en plus de superordinateurs incorporent des GPUs, un gain substantiel peut être également obtenu pour un scientifique faisant l’acquisition d’une nouvelle carte graphique pour son ordinateur de travail : pour obtenir un gain du même ordre avec des CPUs, les factures matérielles et énergétiques seraient bien plus élevées. Il est intéressant d’observer que si actuel-lement la performance d’un GPU outrepasse de loin celle d’un CPU, le nombre de transistors utilisé pour les deux composants est sensiblement le même, ce qui explique la similarité du coût des deux composants. Pour un CPU, une augmentation du nombre de transistors est surtout mise à profit pour diminuer le nombre de cycles d’horloge par instruction, alors qu’elle est utilisée pour multiplier le nombre d’unités de calcul pour un GPU. Comme nous l’avons vu, la première stratégie a atteint ses limites et alors que la seconde continue de bénéficier linéairement de l’aug-mentation des transistors. De plus, l’architecture plus complexe des GPU requiert souvent plus de temps pour s’adapter à une nouvelle finesse de gravure, les GPU actuels (2016) étant gravés en 28nm alors que ceux des CPU le sont en 14nm. Les GPUs ont ainsi plus de potentiel quant aux limites de la gravure lithographique. Enfin, la stratégie de la multiplication du nombre de d’unités de calcul n’entraîne qu’une dépendance linéaire entre augmentation de la puissance de calcul et consommation électrique, et qui est de plus sous linéaire d’une génération à l’autre grâce à la miniaturisation des transistors. Le modèle GPU est donc tout à fait pertinent du point de vue énergétique.
Ainsi, la parallélisation d’un code apparaît comme une étape inévitable pour continuer à béné-ficier du progrès informatique à chaque nouvelle génération de composants. Il demeure cependant quelques restrictions à son usage :
— Le défaut le plus pénalisant est incontestablement la nécessité de devoir coder dans un langage dédié, mais aussi avec un algorithme repensé. De fait, de par l’architecture parallèle des GPU, le calcul sur carte graphique utilise un modèle de programmation spécifique et son apprentissage ne relève pas simplement de la connaissance d’un nouveau langage.
Pour une programmation efficace, il est primordial que l’utilisateur ait conscience de la structure, des limitations et des propriétés du composant qu’il utilise et de leur évolution. Cet investissement d’apprentissage est important du côté des utilisateurs, et de plus le paradigme de programmation plus complexe oblige à consacrer plus de temps que pour une version séquentielle ou parallèle sur CPU. Enfin, pour des problèmes déjà implémentés sur CPU, tout le travail est à refaire.
— La parallélisation du code n’est pas toujours possible. Idéalement, le problème à résoudre doit suivre le modèle SIMD, ce qui constitue la plupart des cas.
— Le nombre d’instructions par donnée chargée doit aussi être suffisamment élevé pour que le temps de charge limité par la bande passante puisse être compensé. Le calcul d’un produit scalaire effectuant une seule addition pour deux données chargées n’offrira pas de bonnes performances. En revanche, pour un produit matriciel N × N, chaque donnée est utilisée N fois, ce qui permet une utilisation totale des unités de calcul si N dépasse une certaine valeur. Un problème de grande taille est par ailleurs toujours souhaitable : si le problème est scindé en un nombre de processus supérieur au nombre d’unités logiques, le GPU peut compenser le temps de chargement de certaines données par des calculs sur les processus dont les données ont déjà été chargées. Cette situation reste la plus courante.
— La performance s’applique surtout à la simple précision, car c’est elle qui est utilisée en quantité par les jeux vidéos, ce qui incite les fabricants de GPU à augmenter ce type d’unité logique. La dernière génération de GPU Nvidia, Maxwell, a fait le sacrifice sur les unités de calcul en double précision : la performance théorique maximale passe de 7 TFLOPS en simple précision à 200 GFLOPS en double précision. Sur CPU, l’écart entre les deux précisions est habituellement d’un demi. Remarquons toutefois que si la double précision est indispensable dans certains cas, elle ne l’est pas le plus souvent, mais est un confort qui évite au programmeur de trop se soucier des contraintes de l’arithmétique finie. En effet, beaucoup de problèmes liés à la précision numérique ont une alternative algorithmique.
Implémentation sur GPU
Dans cette section, nous exposons les détails de l’implémentation de l’algorithme FTIM sur GPU, puis comparons les gains obtenus par rapport à la version CPU. Au delà de la simple va-lorisation du temps passé à l’écriture du code, l’intérêt de cette partie est de fournir un exemple concret pour le lecteur néophyte en matière de programmation GPU et désireux d’en découvrir les aspects primordiaux. L’accent est mis sur la façon dont l’algorithme est pensé pour s’adapter au GPU, puis repensé afin d’éviter les principaux goulots d’étranglement. Pour la bonne com-préhension de cette partie, il est suggéré au lecteur non habitué à la programmation sur GPU de se référer à l’annexe A, où le vocabulaire et les bases du paradigme de programmation sur GPU sont exposées. Pour le calcul des transformées de Fourier, les fonctions prédéfinies dans la librairie cuFFT de Nvidia ont été utilisées. Comme nous l’avons vu, le temps de calcul des trans-formées de Fourier est négligeable dans ce problème. L’objectif de cette partie est de calculer le plus rapidement possible la quantité IF T IM qui s’exprime, dans le cas à 2D par exemple, ainsi : IˆF T IM (i, j) = Nfreq Nelem H(i −p(m −1), j, fl)Sm(fl) Nelem X X X H(i −p(m −1), j, fl)Bscanm(fl) l=1m=1 m=1 (2.17) avec p l’écart en pixels entre deux éléments piézoélectriques, Sm la transformée de Fourier du signal émis par le mème élément, et Bscanm la transformée de Fourier du signal reçu par le mème élément.
L’inversion de la forme d’onde complète
L’inversion de la forme d’onde complète, à laquelle nous référerons par son acronyme anglais, FWI (Full Waveform Inversion), est une technique d’imagerie quantitative issue de la géophysique. En toute généralité, la FWI peut être vue comme une méthode d’assimilation de données, où les données observées sont incorporées avec un traitement mathématique rigoureux à un modèle numérique décrivant la réalité. Mathématiquement, elle repose sur la minimisation itérative de l’écart entre données réelles et données simulées numériquement, en modifiant localement les propriétés physiques du milieu qui influencent la propagation d’onde. Cette méthode a été initiée par les travaux de Claerbout en 1971 [35], qui visaient à reconstruire à partir de données en réflexion les couches des réservoirs de pétrole situés dans le sous-sol. En 1984, Tarantola [36] a exprimé à nouveau le problème en n’impliquant que le calcul de deux champs acoustiques par source et par itération. Cependant, en dépit de cette reformulation, le coût numérique de ce calcul est resté très longtemps prohibitif même à deux dimensions. De plus, l’utilisation de données essentiellement en réflexion n’a pas apporté les résultats escomptés : au lieu de reconstruire une carte de densité et de module de compression, ces méthodes ont montré en pratique qu’elles permettaient de retrouver les interfaces internes du milieu. Ainsi, on réfère à ces deux travaux fondamentaux comme principe d’imagerie, ou technique de migration. A cause du temps de calcul, des modèles de propagation basés sur l’équation eikonale ont d’abord été utilisés et des résultats pertinents sur la structure globale de la Terre à une dimension ont été obtenus, grâce à l’utilisation d’ondes en transmission. Puis, avec l’évolution progressive de la puissance de calcul, l’utilisation de modèles plus réalistes ont permis d’incorporer les effets 2D et 3D, mais aussi de prendre en compte l’aspect fini des fréquences des ondes sismiques. La compréhension de l’intégralité du contenu des sismogrammes est devenue possible et interprétable, et c’est en ce sens que l’on emploie le terme inversion de la forme d’onde complète.
L’objectif de ce chapitre est de présenter la FWI. La formulation mathématique du problème d’optimisation ainsi que les différentes techniques amenant à sa résolution seront exprimées. En-suite, nous verrons les éléments clés de l’inversion, qui influent d’une part sur la capacité à conver-ger vers la solution la plus proche de la réalité en dépit de la forte non linéarité du problème et de la non-unicité possible de la solution, et d’autre part sur la vitesse à laquelle il sera possible de converger vers la solution, qui s’avère en pratique être aussi un facteur critique étant donné la taille du problème et le coût numérique associé. Ces facteurs peuvent être liés à la fois au problème mathématique d’optimisation en lui-même, mais aussi être spécifiques à l’implication de l’équation d’onde. Enfin, comme cette méthode fait intervenir la simulation numérique de la propagation d’onde et que la précision avec laquelle elle est simulée joue également un facteur déterminant, nous décrirons la méthode numérique employée.
Théorie
Mise en équation du problème
Tout comme pour la formulation du problème d’optimisation topologique décrit au début du chapitre précédent, la FWI repose sur la minimisation d’une fonction f(m) qui évalue l’écart entre données réelles pobs(xi, t) et données simulées psyn(xi, t, m), ce qui s’exprime : 1 Nrec T pobs(xi, t) − psyn(xi, t, m) 2dt min f(m) = 2 (3.1)
De même, la modification locale et itérative d’un ou de plusieurs paramètres physiques m du modèle du milieu à imager est le moteur de cette minimisation, et la combinaison finale de chaque paramètre, qui s’apparente à une carte spatiale des paramètres en question est l’image recherchée. Dans le cas de l’optimisation topologique, le paramètre utilisé est la topologie Ω du milieu, et l’image obtenue définit la forme des bords et des trous tels qu’ils sont dans la réalité. Pour la FWI appliquée en géophysique, les paramètres physiques sont des grandeurs continues telles que la densité ρ(x), le tenseur d’anisotropie cijkl(x) ou l’atténuation α (x) du milieu. Dans les deux cas, la propagation d’onde influence les données, ce qui se traduit mathématiquement par un problème d’optimisation contraint par une équation aux dérivées partielles. Dans le cas de l’acoustique linéaire, avec la vitesse de compression c et la masse volumique ρ comme fonctions de la position x, en introduisant le module de compressibilité isotrope κ = ρc2, et en définissant un potentiel scalaire Φ à partir du déplacement s : s = rΦ (3.2)
Cette définition inhabituelle du potentiel à partir du déplacement est héritée de l’origine géo-physique de cette méthode, qui permet une écriture simplifiée des phénomènes de couplage entre milieux fluide et solide.
On définit alors l’équation d’onde suivante, dont la pression est la dérivée seconde : psyn(xr, t) = −∂t2 Φ(xr, t) κ(x)−1∂t2 Φ − r (ρ(x)−1 r Φ) = Ns κ(xs)−1fs(t) (3.3)
Les termes sources fs(t) sont des sources de pression, ce qui correspond par ailleurs à la nature physique des données enregistrées, et montre la pertinence de cette formulation en potentiel. Le cas échéant, une double dérivation numérique des signaux enregistrés aurait été nécessaire, ce qui aurait augmenté le bruit numérique.
On pourra remarquer que l’implication de l’équation d’onde avec des cartes de paramètres c(x), ρ(x) ou κ(x) rend le problème d’optimisation non convexe. De plus, la relative liberté quant au nombre de sources, leur position spatiale et la durée des signaux enregistrés influence largement la possibilité à obtenir la vraie solution. Si la réalisation pratique de l’expérience permet d’affirmer l’existence de la solution en supposant l’absence de bruit expérimental, l’unicité n’est pas garantie non plus, à cause de la non convexité du problème.
Résolution du problème inverse
Choix de la méthode de résolution
De façon générale, le coût numérique d’évaluation de la fonction coût joue un rôle primordial dans le choix de l’algorithme à utiliser. Normalement, pour un problème non convexe comme ici, les algorithmes d’optimisation globale, tels que des méthodes stochastiques ou heuristiques, sont en pratique les plus efficaces quant à la convergence vers l’un des optima globaux. Cependant, ces derniers impliquent un très grand nombre d’évaluations de la fonction coût et/ou de son gra-dient, ce qui les rend inutilisables de façon pratique dès lors que cette évaluation se chiffre en minutes ou plus. Ici, l’évaluation de la fonction coût est très coûteuse car elle n’est pas obtenue de façon analytique ou semi-analytique, comme cela aurait pu être le cas si le milieu de départ était homogène, mais que seul les méthodes numériques permettent d’inclure proprement dans le calcul de la fonction coût l’ensemble des conditions aux limites, ainsi qu’une carte spatiale poten-tiellement variable de l’ensemble des paramètres impliqués dans l’équation d’onde. L’utilisation de ces méthodes numériques pour évaluer la fonction coût, qui dans notre cas est un simulateur de propagation d’onde, est systématiquement coûteuse en temps de calcul, et inévitable ici. Pour limiter le nombre d’estimations de la fonction coût, le seul moyen est de passer par des méthodes d’optimisation locale, de type descente. Le principal inconvénient est lié au caractère local de l’optimisation, qui dans le cas d’un problème non linéaire et non convexe se traduit par le risque fort de converger vers un minimum local et non global. Étant donné la forte non convexité du problème, on peut se représenter la fonction coût comme étant une planète entièrement remplie de montagnes, et l’optimisation locale consiste à descendre vers le lac ou la vallée la plus proche, en fonction de la position initiale choisie, alors que l’on souhaite trouver le point de plus basse altitude existant. Ainsi, sans reformuler plusieurs fois au cours de l’algorithme le problème d’optimisation afin de diminuer le niveau de non-linéarité, la convergence est impossible, et nous verrons no-tamment que la connaissance du problème physique concerné est d’une importance cruciale pour formuler les bonnes hypothèses.
A présent, nous allons voir les principales méthodes de descente, qui permettent la résolution du problème d’optimisation locale. Elles s’expriment autour de la formulation suivante : ( x0 est donné αkdk (3.4)
où xk est la suite d’itérés qui doit converger vers le minimum, αk est le pas de descente et dk la direction de descente.
La différence entre les méthodes que nous allons présenter réside dans le choix de la direction de descente et du pas.
Méthodes de type gradient
La première famille de méthodes consiste à utiliser le gradient rf (x) comme direction de descente. En effet, son action sur la fonctionnelle est optimum à l’ordre 1 dans le sens où : (−rf(x))T rf(x) ≤ dT rf(x) ∀d ∈ Rn tel que ||d|| = ||rf(x)|| (3.5)
Différentes stratégies sont possibles quant au choix du pas α : avec une certaine connaissance du problème, on peut simplement fixer un pas constant pour toutes les itérations. La recherche du pas peut elle même être vue comme un problème d’optimisation visant à résoudre : αk>0 kr min( xk − α f( xk ))) (3.6)
La méthode correspondante, dite de pas optimal, ou de plus grande pente (steepest descent) , peut de prime abord paraître la plus plus intéressante dans la mesure où elle minimise le plus la fonction coût. En pratique, cette méthode est considérée comme mauvaise car elle ne converge que très lentement : par construction, deux directions de descente successives sont orthogonales, et il résulte une trajectoire de descente qui peut s’apparenter à de multiples ’zig-zag’. Pour la recherche du pas de descente, les conditions de Wolfe [37] sont reconnues comme étant les plus efficaces et sont utilisées par la plupart des algorithmes d’optimisation. Elles sont applicables pour tout type d’algorithmes de descente et s’articulent autour des deux conditions suivantes :
— La condition d’Armijo, qui permet d’éviter de choisir un pas trop grand, reformule la condi-tion de décroissance de la fonction coût d’une itération à l’autre ainsi : f(xk + αkdk) ≤ f(xk) + c1αkdkrf(xk) avec 0 < c1 < 1 (3.7)
— La condition de courbure évite de choisir un pas trop petit : dTk rf(xk + αkdk) ≥ c2dTk rf(xk) avec 0 < c2 < 1
Lorsque ces conditions sont utilisées ensemble, on veille aussi à respecter la condition c1 Dans la pratique, on choisit souvent c1 = 10−4 et c2 = 0.99. (3.8)
L’utilisation du gradient comme direction de descente assure la convergence locale vers l’op-timum, mais la vitesse de convergence associée reste souvent assez faible. Dans la pratique, cet algorithme est apprécié pour sa simplicité de mise en œuvre, mais son utilisation est plutôt destinée aux problèmes de petite dimension. Pour des problèmes impliquant une fonction coût quadratique, de type f(x) = xT Ax−bT x, une version modifiée de la descente de gradient existe : la méthode du gradient conjugué. La direction de descente choisie utilise le gradient de l’itération courante, mais aussi ceux des itérations passées, et en utilisant le procédé d’orthogonalisation de Gram-Schmit : k 1, d = f( ) ∀ ≥ k −r xk k,j j k,j h i A j=1 di, di + k−1 α d , avec α = hrf(xk), diiA . (3.9) La force de la méthode est de pouvoir converger vers l’optimum en un nombre d’itérations inférieur
à la dimension du problème. Cependant, pour des problèmes de très grande dimension, et où le calcul du gradient est plus long, cette méthode souffre toujours d’une vitesse de convergence trop faible, qui est linéaire.
|
Table des matières
Introduction
1 Etat de l’art
1.1 Panorama des procédés d’imagerie en acoustique ultrasonore
Le SONAR et l’échographie
Imagerie quantitative : exemple de l’élastographie
1.2 Introduction au calcul scientifique et au GPGPU
Évolution de la puissance de calcul des processeurs
General-Purpose computing on Graphics Processing Units (GPGPU)
2 Le procédé d’imagerie FTIM
2.1 Principe
Origine théorique
Interprétation dans le domaine temporel
Transposition dans le domaine fréquentiel
Avantages et limites de la méthode
Calcul du diagramme de rayonnement pour un milieu homogène semi-infini
2.2 Implémentation sur GPU
Cas bidimensionnel
Cas tridimensionnel
2.3 Application à la PIV
3 L’inversion de la forme d’onde complète
3.1 Théorie
Mise en équation du problème
Résolution du problème inverse
Calcul et considérations sur le gradient
3.2 Principaux facteurs influant sur la convergence et la rapidité de la méthode
Filtrage des données
Fonction coût
Préconditionnement
Régularisation du problème non linéaire
3.3 Simulation numérique de la propagation d’onde
La méthode des éléments finis spectraux
Application à l’équation d’onde
Considérations numériques et parallélisation
4 Application de la FWI à l’échelle ultrasonore : reconstitution de la carte de vitesse d’un milieu inconnu
4.1 Spécificités de l’échelle ultrasonore
4.2 Faibles contrastes de vitesse
4.3 Forts contrastes de vitesse
Page iTable des matières
Modèle avec interfaces
Modèle sans interfaces
4.4 Vers l’inversion de données réelles
Changement de philosophie par rapport au dispositif expérimental échographique
Digitalisation
Bruit
Conclusions et perspectives
Annexe A Paradigme de programmation sur GPU
A.1 Relation avec le CPU
A.2 Architecture matérielle du GPU et code associé
A.3 Les différents types de mémoire .
Télécharger le rapport complet