Code C/C++ le plus rapide pour sélectionner la médiane dans un ensemble de 27 valeurs à virgule flottante

Question

Il s'agit de l'algorithme de sélection bien connu. voir http://en.wikipedia.org/wiki/Selection_algorithm .

J'ai besoin qu'il trouve la valeur médiane d'un ensemble de 3x3x3 valeurs de voxels. Comme le volume est composé d'un milliard de voxels et que l'algorithme est récursif, il a intérêt à être un peu rapide. En général, on peut s'attendre à ce que les valeurs soient relativement proches.

L'algorithme connu le plus rapide que j'ai essayé jusqu'à présent utilise la fonction de partition du tri rapide. J'aimerais savoir s'il en existe un plus rapide.

J'ai "inventé" une solution 20% plus rapide en utilisant deux tas, mais je m'attendais à une solution encore plus rapide en utilisant un hachage. Avant de l'implémenter, j'aimerais savoir si une solution rapide comme l'éclair existe déjà.

Le fait que j'utilise des flottants ne devrait pas avoir d'importance puisqu'ils peuvent être considérés comme des entiers non signés après inversion du bit de signe. L'ordre sera préservé.

EDIT : benchmark et code source déplacés dans une réponse séparée comme suggéré par Davy Landman. Voir ci-dessous pour la réponse de chmike.

EDIT : L'algorithme le plus efficace jusqu'à présent a été référencé ci-dessous par Boojum comme un lien vers la Filtrage médian et bilatéral rapide qui est maintenant la réponse à cette question. La première idée intelligente de cette méthode est d'utiliser le tri radix, la seconde est de combiner la recherche médiane des pixels adjacents qui partagent beaucoup de pixels.

Answer 1

L'algorithme de sélection est en temps linéaire (O(n)). En termes de complexité, vous ne pouvez pas faire mieux que le temps linéaire, car il faut un temps linéaire pour lire toutes les données. Vous n'auriez donc pas pu faire quelque chose qui soit plus rapide du point de vue de la complexité. Peut-être avez-vous quelque chose qui est un facteur constant plus rapide sur certaines entrées ? Je doute que cela fasse une grande différence.

Le C++ comprend déjà l'algorithme de sélection en temps linéaire. Pourquoi ne pas l'utiliser ?

std::vector<YourType>::iterator first = yourContainer.begin();
std::vector<YourType>::iterator last = yourContainer.end();
std::vector<YourType>::iterator middle = first + (last - first) / 2;
std::nth_element(first, middle, last); // can specify comparator as optional 4th arg
YourType median = *middle;

Edit : Techniquement, il s'agit seulement de la médiane pour un conteneur de longueur impaire. Pour un conteneur de longueur paire, il obtiendra la médiane "supérieure". Si vous voulez la définition traditionnelle de la médiane pour une longueur paire, vous devrez peut-être l'exécuter deux fois, une fois pour chacune des deux "médianes" à first + (last - first) / 2 y first + (last - first) / 2 - 1 et en faire la moyenne ou quelque chose comme ça.

Answer 2

EDIT : Je dois m'excuser. Le code ci-dessous était erroné. J'ai le code corrigé, mais j'ai besoin de trouver un fichier de type icc compilateur pour refaire les mesures.

Les résultats de référence des algorithmes considérés jusqu'à présent

Pour le protocole et une brève description des algorithmes, voir ci-dessous. La première valeur est le temps moyen (secondes) sur 200 séquences différentes et la deuxième valeur est l'écart-type.

HeapSort     : 2.287 0.2097
QuickSort    : 2.297 0.2713
QuickMedian1 : 0.967 0.3487
HeapMedian1  : 0.858 0.0908
NthElement   : 0.616 0.1866
QuickMedian2 : 1.178 0.4067
HeapMedian2  : 0.597 0.1050
HeapMedian3  : 0.015 0.0049 <-- best

