Discussion :

[ToTo13]

Bonjour,

je suis amené à faire des calculs de statistiques basiques, telles que moyenne, variance, skewness, kurtosis, mais sur des TRES grandes images, donc sur des TRES grands nombre de valeurs.

Bien entendu, j'atteins immédiatement les limites dues au codage : par exemple, pour calculer la moyenne, il faut théoriquement additionner tous les nombres, puis diviser par le nombre de nombres additionnés.
Le souci c'est qu'en court de route, la somme déborde de la capacité de codage (dans mon cas celle de java) et le calcul est faussé.

Une solution serait de calculer la moyenne à la volée, c'est-à-dire additionner les nombres immédiatement divisés par le nombre total de nombres. Je sais qu'en C, le codage des double offre autant de valeurs entre -1 et 1 que dans tous le reste de la plage de codage, donc le calcul peut très très bien passer comme cela.

Une autre solution, pourrait être de faire cela par bloc de X nombres.

Est ce que quelqu'un connaîtrait LA bonne solution pour traiter ce problème ?


Merci par avance.

[souviron34]

je suis amené à faire des calculs de statistiques basiques, telles que moyenne, variance, skewness, kurtosis, mais sur des TRES grandes images, donc sur des TRES grands nombre de valeurs.

De quels nombres parle-t-on ici ?

Le souci c'est qu'en court de route, la somme déborde de la capacité de codage (dans mon cas celle de java) et le calcul est faussé.

Tu dépasses la valeur des double-précison ??

