HOME | ARCHITMAHBIL |

LE HASARD FAIT BIEN LES CHOSES

(et bien DES choses)

ou comment le hasard nous aide à faire du calcul numérique

Titres:

A télécharger:

Il y a longtemps (une dizaine d'années) après avoir lu dans un bouquin de physique (le Berkeley machin de physique, le tome consacré à l'électricité) sur la manière de calculer numériquement le champs électrique avec des sources arbitraires, j'ai essayé d'appliquer tout ça et écrire un programme dans mon langage préféré: le Pascal (d'ailleurs certains à l'époque m'appelaient monsieur Pascal ce qui est un honneur). Je ne savais pas ce qui m'attendait. Mais non! je n'ai pas eu de difficultés à l'implémenter; voyons! C'est fastoche:

L'équation de Laplace:

L'équation qui décrit le champs électrique et dit en passant le champs gravitationnel, propagation de la chaleur, flexion d'une plaque et encore… cette équation connue sous le doux nom d'équation de Laplace s'écrit d²f/dx²+d²f/dy²=0 (les d c'est des d rond en principe, ne pas confondre avec différentielle totale). Je ne vais pas vous accabler avec les étonnantes propriétés de cette équation car, ou bien vous la connaissez ou bien vous n'avez pas besoin d'en connaître plus (en fait j'ai pas envie d'en dire plus et c'est pas nécessaire pour la suite). Donc, pour la résoudre il faut prendre une boisson, une pizza, s'asseoir confortablement, puis prendre un tableau bi-dimensionnel (ou à n dimensions, peu importe, mais pour le reste c'est bi), des conditions aux limites (des éléments dudit tableau dont les valeurs sont fixés une fois pour toute) puis répéter: pour chaque élément du tableau, remplacer sa valeur par la moyenne de ses quatre voisins[*]; jusqu'à la convergence ou à l'épuisement ou bien à la manifestation suivante du bogue de Windows (si vous utilisez un autre OS c'est pas grave, à chacun ses bogues). Bien sûre c'est pas la seule méthode qui existe ni la meilleure ni la plus rapide, mais je crois que c'est la plus simple.

Le problème des erreurs d'arrondi:

Bien entendu, vous allez utiliser des nombres flottants pour faire tous les calcul. D'ailleurs c'est assez logique (ou plutôt logique tout court). Mais faites attention, si vous utilisez des nombres avec peu de chiffres significatifs vous n'aurez pas le résultat escompté. Ça ne risque pas d'arriver même en utilisant des SINGLES (ou float pour ceux qui préfèrent le C) quand le tableau (réseau) est de petites dimensions, mais il existe bien une limite à cause des erreurs d'arrondi. Plus le réseau est grand et plus le nombre de chiffres significatifs des nombres utilisés doit être grand et dépend de plus, de l'équation à résoudre et de la manière dont elle est résolue. A l'époque j'étais loin de me douter de telles considérations.

Ma première implémentation faisait le calcul sur un tableau de 64x64 qui était affiché à l'aide du mode VGA 256 couleurs. Le programme s'exécutait très rapidement (sur un 486) et l'effet me plaisait bien (l'animation due à l'affichage successif des différentes itérations donne l'impression d'un feu ravageant l'écran, plutôt esthétique:) j'ai voulu donc faire le calcul sur un réseau de 320x200 couvrant tout l'écran. Mais comme j'était sous DOS, la limitation de la taille des données à 64ko m'a mené à essayer d'utiliser juste un octet par cellule du tableau. Il est évident les valeurs calculées ne dépassent jamais les valeurs des conditions aux limites… Et là, la déception. Le calcul était complètement faux à cause des arrondis.


Fig 1: De gauche à droite: calcul fait avec des flottants; entiers et arrondi trunc() (apres division);
arrondi round(); arrondi ceil() et enfin avec notre méthode.
Conditions aux limites: 0 sur les bords et 127 au point milieu. Résultat aprés environ 2000 iterations.

Pour bien faire passer le message, je vais expliquer plus en détails: j'avais forcé les cellules des bords à 0 et mis une "source" (de chaleur ;) au milieu du réseau. C'était ça les conditions aux limites que j'avais données. En principe, au cours du calcul, la "source" doit diffuser jusqu'aux bords comme c'est le cas avec les nombres flottants. Mais il n'en est pas ainsi si on calcule sur des octets: La source diffuse sur une distance finie (quelques cellules) puis s'arrête.