Protocole : générer 27 flottants aléatoires en utilisant des bits aléatoires obtenus par rand(). Appliquer chaque algorithme 5 millions de fois d'affilée (y compris la copie préalable du tableau) et calculer la moyenne et le stdDev sur 200 séquences aléatoires. Code C++ compilé avec icc -S -O3 et exécuté sur Intel E8400 avec 8GB DDR3.

Algorithmes :

HeapSort : tri complet de la séquence en utilisant le tri par tas et en choisissant la valeur centrale. Implémentation naïve utilisant l'accès aux indices.

QuickSort : tri complet en place d'une séquence en utilisant le tri rapide et le choix de la valeur centrale. Implémentation naïve utilisant l'accès aux indices.

QuickMedian1 : algorithme de sélection rapide avec permutation. Implémentation naïve utilisant l'accès aux indices.

HeapMedian1 : méthode de tas équilibré en place avec permutation préalable. Implémentation naïve utilisant l'accès aux indices.

NthElement : utilise l'algorithme STL nth_element. Les données sont copiées dans le vecteur en utilisant memcpy( vct.data(), rndVal, ... ) ;

QuickMedian2 : utilise l'algorithme de sélection rapide avec des pointeurs et la copie dans deux tampons pour éviter le swaping. Basé sur la proposition de MSalters.

HeapMedian2 : variante de l'algorithme que j'ai inventé, utilisant des tas doubles avec des têtes partagées. Le tas de gauche a la plus grande valeur comme tête, le tas de droite a la plus petite valeur comme tête. Initialiser avec la première valeur comme tête commune et la première valeur médiane devinée. Ajouter les valeurs suivantes au tas de gauche si elles sont plus petites que la tête, sinon au tas de droite, jusqu'à ce que l'un des tas soit plein. Il est plein lorsqu'il contient 14 valeurs. Ensuite, on ne considère que le tas plein. S'il s'agit du tas droit, pour toutes les valeurs supérieures à la tête, on ouvre la tête et on insère la valeur. Ignorez toutes les autres valeurs. S'il s'agit du tas de gauche, pour toutes les valeurs inférieures à la tête, on ouvre la tête et on insère la valeur dans le tas. Ignorez toutes les autres valeurs. Lorsque toutes les valeurs ont été traitées, la tête commune est la valeur médiane. Il utilise un indice entier dans le tableau. La version utilisant les pointeurs (64bit) semble être presque deux fois plus lente (~1s).

HeapMedian3 : même algorithme que HeapMedian2 mais optimisé. Il utilise l'index des chars non signés, évite les échanges de valeurs et diverses autres petites choses. Les valeurs de moyenne et de stdDev sont calculées sur 1000 séquences aléatoires. Pour nth_element j'ai mesuré 0.508s et un stdDev de 0.159537 avec les mêmes 1000 séquences aléatoires. HeapMedian3 est donc 33 fois plus rapide que la fonction stl nth_element. Chaque valeur médiane retournée est comparée à la valeur médiane retournée par heapSort et elles correspondent toutes. Je doute qu'une méthode utilisant le hachage puisse être significativement plus rapide.

EDIT 1 : Cet algorithme peut être encore optimisé. La première phase où les éléments sont distribués dans le tas de gauche ou de droite en fonction du résultat de la comparaison n'a pas besoin de tas. Il suffit d'ajouter des éléments à deux séquences non ordonnées. La première phase s'arrête dès qu'une séquence est pleine, c'est-à-dire qu'elle contient 14 éléments (y compris la valeur médiane). La deuxième phase commence par la mise en tas de la séquence complète, puis procède comme décrit dans l'algorithme HeapMedian3. Je fournirai le nouveau code et le benchmark dès que possible.

EDIT 2 : J'ai implémenté et évalué l'algorithme optimisé. Mais il n'y a pas de différence de performance significative par rapport à heapMedian3. Il est même légèrement plus lent en moyenne. Les résultats affichés sont confirmés. Il pourrait y en avoir avec des ensembles beaucoup plus grands. Notez également que je choisis simplement la première valeur comme estimation initiale de la médiane. Comme suggéré, on pourrait profiter du fait que nous recherchons une valeur médiane dans des ensembles de valeurs "chevauchantes". L'utilisation de l'algorithme de la médiane de la médiane permettrait de choisir une bien meilleure estimation de la valeur médiane initiale.

---

Code source de HeapMedian3

// return the median value in a vector of 27 floats pointed to by a
float heapMedian3( float *a )
{
   float left[14], right[14], median, *p;
   unsigned char nLeft, nRight;

   // pick first value as median candidate
   p = a;
   median = *p++;
   nLeft = nRight = 1;

   for(;;)
   {
       // get next value
       float val = *p++;

       // if value is smaller than median, append to left heap
       if( val < median )
       {
           // move biggest value to the heap top
           unsigned char child = nLeft++, parent = (child - 1) / 2;
           while( parent && val > left[parent] )
           {
               left[child] = left[parent];
               child = parent;
               parent = (parent - 1) / 2;
           }
           left[child] = val;

           // if left heap is full
           if( nLeft == 14 )
           {
               // for each remaining value
               for( unsigned char nVal = 27 - (p - a); nVal; --nVal )
               {
                   // get next value
                   val = *p++;

                   // if value is to be inserted in the left heap
                   if( val < median )
                   {
                       child = left[2] > left[1] ? 2 : 1;
                       if( val >= left[child] )
                           median = val;
                       else
                       {
                           median = left[child];
                           parent = child;
                           child = parent*2 + 1;
                           while( child < 14 )
                           {
                               if( child < 13 && left[child+1] > left[child] )
                                   ++child;
                               if( val >= left[child] )
                                   break;
                               left[parent] = left[child];
                               parent = child;
                               child = parent*2 + 1;
                           }
                           left[parent] = val;
                       }
                   }
               }
               return median;
           }
       }

       // else append to right heap
       else
       {
           // move smallest value to the heap top
           unsigned char child = nRight++, parent = (child - 1) / 2;
           while( parent && val < right[parent] )
           {
               right[child] = right[parent];
               child = parent;
               parent = (parent - 1) / 2;
           }
           right[child] = val;

           // if right heap is full
           if( nRight == 14 )
           {
               // for each remaining value
               for( unsigned char nVal = 27 - (p - a); nVal; --nVal )
               {
                   // get next value
                   val = *p++;

                   // if value is to be inserted in the right heap
                   if( val > median )
                   {
                       child = right[2] < right[1] ? 2 : 1;
                       if( val <= right[child] )
                           median = val;
                       else
                       {
                           median = right[child];
                           parent = child;
                           child = parent*2 + 1;
                           while( child < 14 )
                           {
                               if( child < 13 && right[child+1] < right[child] )
                                   ++child;
                               if( val <= right[child] )
                                   break;
                               right[parent] = right[child];
                               parent = child;
                               child = parent*2 + 1;
                           }
                           right[parent] = val;
                       }
                   }
               }
               return median;
           }
       }
   }
}

Answer 3

Puisqu'il semble que vous exécutiez un filtre médian sur un grand tableau de données de volume, vous pourriez vouloir jeter un coup d'œil à l'outil d'analyse des données de volume. Filtrage médian et bilatéral rapide Document du SIGGRAPH 2006. Cet article traite du traitement d'images en 2D, mais vous pourrez peut-être adapter l'algorithme aux volumes en 3D. Cela pourrait vous donner des idées sur la façon de prendre du recul et de considérer le problème d'un point de vue légèrement différent.

Answer 4

Il n'est pas facile de répondre à cette question pour la simple raison que les performances d'un algorithme par rapport à un autre dépendent autant de la combinaison compilateur / processeur / structure de données que de l'algorithme lui-même, comme vous le savez certainement.

Par conséquent, votre approche, qui consiste à en essayer quelques-uns, semble suffisante. Et oui, quicksort devrait être assez rapide. Si vous ne l'avez pas encore fait, vous pourriez essayer insertionsort qui est souvent plus performant sur les petits ensembles de données. Ceci dit, choisissez un algo de tri qui fait le travail assez rapidement. Vous n'obtiendrez généralement pas une vitesse 10 fois supérieure en choisissant le "bon" algo.

