- PVSM.RU - https://www.pvsm.ru -

Точное вычисление средних и ковариаций методом Уэлфорда

Метод Уэлфорда — простой и эффективный способ для вычисления средних, дисперсий, ковариаций и других статистик. Этот метод обладает целым рядом прекрасных свойств:

  • достигает отличных показателей по точности решений;
  • его чрезвычайно просто запомнить и реализовать;
  • это однопроходный онлайн-алгоритм, что крайне полезно в некоторых ситуациях.

Оригинальная статья Уэлфорда была опубликована в 1962 году. Тем не менее, нельзя сказать, что алгоритм сколь-нибудь широко известен в настоящее время. А уж найти математическое доказательство его корректности или экспериментальные сравнения с другими методами и вовсе нетривиально.

Настоящая статья пытается заполнить эти пробелы.

Точное вычисление средних и ковариаций методом Уэлфорда - 1 [1]

Содержание

Введение

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

Сама эта тема оказалась настолько для меня интересной, что я потом даже диссертацию [9] про это написал. Для построения модели из нескольких сотен деревьев приходится решать задачу линейной регрессии десятки тысяч раз, и оказалось сложным достичь хорошего качества во всех случаях, ведь данные для этих моделей весьма разнообразны, и проблемы мультиколлинеарности, регуляризации и вычислительной точности встают в полный рост. А ведь достаточно одной плохой модели в одном листе одного дерева для того, чтобы вся композиция оказалась совершенно непригодной.

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

У статьи следующая структура. В пункте 1 [3] мы рассмотрим простейшую задачу вычисления среднего, на примере которой поймём, что эта задача не так уж и проста, как это кажется на первый взгляд. Пункт 2 [4] вводит используемые в данной статье обозначения, которые пригодятся в разделах 3 [5] и 4 [6], посвящённых выводу формул метода Уэлфорда для, соответственно, вычисления взвешенных средних и взвешенных ковариаций. Если вам неинтересны технические подробности вывода формул, вы можете пропустить эти разделы. Пункт 5 [7] содержит результаты экспериментального сравнения методов, а в заключении [8] находится пример реализации алгоритмов на языке С++.

Модельный код для сравнения методов вычисления ковариации находится в моём проекте на github [10]. Более сложный код, в котором метод Уэлфорда применяется для решения задачи линейной регрессии, находится в другом проекте [11], речь о котором пойдёт в следующих статьях.

1. Погрешности при вычислении средних

Начну c классического примера. Пусть есть простой класс, вычисляющий среднее для набора чисел:

class TDummyMeanCalculator {
private:
    float SumValues = 0.;
    size_t CountValues = 0;
public:
    void Add(const float value) {
        ++CountValues;
        SumValues += value;
    }
    double Mean() const {
        return CountValues ? SumValues / CountValues : 0.;
    }
};

Попробуем применить его на практике следующим способом:

int main() {
    size_t n;
    while (std::cin >> n) {
        TDummyMeanCalculator meanCalculator;
        for (size_t i = 0; i < n; ++i) {
            meanCalculator.Add(1e-3);
        }
        std::cout << meanCalculator.Mean() << std::endl;
    }
    return 0;
}

Что же выведет программа?

Ввод Вывод
10000 0.001000040
1000000 0.000991142
100000000 0.000327680
200000000 0.000163840
300000000 0.000109227

Начиная с некоторого момента, сумма перестаёт меняться после прибавления очередного слагаемого: это происходит, когда SumValues оказывается равным 32768, т.к. для представления результата суммирования типу float просто не хватает разрядности.

Из этой ситуации есть несколько выходов:

  1. Перейти с float на double.
  2. Использовать один из более сложных методов суммирования [12].
  3. Использовать метод Кэхэна [13] для суммирования.
  4. Наконец, можно использовать метод Уэлфорда [14] для непосредственного вычисления среднего.

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

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

