Блог им. b34rcava1ry

Подводные камни программирования или как потерять производительность и не подать виду.

Обычно я просто читаю СМ и не лезу со своим «экспертным» мнением к людям. Но недавно я не выдержал и закусился на форуме, да еще и побил собственный рекорд, ушел в личный бан за 2 сообщения. Суть зарубы  очень проста:  утверждается что найти тренд можно не менее чем за 10млрд операций. За время чтения СМа я видел многое, но это похоже долго еще будет мой топ1 перлов по программированию.

Этот пост я написал для того самого человека, хотя он его и не увидит.

<Disclaimer>
Если ты настоящий программист за деньги, то ничего нового ты тут не найдешь.
Искать тренд мы пока не будем, а будем разбирать похожую задачу.
Для облегчения восприятия весь код написан на python.
</Disclaimer>

Дано: n случайных чисел
Найти: стандартное отклонение для последовательности из m чисел (m < n)

Искать тренд или стандартное отклонение от начала времен не совсем практично. Очевидно, что нужно выбирать некий временной интервал (возможно несколько) для которого будет считаться наша искомая величина, а дальше уже будет некая логика и экзекьюшн, и риск менеджмент, и все такое прочее.
Но вернемся к нашей задаче. Математика нам говорит что стандартное отклонение равно:
Подводные камни программирования или как потерять производительность и не подать виду.
где s:
Подводные камни программирования или как потерять производительность и не подать виду.
Ну что, решим эту задачку в лоб:
from math import sqrt
import random
import time
#n = [random.randint(0, 10) for i in range(10000)]
n = [3, 5, 0, 2, 1, 4, 5, 4, 1, 1, 1, 2, 5, 2, 1, 1, 5, 2, 4, 1]
print(n)
rng = 5
stddev = .0
rng_avr = .0
rng_n2 = .0
rng_elm = []
time_exe = []
tsum = .0
for i in range(len(n)):
    #start_time = time.time()
    rng_elm.append(n[i])
    if i < rng - 1:
        continue
    elif i == rng - 1:
        for j in rng_elm:
            rng_avr += j/rng
        for j in rng_elm:
            rng_n2 += pow((j - rng_avr),2)
        stddev = sqrt(rng_n2/(rng-1))
        print ('\tstddev: %.5f' % stddev)
        #time_exe.append((time.time() - start_time))
    else:
        rng_n2 = .0
        rng_avr = .0
        rng_elm.pop(0)
        for j in rng_elm:
            rng_avr += j/rng
        for j in rng_elm:
            rng_n2 += pow((j - rng_avr),2)
        stddev = sqrt(rng_n2/(rng-1))
        print ('\tstddev: %.5f' % stddev)
        #time_exe.append((time.time() - start_time))
#for i in time_exe:
#    tsum += i/len(time_exe)
#print(tsum)