Pour obtenir des gains de vitesse substantiels, la meilleure solution consiste souvent à utiliser davantage de structure. Quelques idées qui ont fonctionné pour moi dans le passé avec des problèmes à grande échelle :

• Pouvez-vous effectuer un précalcul efficace lors de la création des voxels et stocker 28 flottants au lieu de 27 ?

• Une solution approximative est-elle suffisante ? Si Si c'est le cas, il suffit de regarder la médiane de, disons 9 valeurs, puisque "en général, on peut s'attendre à ce que les valeurs soient relativement proches". Ou vous pouvez le remplacer par la moyenne, tant que les valeurs sont relativement proches.

• Avez-vous vraiment besoin de la médiane pour tous les milliards de voxels ? Vous avez peut-être un test facile pour savoir si vous avez besoin de la médiane, et vous pouvez alors calculer uniquement pour le sous-ensemble concerné.

• Si rien d'autre ne vous aide : regardez le code code asm que le compilateur génère. Vous pouvez peut-être écrire du code asm qui est sensiblement plus rapide (par exemple en effectuant tous les calculs en utilisant des registres).

Edit : Pour ce que cela vaut, j'ai joint le code (partiel) d'insertionsort mentionné dans le commentaire ci-dessous (totalement non testé). Si numbers[] est un tableau de taille N et vous voulez le plus petit P floats triés au début du tableau, appelez partial_insertionsort<N, P, float>(numbers); . Ainsi, si vous appelez partial_insertionsort<27, 13, float>(numbers); , numbers[13] contiendra la médiane. Pour gagner en rapidité, vous devriez également déplier la boucle while. Comme nous l'avons vu plus haut, pour être vraiment rapide, vous devez utiliser vos connaissances sur les données (par exemple, les données sont-elles déjà partiellement triées ? Connaissez-vous les propriétés de la distribution des données ? Je suppose que vous voyez ce que je veux dire).

template <long i> class Tag{};

template<long i, long N, long P, typename T>
inline void partial_insertionsort_for(T a[], Tag<N>, Tag<i>)
{   long j = i <= P+1 ? i : P+1;  // partial sort
    T temp = a[i];
    a[i] = a[j];       // compiler should optimize this away where possible
    while(temp < a[j - 1] && j > 0)
    { a[j] = a[j - 1];
      j--;}
    a[j] = temp;
    partial_insertionsort_for<i+1,N,P,T>(a,Tag<N>(),Tag<i+1>());}

template<long i, long N, long P, typename T>
inline void partial_insertionsort_for(T a[], Tag<N>, Tag<N>){}

template <long N, long P, typename T>
inline void partial_insertionsort(T a[])
 {partial_insertionsort_for<0,N,P,T>(a, Tag<N>(), Tag<0>());}

Answer 5

L'algorithme le plus susceptible d'être utilisé dans votre première tentative est simplement nth_element ; il vous donne à peu près directement ce que vous voulez. Demandez simplement le 14ème élément.

Lors de votre deuxième tentative, l'objectif est de tirer parti de la taille fixe des données. Vous ne voulez pas allouer de mémoire du tout pour votre algorithme. Donc, copiez vos valeurs de voxels dans un tableau pré-alloué de 27 éléments. Choisissez un pivot, et copiez-le au milieu d'un tableau de 53 éléments. Copiez les valeurs restantes de chaque côté du pivot. Ici, vous gardez deux pointeurs ( float* left = base+25, *right=base+27 ). Il y a maintenant trois possibilités : le côté gauche est plus grand, le côté droit est plus grand, ou les deux ont 12 éléments. Le dernier cas est trivial ; votre pivot est la médiane. Sinon, appelez nth_element sur le côté gauche ou le côté droit. La valeur exacte de Nth dépend du nombre de valeurs plus grandes ou plus petites que le pivot. Par exemple, si la division est 12/14, vous avez besoin du plus petit élément plus grand que le pivot, donc Nth=0, et si la division est 14/12, vous avez besoin du plus grand élément plus petit que le pivot, donc Nth=13. Les pires cas sont 26/0 et 0/26, lorsque votre pivot était un extrême, mais ils ne se produisent que dans 2/27e de tous les cas.

La troisième amélioration (ou la première, si vous devez utiliser C et ne disposez pas de nth_element) remplace entièrement nth_element. Vous avez toujours le tableau de 53 éléments, mais cette fois vous le remplissez directement à partir des valeurs des voxels (ce qui vous évite une copie intermédiaire dans un fichier float[27] ). Le pivot dans cette première itération est simplement le voxel[0][0][0]. Pour les itérations suivantes, vous utilisez un deuxième pivot pré-alloué, le voxel[0][0]. float[53] (plus facile si les deux sont de la même taille) et copier les flottants entre les deux. L'étape d'itération de base ici est toujours : copier le pivot au milieu, trier le reste à gauche et à droite. À la fin de chaque étape, vous saurez si la médiane est plus petite ou plus grande que le pivot actuel, ce qui vous permettra de rejeter les flottants plus grands ou plus petits que ce pivot. Par itération, cela élimine entre 1 et 12 éléments, avec une moyenne de 25% des éléments restants.

La dernière itération, si vous avez encore besoin de plus de vitesse, est basée sur l'observation que la plupart de vos voxels se chevauchent de manière significative. Vous pré-calculez pour chaque tranche de 3x3x1 la valeur médiane. Ensuite, lorsque vous avez besoin d'un pivot initial pour votre cube de voxels 3x3x3, vous prenez la médiane des trois. Vous savez a priori qu'il y a 9 voxels plus petits et 9 voxels plus grands que cette médiane des médianes (4+4+1). Ainsi, après la première étape de pivotement, les pires cas sont une répartition 9/17 et 17/9. Il suffirait donc de trouver le 4e ou le 13e élément dans un float[17], au lieu du 12e ou du 14e dans un float[26].

---

Contexte : L'idée de copier d'abord un pivot puis le reste d'un float[N] dans un float[2N-1], en utilisant des pointeurs gauche et droit, est que vous remplissez un subarray float[N] autour du pivot, avec tous les éléments plus petits que le pivot à gauche (indice inférieur) et plus hauts à droite (indice supérieur). Maintenant, si vous voulez le Mième élément, vous pouvez avoir de la chance et avoir M-1 éléments plus petits que le pivot, dans ce cas le pivot est l'élément dont vous avez besoin. S'il y a plus de (M-1) éléments plus petits que le pivot, le Mème élément est parmi eux, donc vous pouvez écarter le pivot et tout ce qui est plus grand que le pivot, et chercher le Mème élément dans toutes les valeurs inférieures. S'il y a moins que (M-1) éléments plus petits que le pivot, vous cherchez une valeur supérieure au pivot. Donc, on élimine le pivot et tout ce qui est plus petit que lui. Soit le nombre d'éléments moins que le pivot, c'est-à-dire à gauche du pivot, soit L. Dans l'itération suivante, vous voulez le (M-L-1)ème élément des (N-L-1)flottants qui sont plus grands que le pivot.

Ce type d'algorithme nth_element est assez efficace car la plupart du travail est consacré à la copie de flottants entre deux petits tableaux, qui seront tous deux dans le cache, et parce que votre état est la plupart du temps représenté par 3 pointeurs (pointeur source, pointeur destination gauche, pointeur destination droit).

Pour montrer le code de base :

float in[27], out[53];
float pivot = out[26] = in[0];     // pivot
float* left = out+25, right = out+27
for(int i = 1; i != 27; ++1)
if((in[i]<pivot)) *left-- = in[i] else *right++ = in[i];
// Post-condition: The range (left+1, right) is initialized.
// There are 25-(left-out) floats <pivot and (right-out)-27 floats >pivot