notacion - El código más rápido C/C++ para seleccionar la mediana en un conjunto de 27 valores de coma flotante

double c++ ejemplo (15)

+1 para todos los que mencionaron nth_element, pero este tipo de código es donde el algoritmo escrito a mano es mejor que STL porque desea generar el código más eficiente para ese compilador que se ejecuta en la única CPU con un conjunto de datos específico. Por ejemplo, para algunas CPU / compilador, la combinación std :: swap (int, int) puede ser más lenta que el swap escrito a mano usando XOR (antes de responder, sé que esto probablemente sea cierto hace 20 años, pero ya no). A veces, el rendimiento se obtiene escribiendo a mano código de ensamblaje específico para su CPU. Si planea aprovechar los procesadores de flujo de la GPU, es posible que deba diseñar su algoritmo en consecuencia.

Usted mencionó usar 2 montones y hacer un seguimiento de la mediana a medida que inserta. Eso es lo que hice hace un tiempo en un proyecto. Cambié la matriz en el lugar y usé solo un montón. No se me ocurrió ningún algoritmo más rápido, pero me gustaría advertirte sobre el uso de la memoria, específicamente la memoria caché de la CPU. Desea tener cuidado con el acceso a la memoria. La memoria caché de la CPU se intercambia por página, por lo que desea que su algoritmo toque la memoria para minimizar la caché de la CPU.

Este es el algoritmo de selección conocido. ver http://en.wikipedia.org/wiki/Selection_algorithm .

Lo necesito para encontrar el valor medio de un conjunto de valores de vóxel de 3x3x3. Como el volumen está compuesto por mil millones de vóxeles y el algoritmo es recursivo, será mejor que sea un poco rápido. En general, se puede esperar que los valores sean relativamente cercanos.

El algoritmo conocido más rápido que he probado hasta ahora usa la función de partición de clasificación rápida. Me gustaría saber si hay uno más rápido.

He "inventado" un 20% más rápido usando dos montones, pero esperaba uno aún más rápido usando un hash. Antes de implementar esto, me gustaría saber si ya existe una solución rápida.

El hecho de que estoy usando flotadores no debería importar, ya que pueden considerarse como un entero sin signo después de invertir el bit de signo. La orden será preservada.

EDITAR: el punto de referencia y el código fuente se movieron a una respuesta separada según lo sugerido por Davy Landman. Vea a continuación la respuesta de chmike.

EDITAR : Boojum hace referencia al algoritmo más eficiente hasta el momento como un enlace al documento de filtrado Fast Median y Bilateral , que ahora es la respuesta a esta pregunta. La primera idea inteligente de este método es utilizar radix sort, la segunda es combinar la búsqueda mediana de píxeles adyacentes que comparten una gran cantidad de píxeles.

Al tener, digamos, un millón de valores diferentes de los que necesita la mediana. ¿Es posible basar su mediana en un subconjunto de esos millones, digamos 10%. Entonces, ¿la mediana está cerca del n-ésimo elemento que divide los valores en 2 subconjuntos iguales (o casi iguales)? Por lo tanto, para encontrar la mediana necesitará menos de O (n) -tiempo (en este caso O (1 / 10n) y de esta manera acercarse a la clasificación óptima con quicksort en O (nlogn)?

Dado que parece que está realizando un filtro de mediana en una gran variedad de datos de volumen, le recomendamos que consulte el documento de filtrado rápido y mediano BANAL de SIGGRAPH 2006. Ese documento trata sobre el procesamiento de imágenes en 2D, pero podría estarlo. capaz de adaptar el algoritmo para volúmenes 3D. Si nada más, puede darte algunas ideas sobre cómo dar un paso atrás y ver el problema desde una perspectiva ligeramente diferente.

El algoritmo de selección es el tiempo lineal (O (n)). En lo que respecta a la complejidad, no se puede hacer mejor que el tiempo lineal, porque lleva tiempo lineal leer en todos los datos. Entonces no podrías haber hecho algo que sea más rápido en cuanto a la complejidad. ¿Quizás tienes algo que es un factor constante más rápido en ciertas entradas? Dudo que haría una gran diferencia.

C ++ ya incluye el algoritmo de selección de tiempo lineal. ¿Por qué no solo usarlo?

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;

Editar: Técnicamente, esa es solo la mediana de un contenedor de longitud impar. Para uno de igual longitud, obtendrá la mediana "superior". Si desea la definición tradicional de mediana para la longitud par, es posible que deba ejecutarla dos veces, una para cada uno de los dos "medios" al first + (last - first) / 2 y first + (last - first) / 2 - 1 y luego promediarlos o algo así.

El algoritmo más probable para usar en su primer intento es solo nth_element; más o menos te da lo que quieres directamente. Solo pide el 14 ° elemento.

En su segundo intento, el objetivo es aprovechar el tamaño de datos fijo. No querrás asignar ningún tipo de memoria a tu algoritmo. Por lo tanto, copie sus valores de vóxel a una matriz preasignada de 27 elementos. Elija un pivote y cópielo en el medio de una matriz de 53 elementos. Copie los valores restantes en cualquier lado del pivote. Aquí tienes dos punteros ( float* left = base+25, *right=base+27 ). Ahora hay tres posibilidades: el lado izquierdo es más grande, el lado derecho es más grande o ambos tienen 12 elementos. El último caso es trivial; tu pivote es la mediana De lo contrario, llame a nth_element en el lado izquierdo o el derecho. El valor exacto de Nth depende de cuántos valores eran más grandes o más pequeños que el pivote. Por ejemplo, si la división es 12/14, necesitas el elemento más pequeño que el pivote, así que Nth = 0, y si la división fue 14/12, necesitas el elemento más grande para el pivote, así que Nth = 13. Los peores casos son 26/0 y 0/26, cuando su pivote fue extremo, pero eso ocurre solo en 2/27 de todos los casos.

La tercera mejora (o la primera, si tiene que usar C y no tiene nth_element) reemplaza nth_element por completo. Aún tiene la matriz de 53 elementos, pero esta vez la completa directamente de los valores de vóxel (guardando una copia provisional en una float[27] ). El pivote en esta primera iteración es solo voxel [0] [0] [0]. Para iteraciones posteriores, utiliza un segundo float[53] preasignado float[53] (más fácil si ambos son del mismo tamaño) y copia floats entre los dos. El paso de iteración básico aquí sigue siendo: copie el pivote en el centro, clasifique el resto a la izquierda y a la derecha. Al final de cada paso, sabrá si la mediana es más pequeña o más grande que el pivote actual, por lo que puede descartar los flotantes más grandes o más pequeños que ese pivote. Por iteración, esto elimina entre 1 y 12 elementos, con un promedio del 25% del resto.

La iteración final, si aún necesita más velocidad, se basa en la observación de que la mayoría de sus vóxeles se superponen significativamente. Precalcula para cada sector 3x3x1 el valor mediano. Luego, cuando necesitas un pivote inicial para tu cubo de vóxel 3x3x3, tomas la mediana de los tres. Sabes a priori que hay 9 vóxeles más pequeños y 9 vóxeles más grandes que esa mediana de medianas (4 + 4 + 1). Entonces, después del primer paso de pivote, los peores casos son una separación de 9/17 y 17/9. Entonces, solo necesitarías encontrar el 4 ° o el 13 ° elemento en un flotador [17], en lugar del 12 ° o 14 ° en un flotador [26].

Trasfondo: la idea de copiar primero un pivote y luego el resto de un flotador [N] en un flotador [2N-1], usando los punteros izquierdo y derecho es que llene un submatriz flotante [N] alrededor del pivote, con todos los elementos más pequeño que el pivote a la izquierda (índice más bajo) y más alto a la derecha (índice más alto). Ahora, si quiere el elemento Mth, puede que tenga suerte y tenga elementos M-1 más pequeños que el pivote, en cuyo caso el pivote es el elemento que necesita. Si hay más de (M-1) elementos más pequeños que el pivote, el M-ésimo elemento está entre ellos, por lo que puedes descartar el pivote y cualquier otro elemento más grande que el pivote, y buscar el elemento Mth en todos los valores inferiores. Si hay menos de (M-1) elementos más pequeños que el pivote, está buscando un valor más alto que el pivote. Entonces, descartarás el pivote y cualquier cosa más pequeña que él. Deje que el número de elementos sea menor que el pivote, es decir, a la izquierda del pivote sea L. En la siguiente iteración, desea el elemento (ML-1) th de los flotantes (NL-1) que son más grandes que el pivote.

Este tipo de algoritmo de nth_element es bastante eficiente porque la mayor parte del trabajo se realiza copiando flotantes entre dos matrices pequeñas, ambas estarán en caché, y porque su estado está representado la mayor parte del tiempo por 3 punteros (puntero de origen, puntero de destino izquierdo , puntero derecho de destino).

Para mostrar el código básico:

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

El nuevo libro de Alex Stepanov, Elementos de la programación, trata un poco sobre cómo encontrar estadísticas de pedidos utilizando el número mínimo de comparaciones promedio, al tiempo que se minimiza la sobrecarga del tiempo de ejecución. Desafortunadamente, se necesita una cantidad considerable de código solo para calcular la mediana de 5 elementos, e incluso entonces da como proyecto la búsqueda de una solución alternativa que utiliza una fracción de una comparación menos en promedio, así que no soñaría con extender esa marco para encontrar la mediana de 27 elementos. Y el libro ni siquiera estará disponible hasta el 15 de junio de 2009. El punto es que, debido a que este es un problema de tamaño fijo, existe un método de comparación directa que es probablemente óptimo.

Además, está el hecho de que este algoritmo no se ejecuta una vez de forma aislada sino varias veces, y entre la mayoría de las ejecuciones, solo cambiarán 9 de los 27 valores. Eso significa que en teoría parte del trabajo ya está hecho. Sin embargo, no he oído hablar de ningún algoritmo de filtrado mediano en el procesamiento de imágenes que aproveche este hecho.

Es posible que desee echar un vistazo al Ejercicio 5.3.3.13 de Knuth. Describe un algoritmo debido a Floyd que encuentra la mediana de n elementos usando (3/2) n + O (n ^ (2/3) log n) comparaciones, y la constante oculta en el O (·) parece no ser demasiado grande en la práctica.

Estoy apostando a que podrías calcularlos sin costo alguno, en un hilo separado mientras se carga desde el disco (o como se generen).

Lo que realmente estoy diciendo es que la "velocidad" no vendrá de un poco de giros, porque 27 valores no son suficientes para que la notación de Big O sea un factor real.

La pregunta no puede ser respondida fácilmente por la simple razón de que el rendimiento de un algoritmo relativo a otro depende tanto de la combinación de compilador / procesador / estructura de datos como del algoritmo en sí, como seguramente sabes

Por lo tanto, su enfoque para probar un par de ellos parece lo suficientemente bueno. Y sí, quicksort debería ser bastante rápido. Si no lo ha hecho, le recomendamos que intente insertionsort, que a menudo rinde mejor en conjuntos de datos pequeños. Dicho esto, basta con resolver un problema de clasificación que hace el trabajo lo suficientemente rápido. Por lo general, no obtendrás 10 veces más rápido simplemente eligiendo el "correcto" algo.

Para obtener aceleraciones sustanciales, la mejor manera es frecuentemente usar más estructura. Algunas ideas que me funcionaron en el pasado con problemas a gran escala:

¿Puedes precalcular eficientemente mientras creas los voxels y almacenar 28 en lugar de 27 flotadores?
¿Es una solución aproximada lo suficientemente buena? Si es así, solo mire la mediana de, digamos 9 valores, ya que "en general se puede esperar que los valores sean relativamente cercanos". O puede reemplazarlo con el promedio siempre que los valores estén relativamente cerca.
¿Realmente necesitas la mediana para todos los miles de millones de vóxeles? Tal vez tengas una prueba fácil, ya sea que necesites la mediana, y solo puedas calcular el subconjunto correspondiente.
Si nada más ayuda: mira el código asm que genera el compilador. Es posible que pueda escribir código asm que sea sustancialmente más rápido (por ejemplo, haciendo todos los calcs usando registros).

Editar: Por lo que vale, he adjuntado el código insertionsort (parcial) mencionado en el comentario a continuación (totalmente no probado). Si los numbers[] son una matriz de tamaño N , y desea los flotantes P más pequeños ordenados al comienzo de la matriz, llame a partial_insertionsort<N, P, float>(numbers); . Por lo tanto, si llama partial_insertionsort<27, 13, float>(numbers); , los numbers[13] contendrán la mediana. Para ganar velocidad adicional, tendrías que desplegar el ciclo while también. Como se discutió anteriormente, para ser realmente rápido, debe usar su conocimiento sobre los datos (por ejemplo, ¿los datos ya están parcialmente ordenados? ¿Conocen las propiedades de la distribución de los datos? Supongo, entienden).

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>());}

