06.12.2023 08:05
SmartEngines 9 1400 Источник

Привет, Хабр! Иногда сидишь, решаешь задачу, и, в процессе решения, чтобы продвинуться на следующий шаг, нужно придумать как сделать что-то очень простое - ну, то что наверняка уже делалось тысячи раз другими людьми. Кинувшись в поисковик перелопачиваешь какое-то количество литературы и вдруг понимаешь что либо ты просто искать не умеешь, либо это действительно никто до тебя не делал, или делал но об этом не писал. В какой-то момент проще просто взять и решить задачу самому…

В этой заметке мы расскажем об одной такой задаче - простой, но которая нам понадобилась для кое-чего другого. Задача - придумать, как при увеличении наблюдаемой выборки быстро пересчитать ее абсолютный центральный момент.

Конкретная задача, которую мы решали, когда наткнулись на эту проблему - не очень существенна - мы пытались ускорить моделирование следующего комбинированного результата распознавания текстовой строки в видеопотоке, которое нам нужно для принятия решение об остановке (мы про это писали раньше на хабре). В процессе разработки метода (описанного, кстати, вот в этом докладе) мы наткнулись на задачу, которая один в один совпадала с задачей быстрого обновления абсолютного центрального момента первого порядка по некоторой выборке. И, к нашему удивлению, не смогли нагуглить, как это обычно делают. На самом деле делать это не трудно, и мы решили рассказать об этом здесь на Хабре, в более общем виде, на случай, если подобная проблема когда-то у кого-нибудь еще возникнет.

Итак, для начала нужно рассказать, что такое моменты выборок и какие они бывают. Для выборки A = \{a_1, ..., a_n\} можно определить два вида моментов k-го порядка:

\text{Начальный:}\quad \nu_k (A) = \frac{1}{n}\sum^n_{i=1} a_i^k,\text{Центральный:}\quad \mu_k (A) = \frac{1}{n}\sum^n_{i=1} \left(a_i - \nu_1 (A) \right)^k.

Для каждого из них так же существуют абсолютные версии:

\text{Абсолютный начальный:}\quad \widetilde{\nu}_k (A) = \frac{1}{n}\sum^n_{i=1} \left|a_i \right|^k,\text{Абсолютный центральный:}\quad \widetilde{\mu}_k (A) = \frac{1}{n}\sum^n_{i=1} \left|a_i - \nu_1 (A) \right|^k.

Первый начальный момент \nu_1(A) - математическое ожидание по выборке (или просто среднее значение по выборке), а 2-й центральный момент \mu_2(A) также называется дисперсией выборки. Очевидно также, что абсолютные моменты четных порядков совпадают с соответствующими обычными моментами тех же порядков.

Задача быстрого обновления моментов выборки формулируется так: пусть наша выборка пополняется по одному элементу за один раз. На n-м шаге процесса существует некоторая выборкаA = \{a_1, ..., a_{n-1}\} размера n-1, и этой выборке добавляется очередной элемент a_n. После добавления a_n в выборку нужно пересчитать текущее значение того или иного момента выборки как можно быстрее, используя какие-нибудь аккумуляторы.

Из четырех типов моментов три легко и быстро обновляются. Действительно, чтобы обновлять начальный или абсолютный начальный моменты порядка k достаточно аккумулировать сумму k-х степеней элементов выборки или их модулей, соответственно. Чтобы понять, как быстро обновлять центральный момент порядка k, можно воспользоваться вот этой формулой, напрямую следующей из разложения формулы центрального момента с помощью бинома Ньютона:

\mu_k(A)=\sum_{j=0}^k C_k^j \left( \frac1n \sum_{i=1}^n a_i^{k-j} \right)(-\nu_1(A))^j.

Для быстрого обновления центрального момента, согласно этой формуле, нам достаточно уметь быстро обновлять первый начальный момент (что мы и так умеем), и аккумулировать суммы всех степеней элементов выборки до k-й включительно.

Что же делать с абсолютным центральным моментом? Если порядок четный, то он совпадает с обычным центральным моментом, это понятно. А если порядок нечетный?..

Давайте мысленно разобьем нашу выборку A на две части: L и R. Часть L=\{a\in A | a < \nu_1(A)\} - это элементы выборки со значениями, меньшими среднего, а часть R=\{a\in A | a \ge \nu_1(A)\} - элементы выборки со значениями, большими и равными среднему. Сумма модулей, возникающая в определении абсолютного центрального момента, может теперь быть раскрыта как две суммы без модулей - отдельно по части L, для которой модули раскрываются "со знаком минус" и по части R, для которой они раскрываются "со знаком плюс":

