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

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

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

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

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

Введение

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

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

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

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

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

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 просто не хватает разрядности.

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

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

Не так просто найти данные, на которых плохо работает решение с использованием 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. Обозначения

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

$display$

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

$S_n^{a,b,...,z} = \sum_{i=1}^n{a_i\cdot b_i\cdot ...\cdot z_i}$

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

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

$S_n^{wx} = \sum_{i=1}^n{w_i\cdot x_i}$

$S_n^{wxy} = \sum_{i=1}^n{w_i\cdot x_i\cdot y_i}$

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

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

$m_n^{wy} = \frac{\sum_{i=1}^n{w_i\cdot 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}$

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

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}\cdot\Big( 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. Экспериментальное сравнение методов

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

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

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

$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$

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

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

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

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

Заключение

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

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.

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

Литература

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