On aurait dit qu'il y a un fond qui absorbe "l'énergie" du réseau (les chiffres significatifs). J'avais donc essayé différent types d'arrondi mais le résultat était tout autant désastreux. C'est là que j'ai eu l'idée d'utiliser les nombres aléatoires pour l'arrondi. En effet, à priori je n'ai aucune information sur la manière dont les informations se perdent au cours des arrondis.

La ou ça devient intéressant :

La méthode utilisée est très simple : au lieu de faire :

tab[i,j] :=(tab[i-1,j]+tab[i+1,j]+tab[i,j-1]+tab[i,j+1]) div 4 ;

il suffit de faire :

tab[i,j] :=((tab[i-1,j]+tab[i+1,j]+tab[i,j-1]+tab[i,j+1])+random(4)) div 4 ;

Il est évident que la somme des quatre voisins peut tenir sur plus d'un octet. Il faut donc faire le calcul intermédiaire sur des mots d'au moins 10bits. D'une manière générale, si nous stockons les valeurs des cellules dans des mots de n bits, il faut faire le calcul sur des mots de n+2 bits.

Ce fut assez incroyable, en exécutant ce code de voir que son comportement approche celui fait avec des nombres flottants. La différence réside dans le fait que les valeurs des cellules " sautent " entre les entiers qui encadrent la valeur calculée avec des nombres flottants même après convergence, ce qui est prévisible. Au fait, je n'ai pas trouvé de test de convergence qui soit pertinent, ça reste une question ouverte.

Plus incroyable encore est le fait que la moyenne de la solution entière (je suis tenté de dire stochastique :) est apparemment égale à la solution réelle (hic). J'ai fait des essais avec différentes conditions aux limites et c'est toujours bon J. C'est magique !

Pourquoi en est-il ainsi :

Bonne question. La vérité est que j'en sait trop rien (comme il a dit lui ;). C'est d'ailleurs l'une des raisons qui m'a amené à publier ce torchon, donc à bon entendeur merci. Je vais quand même tenter une réponse (plutôt quelques éléments de réponse, c'est moins prétentieux) à mes risques et périls.

Tout d'abord tentons quelques expériences. Nous allons faire la comparaison d'une part entre la solution " stochastique " et la solution " réelle " et de l'autre entre la moyenne de la solution stochastique et la solution réelle. Qui dit comparaison dit soustraction donc nous allons étudier en fait les différences R-S et R-M. (Vous avez deviné, j'en suis sure ;)

Stochastique vs réel (R-S):

Fig 2: A gauche, la solution calculée grace à notre méthode (S). A droite, la solution calculée avec des nombres flottants (R).
Fig 3: La différence (R-S) les zones avec des couleurs sombres sont négatives.

Fig 4: Partie decimale de (R). A comparer avec la figure 3.

(du rouge au pourpre les valeurs allant de 0.0 à 1.0)