На закомментированные строчки пока внимание не обращаем.
Сначала немного магии, подключим нужные либы и объявим переменные.
Мы будем перебирать все элементы (n) добавляя новые к концу промежуточного массива (rng_elm), пока не наберем m элементов. После этого мы посчитаем среднее значение (rng_avr), потом квадрат разности (rng_n2). Ну и в конце мы сможем посчитать отклонение и вывести его на экран. Для всех последующих элементов мы будем делать тоже самое, только сначала будем удалять один элемент из начала промежуточного массива.
Запускаем скрипт, должно получится следующее:
[3, 5, 0, 2, 1, 4, 5, 4, 1, 1, 1, 2, 5, 2, 1, 1, 5, 2, 4, 1]
stddev: 1.92354
stddev: 2.07364
stddev: 2.07364
stddev: 1.64317
stddev: 1.87083
stddev: 1.87083
stddev: 1.94936
stddev: 1.30384
stddev: 1.73205
stddev: 1.64317
stddev: 1.64317
stddev: 1.64317
stddev: 2.04939
stddev: 1.64317
stddev: 1.81659
stddev: 1.81659
Ну что ж, все работает. Можно переходить к следующему шагу. Удаляем все #, кроме верхней строчки. Там нужно переписать вот так:
n = [random.randint(0, 10) for i in range(10000)]
#n = [3, 5, 0, 2, 1, 4, 5, 4, 1, 1, 1, 2, 5, 2, 1, 1, 5, 2, 4, 1]
#print(n)
rng = 5000
Запускаем снова и видим что время на обработку одного шага у нас примерно равно 0.0019 секунды. И тут надо бы задать вопрос: А можно ли сделать это быстрее? Конечно можно!
Просто с самого начала нам нужно было немножко поиграть с формулой.
Подводные камни программирования или как потерять производительность и не подать виду.
Прямо из формулы видно, что нам достаточно собирать сумму элементов и их квадратов. Перепишем код.
from math import sqrt
import random
import time
n = [random.randint(0, 10) for i in range(100000)]
#n = [3, 5, 0, 2, 1, 4, 5, 4, 1, 1, 1, 2, 5, 2, 1, 1, 5, 2, 4, 1]
#print(n)
stddev = .0
rng = 50000
rng_sum = 0
rng_n2 = 0
rng_elm = []
time_exe = []
tsum = .0
for i in range(len(n)):
    start_time = time.time()
    rng_sum += n[i]
    rng_n2 += n[i]*n[i]
    rng_elm.append(n[i])
    if i < rng - 1:
        continue
    elif i == rng -1:
        stddev = sqrt((rng_n2 - (rng_sum*rng_sum/rng))/(rng-1))
        print ('\tstddev: %.5f' % stddev)
        time_exe.append((time.time() - start_time))
    else:
        dif = rng_elm.pop(0)
        rng_sum -= dif
        rng_n2 -= dif*dif
        stddev = sqrt((rng_n2 - (rng_sum*rng_sum/rng))/(rng-1))
        print ('\tstddev: %.5f' % stddev)
        time_exe.append((time.time() - start_time))
for i in time_exe:
    tsum += i/len(time_exe)
print(tsum)
Можно убедиться что считается все одинаково закомментировав все как в первом куске кода и выставив rng = 5.
Заметьте что в этот раз я еще больше увеличил интервал, теперь он не 5000, а 50000.
Запускаем и видим что теперь на обработку одного шага тратится 0.00002 (2.3248035863250115e-05) секунды. А все потому что мы избавились от ненужных переборов массивов. В таком виде производительность нашей программы уже не зависит от интервала. Мы с вами магически превратили О(n^2) в О(n).

Зачем я все это написал? Да просто потому что могу.
Какой из этого можно сделать вывод? Ну наверное можно найти тренд и за 10млрд операций, но для этого нужно взять интервал в 1млрд значений. И тут у нормального человека должны появиться сомнения в смысле такой величины,
В комментариях предлагайте темы, возможно я смогу рассказать что-то интересное по ним в следующих сериях рубрики «Программирование для самых маленьких и тупых».




★7
65 комментариев
А вот интересует такая вполне прикладная задача: надо мне посчитать скользящий skewness и kurtosis также за O(N), но без потери точности на длинных рядах (например, минутки или пятиминутки за 10 лет), и чтобы реально без потери точности, т.е. расчеты совпадали с расчетом «в лоб» в том же окне. Причем интересна именно ошибка на маленьких окнах (3-10 наблюдений).
avatar
MadQuant, попробую покрутить в свободное время. Может что и получится.
avatar
b34rcava1ry, среднеквадратичное отклонение так не считают, в лоб — вычислительно неустойчивый алгоритм. Надо по Кнуту — так и быстрее, и точнее, и для окон лучше подходит — можно чиселки как добавлять, так и выбрасывать. http://nabbla1.livejournal.com/83886.html
Zweroboi, дисклеймеры не читай, сразу отвечай)
avatar
Zweroboi, Кнут — это очень хорошо и правильно. Но ссылка почему-то не идет. Впрочем, всегда можно открыть книгу :)
avatar
SergeyJu, https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
MadQuant, а можете дать тестовый датасет?
avatar
А почему «в лоб» не катит? Трудоемкость будет N*расчет одного окна. Нет ни квадратичного никакого иного роста сложности в ростом N.
avatar
SergeyJu, Так O(N * Window) — это уже большая сложность. Например, если минутки, и мне надо посчитать скользящий skewness в окне сутки за 10 лет — это грубо 1440 * (10 * 252 * 1440) ~ 5.2 * 10^9 — как бы уже несколько секунд расчетов, предполагая расчеты за 1 такт и мгновенное обращение к памяти. А там еще сбои кэшей,… Да и контрактов не 1, а 50, например — и вот, мы уже сидим ждем этого счастья несколько минут.

