20 августа 2018 г.

Оценки квантилей распределений потока данных


Последнюю неделю копаю алгоритмы оценок для квантилей выборок. Началось все с префетчинга: есть хранилище данных на диске, часть из них кэшируется в памяти, и надо бы успевать кэшировать предварительно с такой скоростью, чтобы клиентам не приходилось ожидать ввода-вывода. Доступ последовательный, так что все сводится к вопросу "сколько страниц загрузить в текущем раунде упреждающего чтения". И обычно в таких случаях я вижу что-то вроде
int estimatePagesToPrefetch(...){
   return 4;//carefully chosen
}
Иногда (реже) что-то вроде
int estimatePagesToPrefetch(int pagesReadSinceLastPrefetch){
   return ema( pagesReadSinceLastPrefetch ) * 3/2; //150% is enough for everyone
}
И я сам так пишу, но в этот раз душа физика-теоретика захотела бОльшего: можно ли обойтись без с потолка взятых констант? Проблема с константами не только в том, что они непонятно откуда взялись (хотя многие да). Еще одна проблема, что непонятно даже, из каких соображений они взялись. Вот эти 4 страницы, или 150% от средней скорости использования — какого внешнего результата автор от них ожидал? По какому критерию он их выбрал (если вообще выбирал)?

Ну и я спросил себя: а по какому критерию я бы выбирал параметры префетчинга? Мне кажется, префетчинг должен с заданной вероятностью (90%, 95%, 99%...) опережать использование. Это критерий, понятный человеку, он легко пересчитывается в другие latency-oriented критерии. И чтобы его реализовать, нужно оценить соответствующий квантиль распределения величин pagesReadSinceLastPrefetch.

И я начал искать какой-нибудь простой и быстрый способ оценить n-ый квантиль распределения по потоку сэмплов из него. BigData сейчас в моде, так что статей обнаружилось много, но алгоритмы меня не порадовали:
  1. Cложные — не напишешь за 5 минут, если тебе вдруг понадобилось оценить 0.95-й квантиль.
  2. Cами по себе не приспосабливаются к нестационарному распределению: если распределение данных меняется со временем (а в реальных сценариях это почти всегда так), то приходится делать дополнительные усложнения
Но потом я наткнулся статью Frugal Streaming for Estimating Quantiles. И в ней излагаются два совершенно замечательных, очень экономных стохастических алгоритма оценки квантилей потока данных.

Замечательные они вот почему:
  1. Очень бережливо (frugal) относятся к памяти. Frugal1U использует для хранения состояния одну-единственную ячейку памяти (!), Frugal2U использует 2 ячейки + 1 бит. Строго говоря есть еще один конфигурационный параметр, да нужен генератор случайных чисел, у которого есть состояние, и неочевидно, считать ли это все "состоянием алгоритма", или относить к "общей памяти программы". Но даже в самом строгом варианте получается 3-4 ячейки памяти.
  2. Очень простые, как по количеству инструкций, так и когнитивно. Наверное, самое сложное, это поверить, что это вообще может работать, и понять, почему это может довольно хорошо работать. (Я потратил несколько вечеров, чтобы разобраться, почему они работают, и пару недель, пытаясь придумать что-то лучшее — и все еще особо не преуспел).
  3. Автоматически, безо всяких дополнительных усложнений, адаптируются к нестационарностям распределения
  4. Не используют никаких априорных предположений о виде распределения данных в потоке — только само определение квантиля. С одной стороны это недостаток: вкупе с ограничением по используемой памяти получается не очень быстрая сходимость и довольно большая погрешность. А с другой это дает очень высокую устойчивость: "смутить" алгоритм каким-то сложным распределением практически невозможно. Они как в том анекдоте: медленно спускаются с горы, и делают свою работу
Как я понял, авторы разрабатывали эти алгоритмы для чего-то вроде встраиваемых сенсоров, где память и вычислительная мощность очень ограничены, а поток данных может быть большим и зашумленным, и его нужно как-то предобработать, прежде чем выдавать наружу. Но мне кажется, что таким простым и устойчивым алгоритмам найдется применение и много где еще — хотя бы и для построения регуляторов/петель обратной связи.

Frugal1U