class TWelfordMeanCalculator {
private:
    float MeanValues = 0.;
    size_t CountValues = 0;
public:
    void Add(const float value) {
        ++CountValues;
        MeanValues += (value - MeanValues) / CountValues;
    }

    double Mean() const {
        return MeanValues;
    }
};

К вопросу о корректности этого кода мы вернёмся в дальнейшем. Пока же заметим, что реализация выглядит достаточно простой, а на тех данных, что используются в нашей программе, еще и работает с идеальной точностью. Это происходит из-за того, что разность value - MeanValues на первой итерации в точности равна среднему, а на последующих итерациях равна нулю.

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

Ввод Вывод
10000 0.001
1000000 0.001
100000000 0.001
200000000 0.001
300000000 0.001

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

2. Обозначения

Вывод формулы часто можно сделать очень простым и понятным, если выбрать правильные обозначения. Попробуем и мы. Будем считать, что заданы две последовательности вещественных чисел: $x$ и $y$, и последовательность соответствующих им весов $w$:

$x=x_1,x_2,...,x_n,...$

$y=y_1,y_2,...,y_n,...$

$w=w_1,w_2,...,w_n,...$

Поскольку мы хотим вычислять средние и ковариации, нам понадобятся обозначения для взвешенных сумм и сумм произведений. Попытаемся описать их единообразно:

$S_n^{a,b,...,z}=sum_{i=1}^n{a_icdot b_icdot ...cdot z_i}$

Тогда, например, $S_n^{w}$ — это сумма весов первых $n$ элементов, $S_n^{wx}$ — это взвешенная сумма первых $n$ чисел первой последовательности, а $S_n^{wxy}$ — сумма взвешенных произведений:

$S_n^{w}=sum_{i=1}^n{w_i}$

$S_n^{wx}=sum_{i=1}^n{w_icdot x_i}$

$S_n^{wxy}=sum_{i=1}^n{w_icdot x_icdot y_i}$

Также понадобятся средние взвешенные величины:

$m_n^{wx}=frac{sum_{i=1}^n{w_icdot x_i}}{sum_{i=1}^n{w_i}}=frac{S_n^{wx}}{S_n^w}$

$m_n^{wy}=frac{sum_{i=1}^n{w_icdot y_i}}{sum_{i=1}^n{w_i}}=frac{S_n^{wy}}{S_n^w}$

Наконец, введём обозначения для ненормированных "разбросов" $D_{n}^{wxy}$ и нормированных ковариаций $C_{n}^{wxy}$:

$D_{n}^{wxy}=sum_{i=1}^{n} w_i (x_i - m_{n}^{wx})(y_i - m_{n}^{wy})$

$C_{n}^{wxy}=frac{D_{n}^{wxy}}{S_{n}^w}=frac{sum_{i=1}^{n} w_i (x_i - m_{n}^{wx})(y_i - m_{n}^{wy})}{sum_{i=1}^{n}w_{i}}$

3. Вычисление средних

Докажем взвешенный аналог формулы, использованной нами выше для вычисления среднего по методу Уэлфорда. Рассмотрим разность $m_{n+1}^{wx} - m_{n}^{wx}$:

$m_{n+1}^{wx} - m_{n}^{wx}=frac{S_{n+1}^{wx}}{S_{n+1}^{w}} - frac{S_{n}^{wx}}{S_{n}^{w}}=frac{S_{n}^{w}S_{n+1}^{wx} - S_{n}^{wx}S_{n+1}^{w}}{S_{n+1}^{w}S_{n}^{w}}=$

$=frac{S_{n}^{w}S_{n}^{wx} + w_{n+1}x_{n+1}S_n^w - S_{n}^{wx}S_{n}^{w} - w_{n+1}S_n^{wx}}{S_{n+1}^{w}S_{n}^{w}}=frac{w_{n+1}x_{n+1}S_n^w - w_{n+1}S_n^{wx}}{S_{n+1}^{w}S_{n}^{w}}$

