Самый быстрый код C / C ++ для выбора медианы в наборе из 27 значений с плавающей запятой - PullRequest
39 голосов
/ 01 мая 2009

Это хорошо известный алгоритм выбора. см http://en.wikipedia.org/wiki/Selection_algorithm.

Мне нужно, чтобы найти медианное значение набора значений вокселей 3x3x3. Поскольку объем состоит из миллиарда вокселей, а алгоритм рекурсивный, лучше быть немного быстрее. В целом можно ожидать, что значения относительно близки.

Самый быстрый из известных алгоритмов, которые я когда-либо пробовал, использует функцию секционирования быстрой сортировки. Я хотел бы знать, есть ли более быстрый.

Я «изобрел» 20% более быстрый, используя две кучи, но ожидал еще более быстрый, используя хэш. Прежде чем реализовать это, я хотел бы знать, существует ли уже быстрое быстрое решение.

Тот факт, что я использую числа с плавающей запятой, не должен иметь значения, поскольку их можно рассматривать как целое число без знака после инвертирования знакового бита. Заказ будет сохранен.

РЕДАКТИРОВАТЬ: тест и исходный код перенесены в отдельный ответ, как это было предложено Дэви Лэндман. Ниже приведен ответ от chmike.

РЕДАКТИРОВАТЬ : На сегодняшний день Boojum ссылается на наиболее эффективный алгоритм как ссылку на статью Быстрая медиана и двусторонняя фильтрация , которая теперь является ответом на этот вопрос. Первая разумная идея этого метода - использовать радикальную сортировку, вторая - объединить медианный поиск соседних пикселей, которые совместно используют много пикселей.

Ответы [ 15 ]

28 голосов
/ 01 мая 2009

Алгоритм выбора - линейное время (O (n)). Сложность, которую вы не можете сделать лучше, чем линейное время, потому что для считывания всех данных требуется линейное время. Таким образом, вы не могли бы сделать что-то более быстрое в отношении сложности. Возможно, у вас есть что-то, что является постоянным фактором быстрее на определенных входах? Я сомневаюсь, что это будет иметь большое значение.

C ++ уже включает алгоритм линейного выбора времени. Почему бы просто не использовать его?

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;

Редактировать: Технически, это только медиана для контейнера нечетной длины. Для одной четной длины он получит «верхнюю» медиану. Если вам нужно традиционное определение медианы для четной длины, вам, возможно, придется выполнить его дважды, один раз для каждой из двух «середин» на first + (last - first) / 2 и first + (last - first) / 2 - 1, а затем усреднить их или что-то еще.

21 голосов
/ 05 мая 2009

РЕДАКТИРОВАТЬ: я должен извиниться. Код ниже был НЕПРАВИЛЬНЫМ. У меня есть фиксированный код, но мне нужно найти компилятор icc , чтобы повторить измерения.

Результаты тестов алгоритмов, рассмотренных до сих пор

Протокол и краткое описание алгоритмов см. Ниже. Первое значение - это среднее время (в секундах) для 200 различных последовательностей, а второе - 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

Протокол: сгенерируйте 27 случайных чисел с использованием случайных битов, полученных из rand (). Примените каждый алгоритм 5 миллионов раз подряд (включая предыдущее копирование массива) и вычислите среднее значение и stdDev для 200 случайных последовательностей. Код C ++, скомпилированный с помощью icc -S -O3 и работающий на Intel E8400 с 8 ГБ памяти DDR3.

Алгоритмы:

HeapSort: полная сортировка последовательности с использованием сортировки кучи и выбора среднего значения. Наивная реализация с использованием индексного доступа.

Быстрая сортировка: полная последовательность сортировки по месту с использованием быстрой сортировки и выбора среднего значения. Наивная реализация с использованием индексного доступа.

QuickMedian1: алгоритм быстрого выбора с заменой. Наивная реализация с использованием индексного доступа.

HeapMedian1: метод сбалансированной кучи на месте с предыдущей заменой. Наивная реализация с использованием индексного доступа.

NthElement: используется алгоритм STL nth_element. Данные копируются в вектор с использованием memcpy (vct.data (), rndVal, ...);

QuickMedian2: использует алгоритм быстрого выбора с указателями и копирование в два буфера, чтобы избежать подкачки. По предложению MSalters.