Надо строго за O(N)
avatar
MadQuant, в сутках не 3-10 минутных наблюдений. 
и я не понимаю, откуда берется 1440 в квадрате. Что-то некорректно в ваших расчетах.

avatar
SergeyJu, 

1440 — минут в сутках (размеры окна, в котором считаем)
10 — число лет данных
252 — число дней в каждом году данных
1440 — число минут в каждом дне данных

avatar
MadQuant, так окно у Вас 3-10 или 1440?
В принципе, для расчета выборочного к-та асимметрии достаточно получить оценку второго и третьего центрального моментов. Они строятся комбинацией из сумм данных в окне, их квадратов и кубов. Ну, корень еще появится. 
Чтобы это делать в скользящем окне, нужно посчитать эти числа для первого окна, а потом на каждом такте вычесть выбывающий элемент и прибавить входящий. 
Объем вычислений будет const* N, где N=10*252*1440
avatar
SergeyJu, нужно чтобы и для 3-4 быстро работало, и для 1440.
Как сделать быстро — можете не писать, я знаю. Я же и поставил доп. условие — без потери точности.
А вот здесь описанная вами схема
Чтобы это делать в скользящем окне, нужно посчитать эти числа для первого окна, а потом на каждом такте вычесть выбывающий элемент и прибавить входящий.
не работает — происходит накопление ошибок порядка sqrt(машинная точность), которая при возведении в куб и последующего деления на числа с такими ошибками — дают вообще непредказуемые по масштабам ошибки в результате.
avatar
MadQuant, ошибка накапливается не так уж и быстро, можно время от времени перезапускаться. И для маленьких окон вообще можно в целой арифметике скользяшки накапливать, это будет абсолютно точно.
avatar
SergeyJu, хотелось бы красивое решение без припарок
ошибка накапливается не так уж и быстро

Вы не поверите — уже даже в окне 3-4 начинаются косяки из-за машинного округления double'ов, т.е. реализация схемы с прибавлением-вычитанием «в лоб» дает сильно плохой результат при расчете skewnessa / kurtosisa
avatar
MadQuant, а сколько значащих цифр в Ваших данных и какая точность нужна? Может, надо просто убрать константу из всех данных или  перейти к приращениям. 
Если все так плохо, надо подумать о других методах, более устойчивых. Не в смысле наращивания разрядности, а в содержательном смысле. Есть же меры асимметрии, допустим, которые не завязаны на центральные моменты. 

avatar
SergeyJu, задача была довльно абстрактная — сделать быструю библиотечную функцию. Но на практике при реализации «в лоб» описанным вами способом с прибавлением/вычитанием это валилось по точности даже при вычислении skew/kurt от дневных ретурнов.
avatar
MadQuant, значит как-то криво заимплементили
Zweroboi, может быть. Хотя нубом себя не считаю, в финал чемпионата мира по программированию выходил и все такое — но как увидел какие погрешности накапливаются даже алгоритмом Кэхэна при расчете skewness'а в малых окнах — у меня чуть глазки не выпали.
avatar
MadQuant, финал чемпионата мира по программированию это конечно респектос, но в коде где-то косяк был 100% ) даже вероятнее всего не в коде, а в подходе к реализации, то есть алгоритм был заимплеменчен правильно, но он был где-то неправильный сам по себе. Я плотно работал с некоторыми элитными перцами из спортивного программирования в G, и наблюдал такое не раз )
MadQuant, нормально работает такая схема для подсчёта дисперсии, для моментов следующих порядков тоже наверняка формулы можно найти. https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
Zweroboi, от потери точности описанные очевидные схемы не спасают.
К сожалению, даже с алгоритмом Кэхэна не удалось получить удобоваримый результат на малых окнах. В итоге пришлось решать задачу на машинном уровне.
avatar
MadQuant, в той же статье на Вики: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics. Всё отлично работает
MadQuant, https://stackoverflow.com/questions/1058813/on-line-iterator-algorithms-for-estimating-statistical-median-mode-skewnes
Zweroboi, в этой статье все алгоритмы довольно наивные, выводятся на бумажке и реализовывались вот точно так же без всякой википедии. Однако всерьез потери точности никто в этой статье и приведенном вами треде не анализирует — а они для реальных рядов ретурнов получаются катастрофические.
avatar
MadQuant, это нормальные стабильные алгоритмы, которые на даблах практически не дают погрешности. Я не считаю скювнесс и куртосис, но дисперсию считаю по Кнуту как в статье, для каждой секунды за годы, это миллионы точек, и расхождение пренебрежимо мало, может в седьмом знаке было когда я это тестил на точность.
Zweroboi, так погрешность 10^-7 при расчете дисперсии — это дофига! У меня погрешности дисперсии порядка 10^-9 — 10^-10, но это не позволяет в итоге точно посчитать skew & kurt. Имхо, вы недооцениваете проблемы, которые там возникают — вы же учтите, что там вычитание есть, а оно адово убивает точность.
avatar
MadQuant, тут конечно нужно конкретно на код смотреть — как добавляется новая точка данных, и как выбрасывается самая старая в окне — но я гарантирую, что эти все моменты считаются с прекрасной и абсолютно практически достаточной точностью, с константной сложностью на точку данных.
SergeyJu, здесь N это окно по которому считают. И квадратичная сложность присутствует в первом варианте. Мы же там два раза по окну бегаем.
avatar
b34rcava1ry, похоже, мы с Вами о чем-то разном 
avatar
SergeyJu, 
Трудоемкость будет N*расчет одного окна. Нет ни квадратичного никакого иного роста сложности в ростом N.
Я вот об этом. Если считать все целиком то буде O((n-m)*m^2). Но так сложность никто не считает. Нас интересует только сложность расчетов стддев, по этому мы отбрасываем первый множитель.
avatar
Существует весьма большое количество методов ускорения расчетов. Мне кажется, если человек получает образование, связанное с компьютерными вычислениями, он должен знать где что найти и уметь пользоваться методами ускоренных вычислений. 
avatar
SergeyJu, дело же не в том, что катит, а что нет. Результаты то одни и те же, а значит оба решения верные. Но иногда бывают требования к производительности например.
если человек получает образование, связанное с компьютерными вычислениями, он должен знать
как оказалось это совсем не обязательно)
avatar
b34rcava1ry, я отвечал не Вам, а на задачу с маленьким окном. Ваш-то расчет перехода от квадратичной сложности задачи к линейной является правильным, я бы сказал, стандартным приемом. Существует много таких приемов, например, сортировка методом «пузырька» имеет квадратичную сложность, но существует схема сложности N*log(n). Аналогично при  использовании БПФ для расчета дискретного преобразования Фурье или расчетная схема вейвлетов.
avatar
b34rcava1ry, думаю, кто реально парится производительностью(в смысле скорости), смотрит на питон как на говно:)
А то что Вы говорили, что типа он тут дает какое то большое преимущество в читабельности, это навяд ли, цЫклы они и в Африке цЫклы:)
avatar
sortarray sortarray, меньше; {} очень сильно решает. Да и на плюсах далеко не сразу была бы видна разница в перфомансе)
плюс этот код может взять любой и быстро у себя запустить и проверить.
avatar
b34rcava1ry, а что на плюсах мешает быстро запустить и проверить?
avatar
b34rcava1ry, а вы с++ знаете?
avatar
kvazar, им родимым и живу, и иногда verilog'ом.
avatar
b34rcava1ry, 
 иногда verilog'ом