$=frac{w_{n+1}}{S_{n+1}^{w}}Big(x_{n+1} - frac{S_n^{wx}}{S_n^w}Big)=frac{w_{n+1}}{S_{n+1}^{w}}(x_{n+1} - m_n^{wx})$

В частности, если все веса равняются единице, получим, что

$m_{n+1}^{wx}=m_{n}^{wx} + frac{x_{n+1} - m_n^{wx}}{n+1}$

Кстати, формулы для "взвешенного" случая позволяют легко реализовывать операции, отличные от добавления ровно одного очередного числа. Например, удаление числа $x$ из множества, по которому вычисляется среднее — это то же, что добавление числа $-x$ с весом $-1$. Добавление сразу нескольких чисел — то же, что добавление одного из среднего с весом, равным количеству этих чисел, и так далее.

4. Вычисление ковариаций

Покажем, что и для взвешенной задачи верна классическая формула для вычисления ковариации. Сейчас для простоты будем работать с ненормированными величинами.

$D_{n}^{wxy}=sum_{i=1}^{n} w_i (x_i - m_n^{wx})(y_i - m_n^{wy})=S_n^{wxy} - m_n^{wx}S_n^{wy} - S_n^{wx}m_n^{wy}+S_n^wm_n^{wx}m_n^{wy}$

Отсюда уже легко видеть, что

$D_{n}^{wxy}=S_n^{wxy} - frac{S_n^{wx}S_n^{wy}}{S_n^w}$

Эта формула чрезвычайно удобна, в том числе и для онлайн-алгоритма, однако, если величины $S_n^{wxy}$ и ${S_n^{wx}S_n^{wy}}/{S_n^w}$ окажутся близкими и при этом большими по абсолютному значению, её использование приведёт к существенным вычислительным погрешностям.

Давайте попробуем вывести рекуррентную формулу для $D_n^{wxy}$, в каком-то смысле аналогичную формуле Уэлфорда для средних. Итак:

$ D_{n+1}^{wxy}=S_n^{wxy} + w_{n+1}x_{n+1}y_{n+1} - w_{n+1}x_{n+1}frac{S_{n+1}^{wy}}{S_{n+1}^{w}} - frac{S_n^{wx}S_{n+1}^{wy}}{S_{n+1}^w}=$

$=S_n^{wxy} + w_{n+1}x_{n+1}(y_{n+1} - m_{n+1}^{wy}) - frac{S_n^{wx}S_{n+1}^{wy}}{S_{n+1}^w}$

Рассмотрим последнее слагаемое:

$frac{S_n^{wx}S_{n+1}^{wy}}{S_{n+1}^w}=Big(frac{1}{S_n^w} - frac{w_{n+1}}{S_n^wS_{n+1}^w}Big)S_n^{wx}S_{n+1}^{wy}=frac{S_n^{wx}S_{n+1}^{wy}}{S_n^w}-w_{n+1}m_n^{wx}m_{n+1}^{wy}=$

$=frac{S_n^{wx}S_{n}^{wy}}{S_n^w}+w_{n+1}y_{n+1}frac{S_n^{wx}}{S_n^w}-w_{n+1}m_n^{wx}m_{n+1}^{wy}=frac{S_n^{wx}S_{n}^{wy}}{S_n^w}+w_{n+1}m_n^{wx}cdot(y_{n+1}-m_{n+1}^{wy})$

Подставим получившееся в выражение для $D_{n+1}^{wxy}$:

$ D_{n+1}^{wxy}=S_n^{wxy} + w_{n+1}x_{n+1}(y_{n+1} - m_{n+1}^{wy}) - frac{S_n^{wx}S_{n}^{wy}}{S_n^w}-w_{n+1}m_n^{wx}cdot(y_{n+1}-m_{n+1}^{wy})=$

$=Big[S_n^{wxy} - frac{S_n^{wx}S_{n}^{wy}}{S_n^w}Big] + w_{n+1}(x_{n+1}-m_n^{wx})(y_{n+1} - m_{n+1}^{wy})=$

$=D_{n}^{wxy} + w_{n+1}(x_{n+1}-m_n^{wx})(y_{n+1} - m_{n+1}^{wy})$

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

double WelfordCovariation(const std::vector<double>& x, const std::vector<double>& y) {
    double sumProducts = 0.;
    double xMean = 0.;
    double yMean = 0.;

    for (size_t i = 0; i < x.size(); ++i) {
        xMean += (x[i] - xMean) / (i + 1);
        sumProducts += (x[i] - xMean) * (y[i] - yMean);
        yMean += (y[i] - yMean) / (i + 1);
    }

    return sumProducts / x.size();
}

Также интересным является вопрос об обновлении величины собственно ковариации, т.е. нормированной величины:

$C_{n+1}^{wxy}=frac{sum_{i=1}^{n+1} w_i (x_i - m_{n+1}^{wx})(y_i - m_{n+1}^{wy})}{sum_{i=1}^{n+1}w_{i}}=frac{D_{n+1}^{wxy}}{S_{n+1}^w}=$

$=frac{1}{S_{n+1}^w}cdotBig( D_n^{wxy} + w_{n+1}(x_{n+1}-m_n^{wx})(y_{n+1} - m_{n+1}^{wy})Big)=$

$=frac{D_n^{wxy}}{S_{n+1}^w}+frac{w_{n+1}}{S_{n+1}^w}(x_{n+1}-m_n^{wx})(y_{n+1} - m_{n+1}^{wy})$

Рассмотрим первое слагаемое:

$frac{D_n^{wxy}}{S_{n+1}^w}=Big(frac{1}{S_n^w}-frac{w_{n+1}}{S_{n+1}^wS_n^w}Big)D_n^{wxy}=frac{D_n^{wxy}}{S_n^w}Big(1 - frac{w_{n+1}}{S_{n+1}^w}Big)=C_n^{wxy}Big(1-frac{w_{n+1}}{S_{n+1}^w}Big)$

Вернёмся теперь к рассмотрению $C_{n+1}^{wxy}$:

$C_{n+1}^{wxy}=C_n^{wxy}Big(1-frac{w_{n+1}}{S_{n+1}^w}Big) + frac{w_{n+1}}{S_{n+1}^w}(x_{n+1}-m_n^{wx})(y_{n+1} - m_{n+1}^{wy})$

Это можно переписать, например, так:

$C_{n+1}^{wxy}=C_n^{wxy} + frac{w_{n+1}}{S_{n+1}^w}Big((x_{n+1}-m_n^{wx})(y_{n+1} - m_{n+1}^{wy}) - C_n^{wxy}Big)$

Получилась формула, действительно обладающая удивительным сходством с формулой для обновления среднего!

5. Экспериментальное сравнение методов

Я написал небольшую программу [10], в которой реализуются три способа вычисления ковариации:

  1. "Стандартный" алгоритм, напрямую вычисляющий величины $S_n^{wxy}$, $S_n^{wx}$ и $S_n^{wy}$.
  2. Алгоритм, использующий метод Кэхэна.
  3. Алгоритм, использующий метод Уэлфорда.

Данные в задаче формируются следующим образом: выбираются два числа, $m_x$ и $m_y$ — средние двух выборок. Затем выбираются еще два числа, $d_x$ и $d_y$ — соответственно, отклонения. На вход алгоритмам подаются последовательности чисел вида

$x_i=m_x pm d_x,$

$y_i=m_y pm d_y,$

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

$C_{n}^{xy}=frac{sum_{i=1}^{n} (x_i - m_x)(y_i - m_y)}{n}=d_x cdot d_y $

Истинная ковариация константна и не зависит от числа слагаемых, поэтому мы можем вычислять относительную погрешность вычисления для каждого метода на каждой итерации. Так, в текущей реализации $d_x=d_y=1$, а средние принимают значения 100000 и 1000000.

На первом графике изображена зависимость величины относительной погрешности при использовании наивного метода для среднего, равного 100 тысячам. Этот график демонстрирует несостоятельность наивного метода: начиная с какого-то момента погрешность начинает быстро расти, достигая совершенно неприемлемых величин. На тех же данных методы Кэхэна и Уэлфорда не допускают существенных ошибок.