HeapMedian2: вариант моего изобретенного алгоритма, использующего двойные кучи с общими головками. Левая куча имеет наибольшее значение как голова, правая имеет наименьшее значение как голова. Инициализируйте с первым значением в качестве общей головы и первым медианным значением. Добавьте последующие значения в левую кучу, если она меньше, чем head, в противном случае в правую кучу, пока одна из куч не заполнится. Он полон, когда содержит 14 значений. Тогда рассмотрим только полную кучу. Если это правильная куча, для всех значений больше, чем голова, выдвиньте голову и вставьте значение. Игнорируйте все другие значения. Если это левая куча, для всех значений, меньших, чем голова, выдвиньте голову и вставьте ее в кучу. Игнорируйте все другие значения. Когда все значения получены, общая голова является медианным значением. Он использует целочисленный индекс в массиве. Версия с использованием указателей (64 бита) оказалась почти в два раза медленнее (~ 1 с).

HeapMedian3: тот же алгоритм, что и в HeapMedian2, но оптимизированный. Он использует беззнаковый индекс, избегает обмена значениями и других мелочей. Средние значения и значения stdDev вычисляются для более 1000 случайных последовательностей. Для nth_element я измерил 0,508 с и stdDev 0,159537 с теми же 1000 случайными последовательностями. Таким образом, HeapMedian3 в 33 раза быстрее, чем функция nth_element stl. Каждое возвращаемое медианное значение проверяется по медианному значению, возвращенному heapSort, и все они совпадают. Я сомневаюсь, что метод, использующий хэш, может быть значительно быстрее.

EDIT 1: этот алгоритм может быть дополнительно оптимизирован. Первая фаза, где элементы отправляются в левой или правой куче на основе результата сравнения, не требует куч. Достаточно просто добавить элементы к двум неупорядоченным последовательностям. Фаза 1 останавливается, как только одна последовательность заполнена, что означает, что она содержит 14 элементов (включая медианное значение). Второй этап начинается с накапливания полной последовательности, а затем продолжается, как описано в алгоритме HeapMedian3. Я предоставлю новый код и тест как можно скорее.

РЕДАКТИРОВАТЬ 2: я реализовал и протестировал оптимизированный алгоритм. Но нет существенной разницы в производительности по сравнению с heapMedian3. Это даже немного медленнее в среднем. Показанные результаты подтверждены. Там может быть с гораздо большими наборами. Также обратите внимание, что я просто выбираю первое значение в качестве начального среднего значения. Как и предполагалось, можно извлечь пользу из того факта, что мы ищем медианное значение в «перекрывающихся» наборах значений. Использование алгоритма медианы медианы помогло бы выбрать намного лучшее начальное предположение о медиане.


Исходный код 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;
           }
       }
   }
} 
14 голосов
/ 12 сентября 2009

Поскольку звучит так, будто вы выполняете медианный фильтр для большого массива объемных данных, вы можете взглянуть на статью Быстрая медиана и двусторонняя фильтрация из SIGGRAPH 2006. Эта статья посвящена с обработкой 2D-изображений, но вы можете адаптировать алгоритм для 3D-объемов. Если ничего другого, это может дать вам некоторые идеи о том, как сделать шаг назад и посмотреть на проблему с несколько иной точки зрения.

13 голосов
/ 01 мая 2009

На этот вопрос нелегко ответить по той простой причине, что производительность одного алгоритма относительно другого зависит в такой же степени от комбинации компилятора / процессора / структуры данных, как и от самого алгоритма, как вы наверняка знаете

Поэтому ваш подход попробовать пару из них кажется достаточно хорошим. И да, быстрая сортировка должна быть довольно быстрой. Если вы этого не сделали, вы можете попробовать сортировку вставок, которая часто работает лучше для небольших наборов данных. Это сказало, просто выбирайте алгоритм сортировки, который делает работу достаточно быстро. Как правило, вы не получите в 10 раз быстрее, просто выбрав «правильный» алгоритм.

Чтобы добиться существенного ускорения, лучше часто использовать больше структур. Некоторые идеи, которые работали для меня в прошлом с крупномасштабными проблемами:

  • Можете ли вы эффективно предварительно рассчитать при создании вокселей и сохранить 28 вместо 27 поплавков?

  • Достаточно ли приблизительное решение? Если так что, просто посмотрите на медиану, скажем, 9 значения, так как "в целом это может быть Ожидается, что значения относительно закрыть. "Или вы можете заменить его среднее до тех пор, пока значения относительно близко.

  • Вам действительно нужна медиана для всех миллиарды вокселей? Может быть, у вас есть легко проверить, нужно ли вам Медиана, а затем может только рассчитать для соответствующего подмножества.

  • Если больше ничего не помогает: посмотрите на ASM-код, который генерирует компилятор. Вы могли бы написать код asm, который существенно быстрее (например, делая все вызовы с использованием регистров).

Редактировать: Для чего бы то ни было, я прикрепил (частично) код вставки, упомянутый в комментарии ниже (полностью не проверен). Если numbers[] - это массив размером N, и вы хотите, чтобы самые маленькие значения P были отсортированы в начале массива, вызовите partial_insertionsort<N, P, float>(numbers);. Следовательно, если вы позвоните partial_insertionsort<27, 13, float>(numbers);, numbers[13] будет содержать медиану. Чтобы получить дополнительную скорость, вам также придется развернуть цикл while. Как уже говорилось выше, чтобы получить действительно быстрый результат, вы должны использовать свои знания о данных (например, данные уже частично отсортированы? Знаете ли вы свойства распределения данных? Я полагаю, вы получаете дрейф).

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>());}
6 голосов
/ 01 мая 2009

Наиболее вероятный алгоритм для использования в первой попытке - просто nth_element; это в значительной степени дает вам то, что вы хотите напрямую. Просто спросите 14-й элемент.

С вашей второй попытки цель состоит в том, чтобы воспользоваться фиксированным размером данных. Вы не хотите выделять какую-либо память вообще из-за вашего алгоритма. Итак, скопируйте ваши значения вокселей в предварительно выделенный массив из 27 элементов. Выберите стержень и скопируйте его в середину массива из 53 элементов. Скопируйте оставшиеся значения по обе стороны от оси. Здесь вы держите два указателя (float* left = base+25, *right=base+27). Теперь есть три варианта: левая сторона больше, правая сторона больше или у обоих по 12 элементов. Последний случай тривиален; ваш стержень является медианой. В противном случае вызовите nth_element с левой или правой стороны. Точное значение Nth зависит от того, сколько значений было больше или меньше, чем опорная точка. Например, если деление 12/14, вам нужен наименьший элемент больше, чем стержень, поэтому Nth = 0, а если деление было 14/12, вам нужен самый большой элемент, меньший стержень, поэтому Nth = 13. Худшие случаи - это 26/0 и 0/26, когда ваш пивот был экстремальным, но это случается только в 2/27 от всех случаев.

Третье улучшение (или первое, если вам нужно использовать C и у вас нет nth_element) полностью заменяет nth_element. У вас все еще есть массив из 53 элементов, но на этот раз вы заполняете его непосредственно из значений вокселей (сохраняя промежуточную копию в float[27]). Поворот в этой первой итерации - просто воксель [0] [0] [0]. Для последующих итераций вы используете второй предварительно выделенный float[53] (проще, если оба имеют одинаковый размер) и копирование с плавающей точкой между ними. Основной шаг итерации здесь все еще: скопируйте ось в середину, отсортируйте остальное влево и вправо. В конце каждого шага вы будете знать, является ли медиана меньше или больше, чем текущий пивот, так что вы можете сбросить поплавки больше или меньше, чем этот пивот. За одну итерацию это исключает от 1 до 12 элементов, в среднем 25% от оставшихся.

Последняя итерация, если вам все еще нужна большая скорость, основана на наблюдении, что большинство ваших вокселей значительно перекрываются. Для каждого среза 3x3x1 вы предварительно рассчитываете медианное значение. Затем, когда вам нужен начальный шарнир для воксельного куба 3x3x3, вы берете медиану из трех. Априори вы знаете, что на 9 вокселей меньше и на 9 вокселей больше, чем медиана медиан (4 + 4 + 1). Таким образом, после первого шага разворота наихудшими являются деление 9/17 и 17/9. Таким образом, вам нужно будет найти только 4-й или 13-й элемент в поплавке [17] вместо 12-го или 14-го в поплавке [26].


Справочная информация. Идея скопировать сначала сводную, а затем оставшуюся часть поплавка [N] в поплавок [2N-1], используя левый и правый указатели, заключается в том, что вы заполняете подрешетку поплавка [N] вокруг этой оси, со всеми элементами меньше, чем стержень слева (нижний индекс) и выше вправо (более высокий индекс). Теперь, если вам нужен элемент Mth, вам может повезти, и элементы M-1 будут меньше, чем стержень, и в этом случае стержень - это элемент, который вам нужен. Если имеется больше, чем (M-1) элементов, меньших, чем центр, среди них есть M-элемент, так что вы можете отбросить элемент и все, что больше, чем элемент Pivot, и seacrh для элемента Mth во всех более низких значениях. Если на меньше , чем (M-1) элементов меньше, чем пивот, вы ищете значение выше, чем пивот. Таким образом, вы откажетесь от оси и всего, что меньше ее. Пусть число элементов на меньше , чем в сводной точке, то есть слева от этой оси будет L. На следующей итерации вы хотите, чтобы (ML-1) -й элемент (NL-1) с плавающей точкой больше оси.

Этот тип алгоритма nth_element довольно эффективен, поскольку большая часть работы проводится при копировании поплавков между двумя небольшими массивами, оба из которых будут в кеше, а также потому, что ваше состояние в большинстве случаев представлено 3 указателями (указатель источника, левый указатель назначения, правый указатель назначения).

Показать основной код:

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
4 голосов
/ 11 сентября 2009

Сортировочная сеть, созданная с использованием алгоритма Бозе-Нельсона, найдет медиану напрямую без циклов / рекурсии, используя 173 сравнения. Если у вас есть возможность выполнять параллельные сравнения, такие как использование векторно-арифметических инструкций, вы можете сгруппировать сравнения в всего 28 параллельных операций.

Если вы уверены, что числа с плавающей запятой нормализованы, а не числа (qs), то вы можете использовать целочисленные операции для сравнения операций с плавающей запятой IEEE-754, которые могут работать более эффективно на некоторых процессорах.

Прямое преобразование этой сортировочной сети в C (gcc 4.2) дает наихудший случай из 388 тактов на моем Core i7.

Сортировка сетей

3 голосов
/ 01 мая 2009

Полагаю, вам лучше всего выбрать существующий алгоритм сортировки и попытаться выяснить, можете ли вы его адаптировать так, чтобы набор не нуждался в полной сортировке. Для определения медианы вам потребуется не более половины отсортированных значений, достаточно будет либо нижней, либо верхней половины:

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 |

Другая половина будет набором несортированных значений, которые просто имеют свойство быть больше / меньше или равны наибольшему / наименьшему отсортированному значению.

Но у меня нет готового алгоритма для этого, это просто представление о том, как вы могли бы сократить путь при сортировке.

2 голосов
/ 01 мая 2009

+ 1 для всех, кто упомянул nth_element, но этот вид кода лучше, чем рукописный алгоритм лучше, чем STL, потому что вы хотите сгенерировать наиболее эффективный код для этого одного компилятора, работающего на одном CPU с конкретным набором данных. Например, для некоторой комбинации процессор / компилятор std :: swap (int, int) может быть медленнее, чем рукописный своп с использованием XOR (прежде чем ответить, я знаю, что это, вероятно, верно 20 лет назад, но не больше). Иногда производительность достигается за счет написания кода ассемблера, специфичного для вашего процессора. Если вы планируете использовать потоковые процессоры GPU, вам, возможно, придется разработать свой алгоритм соответствующим образом.

Вы упомянули, что используете 2 кучи и следите за медианой при вставке. Это то, что я сделал некоторое время назад в проекте. Я изменил массив на месте и использовал только одну кучу. Я не мог придумать более быстрого алгоритма, но я хотел бы предупредить вас об использовании памяти, в частности, кэш-памяти процессора. Вы хотите быть осторожным с доступом к памяти. Кэш-память ЦП переставляется на страницу, поэтому вы хотите, чтобы ваш алгоритм касался памяти, которая находится близко друг к другу, чтобы минимизировать потерю кэш-памяти ЦП.

2 голосов
/ 01 мая 2009

Новая книга Алексея Степанова Элементы программирования подробно рассказывает о поиске статистики заказов с использованием минимального числа средних сравнений при минимизации накладных расходов времени выполнения. К сожалению, для вычисления медианы из 5 элементов требуется значительный объем кода, и даже в этом случае он дает в качестве проекта поиск альтернативного решения, которое использует в среднем долю сравнения меньше, поэтому я бы не стал мечтать о его расширении. рамки для нахождения медианы 27 элементов. И книга даже не будет доступна до 15 июня 2009 года. Дело в том, что, поскольку это проблема фиксированного размера, существует метод прямого сравнения, который доказуемо оптимален.

Кроме того, существует тот факт, что этот алгоритм запускается не один раз в отдельности, а довольно много раз, и между большинством прогонов изменится только 9 из 27 значений. Это означает, что в теории часть работы уже сделана. Однако я не слышал о каких-либо алгоритмах медианной фильтрации при обработке изображений, которые бы использовали этот факт.

1 голос
/ 06 октября 2010

При наличии, скажем, миллиона различных значений, из которых вам нужна медиана. Можно ли основать вашу медиану на подмножестве этих миллионов, скажем, 10%. Так что медиана близка к n-му элементу, который делит значения на 2 равных (или почти равных) подмножества? Поэтому для нахождения медианы вам понадобится меньше, чем за O (n) раз (в данном случае O (1/10n), и тем самым приблизитесь к оптимальной сортировке с быстрой сортировкой по O (nlogn)?

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...