b34rcava1ry
b34rcava1ry личный блог
09 августа 2017, 17:05

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

Обычно я просто читаю СМ и не лезу со своим «экспертным» мнением к людям. Но недавно я не выдержал и закусился на форуме, да еще и побил собственный рекорд, ушел в личный бан за 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млрд значений. И тут у нормального человека должны появиться сомнения в смысле такой величины,
В комментариях предлагайте темы, возможно я смогу рассказать что-то интересное по ним в следующих сериях рубрики «Программирование для самых маленьких и тупых».




65 Комментариев
  • MadQuant
    09 августа 2017, 17:14
    А вот интересует такая вполне прикладная задача: надо мне посчитать скользящий skewness и kurtosis также за O(N), но без потери точности на длинных рядах (например, минутки или пятиминутки за 10 лет), и чтобы реально без потери точности, т.е. расчеты совпадали с расчетом «в лоб» в том же окне. Причем интересна именно ошибка на маленьких окнах (3-10 наблюдений).
      • Пафос Респектыч
        09 августа 2017, 17:46
        b34rcava1ry, среднеквадратичное отклонение так не считают, в лоб — вычислительно неустойчивый алгоритм. Надо по Кнуту — так и быстрее, и точнее, и для окон лучше подходит — можно чиселки как добавлять, так и выбрасывать. http://nabbla1.livejournal.com/83886.html
        • SergeyJu
          09 августа 2017, 18:09
          Zweroboi, Кнут — это очень хорошо и правильно. Но ссылка почему-то не идет. Впрочем, всегда можно открыть книгу :)
  • SergeyJu
    09 августа 2017, 17:21
    А почему «в лоб» не катит? Трудоемкость будет N*расчет одного окна. Нет ни квадратичного никакого иного роста сложности в ростом N.
    • MadQuant
      09 августа 2017, 17:35
      SergeyJu, Так O(N * Window) — это уже большая сложность. Например, если минутки, и мне надо посчитать скользящий skewness в окне сутки за 10 лет — это грубо 1440 * (10 * 252 * 1440) ~ 5.2 * 10^9 — как бы уже несколько секунд расчетов, предполагая расчеты за 1 такт и мгновенное обращение к памяти. А там еще сбои кэшей,… Да и контрактов не 1, а 50, например — и вот, мы уже сидим ждем этого счастья несколько минут.

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

        • MadQuant
          09 августа 2017, 17:48
          SergeyJu, 

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

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

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

                    • MadQuant
                      09 августа 2017, 18:31
                      SergeyJu, задача была довльно абстрактная — сделать быструю библиотечную функцию. Но на практике при реализации «в лоб» описанным вами способом с прибавлением/вычитанием это валилось по точности даже при вычислении skew/kurt от дневных ретурнов.
                      • Пафос Респектыч
                        09 августа 2017, 19:35
                        MadQuant, значит как-то криво заимплементили
                        • MadQuant
                          09 августа 2017, 19:48
                          Zweroboi, может быть. Хотя нубом себя не считаю, в финал чемпионата мира по программированию выходил и все такое — но как увидел какие погрешности накапливаются даже алгоритмом Кэхэна при расчете skewness'а в малых окнах — у меня чуть глазки не выпали.
                          • Пафос Респектыч
                            09 августа 2017, 20:07
                            MadQuant, финал чемпионата мира по программированию это конечно респектос, но в коде где-то косяк был 100% ) даже вероятнее всего не в коде, а в подходе к реализации, то есть алгоритм был заимплеменчен правильно, но он был где-то неправильный сам по себе. Я плотно работал с некоторыми элитными перцами из спортивного программирования в G, и наблюдал такое не раз )
              • Пафос Респектыч
                09 августа 2017, 18:19
                MadQuant, нормально работает такая схема для подсчёта дисперсии, для моментов следующих порядков тоже наверняка формулы можно найти. https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
                • MadQuant
                  09 августа 2017, 18:26
                  Zweroboi, от потери точности описанные очевидные схемы не спасают.
                  К сожалению, даже с алгоритмом Кэхэна не удалось получить удобоваримый результат на малых окнах. В итоге пришлось решать задачу на машинном уровне.
                  • Пафос Респектыч
                    09 августа 2017, 19:32
                    MadQuant, в той же статье на Вики: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics. Всё отлично работает
                  • Пафос Респектыч
                    09 августа 2017, 19:34
                    MadQuant, https://stackoverflow.com/questions/1058813/on-line-iterator-algorithms-for-estimating-statistical-median-mode-skewnes
                    • MadQuant
                      09 августа 2017, 19:46
                      Zweroboi, в этой статье все алгоритмы довольно наивные, выводятся на бумажке и реализовывались вот точно так же без всякой википедии. Однако всерьез потери точности никто в этой статье и приведенном вами треде не анализирует — а они для реальных рядов ретурнов получаются катастрофические.
                      • Пафос Респектыч
                        09 августа 2017, 19:51
                        MadQuant, это нормальные стабильные алгоритмы, которые на даблах практически не дают погрешности. Я не считаю скювнесс и куртосис, но дисперсию считаю по Кнуту как в статье, для каждой секунды за годы, это миллионы точек, и расхождение пренебрежимо мало, может в седьмом знаке было когда я это тестил на точность.
                        • MadQuant
                          09 августа 2017, 21:27
                          Zweroboi, так погрешность 10^-7 при расчете дисперсии — это дофига! У меня погрешности дисперсии порядка 10^-9 — 10^-10, но это не позволяет в итоге точно посчитать skew & kurt. Имхо, вы недооцениваете проблемы, которые там возникают — вы же учтите, что там вычитание есть, а оно адово убивает точность.
                      • Пафос Респектыч
                        09 августа 2017, 19:56
                        MadQuant, тут конечно нужно конкретно на код смотреть — как добавляется новая точка данных, и как выбрасывается самая старая в окне — но я гарантирую, что эти все моменты считаются с прекрасной и абсолютно практически достаточной точностью, с константной сложностью на точку данных.
      • SergeyJu
        09 августа 2017, 18:17
        b34rcava1ry, похоже, мы с Вами о чем-то разном 
  • SergeyJu
    09 августа 2017, 17:27
    Существует весьма большое количество методов ускорения расчетов. Мне кажется, если человек получает образование, связанное с компьютерными вычислениями, он должен знать где что найти и уметь пользоваться методами ускоренных вычислений. 
      • SergeyJu
        09 августа 2017, 17:37
        b34rcava1ry, я отвечал не Вам, а на задачу с маленьким окном. Ваш-то расчет перехода от квадратичной сложности задачи к линейной является правильным, я бы сказал, стандартным приемом. Существует много таких приемов, например, сортировка методом «пузырька» имеет квадратичную сложность, но существует схема сложности N*log(n). Аналогично при  использовании БПФ для расчета дискретного преобразования Фурье или расчетная схема вейвлетов.
      • sortarray sortarray
        09 августа 2017, 17:42
        b34rcava1ry, думаю, кто реально парится производительностью(в смысле скорости), смотрит на питон как на говно:)
        А то что Вы говорили, что типа он тут дает какое то большое преимущество в читабельности, это навяд ли, цЫклы они и в Африке цЫклы:)
          • Displacer
            09 августа 2017, 20:26
            b34rcava1ry, а что на плюсах мешает быстро запустить и проверить?
      • Трейдер Вася
        09 августа 2017, 17:50
        b34rcava1ry, а вы с++ знаете?
          • Андрей К
            09 августа 2017, 23:34
            b34rcava1ry, 
             иногда verilog'ом
            на какой железяке?
              • Андрей К
                10 августа 2017, 10:21
                b34rcava1ry, спасибо, жаль без сетевых портов
              • Андрей К
                10 августа 2017, 10:23
                b34rcava1ry, а вы на нее не навешивали расширение с сетевыми портами? там по ссылке есть такое
                  • Андрей К
                    10 августа 2017, 10:32
                    b34rcava1ry, спасибо
                      • Андрей К
                        10 августа 2017, 10:59
                        b34rcava1ry, да я тут мало что понял, просто интересны возможности верилога
                  • Андрей К
                    10 августа 2017, 10:34
                    b34rcava1ry, перечитал еще раз вашу беседу. Это вы все вышеизложенное на verilog переложили?
                      • Андрей К
                        10 августа 2017, 11:06
                        b34rcava1ry, то есть все таки основную математику вы вынесли все таки во внешнее приложение? раз используете fix/fast + bypass. То есть такие карты не тянут такие мощности?
                          • Андрей К
                            10 августа 2017, 11:17
                            b34rcava1ry, спасибо.
                              • Андрей К
                                10 августа 2017, 11:24
                                b34rcava1ry, ну да, логично
  • Alexey
    09 августа 2017, 18:52
    Что это было?
  • Чарльз Маккей
    09 августа 2017, 20:04
    Есть один недостаток у программистов — они считают себя Богами.
  • Jkrsss
    09 августа 2017, 20:28
    Честно я не понял. Зачем среднеквадратичное отклонение или ошибка нужна для нахождения тренда это типа по леме 3 сигм?

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

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

  • Андрей К
    09 августа 2017, 23:09
    Почитал вашу беседу. Если не секрет, что нужно торговать, чтобы применить такие знания?
    • SECRET
      09 августа 2017, 23:24
      Андрей К, думаю синусоиду или белый шум :)
      • Андрей К
        09 августа 2017, 23:33
        SECRET, да это был диалог квантов. Примерно наверное это и торгуют.
        MadQuant иногда пишет. Я так понял корзины инструментов прогнанные на истории в 15 лет
        • Евгений
          10 августа 2017, 02:29
          Андрей К, это не диалог квантов. Это диалог программистов, ошибочно считающих себя квантами :-)
  • SECRET
    09 августа 2017, 23:10
    Вдумайтесь!!! какая потеря точности? какие объемы и сложности вычислений? вы далеко не число пи считаете до 57854 знака после запятой. Вы считаете слабоприменимые и низкоэффективные для трейдинга величины по громоздким ресурсоемким формулам. Может стоит переосмыслить задачу и появится простое красивое решение? ;)
    • SergeyJu
      10 августа 2017, 11:44
      SECRET, мне кажется, что обсуждается деятельность не практического спекулянта, а программиста методов вычислений. Видит Бог, ни точность до 10 десятичного знака, ни моменты порядка выше 2 для торговли нах не нужны. Но входят в инструментарий просто потому что есть в наличии. 

Активные форумы
Что сейчас обсуждают

Старый дизайн
Старый
дизайн