\widetilde{\mu}_k(A)=\frac1n \left( \sum_{l\in L} (\nu_1(A)-l)^k + \sum_{r\in R} (r-\nu_1(A))^k \right).

Раскроем эту формулу с помощью бинома Ньютона:

\widetilde{\mu}_k(A)=\frac1n \left( \sum_{l\in L} \sum_{j=0}^k (\nu_1(A))^{k-j}(-1)^kl^jC_k^j + \sum_{r\in R}\sum_{j=0}^kr^{k-j}(-1)^j(\nu_1(A))^jC_k^j \right).

И, аналогично формуле для обычного центрального момента, вынесем за скобки и за знаки суммы все, что можем:

\widetilde{\mu}_k(A)=\frac1n \left( \sum_{j=0}^k(-1)^jC_k^j \left( (\nu_1(A))^j\sum_{r\in R}r^{k-j}+(\nu_1(A))^{k-j}\sum_{l\in L} l^j \right) \right).

Как теперь видно, чтобы быстро пересчитывать абсолютный центральный момент, нужно уметь быстро пересчитывать первый начальный (это мы умеем), и аккумулировать суммы всех степеней элементов выборки до k-й включительно, но не по всей выборке, а отдельно по подмножествам L и R. Учитывая что мы и так умеем аккумулировать все степени всех элементов всей выборки A, то достаточно только придумать, как аккумулировать суммы по подмножеству L, а оставшееся вычислять простым вычитанием из общей суммы по выборке.

Аккумулировать суммы всех нужных нам степеней элементов множества L легко, если хранить выборку A в виде сбалансированного дерева поиска (к примеру, АВЛ-дерева, декартова дерева, или любого другого). Пусть каждому элементу выборки A соответствует вершина дерева поиска, а также в вершине содержатся суммы всех нужных нам степеней всех элементов в поддереве с корнем в этой вершине. Тогда в таком дереве можно реализовать ответ на вопрос вида "какова сумма j-х степеней всех элементов выборки A, меньших заданного значения x?" (в качестве x мы всегда будем использовать \nu_1(A)). Действительно, чтобы ответить на такой вопрос, достаточно спуститься по дереву в поиске значения x и каждый раз при переходе в правое поддерево добавлять к ответу j-ю степень текущей вершины и ее левого поддерева. Трудоемкость операции - O(\log n).

Трудоемкость добавления нового элемента в сбалансированное дерево будет составлять O(k\log n) - логарифм на перестроение структуры, и мультипликатор k, поскольку при каждой операции над деревом придется пересчитывать k аккумуляторов в нескольких вершинах. За эту же трудоемкость мы можем получить значения всех аккумуляторов, которые используются в нашей формуле абсолютного центрального момента, а трудоемкость всех остальных операций для ее вычисления не превышаетO(k).

Таким образом, у нас получилось обновить k-й абсолютный центральный момент выборки размера n за время O(k\log n), хотя для этого нам пришлось содержать сбалансированное дерево с n вершинами и с аккумулятором размера O(k) в каждой вершине.

Пару слов про то, зачем это нам вообще понадобилось. Ранее мы писали на хабре, что нас интересует задача поиска оптимального момента остановки - когда принимать решение о том, что захват кадров при распознавании в видеопотоке можно прекратить, и текущий результат принять как окончательный. Алгоритм, которым мы пользовались, состоял в моделировании возможного следующего наблюдения (т.е. следующего комбинированного результата распознавания) и в оценке расстояния от текущего результата до возможного следующего. Однако моделирование следующего результата, которое мы использовали, было достаточно накладным - оно нас устраивало теоретически, но практически требовало слишком много вычислений, да еще и это количество вычислений линейно зависело от количества уже просмотренных кадров. Мы придумали (и описали в этом докладе) более дешевый способ такого моделирования, и для обновления значения расстояния от текущего результата распознавания до моделируемого нам пришлось решать задачу, один в один совпадающую с задачей обновления 1-го абсолютного центрального момента. Решив ее так, как мы описали выше, с использованием декартовых деревьев в качестве поисковой структуры данных, мы получили способ принятия решения об остановке процесса распознавания, который вычисляет свое решение в десятки раз быстрее, чем то, что мы использовали до этого.

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