Идея алгоритма настолько проста, что кажется, тебя разыгрывают: пусть нужно оценить 50% квантиль (медиану). Медиана — это такое X, что элементы (сэмплы) из потока будут равновероятно оказываться как справа, так и слева от X. Если относительно X сэмплы с одной стороны выпадают чаще, чем с другой — значит X это не медиана, а настоящая медиана лежит как раз с той стороны, куда сэмплы попадают чаще. Так возьмем любое X за начальное приближение, выберем небольшой шаг step, и будем на каждый приходящий сэмпл обновлять X[i+1] -> X[i] +/- step, смотря в какую сторону от X попал текущий сэмпл. Статистически мы будем чаще делать шаг к медиане, чем от нее, так что рано или поздно мы до нее дошагаем, и дальше будем топтаться около.
step = 1
X    = 0   # initial estimation of median
for x in stream:
    if x > X:      
        X += step
    elif x < X: 
        X -= step
Ок, с медианой понятно (понятно же?). А если я хочу 0.9-квантиль? Значение 0.9 квантиля, X0.9, это такая величина, что элементы из потока данных будут с вероятностью 0.9 падать слева от него, и с вероятностью 0.1 справа. Тут авторы предлагают буквально "...выльем чайник на плиту и сведем задачу к предыдущей": давайте сгенерируем дополнительную булевую случайную величину, с "обратным" распределением (0.1 к 0.9), и скомбинируем ее с вероятностью попадания справа-слева (через AND). Для этой комбинированной случайной величины значение X0.9 будет как-бы-медианой, а медиану мы уже искать умеем:
step = 1
Q    = 0.9 # quantile we're about to estimate
Xq   = 0   # initial estimation of Q-th quantile
for x in stream:
    r = random()                   # r in [0,1]
    if x > Xq and r < Q:        # Prob(x > Xq)=0.1, Prob(r < Q) = 0.9
        Xq += step
    elif val < Xq and r > Q:    # Prob(x < Xq)=0.9, Prob(r > Q) = 0.1
        Xq -= step

(Без генератора случайных чисел можно и обойтись: достаточно заметить, что матожидание величины шага вправо будет 0.1*step, а влево 0.9*step. Если заменить "шаг step с вероятностью 0.1" на "шаг 0.1*step всегда", то алгоритм станет проще и ГСЧ будет не нужен. Странно, что авторы статьи не отметили такой возможности: в тестах, что я делал, оба алгоритма ведут себя практически неразличимо, и только при увеличении можно различить чуть разное поведение на уровне отдельных шагов. Мое предположение, что просто анализ такого алгоритма сложнее, и авторы не стали заморачиваться)


Алгоритм выглядит подозрительно просто, у меня сразу возникли два вопроса:
1. Как быстро он сходится?
2. Насколько стабильно?

Про скорость сходимости можно сразу заметить, что алгоритм двигает текущую оценку не более чем на step за итерацию. Так что если между текущей оценкой и истинным значением квантиля дистанция D, то алгоритм сойдется точно не раньше, чем за N=D/step шагов. Более аккуратный анализ (см статью) дает такую оценку: за не более чем $$N\frac{|ln(1-P)|}{e}$$ шагов алгоритм с вероятностью P хотя бы раз попадет в e-окрестность квантиля. Оценка отличается от N только на константу $$\frac{|ln(1-P)|}{e}$$, но константа эта довольно не маленькая: если хочется с вероятностью 90% попасть в 5% окрестность квантиля, то ждать придется "всего лишь" не больше чем N*|ln(1-0.9)|/0.05 =~ 46*N шагов.

Почему так медленно: пусть истинное значение квантиля Q это точка X. Алгоритм сейчас находится в точке X' которая является квантилем Q' (т.е $$Pr(x<X') = Q'$$). Тогда матожидание смещения на текущем шаге будет просто step*(Q'-Q). Конкретный закон эволюции матожидания отсюда не получишь — нужно знать вид функции распределения P(X) — но и по общему виду понятно, что приближение (матожидания) к истинному значению квантиля будет асимптотическим: чем ближе подходим, тем медленнее идем. Если бы надо было "дошагать" матожиданием, то ждать пришлось бы вообще бесконечно, но ближе к концу существенную роль начинает играть накопившийся статистический разброс вокруг среднего: за счет него ждать бесконечно долго нам почти наверняка не придется.

Теперь про стабильность: поскольку алгоритм статистический, то даже сойдясь к значению квантиля он не остановится, а будет продолжать топтаться вокруг. Насколько далеко вокруг?

В первом приближении в окрестности заданного квантиля алгоритм делает шаги вправо и влево с одинаковой вероятностью. Это похоже на одномерное броуновское движение, про которое известно, что среднее расстояние от начала пропорционально корню из количества шагов. Что звучит довольно хреново, потому что означает, что со временем наша оценка уплывет сколь угодно далеко. Но это верно только вблизи целевого квантиля: чем дальше алгоритм от него отходит, тем больше будут отличаться вероятности шагнуть "от" и "к" — и вероятность шагнуть ближе к квантилю будет возрастать. В статье оценивают это более аккуратно, и получается такая оценка:

Пусть δ это максимальное значение (дискретной) функции распределения элементов из потока (т.е. максимальная вероятность выпадения какого-то конкретного значения). Тогда вероятность того, что за N шагов алгоритм отклонится от квантиля Q больше, чем на ΔQ оценивается так:
$$Pr(|Q(N)-Q| > \Delta Q) <= N e^-\frac{(2\Delta Q)^2}{\delta}$$


Например: если в потоке 1000 разных равновероятных значений, то δ=1/1000, и если мы хотим ограничиться 5% отклонением от квантиля, то можем быть спокойны: за пределы 5% алгоритм вылезет не чаще, за 10000 шагов. А вот если потребовать 4% коридор, то за него мы будем вылезать уже за 300-400 шагов.

Мне обе оценки показались не очень интуитивны (особенно вторая). Так что я посчитал и численно смоделировал эволюцию распределения вероятностей результатов алгоритма по шагам. Для стабильного распределения входных данных алгоритм сходится к стабильному же распределению вокруг истинной медианы — т.е. дисперсия вырастет до какой-то величины, и дальше расти не будет. На графике результат симуляции для самого простого, равномерного распределения входных данных U[0..128), и шага в 1. Красным показана настоящая медиана (=64), синим — наиболее вероятный результат алгоритма, а пунктир очерчивает интервал, в который алгоритм попадет с 90% вероятностью.


Видно, что наиболее вероятное значение даже после 400 шагов все еще отличается от настоящей медианы, а вот верхняя граница 90% интервала пересекла медиану уже на 120 шаге. На следующем графике — вероятность достижения медианы к шагу N:


Т.е. минимально мы могли бы дошагать до медианы за 64 шага, но с вероятностью 20% дошагаем только к 120-му шагу, с вероятностью 60% к 150, а к 200-му шагу уже почти наверняка дошагаем. Сильно лучше, чем показывает оценка из статьи (там нам обещали не более чем за 40*N), ну так равномерное распределение и не самый сложный случай. С другой стороны разброс вокруг медианы довольно приличный: с вероятностью 10% мы обнаружим алгоритм за пределами интервала (35%, 65%).

А если пересчитать все то же самое с разбросом входных данных в 8 раз больше (0..1024) и тем же шагом 1, то время 90% достижения медианы увеличится до 1700 шагов, зато 90% интервал вокруг медианы сожмется до (45%, 55%)

Виден компромисс между скоростью сходимости и стабильностью. В первом случае шаг составлял 1/64 от дистанции до медианы, и алгоритм почти сошелся за 200 шагов, и держит медиану с точностью +/- 15%. Во втором случае шаг был 1/512 от дистанции до медианы, и алгоритм почти сошелся за 1700 шагов, и держит медиану с точностью +/-5%.

Тут же возникает идея подстраивать шаг: начинать с большого, и постепенно уменьшать, как это делают во алгоритмах имитации отжига. Но это подходит только если мы знаем, что распределение стационарно, если же оно меняется со временем, то уменьшая шаг мы будем одновременно замедлять подстройку к изменениям в распределении. Хорошо бы менять шаг в зависимости от того, насколько мы далеко от цели — но как оценить, насколько мы далеко? Тут авторы достают свой второй козырь:

Frugal2U

Этот алгоритм расширяет предыдущий: основная логика та же, но к ней добавляется подстройка величины шага. Чтобы организовать подстройку, требуется еще одна ячейка памяти для хранения текущего шага, currentStep, и еще один бит sign для хранения направления предыдущего шага. И идея, снова, очень простая: давайте запоминать, в каком направлении мы шагнули (sign), и увеличивать текущий шаг каждый раз, как мы делаем следующий шаг в одном и том же направлении, что и предыдущий:
step = 1
currentStep=step
sign = 1
Q    = 0.9       # quantile we're about to estimate
Xq   = 0         # initial estimation of Q-th quantile
for x in stream:
  r = np.random.random()
  if x < Xq and r > Q:
      currentStep -= sign * step
      Xq -= max(currentStep, step)
      if value > Xq:  # too much of adjustment
          currentStep += Xq - value
          Xq = value
      if sign > 0 and currentStep > step:          
          currentStep = step              # if sign was changed -> reset currentStep to minimum
      sign = -1
  elif x > Xq and r < Q:
      currentStep += sign * step
      Xq += max(currentStep, step)
      if value < Xq:  # too much of adjustment
          currentStep += value - Xq
          Xq = value
      if sign < 0 and currentStep > step:          
          currentStep = step             # if sign was changed -> reset currentStep to minimum
      sign = 1
(код выглядит громоздко, потому что довольно много места занимает обработка выходов за границы разумности)

Почему это работает? На пальцах: если вероятность пойти в какую-то сторону один раз равна P, то вероятность пойти туда же N раз подряд — P^N, и с такой вероятностью шаг вырастет до N*step. Чем больше P, тем больше вероятность достичь бОльшего шага. Но в нашем алгоритме вероятность P зависит от расстояния до истинного значения квантиля: чем мы дальше от цели, тем больше вероятность шагнуть именно в ее сторону.

Более количественно: пусть мы сделали N шагов в одну сторону подряд — вероятность этого P^N, и шаг вырос до максимума N*step, т.е. средний шаг step*(N+1)/2. То есть с вероятностью P^N мы шли со средним шагом step*(N+1)/2. Просуммируем по всем N (спонсор вычисления бесконечных сумм wolfram alpha), и получим матожидание шага:
$$E[D(P)] = step\cdot\sum_{i=1}^{\infty}P^i \frac{(i+1)}{2} = step\frac{(2-P)*P}{2(1-P)^2}$$
Переводя на русский: схема "увеличиваем шаг пока идем в одну сторону" неявно задает вот такую зависимость среднего размера шага от дистанции до цели (целью у нас выступает медиана, P=0.5).

Я построил на одном графике сразу две зависимости: E[D(P)] и E[D(1-P)], потому что мы ведь из каждой точки можем пойти в обе стороны, и удобно сразу видеть, как меняется матожидание шага в одну, и в другую сторону


Оба графика пересекаются в P=0.5 — это как раз стационарная точка, к которой алгоритм сходится. А вот по мере отклонения от 0.5 графики расходятся, и очень быстро: уже на расстоянии в 0.3 от равновесия средний шаг в направлении равновесия будет >10 единиц. То есть если мы на расстоянии 0.3 от медианы, то приближаться мы будем примерно десятикратными шагами. Похоже, несмотря на простоту, такая схема управления размером шага может вполне хорошо работать!

(По мере приближения к 1 формула дает шаг равный бесконечности. Но сильно переживать по этому поводу не стоит: ведь получена в предположении, что P постоянно, т.е. мы стоим на месте. А чтобы шаг действительно рос, нам нужно двигаться — ведь шаг увеличивается только на 1 за итерацию. И как только мы начнем двигаться, то мы начнем и приближаться к истинному квантилю, а значит и P начнет отходить от 1, и приближаться к 0.5. Т.е. порвать вселенную нам не грозит, но набрать достаточно большой шаг можно — поэтому в алгоритме выше и есть ветки, проверяющие перехлест)

Под конец пару примеров, как алгоритмы ведут себя на тестовых данных.
Пример 1: оценка 80% квантиля на выборке из равномерного распределения. Красным показан "эталон" для сравнения — сглаженный "честный" квантиля на истории из последних 400 элементов:

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

Пример 2: ступеньки + равномерное распределение. Здесь можно лучше рассмотреть скорости перехода — видно, насколько быстро перестраивается Frugal2U, насколько сильно запаздывает Frugal1U.

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

Комментариев нет:

Отправить комментарий