C'est pas dans l'algo en tant que tel ? (l'image ou son exploration ?)

Une solution serait de calculer la moyenne à la volée, c'est-à-dire additionner les nombres immédiatement divisés par le nombre total de nombres.

Je ne suis pas certain qu'on obtienne le même résultat...

1 + 10 + 2 + 20 => 8.25
1+10 => 5.25
5.25 + 2 => 3.625
3.625 + 20 => 11.8125

[gpcbitnik38]

comme tu le dit faire ça en x fois me semble le plus simple

tu cherche ta valeur max ( peu etre que tu la connait a l'avance meme)
tu découte ton image en n zone de tel sort que meme si une zone comptenait que le max la somme ne deborderait pas.
tu calcule tes moyenne sur chaque zone puis ta moyenne global en faisant la moyenne des moyennes ( avec toute un coef de 1 puisque les zones on le meme nombre de pixels)
l'inconveniant est quetu a plus de variable :s


tu peux envisager un découpage "optimal" en réalisant la somme d'autant de pixel que nécéssaire jusqu'au debordement, au moment ou tu deborde tu stoque le nombre de pixel deja compter tu fait ta moyenne de ta zone et tu passe a la suivante.
a la fin tu fait ta moyenne générale en oubliant pas de pondérer chaque moyenne par sont nombre de pixel sur le nombre de pixel total.


après je sais pas si c'est adapté pour les moment d'ordre superieur a 2 je me souvient plus de leur definition


faire la moyenne "à la volé" risque de te prendre un temps fou (2 divisions une multiplication et une somme a chaque etape ) et je suis pas sur que cela règle ton probleme.
si m(n) est ta moyenne n ton nombre d'élement deja pris en compte et x l'element courant
m(n+1)=[n/(n+1)]*m(n)+x/(n+1) ; je sais pas exactement si c'est de ça dont tu parlais mais si c'est le cas tu risque d'avoir exactement le même probleme sur la précision de n/(n+1) et de 1/(n+1)

[Alexis.M]

Je ne sais pas si c'est LA bonne solution mais l'approche par bloc me parait raisonnable.

Avec l'approche où tu divises tout le monde d'abord par le nombre total d'éléments:
<x> = \sum_i x_i / N
tu risques de propager des erreurs d'arrondis si N est très grand par rapport a x_i.

L'approche par bloc est en quelque sorte intermédiaire:
<x> = \sum_i (S_i / N)
où S_i est la somme calculée pour le bloc i. S_i peut alors être raisonnablement grand.

Le gros avantage de cette méthode est surtout qu'elle parallélisable, tu peux calculer les sommes en parallèle et ensuite additionner les résultats (du map-reduce en quelque sort).

Edit:
Une petite expérience avec numpy:

Code python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
 
>>> import numpy as np
>>> X = np.ones((100000,1000))
>>> print (X/X.size).sum() - 1.0 # somme des élements individuels
2.28986718476e-09
>>> print (X.sum(1)/X.size).sum() - 1.0 # blocks de taille 1000
-1.9162449405e-12
>>> print (X.sum(0)/X.size).sum() - 1.0 # blocks de taille 100000
6.66133814775e-16

[souviron34]

je suis d'accord, mais j'ai quand même du mal à voir les limites du calcul normal..

valeur max d'un double-précision sur 32 bits 10^308

ça dépasse ça ??????????

[Alexis.M]

je suis d'accord, mais j'ai quand même du mal à voir les limites du calcul normal..

valeur max d'un double-précision sur 32 bits 10^308

ça dépasse ça ??????????

Je ne suis pas sûr que ce soit ce que voulais dire Toto13, mais ce qui est sûr c'est que si tu fais:

Code python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
>>> 1e16 + 1 == 1e16
True

Et avec de très grandes images (certaines acquisition par microscopes peuvent atteindre plusieurs Giga pixels), tu dois pouvoir arriver facilement à ce type de cas.

[souviron34]

Je ne suis pas sûr que ce soit ce que voulais dire Toto13, mais ce qui est sûr c'est que si tu fais:

Code python :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
>>> 1e16 + 1 == 1e16
True

Et avec de très grandes images (certaines acquisition par microscopes peuvent atteindre plusieurs Giga pixels), tu dois pouvoir arriver facilement à ce type de cas.

ben si ton image est sur 8 bits (255 max d'une valeur) il te faut 10^14 points pour y arriver, soit une image de 10^7 pixels de côté

(on est au dessus du téra-pixel, très très largement au dessus du Giga pixel)

[ToTo13]

je suis d'accord, mais j'ai quand même du mal à voir les limites du calcul normal..

valeur max d'un double-précision sur 32 bits 10^308

ça dépasse ça ??????????

Ben oui :-(
Et du coup en Java la valeur repars à Double.MIN_VALUE.
Ca arrive facilement pour les Kurtosis et autres Skewness.

Je suis d'accord que faire la division à la volée est très coûteux en temps de calculs.
Comme je l'avais mentionné, il y a un intermédiaire avec la méthode par bloc et il semblerait que les solutions proposées tournent autour de cela.

Mais il persiste toujours des soucis (sans doute irréductibles):
- 1 - Travailler par bloc réduit considérablement les risques de débordement, mais ne les supprime pas.
- 2 - On peut tester si on risque le débordement, mais cela nécessite un test par itération/nombre, donc on retombe dans des problèmes de performances.

[souviron34]

attends, il y a un truc que je comprend pas trop..

Tu travailles sur des images, c'est bien ça ?

Quelle profondeur ? 8 16, 24 bits ? (en gris, en R,G,B (3 canaux)) ?? (ta valeur max et min pour un pixel c'est quoi . -127/+128 ? -32767/+32768 ?)

Tes images font quelles dimensions ? parce que j'ai vraiment du mal à voir comment pour la somme ou la varaince tu peux dépasser 10^308....

Un moyen simple de réduire la dynamique de la somme est de diviser (en réel) la valeur du pixel par la plage, et de n'avoir du coup forcément que des valeurs entre 0 et 1 (ou -1 et 1)..

Puis tu peux redresser une fois la moyenne/variance/.... faite en remultipliant par la plage...

Exemple : 8 bits

0 <= pixel <= 255

image originale => image (tableau de réels) = image ori / 255.0

calcul de la moyenne : entre 0 et 1

Vraie valeur = moyenne * 255.0

[ToTo13]

C'est effectivement une possibilité, merci.

[souviron34]

en fait tu peux même le faire en dynamique, sans passer par l'image intermédiaire

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
 
s = 0
 
pour tout i
  pour tout j
     s = s +pix(ij) / plage
  fin pour
fin pour
 
s = s * plage / (N*M)

[ToTo13]

J'ai bien pensé à cela.
La seule chose qui me dérange dans la méthode (même si elle résout le problème) c'est qu'elle introduit N divisions supplémentaire, (ce qui nuit sérieusement aux performances).

[souviron34]

pas vraiment

il suffit de calculer la division avant, et donc de ne faire qu'une multiplication (ce qui va beaucoup plus vite)

[souviron34]

et tu peux accélérer en remplaçant

pix[i,j] par image++

c'est à dire en remplaçant l'exploration "indicielle" par l'exploration par pointeur

tu économises 2 multiplications et 2 additions par point, donc en fait tu vas 2 fois plus vite que juste référencer le point sans riien lui faire

[Aleph69]

Salut,

peut-être pourrais-tu essayer de diviser ton vecteur par une puissance de 2 de manière à te ramener à des nombres entre 0 et 1 sans faire aucune erreur d'arrondi. La capacité limitée de ta mémoire vive devrait suffire à assurer que la somme des nombres obtenus, qui est inférieure au nombre d'éléments, ne dépasse pas la valeur maximale de ton type de données. Ainsi, tu peux calculer ta moyenne comme d'habitude puis multiplier pas l'inverse de ta puissance de 2 pour retrouver le résultat initial souhaité.