на какой железяке?
avatar
Андрей К, вот на такой.
avatar
b34rcava1ry, спасибо, жаль без сетевых портов
avatar
b34rcava1ry, а вы на нее не навешивали расширение с сетевыми портами? там по ссылке есть такое
avatar
Андрей К, 4х портовая у меня)
avatar
b34rcava1ry, спасибо
avatar
Андрей К, тоже хотите свою сетевуху с плавующей средней и Герчиком?
avatar
b34rcava1ry, да я тут мало что понял, просто интересны возможности верилога
avatar
Андрей К, ну у него безграничные возможности, правда нужно знать цифровую и аналоговую схемоту.
avatar
b34rcava1ry, перечитал еще раз вашу беседу. Это вы все вышеизложенное на verilog переложили?
avatar
Андрей К, почти. На этой карте делается сетевуха с потоковой обработкой FIX/FAST и kernel bypass драйверами.
avatar
b34rcava1ry, то есть все таки основную математику вы вынесли все таки во внешнее приложение? раз используете fix/fast + bypass. То есть такие карты не тянут такие мощности?
avatar
Андрей К, там есть такой момент, что чип то не резиновый. Это довольно старая карта и на ней логика сетевухи и проца для обработки протокола заняли 2/3 места. С одной стороны там еще многое можно уместить, с другой стороны туда не всякую математику стоит переносить в виду того что эта FPGA в лучшем случае работает на 300mhz, а CPU сейчас 3ghz.
avatar
b34rcava1ry, спасибо.
avatar
Андрей К, фактически чтобы получилось действительно быстрее, нужно чтобы алгоритм, который переносится на ПЛИС отрабатывал на ней быстрее чем время 2ух транзакций по PCI. Как-то так.
avatar
b34rcava1ry, ну да, логично
avatar
Что это было?
avatar
Есть один недостаток у программистов — они считают себя Богами.
avatar
Честно я не понял. Зачем среднеквадратичное отклонение или ошибка нужна для нахождения тренда это типа по леме 3 сигм?

И по первой формуле среднеквадратичного отклонения она работает только при n<40 при n>40 разность квадратов. 

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

avatar
Jkrsss, да просто похожая задачка. Цель была показать возможности оптимизации.
avatar
Почитал вашу беседу. Если не секрет, что нужно торговать, чтобы применить такие знания?
avatar
Андрей К, думаю синусоиду или белый шум :)
avatar
SECRET, да это был диалог квантов. Примерно наверное это и торгуют.
MadQuant иногда пишет. Я так понял корзины инструментов прогнанные на истории в 15 лет
avatar
Андрей К, это не диалог квантов. Это диалог программистов, ошибочно считающих себя квантами :-)
avatar
Вдумайтесь!!! какая потеря точности? какие объемы и сложности вычислений? вы далеко не число пи считаете до 57854 знака после запятой. Вы считаете слабоприменимые и низкоэффективные для трейдинга величины по громоздким ресурсоемким формулам. Может стоит переосмыслить задачу и появится простое красивое решение? ;)
avatar
SECRET, мне кажется, что обсуждается деятельность не практического спекулянта, а программиста методов вычислений. Видит Бог, ни точность до 10 десятичного знака, ни моменты порядка выше 2 для торговли нах не нужны. Но входят в инструментарий просто потому что есть в наличии. 
avatar

теги блога b34rcava1ry

....все тэги



UPDONW
Новый дизайн