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

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


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


uynd2biz-sh2tfwo_fyocftj2cw.jpeg


Линейные модели часто используют благодаря нескольким очевидным преимуществам:


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


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


Содержание


  • 1. Многомерная линейная регрессия
  • 2. Общее описание решения
  • 3. Применение метода Уэлфорда
  • 4. Реализация
  • 5. Экспериментальное сравнение методов
  • Заключение
  • Литература


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


В задаче многомерной линейной регрессии требуется восстановить неизвестную зависимость наблюдаемой вещественной величины (которая является значением неизвестной целевой функции) $y^*$ от набора вещественных же признаков $f_1, f_2,..., f_m$.


Обычно для каждого отдельного наблюдения записываются значения признаков и целевой функции. Тогда из значений признаков для всех наблюдений можно составить матрицу признаков $X$, а из значений целевой функции — вектор её значений $y$:


$X_{ij} = f_j(x_i), 1 \leq i \leq n, 1 \leq j \leq m$


$y_i = y^*(x_i), 1 \leq i \leq n, 1 \leq j \leq m$


Здесь $x_i$ — некоторое конкретное наблюдение. Таким образом, строчки нашей матрицы соответствуют наблюдениям, а столбцы — признакам, а общее количество наблюдений и, следовательно, строк в матрице, — $n$.


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


Найти решение задачи означает построить некоторую решающую функцию $a : \mathbb{R^m} \rightarrow \mathbb{R}$. Сейчас мы говорим о задаче линейной регрессии, поэтому мы будем считать, что решающая функция является линейной. На самом деле, это значит, что решающая функция представляет из себя просто скалярное произведение вектора признаков на некоторый вектор весов:


$a_{w}(x_i) = \sum_{j=1}^{m} w_j \cdot f_j(x_i)$


Стало быть, решением нашей задачи является вектор весов признаков $w = \langle w_1, w_2, ..., w_m \rangle$.


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


zvb17oyksje3b8awn3qj47odl-o.png


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


$Q(a,X,y)= \sqrt{\frac{1}{n} \sum_{i=1}^{n} (a(x_i) - y_i)^2}$


Итак, задачей многомерной линейной регрессии является нахождение набора весов $w$ такого, что значение функционала потерь $Q(a_{w},A,y)$ достигает на нём своего минимума.


Для анализа обычно удобно отказаться от радикала и усреднения в функционале потерь:


$ Q_1(a,X,y)= \sum_{i=1}^{n} (a(x_i) - y_i)^2 $


Для дальнейшего изложения будет отождествить наблюдения с наборами соответствующих им факторов:


$ x_{ij} = f_j(x_i) $


Нетрудно видеть, что функционал потерь теперь записывается просто как скалярный квадрат: вектор предсказаний можно записать как $Xw$, а вектор отклонений от правильных ответов — как $Xw-y$. Тогда:


$ Q_1(a,X,y)= \langle Xw-y,Xw-y \rangle $


2. Общее описание решения


Из регрессионного анализа хорошо известно, что вектор-решение в задаче многомерной линейной регрессии находится как решение системы линейных алгебраических уравнений:


$(X^{T}X)w=X^{T}y$


Получить доказательство этому несложно. Действительно, давайте рассмотрим частную производную функционала $Q_1$ по весу $j$-го признака:


$ \frac{\partial Q_1}{\partial w_j} = 2 \cdot \sum_{i=1}^{n} (\langle w, x_i \rangle - y_i) \cdot x_{ij} $


Если приравнять эту производную нулю, получим:


$ \sum_{i=1}^{n} (\sum_{k=1}^{m} w_k \cdot x_{ik} - y_i) \cdot x_{ij} = 0 $


$ \sum_{i=1}^{n} \Big(\sum_{k=1}^{m} w_k \cdot x_{ik}\cdot x_{ij}\Big) = \sum_{i=1}^{n}x_{ij}y_i $


Теперь можно поменять порядок суммирования:


$ \sum_{k=1}^{m} w_k \sum_{i=1}^{n} \Big(x_{ik} \cdot x_{ij} \Big) = \sum_{i=1}^{n}x_{ij}y_i $


Теперь понятно, что коэффициент при $w_k$ слева — это соответствующий элемент матрицы $X^T X$:


$\sum_{i=1}^{n} (x_{ik} x_{ij} )= (X^{T} X)_{jk} $


В свою очередь, правая часть построенного уравнения — это элемент вектора $X^T y$:


$ \sum_{i=1}^{n}x_{ij}y_i = (X^T y)_j $


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


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


При решении СЛАУ $(X^{T}X)w=X^{T}y$ возникают различные проблемы, например, мультиколлинеарность признаков, которая приводит к вырожденности матрицы $(X^{T}X)$. Кроме того, необходимо выбрать метод для решения указанной СЛАУ.


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


Однако сейчас нас будет интересовать намного более приземлённая проблема: как сформировать уравнения для СЛАУ? Если эти уравнения сформированы корректно, а вычислительные погрешности невелики, можно надеяться, что выбранный метод приведёт к хорошему решению. Если же ошибки появились уже на этапе формирования матрицы системы и вектора в правой части, никакой метод решения СЛАУ не достигнет успеха.


Из сказанного выше ясно, что для составления СЛАУ необходимо вычислять коэффициенты двух типов:


$ (X^TX)_{kj} = \sum_{i=1}^{n} (x_{ik} x_{ij} )$


$ (X^Ty)_{j} = \sum_{i=1}^{n} (x_{ij}y_i)$


Если бы выборка была центрированной, это были бы выражения для ковариаций, о которых говорили в самой первой статье!


С другой стороны, выборку можно центрировать самостоятельно. Пусть $X'$ — матрица центрированной выборки, а $y'$ — вектор ответов для центрированной выборки:


$x'_{ij} = x_{ij} - \frac{\sum_{k=1}^{n}x_{kj}}{n} $


$y'_{j} = y_{j} - \frac{\sum_{k=1}^{n}y_k}{n} $


Обозначим через $\bar{x}_j$ среднее значение $j$-го признака, а через $\bar{y}$ — среднее значение целевой функции:


$\bar{x}_j = \frac{1}{n}\sum_{k=1}^{n}x_{kj}$


$\bar{y} = \frac{1}{n}\sum_{k=1}^{n}y_{k}$


Пусть решением центрированной задачи является вектор коэффициентов $w'$. Тогда для $i$-го события величина предсказания будет вычисляться следующим образом:


$ a'(x') = \sum_{j=1}^m{x'_j \cdot w'_j} = \sum_{j=1}^m{(x_j - \bar{x}_j)\cdot w'_j} = \sum_{j=1}^m{x_j \cdot w'_j} - \sum_{j=1}^m{\bar{x}_j\cdot w'_j} $


При этом $a'(x')$ приближает центрированную величину. Отсюда можем окончательно записать решение исходной задачи:


$ a(x) = \sum_{j=1}^m{x_j \cdot w'_j} - \sum_{j=1}^m{\bar{x}_j\cdot w'_j} + \bar{y}$


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


$\bar{y} - \sum_{j=1}^m{\bar{x}_j\cdot w'_j}$


4. Реализация


Реализацию методов многомерной линейной регрессии я поместил в ту же программу, о которой рассказывал в прошлый раз.


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


Решатель задачи многомерной линейной регрессии по методу Уэлфорда реализован в виде класса TWelfordLRSolver: 


class TWelfordLRSolver {
protected:
    double GoalsMean = 0.;
    double GoalsDeviation = 0.;

    std::vector FeatureMeans;
    std::vector FeatureWeightedDeviationFromLastMean;
    std::vector FeatureDeviationFromNewMean;
    std::vector LinearizedOLSMatrix;

    std::vector OLSVector;

    TKahanAccumulator SumWeights;

public:
    void Add(const std::vector& features, const double goal, const double weight = 1.);
    TLinearModel Solve() const;
    double SumSquaredErrors() const;

    static const std::string Name() {
        return "Welford LR";
    }

protected:
    bool PrepareMeans(const std::vector& features, const double weight);
};


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


При добавлении элемента первым делом вычисляются новый вектор средних значений признаков и величины отклонений.


erpatgqlj6o8uccfrzwmnypt-s8.png


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


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


npwau8b0rkishxdcojg2ujhmisi.png


Экономия на арифметических операциях здесь более чем оправдана, т.к. эта процедура обновления матрицы является наиболее затратной с т.з. ресурсов CPU с асимптотикой $O(m^2)$. Остальные операции, в т.ч. обновление вектора правых частей, асимптотически линейны:


zsghwb7kt0xw246t27masvbrdg4.png


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


TLinearModel TWelfordLRSolver::Solve() const {
    TLinearModel model;
    model.Coefficients = NLinearRegressionInner::Solve(LinearizedOLSMatrix, OLSVector);
    model.Intercept = GoalsMean;

    const size_t featuresCount = OLSVector.size();
    for (size_t featureNumber = 0; featureNumber < featuresCount; ++featureNumber) {
        model.Intercept -= FeatureMeans[featureNumber] * model.Coefficients[featureNumber];
    }

    return model;
}


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


Так же, как и в прошлый раз, реализации на основе метода Уэлфорда (welford_lr и normalized_welford_lr) сравниваются с «наивным» алгоритмом, в котором матрица и правая часть системы вычисляются по стандартным формулам (fast_lr).


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


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


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


linear_regression research-lr --features kin8nm.features
injure factor: 1
injure offset: 1
   fast_lr                                        time: 0.002304    R^2: 0.41231
   welford_lr                                     time: 0.003248    R^2: 0.41231
   normalized_welford_lr                          time: 0.003958    R^2: 0.41231

injure factor: 0.1
injure offset: 10
   fast_lr                                        time: 0.00217    R^2: 0.41231
   welford_lr                                     time: 0.0032    R^2: 0.41231
   normalized_welford_lr                          time: 0.00398    R^2: 0.41231

injure factor: 0.01
injure offset: 100
   fast_lr                                        time: 0.002164    R^2: -0.39518
   welford_lr                                     time: 0.003296    R^2: 0.41231
   normalized_welford_lr                          time: 0.003846    R^2: 0.41231

injure factor: 0.001
injure offset: 1000
   fast_lr                                        time: 0.002166    R^2: -1.8496
   welford_lr                                     time: 0.003114    R^2: 0.41231
   normalized_welford_lr                          time: 0.003978    R^2: 0.058884

injure factor: 0.0001
injure offset: 10000
   fast_lr                                        time: 0.002165    R^2: -521.47
   welford_lr                                     time: 0.003177    R^2: 0.41231
   normalized_welford_lr                          time: 0.003931    R^2: 0.00013929

full learning time:
   fast_lr                                        0.010969s
   welford_lr                                     0.016035s
   normalized_welford_lr                          0.019693s


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


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


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


Заключение


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


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


Литература


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

© Habrahabr.ru