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

Ответ 1

Так как кажется, что вы выполняете срединный фильтр на большом массиве данных тома, вы можете взглянуть на Fast Median и Bilateral Filtering от SIGGRAPH 2006. В этом документе рассматривается обработка 2D-изображений, но вы можете адаптировать алгоритм для трехмерных томов. Если ничего другого, это может дать вам некоторые идеи о том, как отступить и посмотреть на проблему с немного другой точки зрения.

Ответ 2

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

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

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, а затем усреднить их или что-то в этом роде.

Ответ 3

РЕДАКТОР: Я должен извиниться. Код ниже был НЕПРАВИЛЬНЫМ. У меня есть фиксированный код, но вам нужно найти компилятор 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 случайных последовательностей. С++, скомпилированный с icc -S-O3 и запущенный на Intel E8400 с 8 ГБ DDR3.

Алгоритмы:

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

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

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

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

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

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

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

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

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

EDIT 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;
           }
       }
   }
} 

Ответ 4

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

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

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

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

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

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

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

Изменить: Для чего это стоит, я добавил (частично) код insertionsort, упомянутый в комментарии ниже (полностью непроверенный). Если 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>());}

Ответ 5

Наиболее вероятным алгоритмом для использования в вашей первой попытке является только 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 медианное значение. Затем, когда вам нужен начальный стержень для кубка voxel 3x3x3, вы берете медианную из трех. Вы знаете априори, что 9 вокселей меньше и 9 вокселей больше, чем медиана медиан (4 + 4 + 1). Таким образом, после первого этапа поворота наихудшими являются 9/17 и 17/9 раскол. Таким образом, вам нужно будет найти только 4-й или 13-й элемент в поплавке [17] вместо 12-го или 14-го поплавка [26].

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

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

Чтобы показать базовый код:

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

Ответ 6

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

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 |

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

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

Ответ 7

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

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

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

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

Ответ 8

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

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

Ответ 9

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

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

Ответ 10

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

Ответ 11

Если вы хотите, чтобы алгоритмы искали книги Дональда Кнута.

PS. Если вы думаете, что придумали что-то лучше, то вы должны показать, что сложность аналогична или более сложна для известных алгоритмов. С другой стороны, для вариантов, основанных на ведро и радиксе, есть O (n), а quick-sort - только O (n.log(n)). Метод, который на 20% быстрее, по-прежнему равен O (n.log(n)), пока вы не сможете показать алгоритм: -)

Ответ 12

Я уверен, что вы могли бы рассчитать их за нулевую стоимость - в отдельном потоке при загрузке с диска (или, тем не менее, сгенерированы).

То, что я на самом деле говорю, заключается в том, что "скорость" не будет исходить из бит-трюка, потому что 27 значений недостаточно для обозначения Big O, чтобы быть реальным фактором.

Ответ 13

Возможно, вам стоит взглянуть на упражнение Кнута 5.3.3.13. Он описывает алгоритм Флойда, который находит медиану n элементов, используя (3/2) n + O (n ^ (2/3) log n) сравнения, а константа, скрытая в O (·), кажется, не слишком большой на практике.

Ответ 14

Если существует 3x3x3 = 27 возможных значений (если да, то почему float?), можете ли вы создать массив из 27 элементов и подсчитать каждую возможность за один проход через данные?

Ответ 15

Мой супер быстрый алгоритм вычисления медианы одномерного набора данных выполняет задание за три прохода и не требует сортировки (!!!) набора данных.

Очень общее описание выглядит следующим образом:

• Пасс 1: сканирует 1-D набор данных и собирает некоторую статистическую информацию из набора данных
• Pass 2: использует статистическую информацию набора данных и применяет некоторый интеллектуальный анализ данных для создания промежуточного (вспомогательного) массива
• Pass 3: сканирует промежуточный (вспомогательный) массив, чтобы найти медиану

Алгоритм предназначен для поиска медианов чрезвычайно больших одномерных наборов данных, превышающих 8GE (элементы гига) значений одноточечной плавающей точки (на настольной системе с 32 ГБ физической памяти и 128 ГБ виртуальной памяти) или для поиска медианов небольших наборов данных в жесткой среде реального времени.

Алгоритм:

• быстрее, чем классический алгоритм, основанный на алгоритме сортировки кучи или слияния, в ~ 60 - ~ 75 раз
• реализован на чистом языке C
• не использует функции встроенных функций Intel.
• не использует встроенные инструкции ассемблера
• абсолютно переносимый между компиляторами C/С++, такими как MS, Intel, MinGW, Borland, Turbo и Watcom.
• абсолютно переносимый между платформами

С уважением, Сергей Костров