Mi algoritmo súper rápido para el cálculo de la mediana de un conjunto de datos 1-D hace el trabajo en tres pasos y no necesita ordenar (!!!) el conjunto de datos.

Una descripción muy genérica es la siguiente:

Pase 1: escanea el conjunto de datos 1-D y recopila cierta información estadística del conjunto de datos
Pase 2: utiliza la información estadística del conjunto de datos y aplica algunos datos de minería de datos para crear una matriz intermedia (auxiliar)
Pase 3: escanea la matriz intermedia (ayudante) para encontrar la mediana

El algoritmo está diseñado para encontrar medianas de conjuntos de datos 1-D extremadamente grandes mayores que 8GE (elementos giga) de valores de punto flotante de precisión simple (en un sistema de escritorio con 32 GB de memoria física y 128 GB de memoria virtual) o para encontrar medianas de pequeños conjuntos de datos en un entorno duro en tiempo real.

El algoritmo es:

más rápido que el algoritmo clásico, basado en el algoritmo de clasificación Heap o Merge, en ~ 60 - ~ 75 veces
implementado en lenguaje C puro
no usa ninguna función intrínseca de Intel
no usa ninguna instrucción de ensamblador en línea
absolutamente portátil entre compiladores C / C ++, como MS, Intel, MinGW, Borland, Turbo y Watcom
absolutamente portátil entre plataformas

Saludos cordiales, Sergey Kostrov

Si desea ver algoritmos, busque los libros de Donald E. Knuth.

PD. Si crees que has inventado algo mejor, entonces deberías poder demostrar que la complejidad es similar o mejor a la complejidad de los algoritmos conocidos. Que para las variaciones basadas en cubo y raíz es O (n), por otro lado, ordenar rápidamente es solo O (n.log (n)). Un método que es un 20% más rápido sigue siendo O (n.log (n)) hasta que pueda mostrar el algoritmo :-)

Si hay 3x3x3 = 27 valores posibles (en caso afirmativo, ¿por qué los flotadores?), ¿Puede crear una matriz de 27 elementos y contar cada posibilidad en una sola pasada a través de los datos?

Supongo que tu mejor opción es tomar un algoritmo de clasificación existente y tratar de averiguar si puedes adaptarlo para que el conjunto no tenga que estar totalmente ordenado. Para determinar la mediana, necesita como máximo la mitad de los valores ordenados, ya sea la mitad inferior o superior sería suficiente:

original: | 5 | 1 | 9 | 3 | 3 | sorted: | 1 | 3 | 3 | 5 | 9 | lower half sorted: | 1 | 3 | 3 | 9 | 5 | higher half sorted: | 3 | 1 | 3 | 5 | 9 |

La otra mitad sería un cubo de valores sin clasificar que simplemente comparten la propiedad de ser más grande / más pequeño o igual al valor ordenado más grande / más pequeño.

Pero no tengo un algoritmo listo para eso, es solo una idea de cómo puedes tomar un atajo en tu clasificación.

Una red de clasificación generada utilizando el algoritmo de Bose-Nelson encontrará la mediana directamente sin bucles / recursión utilizando 173 comparaciones. Si tiene la posibilidad de realizar comparaciones en paralelo, como el uso de las instrucciones vector aritméticas, entonces puede agrupar las comparaciones en tan solo 28 operaciones paralelas.

Si está seguro de que los flotadores están normalizados y no (cs) NaN, entonces puede usar operaciones enteras para comparar flotantes IEEE-754 que pueden funcionar de manera más favorable en algunas CPU.

Una conversión directa de esta red de clasificación a C (gcc 4.2) produce un peor caso de 388 ciclos de reloj en mi Core i7.

Clasificación de redes

EDITAR: Tengo que disculparme. El código a continuación era INCORRECTO. Tengo el código fijo, pero necesito encontrar un compilador icc para rehacer las medidas.

Los resultados de referencia de los algoritmos considerados hasta ahora

Para el protocolo y breve descripción de los algoritmos, ver a continuación. El primer valor es el tiempo medio (segundos) sobre 200 secuencias diferentes y el segundo valor es stdDev.

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

Protocolo: genera 27 flotantes aleatorias utilizando bits aleatorios obtenidos de rand (). Aplique cada algoritmo 5 millones de veces seguidas (incluida la copia de matriz anterior) y calcule el promedio y el stdDev en 200 secuencias aleatorias. Código de C ++ compilado con icc -S-O3 y ejecutado en Intel E8400 con 8GB DDR3.

Algoritmos:

HeapSort: tipo completo de secuencia que utiliza la clasificación de montón y elige el valor medio. Implementación ingenua utilizando acceso de subíndice.

QuickSort: tipo de secuencia completa en el lugar usando clasificación rápida y selección del valor medio. Implementación ingenua utilizando acceso de subíndice.

QuickMedian1: algoritmo de selección rápida con intercambio. Implementación ingenua utilizando acceso de subíndice.

HeapMedian1: en lugar de método de pila equilibrada con intercambio previo. Implementación ingenua utilizando acceso de subíndice.

NthElement: utiliza el algoritmo STL nth_element. Los datos se copian en el vector utilizando memcpy (vct.data (), rndVal, ...);

QuickMedian2: utiliza un algoritmo de selección rápida con punteros y copia en dos almacenamientos intermedios para evitar el cambio. Basado en la propuesta de MSalters.

HeapMedian2: variante de mi algoritmo inventado usando montones dobles con cabezas compartidas. El montón izquierdo tiene el mayor valor como cabeza, el derecho tiene el valor más pequeño como cabeza. Inicialice con el primer valor como conjetura de cabeza común y primer valor mediano. Agregue valores posteriores al montón izquierdo si es más pequeño que la cabeza, de lo contrario al montón derecho, hasta que uno del montón esté lleno. Está lleno cuando contiene 14 valores. Luego considera solo el montón completo. Si es el montón correcto, para todos los valores más grandes que la cabeza, pop head e insert value. Ignora todos los otros valores. Si es el montón izquierdo, para todos los valores más pequeños que la cabeza, separe el encabezado e insértelo en el montón. Ignora todos los otros valores. Cuando todos los valores han sido procesados, la cabeza común es el valor mediano. Utiliza un índice entero en una matriz. La versión que usa punteros (64 bits) parece ser casi dos veces más lenta (~ 1s).

HeapMedian3: el mismo algoritmo que HeapMedian2 pero optimizado. Utiliza el índice de caracteres sin signo, evita el intercambio de valores y otras cosas pequeñas. Los valores de media y stdDev se calculan a través de 1000 secuencias aleatorias. Para nth_element medí 0.508s y un stdDev de 0.159537 con las mismas 1000 secuencias aleatorias. HeapMedian3 es por lo tanto 33 veces más rápido que la función nth_element stl. Cada valor mediano devuelto se compara con el valor mediano devuelto por heapSort y todos coinciden. Dudo que un método que use hash pueda ser significativamente más rápido.

EDIT 1: Este algoritmo se puede optimizar aún más. La primera fase donde los elementos se envían en el montón izquierdo o derecho según el resultado de la comparación no necesita montones. Es suficiente simplemente agregar elementos a dos secuencias desordenadas. La fase uno se detiene tan pronto como una secuencia está llena, lo que significa que contiene 14 elementos (incluido el valor mediano). La segunda fase comienza al heapificar la secuencia completa y luego proceder como se describe en el algoritmo HeapMedian3. Proporcionaré el nuevo código y punto de referencia lo antes posible.

EDIT 2: Implementé y comparé el algoritmo optimizado. Pero no hay una diferencia de rendimiento significativa comparada con heapMedian3. Es incluso un poco más lento en promedio. Se muestran los resultados confirmados. Puede haber con conjuntos mucho más grandes. Tenga en cuenta también que simplemente escojo el primer valor como conjetura mediana inicial. Como se sugirió, uno podría beneficiarse del hecho de que buscamos un valor mediano en conjuntos de valores "superpuestos". Usar la mediana del algoritmo mediano ayudaría a obtener una mejor estimación inicial del valor de la mediana.

Código fuente 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; } } } }