Комментарии (9)


  1. adeshere
    06.12.2023 08:53
    #26235722
    +4

    Спасибо за подробный разбор алгоритма! Но вообще-то почти в точности такая же технология используется уже несколько десятилетий, и называется "обработка в скользящем окне". Например, мы у себя (при работе с временными рядами) используем подобный подход с 1992 года. При добавлении нового элемента в выборку все суммы инкрементируются, при удалении элемента - декрементируются. При этом как поступающий в выборку, так и удаляемый из нее элемент в нашем случае может быть так называемым пропуском, который кодируется Nan или другим особым значением. Поэтому перед инкрементом/декрементом сумм выполняется дополнительная проверка на Nan.

    А вот деревья мы не используем, так как абсолютные моменты бывают нужны

    очень редко,

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

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

    Ну и еще одно отличие от описанного в статье решения состоит в том, что нам приходится учитывать ограниченную точность FP-вычислений, которая становится особенно критичной при расчете высших моментов. А ряды длинные (миллиарды значений данных). Поэтому нам приходится принимать специальные меры для предотвращения потери точности:

    Во-первых, из данных перед расчетом сумм удаляется среднее. Заранее мы его, конечно, не знаем, но хотя бы примерно оценить можем. В нашем случае ожидаемое среднее берется

    из паспорта сигнала

    тут нам повезло, у нас в программе система анализа данных интегрирована с СУБД. Но думаю, что в большинстве практических случаев какие-то ожидания по поводу среднего значения сделать обычно можно

    Это позволяет избежать вычитания двух больших близких чисел, и очень существенно повышает точность.

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

    Полные тексты статей

    которые упомянуты выше, можно найти вот тут.

    а сама программа

    в виде exe-шников доступна для скачивания вот здесь,

    А ее исходники и файлы проектов - вот тут.


    1. SmartEngines Автор
      06.12.2023 08:53
      #26236320
      +1

      Спасибо за ваш интерес!

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

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

      Кстати, по поводу вычислительной стабильности, работа 2016-го года, на которую мы ссылаемся в тексте (https://link.springer.com/article/10.1007/s00180-015-0637-z) о быстром обновлении центральных моментов (не абсолютных) почти целиком посвящена проблемам накопления ошибки вычислений с плавающей точкой. Этот вопрос там хорошо и подробно рассмотрен, рекомендуем.


    1. avdx
      06.12.2023 08:53
      #26236646
      +1

      А откуда берутся значения с плавающей точкой? Если речь о данных из "реального мира", то обычно они приходят с АЦП, т.е. по факту целочисленны. В этом случае можно все накопление/обновление сумм выполнять в целочисленной арифметике и никаких проблем с точностью не будет.


      1. adeshere
        06.12.2023 08:53
        #26236962

        А откуда берутся значения с плавающей точкой? 

        Это от задачи зависит.

        Например, у нас в геофизическом мониторинге основная задача - это изучение природных явлений и закономерностей, а не

        управление ими

        То есть, речь чаще всего идет о поиске каких-то эффектов или проверке гипотез по ретроспективным данным. Обычно для этого анализируются сигналы, накопленные за много лет. А они копятся в физических единицах, естественно - так как разряды АЦП после каждой перенастройки аппаратуры разные, а ряд должен быть однородным. Не говоря уже о том, что для хранения коэффициентов перевода (если в базе лежит выходной сигнал АЦП) потребуется отдельная нехилая БД...

        Потенциально задачи реального времени в геофизике тоже есть. Например, чтобы спрогнозировать землетрясение в реальном времени, нужен онлайн-мониторинг (измерения в реальном времени) + четкий алгоритм прогноза. Но у нас пока что такого алгоритма готового нет. Мы его только лишь ищем и пытаемся построить, используя архивные ряды данных. Фактически на данный момент более-менее надежные алгоритмы есть только для оценки сейсмической опасности (фактически это прогноз на бесконечное время). Долгосрочный (на многие годы) или среднесрочный прогноз - это все еще полуэкспериментальный режим, и там обязательно на каком-то этапе экспертная оценка используется. Поэтому обработка данных для таких целей идет не в реальном времени, а в пакетном режиме (порциями), и с некоторой задержкой всегда. А это уже совсем другой стиль работы с сигналами.

        Конечно, бывают особые случаи вроде прогноза цунами по факту возникновения сильного землетрясения... но это не моя область.

        И в таких задачах без человека не обойтись. А везде, где участвует человек, гораздо удобнее работать в физических единицах - градусах там или нанострейнах, а не в условных "отсчетах АЦП"

        Ну и второй нюанс - наши измерительные системы обычно стоят в труднодоступных местах, и данные оттуда скачиваются порциями. Обычно это раз в сутки или даже реже. Потом они проходят кучу этапов первичной обработки, в сигнал вводятся всякие коэффициенты и поправки, включая пересчет в физические единицы. И уже в таком виде заносятся в базу данных. Чтобы абстрагировать все эти уровни от конечного пользователя данных.

        Соответственно, данные для всяких расчетов берутся не с АЦП, а из БД, где они лежат уже в виде значений с плавающей точкой.

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


        1. avdx
          06.12.2023 08:53
          #26237114
          +1

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


          1. adeshere
            06.12.2023 08:53
            #26237414
            +1

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

            Это понятная схема, но она не всегда применима на практике, так как плотность real-шкалы неравномерна, в отличие от integer. И очень часто это существенно. Так что проблема вовсе не сводится к ошибке квантования (в обычном ее понимании).

            Например, в экспериментальных рядах очень часто есть выбросы, которые можно сохранять с худшей точностью, что прекрасно реализует real-формат. Если же разбить весь диапазон на равноотстоящие уровни (что и произойдет при преобразовании в integer), то точность записи основной массы значений резко упадет. Проблема тут в том, что эти выбросы обычно имеют самую разную амплитуду: от гигантских и до пограничных с "нормой", что не позволяет отбраковать их все

            еще на предварительной стадии (до записи в базу).

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

            под каким углом мы его изучаем

            У нас же почти каждый сигнал - это плод многолетних усилий целого коллектива. Он измеряется один раз, и повторить эти измерения уже невозможно. Поэтому мы потом эти ряды рассматриваем буквально "под микроскопом", многократно возвращаясь к одним и тем же данным

            Иногда для идентификации выбросов достаточно тонкие алгоритмы приходится применять.

            И еще один ключевой момент.

            Если Вы продумали, что в результате хитрых трюков первичной обработки у нас в итоге все равно сохраняется примерно однозначное соответствие между исходными int-значениями с выхода АЦП и "физическими" real-значениями, то

            это совершенно не так

            При геофизическом мониторинге длина рядов может достигать десятков лет. За это время аппаратура и установка могут меняться неоднократно. А уж замена/регулировка усилителя (АЦП) - это и вовсе рутинная процедура. Например, когда наклон в результате дрейфа уплывает за границу диапазона, мы меняем смещение (электронно или механически), возвращая сигнал обратно к центру шкалы. После этого все формулы пересчета в физ. единицы меняются, и тот же самый отсчет АЦП будет соответствовать совершенно другой физической величине. Другая типичная ситуация - это коррекция коэффициента усиления. Мы же не знаем заранее, как себя поведет сигнал с уникального датчика в уникальных условиях. Поэтому время от времени может возникнуть необходимость расширить динамический диапазон (если сигналу его не хватает) или, наоборот, сузить (если сигнал упорно занимает лишь небольшую часть полной шкалы). Это тоже приводит к пересчету всех формул.

            Из-за этого существенная потеря информации при ре-квантовании в общем случае неизбежна. Поэтому даже в тех случаях, когда вроде бы можно было бы вернуться к целым значениям,

            мы этого не делаем

            И причина тут даже не одна, а все две:

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

            2) Ну и во-вторых, мы свои алгоритмы модифицируем иногда. И сделать это гораздо проще, когда алгоритм существует в одном варианте, чем когда существует зоопарк версий...

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

            В общем, спасибо Вам за желание подкинуть конструктивный совет - кому-то он здесь наверняка пригодится (ведь случайные люди такую статью не полезут читать ;-).

            Но вот именно нам - увы...


  1. wataru
    06.12.2023 08:53
    #26237456
    +1

    Таким образом, у нас получилось обновить k-й абсолютный центральный момент выборки размера n за время O(k\log n), хотя для этого нам пришлось содержать сбалансированное дерево с n вершинами и с аккумулятором размера O(k) в каждой вершине.

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

    Идея в том, чтобы поддерживать 2 дерева/множества - одно для чисел меньше-равных среднего, второе - для чисел больше. Все операции разбиваются на 2 части: сначала добавляем новые числа и пересчитываем среднее, потом изменяем деревья, чтобы числа были на своих местах.

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

    O(k) аккумуляторов тут надо всего 2 - для левого и правого множества целиком. Пересчитываем их при перемещении элементов между множествами.

    Недостаток этого метода, что тут не так просто удалять числа из набора. Ибо они могут быть в середине хипа, а не на вершине. Но это можно реализовать, если хранить не heap, а балансированное дерево, но опять же, стандартное - без аккумуляторов.


    1. amishaa
      06.12.2023 08:53
      #26237926

      Только может ттак получиться, что много элементов каждый раз будут менять heap. Например, если данные вот такие: 1, 0, 0, -1.1, 1.1, 0, 0, -1.1, 1.1, .... то среднее всё время скачет от "чуть больше нуля" до "чуть меньше нуля" и все нули нужно перекладывать между хипами.


      1. wataru
        06.12.2023 08:53
        #26237944

        Да, действительно.