Re: essayer Julia ?

Tu détournes habilement le sujet, mais je pense que tu as tort sur le fond.

Je ne pense pas détourner le sujet, et je pense fondamentalement avoir raison sur le fond. ;-)

quand on code dans la plupart des langages (qui sont destinés à produire des logiciels et pas des démonstrations mathématiques), on doit avant tout penser aux aspects informatiques des problèmes (efficacité de l'algorithme en temps et en mémoire, gestion des arrondis, évolutivité, modularité et clarté du code). Dans ce cadre, les maths sont un outil, et pas un état d'esprit.

Et les concepts que tu mets en branle, dans ton esprit, quand tu penses aux aspects informatiques des problèmes, ils relèvent de quelle science à ton avis ? ;-) Je t'ai donné l'exemple d'une équipe membre de l'Institut de Recherche en Informatique Fondamentale qui s'associe avec une autre équipe de l'Institut Mathématique de Jussieu pour dispenser une formation intituler Logique Mathématique et Fondements de l'Informatique, et tu ne vois toujours pas le rapport ?

Tu me fais penser à M. Jourdain : il faisait de la prose sans le savoir mais, quand on lui a expliqué ce qu'était la prose, il a au moins reconnu qu'il en faisait. Toi c'est un peu différent, tu fais des mathématiques sans le savoir (sans doute par ce que tu ignores ce que sont les mathématiques et que tu n'en reconnais pas toujours quand tu en vois), je t'expliques qu'en réalité tu en fais quand tu programmes, mais tu restes dans le déni et prétends que tu n'en fais pas.

Je vais le dire autrement avec le théorème de Pythagore. Voilà un théorème qui dit : donne moi un triangle, je te construirais trois carrés dont la surface de l'un et la somme de la surface des autres. Autrement dit c'est une fonction qui prend en entrée un triangle et retourne un triplet de carré. Alors assurément, comme tout théorème, il a plus d'une démonstration mais elles font toutes la même chose. Cela étant, dans toutes ces démonstrations, il y en a qui sont plus efficaces que d'autres pour produire la sortie. C'est pareil pour toutes les fonctions que tu codes : ce sont des preuves de théorèmes mais certaines sont plus efficaces que d'autres. Que tu l'ignores ou que tu ne le vois pas, c'est une chose; que ce soit faux, s'en est une autre. ;-)

Pour revenir au débat d'origine avec aurelienpierre :

Mais je ne sais pas si tu as fait semblant de mal comprendre là où je voulais en venir, ou si je ne m'étais pas exprimé clairement. La question qu'on discutait, c'était de dire que la vectorisation était une manière intuitive en mathématique d'aborder un problème, et que les gens qui considéraient une boucle FOR plus intuitive qu'un calcul vectorisé avaient, en gros, un problème de formation.

Je ne sais pas trop ce qu'il faut entendre dans votre discussion par le terme vectorisation. S'agit-il des instructions SIMD des CPU ou de manipuler des structures de données abstraites représentant le concept mathématique de vecteur que l'on trouve en algèbre linéaire ?

Pour ce qui est des idiomes des langages, en python on utilisera volontiers des itérateurs plutôt que des boucles FOR (en C++ aussi, il me semble qu'il y a des itérateurs dans la STL). Le produit scalaire entre deux vecteurs se définira ainsi :

import operator
def dotprod(vec1, vec2):
 return sum(map(operator.mul, vec1, vec2))

et non avec une boucle FOR. Il me semble que c'était déjà, là, une des choses que voulait faire remarquer aurelienpierre. Dans un langage comme le C, assurément on fera la même chose avec une boucle FOR mais parce c'est là l'idiome du langage pour faire ce genre de calcul.

Revenons au calcul du produit matriciel et à la quette d'optimisation. En C, la traduction naïve de la chose donnerait :

for (i = 0; i < N; ++i)
 for (j = 0; j < N; ++j)
 for (k = 0; k < N; ++k)
 res[i][j] += mul1[i][k] * mul2[k][j];

Ici comme on parcourt la deuxième matrice colonne par colonne, sur de grosses matrices on a du cache miss. En la transposant d'abord on a :

double tmp[N][N];
 for (i = 0; i < N; ++i)
 for (j = 0; j < N; ++j)
 tmp[i][j] = mul2[j][i];
 for (i = 0; i < N; ++i)
 for (j = 0; j < N; ++j)
 for (k = 0; k < N; ++k)
 res[i][j] += mul1[i][k] * tmp[j][k];

Les exemples de code sont issus de Memory part 5: What programmers can do, une série d'articles sur LWN par Ulrich Drepper au sujet du fonctionnement de la mémoire et des caches. Rien que là, dans ses benchmarks, il a un gain de 76.6%.

Néanmoins, il faut allouer une matrice temporaire : c'est lourd et on a pas toujours l'envie ni la place de faire. Il propose alors mieux :

#define SM (CLS / sizeof (double))
 for (i = 0; i < N; i += SM)
 for (j = 0; j < N; j += SM)
 for (k = 0; k < N; k += SM)
 for (i2 = 0, rres = &res[i][j],
 rmul1 = &mul1[i][k]; i2 < SM;
 ++i2, rres += N, rmul1 += N)
 for (k2 = 0, rmul2 = &mul2[k][j];
 k2 < SM; ++k2, rmul2 += N)
 for (j2 = 0; j2 < SM; ++j2)
 rres[j2] += rmul1[k2] * rmul2[j2];

et il compile le code avec gcc -DCLS=$(getconf LEVEL1_DCACHE_LINESIZE) pour optimiser le code pour la machine sur lequel il est compilé : CLS représente la taille d'une ligne de cache de niveau 1 sur la machine. Et là ce qu'il fait, avec des boucles FOR parce que tel est l'idome du C, c'est suivre la courbe en Z de Lebesgue (cf. mon premier commentaire) en adaptant la taille des zigzag à celui de la ligne de cache.

Il évite ainsi d'allouer une matrice temporaire pour calculer la transposer et il gagne 6.1% de plus qu'avec le code précédent. Mais au fond ce qu'il vient d'écrire ce n'est que la traduction dans le langage formel qu'est le C d'une pensée qui est mathématique de part en part.

Il conclue, enfin, en disant que l'on peut aller encore plus loin en vectorisant (instruction SIMD) le code et gagner encore 7.3%. À l'arrivée, il a un code qui va 10% plus vite que la boucle FOR naïve.

Ceci étant, les compilateurs appliquent déjà des optimisations de ce genre (pas forcément sur ces problèmes, mais sur d'autres) mais pour ce faire leurs auteurs, eux, connaissent l'outillage conceptuel mathématique nécessaire et il vaut mieux les laisser faire que d'essayer de le faire soi-même (ce que tu as reconnu ;-).

Encore un autre exemple, si tu n'est toujours pas convaincu. Voici une liste chaînée :

: -> : -> : -> : -> []
| | | |
1 2 3 4

elle a sa petite sœur, la liste doublement chaînée :

] <- : <- : <- FOCUS -> : -> []
 | | | |
 1 2 3 4

en programmation fonctionnelle on appelle cela le zipper sur une liste. Et bien le zipper (ou liste doublement chaînée) et le type dérivé du type des listes chaînées. Explication ici : The algebra (and calculus!) of algebraic data types, on tu verras du développement en série entières et du calcul différentiel sur des types. ;-)

Sapere aude ! Aie le courage de te servir de ton propre entendement. Voilà la devise des Lumières.