О современных процессорах и несовременных алгоритмах
Профайлю потихоньку RawDigger и после серии очевидных улучшений - верхнюю строчку в профиле занимает подсчет статистики всего изображения. Среднее/минимум/максимум/дисперсия по каналам.
У меня там используется классический однопроходный алгоритм (я запутался в цитатнике википедийном, назвать автора не могу, сами смотрите). Чтобы не всухую, скопирую его сюда прямо из википедии:
n = n + 1
delta = x - mean
mean = mean + delta/n
M2 = M2 + delta*(x - mean)
variance = M2/(n - 1)
return variance
Посмотрел я на это дело - и спрограммировал в два прохода. Прямо вот по классике, одним проходом считаем среднее, вторым проходом - сумму квадратов отклонений. По компенсированной (в смысле точности) формуле, хотя на практике разницы не вижу.
Прикол в том, что скорость работы у скалярного варианта ровно такая же в 2.2 раза лучше (см. update), хотя количество проходов по памяти - вроде стало вдвое больше, а количество инструкции - ну примерно одного порядка.
Я собственно к тому клоню, что два прохода по памяти - могут оказаться лучше одного прохода, но с плохим (для расчетов, векторизации) телом цикла. Или, как минимум, не хуже.
Update: в первой версии этого поста было написано, что скорость работы 1- и 2-проходного варианта - одинакова. Это действительно было так пока я не заметил слона - порядок обхода памяти был неудачный. Замена на удачный - ускорила однопроходный вариант (с делением и зависимостью команд) в ~1.3 раза, а двухпроходный - втрое.
Comments
Есть удобная формула для
Есть удобная формула для расчета variance - используется, например, когда заранее неизвестно количество данных.
Нужно посчитать две суммы: сумму значений и сумму их квадратов. Затем
variance = (сумма квадратов - сумма^2/N)/(N-1)
В цикле нужно только складывать и умножать, а делить - только один раз в самом конце.
Взял отсюда: http://www.dspguide.com/ch2/2.htm
Проблема этой формулы, что
Проблема этой формулы, что она в условиях реальной арифметики ограниченной точности - может очень сильно провраться.
На эти грабли я уже успел наступить. Кажется, не в публичных версиях RawDigger, уже хорошо.
Вот отголоски этой темы,
Вот отголоски этой темы, когда я налетел на ошибку в округлении при *сложении*: http://blog.lexa.ru/2011/03/14/o_plavayushchei_tochke.html
тА чем "наивный" алгоритм не устраивает? Неужели на этих дан
тА чем "наивный" алгоритм не устраивает? Неужели на этих данных он дает заметные потери точности?
Да, там на миллионах (зачастую близких) значений - полная ка
Да, там на миллионах (зачастую близких) значений - полная катастрофа. Вообще ни в какие ворота, ошибки очень большие (бывают).
У меня из той же степи, но другое было клево (тоже потеря то
У меня из той же степи, но другое было клево (тоже потеря точности при суммировании целых в float):
минимум - 59259
максимум - 59259
дисперсия - 0
среднее - 59259.7
Когда префетчер вовремя тащит данные в L1D через L2 и LLC, к
Когда префетчер вовремя тащит данные в L1D через L2 и LLC, как в этом случае правильного порядка обхода, действительно не такая большая разница 2 раза ходить по памяти или 1, если с 2 CPI лучше.
Раз все берется из L1D, для сравнения может помочь software.intel.com/en-us/articles/intel-architecture-code-analyzer
Я мерял на ~70M данных (36Mpix файл, 2 байта на пиксель). Т.
Я мерял на ~70M данных (36Mpix файл, 2 байта на пиксель). Т.е. разница должна быть заметной.
Но CPI в данном случае бьет transfer.
да, деление всё портит. проще прочитать 2 раза, но без делен
да, деление всё портит. проще прочитать 2 раза, но без деления.
Почему бы не использовать <a href="http://www.strchr.com/sta
Почему бы не использовать более простой однопроходный алгоритм? Если важна точность, можно считать сумму квадратов в double, всё равно будет быстрее.
Вы не поверите - но для миллионов пикселей даже double не хв
Вы не поверите - но для миллионов пикселей даже double не хватает, оно сильно провирается.
Т.е. у меня этот алгоритм в одной версии был - и был срочно устранен.
А если добавить <a href="http://en.wikipedia.org/wiki/Kahan_
А если добавить Kahan's summation algorithm?
Да, возможно что-то скомпенсированное помогло бы. Но в реаль
Да, возможно что-то скомпенсированное помогло бы. Но в реальности у меня сейчас (многопоточный) двухпроходный вариант работает что-то в районе 40мс на 36-мегапиксельном файле и это пока устраивает, уже не буду переделывать.
А про порядок обхода памяти можно подробнее?
А про порядок обхода памяти можно подробнее?
Ну вот для изображения Хорошо: for(row=0...) for(col=0..)
Ну вот для изображения
Хорошо:
for(row=0...)
for(col=0..)
image[row*rowpitch+col] =...
А плохо, соответственно, если не подряд:
for(col=0..)
for(row=0...)
image[row*rowpitch+col] =...
Чёрт, я хотел предположить, что обход по столбцам сделан вме
Чёрт, я хотел предположить, что обход по столбцам сделан вместо строчек, но побоялся ))
Спасибо за ответ.
Прикол в том, что для варианта с делением (однопроходного) р
Прикол в том, что для варианта с делением (однопроходного) разница была не настолько большая.
Ну, видимо трудоёмкость операции перевешивала. Хотя забавно,
Ну, видимо трудоёмкость операции перевешивала. Хотя забавно, да.
Надо суммы считать по частям.
Надо суммы считать по частям. Например, сначала по каждой строке отдельно, а потом просуммировать суммы строк.
И получится двухпроходный
И получится двухпроходный алгоритм. О чем и речь.
Ну и что? Само количество
Ну и что? Само увеличение количества проходов ничего не ухудшает.
Я в первую очередь имел в
Я в первую очередь имел в виду проблему точности.
Для точности - есть
Для точности - есть "компенсированный алгоритм".
А некомпенсированный может провираться даже на очень небольших данных (вроде строки)
"А некомпенсированный может
"А некомпенсированный может провираться даже на очень небольших данных (вроде строки)"
Вот этого я как-то не понимаю. Даже страшно жить :) Уже сумму пару тысяч чисел подсчитать невозможно...
Ага, при некотором невезении
Ага, при некотором невезении - жить страшно.
Вот я про это писал, как выясняется полтора года назад (а как вчера было): http://blog.lexa.ru/2011/03/14/o_plavayushchei_tochke.html
Везенье не при чем. Все
Везенье не при чем. Все оценивается, мне кажется. Если, скажем, у вас значения 16-битные (судя по 59 тысячам), а мантисса у float 24 бита, то сумма 2^(24-16) = 2^8 = 256 слагаемых будет вычисляться точно, а дальше... Так что, лучше суммировать вообще в 32-разрядных целых, тогда можно сложить точно до 2^(31-16) = 2^15 = 31768 значений - для одной строки как раз достаточно. Ну, там еще квадраты надо складывать... double для одной строки должно хватать.
Да, double хватает для
Да, double хватает для строки, но если потом сложить тысячи этих double (по числу строк) - то что будет?
Для реальной фотографии все
Для реальной фотографии все будет отлично. Вырожденные случаи... надо посмотреть...
Почитал внимательней исходный
Почитал внимательней исходный Ваш пост, проблема то, как понял, снята...
Проблема была снята double +
Проблема была снята double + компенсация. Прямо по википедии. Чтобы вообще не думать о вырожденных случаях.