Pour chaque cellule cette différence est le plus souvent inférieure à 1 mais elle peut être décalée de 2 ou même de 3 unités mais beaucoup plus rarement. Je présume que les valeurs (entières des cellules) sont " choisies " de manière à se rapprocher le plus possible de la valeur réelle. La probabilité qu'une cellule prenne la valeur la plus proche de la solution réelle est maximale puis diminue brutalement avec l'écart à cette solution. On peut supposer qu'il s'agit de quelque chose comme une distribution gaussienne, mais il ne faut présumer de rien (bien que je l'ai déjà fait J).

Vous pouvez remarquer la structure qui apparaît sur l'image (R-S), c'est tout bonnement la partie fractionnaire de la solution réelle.

Moyenne vs réel (R-M) :

Fig 5: La Moyenne de la solution calculée grace à notre méthode (S)->(M) aprés quelques dizaines d'itérations.
Fig 5: La différence entre la solution dite réelle et la Moyenne (celle de la fig 5) (R-M). Les couleurs sombres correspondent à des valeurs négatives. (R-M) dépend du nombre d'itérations.

Tout d'abord, il faut bien entendu commencer à calculer la moyenne quand nos solutions réelle et stochastique convergent. A partir de cet instant, la différence (R-M) prend l'aspect d'un terrain fractal, une texture de type perlin. Notez que les valeurs sont toutes inférieures à un et diminuent graduellement. De plus les petites structures disparaissent plus rapidement que les grandes. On peut montrer facilement que cette différence tend vers zéro. En effet, soit R la solution réelle et Mi la moyenne à l'itération i après convergence. Nous avons :

Mi = somme(Si,i)/i ; Si est la solution stochastique à l'itération i après convergence.

R - Mi = R - somme(Si,i)/i = somme(R-Si,i)/i ;

Nous savons (même si ce n'est pas encore démontré) que la somme de tous les R-Si est nulle. Donc R-Mi tend vers zéro quand i tend vers l'infini. Cette convergence, vous vous en doutez est très lente mais effective.

Vous pouvez toujours utiliser R-M pour générer des terrains fractals ou des textures. Mais vous pouvez aussi utiliser ce programme à la place. Celui-ci simule R-M sans avoir à calculer R, S et M. vous pouvez étudier son fonctionnement dans ces sources (voir le programme lapnoise.exe).

Jusqu'à maintenant nous ne savons toujours pas pourquoi nous obtenons ces résultats. En essayant de monter tout ça j'ai eu d'affreuses migraines, les calculs étant trop compliqués pour moi. Mais je vais quand même aborder le sujet en résonnant sur un cas simple.

Supposons que les nombres que nous utilisons n'ont qu'un seul bit (pas d'arrières pensées hein ;). Je voulais dire par là que les valeurs des cellules ne peuvent prendre que les valeurs 0 ou 1 (mais si même dans ce cas ça marche J). Si une valeur approche un réel r compris entre 0 et 1 sa probabilité d'être à 1 est r et (1-r) pour 0. Prenons une cellule quelconque. Chacune de ses voisines à une probabilité p d'être égale à 1 et (1-p) d'être à zéro. Soit p1 p2 p3 p4 ces probabilités d'être à 1. Ici on fait comme pour calculer la probabilité d'avoir telle ou telle somme avec un jeu de quatre pièces (n'oubliez pas que la formule reste centre = moyenne des voisines).

Pour avoir une somme de 4 nous avons la probabilité :

P4 = p1*p2*p3*p4.

Avant de continuer, remarquer qu'on prend le produit car les variables sont indépendantes. Ça marche ainsi, mais il reste à prouver qu'elles le sont (indépendantes) et je ne sais pas comment.

Pour avoir une somme de 3 la probabilité est :

P3= (1-p1)*p2*p3*p4+p1*(1-p2)*p3*p4+…(devinez la suite)

P3=p1*p2*p3+p1*p2*p4+p1*p3*p4+p2*p3*p4- 4*p1*p2*p3*p4.

Pour avoir une somme de 2 la brobabilité est:

P2= (1-p1)*(1-p2)*p3*p4+(1-p1)*p2*(1-p3)*p4+...

P2=6*p1*p2*p3*p4-3*p1*p2*p3-3*p1*p2*p4-...-3*p2*p3*p4+p1*p2+...+p3*p4.

Pour 1:

P1=(1-p1)*(1-p2)*(1-p3)*p4+...

P1= -4*p1*p2*p3*p4+3*p1*p2*p3+...+3*p2*p3*p4-2*p1*p2-...-2*p2*p4+p1+p2+p3+p4.

et pour 0:

P0=(1-p1)*(1-p2)*(1-p3)*(1-p4).

P0=p1*p2*p3*p4-p1*p2*p3-...-p2*p3*p4+p1*p2+...+p2*p4-(p1+p2+p3+p4)+1.

Puisque nous faisons la moyenne et que celle ci ne donne un nombre entier que pour 4 (->1) et 0 (->0) les valeurs intermédiaires nous donnent la probabilité d'être à 0 ou à 1 donc il faut multiplier P3 par 3 / 4 …etc. Donc, la probabilité que la cellule centrale soit à 1 est :

PP1=P4+3/4*P3+1/2*P2+1/4*P1.

Et à 0 la probabilité est :

PP0=P0+3/4*P1+1/2*P2+1/4*P3.

Après simplification on trouve :

PP1=(p1+p2+p3+p4)/4.

PP0=1-(p1+p2+p3+p4)/4=1-PP1.

Ca veut dire tout simplement que les probabilités suivent la même équation. La probabilité pour la cellule centrale d'avoir pour valeur 1 est la moyenne des probabilités de ses voisines d'avoir pour valeur 1. Magique !

Vous avez remarqué qu'il s'agit de variables aléatoires prenant des valeurs dans [0..1] (sinon c'est pas grave). Pour le calcul des nombres entiers sur n bits ou plus généralement dans des intervalles entiers [0..m], il faut raisonner sur des variables aléatoires prenant des valeurs sur ce même intervalle. Il va sans dire que c'est beaucoup plus compliqué. C'est pour ça que je laisse tomber.

Aller plus loin :

Nous avons vu qu'avec l'arrondi aleatoire les erreurs d'arrondi dans le cas de la resolution de l'equation de Laplace ne s'accumulent plus au dela d'un certain niveau. Ce niveau reste a definir. Neanmoins on peut poser l'hypothese que l'erreure d'arrondi suit dans notre cas une loi normale. C'est peut etre du a l'ajout continuel du terme aleatoire. Mais ce qui est surprenant, c'est que la deviation par rapport a la valeur exacte est petit (rarement superieure a 1): la solution fluctue autour de la valeur exacte et l'encadre remarquablement! Tout cela reste empirique: Qui peut demontrer! (a moins que je ne me trompe)

Nous avons expérimenté l'arrondi stochastique (c'est comme ça que je le baptise... c'est certainement un abus car le terme et le concept exite deja: J.M. Chesneaux a developpe une bibliothque pour utiliser l'arithmetique stochastique due a J.Vignes et son equipe: C'est une methode de controle des erreurs d'arrondi. Voir aussi le travail de Stott Parker sur l'arithmetique de Monte-Carlo et ce toujours pour le controle des erreurs d'arrondi et tout ce qui s'ensuit) pour la résolution de l'équation de Laplace, mais rien ne nous empêche de le faire pour d'autres et pourquoi pas généraliser son utilisation dans tous les domaines de calcul numérique.

Je l'ai essayé avec le fameux effet WATER. Cette méthode permet d'avoir un effet régulier même quand on fait le calcul sur des nombres avec peu de chiffres significatifs. Cet effet est en realite moin ensible aux erreurs d'arrondi: avec un nombre de chiffre significatifs donne, la taille du reseau peut etre pluss grande tout en ayant un effet correct (pour 7bits on a environ une grille de 100x100) mais au dela l'influence de la topologie du reseau devient de plus en plus importante. L'utilisation des nombres (stochastiques) permet d'eviter cette influence (plus de precissions prochainement).

On pourrait définir un nouveau type de nombres flottants utilisant l'arrondi stochastique (Les recherces font bon train :). En effet au cours d'un calcul avec les nombres flottants classiques il y a toujours des chiffres significatifs qui se perdent surtout lors des multiplications et des divisions et aussi lors de l'addition de nombres dont les magnitudes sont très différentes. La méthode est en principe simple : on ajoute au résultat un nombre (pseudo) aléatoire fonction de la partie éliminée lors de l'arrondi. ceci permettrait de mieux contrôler les erreurs d'arrondi et les problèmes qui s'en suivent. Je n'ai pas abordé sérieusement ce dernier point .

Avant d'en finir, énumérons les questions en suspens :

Un peu de tmah'bil :

Je ne sais pas pourquoi, mais tout ça me rappelle les trucs de la physique quantique. Ca me donne toujours une impression d'irréalité. Je n'arrive toujours pas à COMPRENDRE comment ça marche. je vous l'ai dit, c'est magique.

Et si le monde ne connaissait que les nombres entiers et le hasard? ça c'est une autre histoire...

Contact : nait_merzouk_abdelaziz@yahoo.fr

Copyright © ;) Naït-Merzouk Abdelaziz. Janvier 2001... et quelques notes depuis octobre de la meme annee (abscence d'accents vous voyez!)

 

Note: Cet "article" pourrait bien entrer dans la categorie "Plagiat par anticipation" ;)) et ce, du fait que je n'ai pu trouver aucune information sur le sujet a part l'arithmetique stochastique (merci J.M.M.) et l'arithmetique de Monte-Carlo. En effet celles-ci s'occupent de controle des erreurs d'arrondi.

1