Точное вычисление средних и ковариаций методом Уэлфорда - 62

Второй график построен для метода Кэхэна при среднем, равном одному миллиону. Ошибка не растёт с увеличением количества слагаемых, но, хотя она и существенно ниже, чем ошибка "наивного" метода, всё равно она слишком велика для практических применений.

Точное вычисление средних и ковариаций методом Уэлфорда - 63

Метод Уэлфорда, в свою очередь, и на этих данных демонстрирует идеальную точность!

Заключение

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

class TWelfordMeanCalculator {
private:
    double Mean = 0.;
    size_t Count = 0;
public:
    void Add(const float value) {
        ++Count;
        Mean += (value - Mean) / Count;
    }
    double GetMean() const {
        return Mean;
    }
};

class TWelfordCovariationCalculator {
private:
    size_t Count = 0;
    double MeanX = 0.;
    double MeanY = 0.;
    double SumProducts = 0.;
public:
    void Add(const double x, const double y) {
        ++Count;
        MeanX += (x - MeanX) / Count;
        SumProducts += (x - MeanX) * (y - MeanY);
        MeanY += (y - MeanY) / Count;
    }
    double Covariation() const {
        return SumProducts / Count;
    }
};

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

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

Литература

  1. github.com: Different covariation calculation methods [10]
  2. machinelearning.ru: Сложение большого множества чисел, существенно отличающихся по величине [12]
  3. ru.wikipedia.org: Алгоритм Кэхэна [13]
  4. en.wikipedia.org: Algorithms for calculating variance: Online algorithm [14]

Автор: ashagraev

Источник [1]


Сайт-источник PVSM.RU: https://www.pvsm.ru

Путь до страницы источника: https://www.pvsm.ru/programmirovanie/261056

Ссылки в тексте:

[1] Image: https://habrahabr.ru/post/333426/

[2] Введение: #vvedenie

[3] 1. Погрешности при вычислении средних: ##1-pogreshnosti-pri-vychislenii-srednih

[4] 2. Обозначения: #2-oboznacheniya

[5] 3. Вычисление средних: #3-vychislenie-srednih

[6] 4. Вычисление ковариаций: #4-vychislenie-kovariaciy

[7] 5. Экспериментальное сравнение методов: #5-eksperimentalnoe-sravnenie-metodov

[8] Заключение: #zaklyuchenie

[9] диссертацию: https://yadi.sk/d/kKLS-3eNbqNXR

[10] в моём проекте на github: https://github.com/ashagraev/covariation

[11] другом проекте: https://github.com/ashagraev/linear_regression

[12] методов суммирования: http://www.machinelearning.ru/wiki/index.php?title=%D0%A1%D0%BB%D0%BE%D0%B6%D0%B5%D0%BD%D0%B8%D0%B5_%D0%B1%D0%BE%D0%BB%D1%8C%D1%88%D0%BE%D0%B3%D0%BE_%D0%BC%D0%BD%D0%BE%D0%B6%D0%B5%D1%81%D1%82%D0%B2%D0%B0_%D1%87%D0%B8%D1%81%D0%B5%D0%BB%2C_%D1%81%D1%83%D1%89%D0%B5%D1%81%D1%82%D0%B2%D0%B5%D0%BD%D0%BD%D0%BE_%D0%BE%D1%82%D0%BB%D0%B8%D1%87%D0%B0%D1%8E%D1%89%D0%B8%D1%85%D1%81%D1%8F_%D0%BF%D0%BE_%D0%B2%D0%B5%D0%BB%D0%B8%D1%87%D0%B8%D0%BD%D0%B5

[13] метод Кэхэна: https://ru.wikipedia.org/wiki/%D0%90%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_%D0%9A%D1%8D%D1%85%D1%8D%D0%BD%D0%B0

[14] метод Уэлфорда: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm