ashagraev Aug 16 2017 at 11:23

Метод Уэлфорда и одномерная линейная регрессия

7 min

11K

Programming*Algorithms*Mathematics*Machine learning*

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

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

1. Одномерная линейная регрессия

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

Нашей задачей будет восстановить линейную зависимость между признаками и ответами, то есть, построить линейную решающую функцию $f : \mathbb{R} \rightarrow \mathbb{R}$ :

$f(x_i) = \alpha \cdot x_i + \beta$

При этом минимизируется среднеквадратичный функционал потерь:

$Q(f,x,y,w)= \sqrt{\frac{1}{\sum_{i=1}^{n} w_i} \sum_{i=1}^{n} w_i \cdot (f(x_i) - y_i)^2}$

Для анализа проще работать с формулой, свободной от радикала и нормировки:

$Q_1(f,x,y,w)=\sum_{i=1}^{n}w_i\cdot(f(x_i)-y_i)^2=\sum_{i=1}^{n} w_i \cdot(\alpha x_i+\beta-y_i)^2$

Поскольку точки минимума для функционалов $inline$ и $inline$ совпадают, такая замена корректна.

2. Решение для центрированных выборок

Для функционала потерь $inline$ легко записать производные по $\alpha$ и $\beta$ :

$\frac{\partial Q_1(f,x,y,w)}{\partial \alpha} = 2 \cdot \sum_{i=1}^{n} w_i x_i (\alpha x_i + \beta - y_i)$

$\frac{\partial Q_1(f,x,y,w)}{\partial \beta} = 2 \cdot \sum_{i=1}^{n} w_i (\alpha x_i + \beta - y_i)$

Если приравнять их к нулю, получим:

$\alpha \cdot \sum_{i=1}^{n} w_i x^2_i + \beta \cdot \sum_{i=1}^{n} w_i x_i - \sum_{i=1}^{n} w_i x_i y_i = 0$

$\alpha = \frac{\sum_{i=1}^{n} w_i x_i y_i - \beta \cdot \sum_{i=1}^{n} w_i x_i}{\sum_{i=1}^{n} w_i x^2_i}$

$\beta = \sum_{i=1}^{n} w_i y_i - \alpha \cdot \sum_{i=1}^{n} w_i x_i$

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

Если бы оптимальный параметр $\beta$ равнялся нулю, найти решение не составило бы труда. Можно заметить, что центрирование, стандартный для машинного обучения способ предобработки выборки, приводит именно к этому эффекту. Действительно, рассмотрим задачу для центрированных переменных:

$x'_i = x_i - \frac{\sum_{i=1}^n w_i x_i}{\sum_{i=1}^n w_i}$

$y'_i = y_i - \frac{\sum_{i=1}^n w_i y_i}{\sum_{i=1}^n w_i}$

Сумма взвешенных признаков теперь равняется нулю, так же, как и сумма взвешенных ответов:

$\sum_{k=1}^{n} w_k x_k' = \sum_{k=1}^{n} w_k \cdot \Big(x_k - \frac{\sum_{i=1}^n w_i x_i}{\sum_{i=1}^n w_i}\Big) = \sum_{k=1}^{n} w_k x_k - \sum_{i=1}^n w_i x_i = 0$

$\sum_{k=1}^{n} w_k y_k' = \sum_{k=1}^{n} w_k \cdot \Big(y_k - \frac{\sum_{i=1}^n w_i y_i}{\sum_{i=1}^n w_i}\Big) = \sum_{k=1}^{n} w_k y_k - \sum_{i=1}^n w_i y_i = 0$

Тогда оптимальное значение свободного параметра равняется нулю:

$\beta' = \sum_{i=1}^{n} w_i y_i' - \alpha' \cdot \sum_{i=1}^{n} w_i x'_i = 0$

А это значит, что и оптимальное значение параметра $\alpha'$ найти легко:

$\alpha' = \frac{\sum_{i=1}^{n} w_i x'_i y'_i}{\sum_{i=1}^{n} w_i x'^2_i}$

3. Решение в общем случае

Теперь давайте попробуем вернуться к общему случаю нецентрированных данных. Если $inline$ — решающая функция для центрированного случая, значения которой определяются формулой

$f'(x'_k) = \alpha' \cdot x'_k$

и приближают величину $y'_k=y_k - \frac{\sum_{i=1}^n y_i x_i}{\sum_{i=1}^n w_i}$ , то следующая решающая функция аппроксимирует уже величину $inline$ :

$f(x_k) = \alpha' \cdot \Big(x_k - \frac{\sum_{i=1}^n w_i x_i}{\sum_{i=1}^n w_i}\Big) + \frac{\sum_{i=1}^n y_i x_i}{\sum_{i=1}^n w_i} = \alpha' \cdot x_k + \Big(\frac{\sum_{i=1}^n y_i x_i}{\sum_{i=1}^n w_i} - \alpha' \cdot \frac{\sum_{i=1}^n w_i x_i}{\sum_{i=1}^n w_i}\Big)$

Поэтому решение изначальной задачи одномерной линейной регрессии можно записать следующим образом:

$\alpha = \frac{\sum_{i=1}^{n} w_i (x_i - m_n^{wx})(y_i - m_n^{wy})}{\sum_{i=1}^{n} w_i (x_i - m_n^{wx})(x_i - m_n^{wx})}$

$\beta = m_n^{wy} - \alpha \cdot m_n^{wx}$

Здесь мы используем введенное в прошлой статье обозначение для взвешенного среднего:

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

$m_n^{wy} = \frac{\sum_{i=1}^n w_i y_i}{\sum_{i=1}^n w_i}$

Можно понять, что такой переход корректен, еще одним способом. Если решение для центрированных данных оптимально, то параметры $\alpha'$ и $\beta'$ доставляют минимум функционалу потерь $inline$ :

$Q_1(f',x',y',w)=\sum_{i=1}^{n}w_i\cdot(\alpha' \cdot x'_i + \beta' - y'_i)^2$

Произведём теперь замену переменных, вернувшись к нецентрированным данным:

$Q_1(f',x',y',w)=\sum_{i=1}^{n}w_i\cdot\Big(\alpha'\cdot(x_i - m_n^{wx}) - y_i + m_n^{wy}\Big)^2=$

$=\sum_{i=1}^{n}w_i\cdot\Big(\alpha'\cdot x_i + (m_n^{wy} - \alpha' \cdot m_n^{wx}) - y_i \Big)^2$

Получившееся выражение описывает значение функционала потерь $inline$ для несмещённых данных в соответствии с формулами для $\alpha$ и $\beta$ , которые мы получили выше. Значение функционала при этом достигает минимума, следовательно, задача решена корректно!

4. Применение метода Уэлфорда

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

$\alpha = \frac{D_{n}^{wxy}}{D_{n}^{wxx}} = \frac{C_{n}^{wxy}}{C_{n}^{wxx}}$

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

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

void TWelfordSLRSolver::Add(const double feature, const double goal, const double weight) {
    SumWeights += weight;
    if (!SumWeights) {
        return;
    }
    const double weightedFeatureDiff = weight * (feature - FeaturesMean);
    const double weightedGoalDiff = weight * (goal - GoalsMean);
    FeaturesMean += weightedFeatureDiff / SumWeights;
    FeaturesDeviation += weightedFeatureDiff * (feature - FeaturesMean);
    GoalsMean += weightedGoalDiff / SumWeights;
    GoalsDeviation += weightedGoalDiff * (goal - GoalsMean);
    Covariation += weightedFeatureDiff * (goal - GoalsMean);
}

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

template <typename TFloatType>
void TWelfordSLRSolver::Solve(TFloatType& factor, TFloatType& intercept, const double regularizationParameter = 0.1) const {
    if (!FeaturesDeviation) {
        factor = 0.;
        intercept = GoalsMean;
        return;
    }
    factor = Covariation / (FeaturesDeviation + regularizationParameter);
    intercept = GoalsMean - factor * FeaturesMean;
}

Величина GoalsDeviation здесь не используется, но она потребуется нам в дальнейших статьях.

Объединение всех вычислений в рамках одного класса позволяет избежать некоторых накладных расходов. Скажем, если бы в реализации использовались два объекта для хранения средних и три объекта для хранения ковариаций (признаки с признаками, ответы с ответами, признаки с ответами), то сумма весов обновлялась бы пять раз для каждого примера из выборки.

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

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

Сравнивать будем, как обычно, "наивные" методы, методы, основанные на суммировании методом Кэхэна, и методы, основанные на методе Уэлфорда.

"Наивные" методы непосредственно применяют формулы для вычисления ковариаций:

void Add(const double feature, const double goal, const double weight = 1.) {
    SumFeatures += feature * weight;
    SumSquaredFeatures += feature * feature * weight;

    SumGoals += goal * weight;
    SumSquaredGoals += goal * goal * weight;

    SumProducts += goal * feature * weight;

    SumWeights += weight;
}

Класс является шаблонным и имеет специализации со счетчиками типа double и типа TKahanAccumulator.

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

Для проверки разработанных методов воспользуемся модельными данными из коллекции LIAC. Некоторые из наборов данных для удобства помещены в директорию data в формате, который понимает написанная программа.

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

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

Например, для выборки kin8nm результаты работы получаются следующими:

injure factor: 1
injure offset: 1
   fast_bslr                                      time: 0.001322    R^2: 0.27359
   kahan_bslr                                     time: 0.002999    R^2: 0.27359
   welford_bslr                                   time: 0.00432     R^2: 0.27359
   normalized_welford_bslr                        time: 0.004288    R^2: 0.27359

injure factor: 0.1
injure offset: 10
   fast_bslr                                      time: 0.001256    R^2: 0.27359
   kahan_bslr                                     time: 0.002948    R^2: 0.27359
   welford_bslr                                   time: 0.004303    R^2: 0.27359
   normalized_welford_bslr                        time: 0.004275    R^2: 0.27359

injure factor: 0.01
injure offset: 100
   fast_bslr                                      time: 0.001283    R^2: 0.27359
   kahan_bslr                                     time: 0.003015    R^2: 0.27359
   welford_bslr                                   time: 0.004304    R^2: 0.27359
   normalized_welford_bslr                        time: 0.004285    R^2: 0.27359

injure factor: 0.001
injure offset: 1000
   fast_bslr                                      time: 0.001262    R^2: 0.27324
   kahan_bslr                                     time: 0.002977    R^2: 0.27359
   welford_bslr                                   time: 0.004329    R^2: 0.27359
   normalized_welford_bslr                        time: 0.00428     R^2: 0.27359

injure factor: 0.0001
injure offset: 10000
   fast_bslr                                      time: 0.00128     R^2: -59.271
   kahan_bslr                                     time: 0.003009    R^2: -0.0005269
   welford_bslr                                   time: 0.004304    R^2: 0.27359
   normalized_welford_bslr                        time: 0.00428     R^2: 0.27359

full learning time:
   fast_bslr                                      0.006403s
   kahan_bslr                                     0.014948s
   welford_bslr                                   0.02156s
   normalized_welford_bslr                        0.021408s

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

Заключение

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

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

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

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

Литература

habrahabr.ru: Точное вычисление средних и ковариаций методом Уэлфорда
github.com: Linear regression problem solvers with different computation methods
github.com: The collection of ARFF datasets of the Connectionist Artificial Intelligence Laboratory (LIAC)
machinelearning.ru: Одномерная линейная регрессия

Tags:

Hubs: