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

машинное обучение

Формулы лучше отображаются на моем сайте,Мистер Лу.info/machine- сейчас..., Добро пожаловать в гости.

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

Normal Equation

Для функции потерь вы можете сделать ее производную равной нулю и найти экстремум функции потерь.

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

Предполагая, что наша модель имеет только одномерные данные, модель представляет собой прямую линиюf(x) = ax + b, мы делимсяmобучающие данные, а функция потерь представляет собой среднее значение суммы квадратов ошибок:

L(a, b) =\frac{1}{m} \sum_{i=1}^m[(ax_i + b) - y_i]^2

даaиbПроизводные получаются отдельно, когда производная равна 0, функция потерь наименьшая.

\begin{aligned} \frac{\partial}{\partial a}L(a, b) &= \frac{2}{m}\sum_{i=1}^mx_i(ax_i + b -y_i) \\ &=\frac{2a}{m}\sum_{i=1}^m{x_i}^2 + \frac{2b}{m}\sum_{i=1}^mx_i-\frac{2}{m}\sum_{i=1}^m{x_iy_i}  \\ &= \frac{2}{m}(a\sum_{i=1}^m{x_i}^2 + b\sum_{i=1}^mx_i - \sum_{i=1}^m{x_iy_i}) \\ &= \frac{2}{m}(a\sum_{i}^m{x_i}^2-\sum_{i=1}^m(y_i-b)x_i) \end{aligned}
\begin{aligned} \frac{\partial}{\partial b}L(a, b) &= \frac{2}{m}\sum_{i=1}^max_i+b-y_i \\ &= \frac{2a}{m}\sum_{i=1}^m{x_i} + 2b - \frac{2}{m}\sum_{i=1}^m{y_i} \\ &= \frac{2}{m}(\sum_{i=1}^max_i + mb - \sum_{i=1}^my_i) \\ &= \frac{2}{m}(mb - \sum_{i=1}^m(y_i -ax_i)) \end{aligned}

Две приведенные выше формулы являются функциями потерьLправильноaиbПроведите смещение. Когда производная равна 0, можно получить минимальное значение функции потерь, то есть оптимальное решение можно получить из двух приведенных выше формул.a^*иb^*.

\frac{\partial}{\partial a}L(a, b) = 0 \Rightarrow \frac{2}{m}(a\sum_{i=1}^m{x_i}^2-\sum_{i=1}^m(y_i-b)x_i) = 0
\frac{\partial}{\partial b}L(a, b) = 0 \Rightarrow \frac{2}{m}(mb - \sum_{i=1}^m(y_i -ax_i)) = 0

Оптимальное решение:

a^*=\frac{\sum_{i=1}^my_i(x_i-\bar{x})}{\sum_{i=1}^mx_{i}^2-\frac{1}{m}(\sum_{i=1}^mx_i)^2}
b^*=\frac{1}{m}\sum_{i=1}^m(y_i-ax_i)

в,\bar{x}=\frac{1}{m}\sum_{i=1}^mx_i, который\boldsymbol{x}среднее значение .

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

Множественная линейная регрессия

В более общем случае признаки многомерны:

f(\boldsymbol{x_i}) = \boldsymbol{w}^T \boldsymbol{x_i}

В приведенной выше формуле мы фактически используемb=w_{0}для представления параметраb,будетbдобавленный к вектору признаков,nРазмерная функция расширяется доn+1измерение, то есть в каждом\boldsymbol{x_i}Добавляет элемент со значением 1 в . Формы отдельных векторов и матриц показаны в формулах ниже.

\boldsymbol{w} = \left[  \begin{array}{c} w_0\\ w_1\\ w_2\\ \vdots\\ w_n \\ \end{array} \right]  \quad \boldsymbol{X} = \left[  \begin{array}{c} x_1\\ x_2\\ \vdots\\ x_m \end{array} \right] =  \begin{pmatrix}         1 & x_{1,1} & x_{1,2} & \cdots & x_{1,n} \\         1 & x_{2,1} & x_{2,2} & \cdots & x_{2,n}\\         \vdots & \vdots & \ddots & \vdots & \vdots \\         1 & x_{m,1} & x_{m,2} & \cdots & x_{m,n} \\ \end{pmatrix} \quad \boldsymbol{y} = \left[  \begin{array}{c} y_1\\ y_2\\ \vdots\\ y_m \end{array} \right]

где вектор\boldsymbol{w}Представляет вес каждой функции в модели; матрица\boldsymbol{X}Каждая строка является образцом, и каждый образец имеетn+1собственные значения, которые являются значениями выборки в разных измерениях; вектор\boldsymbol{y}является истинным значением. Суммарный член может быть выражен в виде скалярного произведения:\sum_{j=0}^{n}w_jx_{i,j} =\boldsymbol{w}^\top \boldsymbol{x_i}. В виде матричного умножения это можно выразить так:\boldsymbol{y} = \boldsymbol{X}\boldsymbol{w}.

y_i = \sum_{j=1}^{j+1}w_jx_{i,j} =\boldsymbol{w}^\top \boldsymbol{x_i} \Rightarrow \boldsymbol{y} = \boldsymbol{X}\boldsymbol{w}

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

Обратите внимание, что выделенные жирным шрифтом символы в формуле представляют скаляры, а жирные — векторы или матрицы.

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

多元线性回归一般寻找最优超平面
Множественная линейная регрессия обычно находит оптимальную гиперплоскость

L(\boldsymbol{w}) = (\boldsymbol{X}\boldsymbol{w}-\boldsymbol{y})^\top (\boldsymbol{X}\boldsymbol{w}-\boldsymbol{y}) = ||\boldsymbol{X}\boldsymbol{w}-\boldsymbol{y}||_2^2

Функция потерь множественной линейной регрессии по-прежнему рассчитывается с использованием квадрата «прогнозируемое значение — истинное значение», а приведенная выше формула представляет собой векторное представление функции потерь всей модели. Здесь появляется вертикальная линия, которую называют квадратом нормы L2. Норма обычно используется для векторов, а также является математическим символом, часто используемым в области машинного обучения.Следующая формула показывает вектор\boldsymbol{x}Квадрат L2-нормы и ее производные.

\\ ||\boldsymbol{x}||_2^2 = ((\sum_{i=1}^Nx_i^2)^{\frac{1}{2}})^2 \qquad \qquad \qquad \qquad \nabla||\boldsymbol{x}||_2^2=2\boldsymbol{x}

Для вектора в формуле функции потерь линейной регрессии\boldsymbol{w}Возьмем производную, приравняв производную к нулю:

\frac{\partial}{\partial \boldsymbol{w}}L(\boldsymbol{w}) = 2\boldsymbol{X}^\top (\boldsymbol{X}\boldsymbol{w} - \boldsymbol{y}) = 0
\boldsymbol{X}^\top\boldsymbol{X}\boldsymbol{w} = \boldsymbol{X}^\top \boldsymbol{y} \Rightarrow (\boldsymbol{X}^\top\boldsymbol{X})^{-1} \boldsymbol{X}^\top\boldsymbol{X}\boldsymbol{w}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}
\boldsymbol{w}^* = (\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top \boldsymbol{y}

Приведенная выше формула представляет собой вектор\boldsymbol{w}, которое является матричным уравнением. Использование этого метода для поиска оптимального решения фактически является решением матричного уравнения, которое на английском языке называется Normal Equation.

У этого метода есть проблема, которая должна была быть упомянута в курсе линейной алгебры,\boldsymbol{X}^\top\boldsymbol{X}Система уравнений может быть решена только в том случае, если она является полноранговой или положительно определенной. Что именно означает «полный ранг» или «положительная концентрация»? С точки зрения непрофессионала, данные в выборке должны быть достаточно богатыми и репрезентативными, чтобы иметь единственное решение матричного уравнения, иначе матричное уравнение будет иметь несколько наборов решений. Если функции имеют десятки тысяч измерений, но для обучения используются только десятки выборок, нам трудно получить удовлетворительное оптимальное решение.

У вышеуказанного метода также есть проблема: количество вычислений инверсии матрицы в формуле относительно велико, а сложностьO(n^3)уровень. Когда размерность признака достигает более миллиона или количество отсчетов чрезвычайно велико, время расчета очень велико, а память одного компьютера не может даже хранить эти параметры, и метод решения матричного уравнения нереалистичен.

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

Градиентный спуск

При решении задачи минимизации функции потерь, или решении задачи оптимизации, минимизирующей функцию потерь, часто используется метод перебора. В частности, в качестве начальной точки выбирается начальная точка, затем запускается непрерывный поиск, и функция потерь постепенно уменьшается, и при достижении конечного условия итерации поиска эта позиция является конечным результатом алгоритма поиска. Давайте случайным образом угадаем\boldsymbol{w}, тогда правильно\boldsymbol{w}значение постоянно корректируется, чтобы позволитьL(\boldsymbol{w})постепенно становится меньше, лучше всего найти такой, чтобыL(\boldsymbol{w})самый маленький\boldsymbol{w}.

В частности, мы можем рассмотреть использование метода градиентного спуска (Gradient Descent), этот метод из определенного\boldsymbol{w}Начните с начального значения , а затем постепенно обновляйте веса, или каждый раз перезаписывайте исходные значения вновь вычисленными значениями:

w_j := w_j - \alpha \frac \partial {\partial{w_j}}L(\boldsymbol{w})

здесь\alphaТакже известен как скорость обучения,\frac \partial {\partial{w_j}}L(\boldsymbol{w})это Градиент. В классе исчисления упоминалось, что в определенной точке функция быстрее всего изменяется по направлению градиента. Потому что мы хотим минимизировать функцию потерьL, поэтому каждый раз спускаемся по градиенту,LОпустите самое быстрое направление движения.

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

损失函数沿梯度下降的过程

вернуться к скорости обучения\alphaначальство,\alphaПредставляет, насколько мы уверены в градиенте в определенной точке. В целом,0 < \alpha < 1.\alphaЧем больше значение, тем быстрее мы хотим, чтобы функция потерь падала,\alphaЧем он меньше, тем медленнее мы хотим, чтобы функция потерь падала. если\alphaЕсли он установлен неправильно, каждый раз, когда шаг слишком велик, функция потерь может не быстро сходиться к минимальному значению. Как показано ниже, функция потерь также изо всех сил пытается сходиться к минимальному значению после длительного времени. В практических приложениях\alphaОн часто меняется в зависимости от количества итераций, например, в начале он больше, а затем постепенно становится меньше.

收敛速度过慢的情形

Мы упоминали ранее,\boldsymbol{w}является вектором, предполагая, что онnобъемный, обновление\boldsymbol{w}, мы одновременноnвсе размеры\boldsymbol{w}значение обновляется, гдеjИзмерение должно использовать приведенную выше формулу обновления веса.

Далее мы просто выводим формулу градиента.Во-первых, учтите, что есть только одна обучающая выборка.(\boldsymbol{x_i}, y_j)Случай. Зависит отL(w)=\frac{1}{2}(f(\boldsymbol{x_i})-y_i)^2,в,\frac{1}{2}постоянный член, который не влияет на значение оптимального решения, в основном для удобства вывода. Вы можете получить:

\begin{aligned} \frac \partial {\partial w_j}L(\boldsymbol{w}) & = \frac \partial {\partial w_j} \frac  12(f(\boldsymbol{x_i})-y_i)^2\\ & = 2 \cdot\frac 12(f(\boldsymbol{x_i})-y_i)\cdot \frac \partial {\partial w_j}  (f(\boldsymbol{x_i})-y_i) \\ & = (f(\boldsymbol{x_i})-y_i)\cdot \frac \partial {\partial w_j}(\sum^n_{j=0} w_jx_{i,j}-y_i) \\ & = (f(\boldsymbol{x_i})-y_i) x_{i,j} \end{aligned}

Для одной обучающей выборки правила обновления для каждого градиента следующие:

w_j := w_j - \alpha (f(\boldsymbol{x_i})-y_i) x_{i,j}

У этого правила есть несколько свойств, которые кажутся естественными и интуитивно понятными:

  • обновленный размер с(f(\boldsymbol{x_i})-y_i)пропорциональный.
  • Когда мы сталкиваемся с предсказанным значением обучающей выборкиf(\boldsymbol{x_i})иy_iКогда истинное значение очень близко, вы обнаружите, что в основном нет необходимости изменять параметры; наоборот, если наше предсказанное значениеf(\boldsymbol{x_i})иy_iЕсли истинное значение имеет большую погрешность (например, расстояние очень большое), то параметры нужно больше подгонять. Это также соответствует динамическому графику градиентного спуска, показанному ранее.

пакетный градиентный спуск

Когда имеется только одна обучающая выборка, мы выводим правило LMS. Когда тренировочный наборmПри использовании обучающей выборкиL(\boldsymbol{w}) = \frac{1}{2}\sum_{i=1}^m(f(\boldsymbol{x_i}) - y_i)^2. При выводе необходимо добавить только данные нескольких обучающих выборок.

L(\boldsymbol{w}) = \frac{1}{2} \lbrace (f(\boldsymbol{x_1}) - y_i)^2 + \cdots + (f(\boldsymbol{x_m}) - y_m)^2 \rbrace
\begin{aligned} \frac \partial {\partial w_j}L(w) & = \frac \partial {\partial w_j} \frac  12\lbrace (f(\boldsymbol{x_1}) - y_i)^2 + \cdots + (f(\boldsymbol{x_m}) - y_m)^2 \rbrace\\ & = 2 \cdot\frac 12(f(\boldsymbol{x_1})-y_1)\cdot \frac \partial {\partial w_j}  (f(\boldsymbol{x_1})-y_1) + \cdots\\ & = (f(\boldsymbol{x_1})-y_1)\cdot \frac \partial {\partial w_j}(\sum^n_{j=0} w_jx_{1,j}-y_1) + \cdots \\ & = (f(\boldsymbol{x_1})-y_1) x_{1,j} + \cdots \\ & = \sum_{i=1}^m(f(\boldsymbol{x_i}) - y_i)x_{i,j} \end{aligned}

Таким образом, можно сделать вывод, что каждыйw_jПроизводное от :

w_j := w_j - \alpha \sum^m_{i=1}(f(\boldsymbol{x_{i}})-y_i)x_{i,j}

В частности, алгоритм такой:

\begin{aligned} &\quad Repeat \ k \ iterations \quad \{ \\ & \quad\quad for \ every \ \ w_j \ in \ \boldsymbol{w} \ w_j := w_j - \alpha \sum^m_{i=1}(f(\boldsymbol{x_{i}})-y_i)x_{i,j}\quad \\ &\quad\} \end{aligned}

В этом методе используются все образцы во всем обучающем наборе для обновления параметров на каждой итерации, также известный как пакетный градиентный спуск (BGD). Функция потерь для линейной регрессииLЭто выпуклая квадратичная функция (Convex Quadratic Function), локальное минимальное значение выпуклой функции является глобальным минимальным значением, а оптимизационная задача линейной регрессии имеет только одно глобальное решение. То есть, если предположить, что скорость обучения не\alphaЕсли параметр слишком велик, а количество итераций достаточно велико, метод градиентного спуска всегда будет сходиться к глобальному минимуму.

Стохастический градиентный спуск

Пакетный градиентный спуск учитывает все выборки при обновлении параметров. Когда объем данных большой и признаков много, использовать весь объем данных в каждой итерации нереально, а сам полный объем данных содержит много избыточной информации, Чем больше объем данных, тем больше избыточная информация При поиске оптимального решения избыточная информация мало помогала. Компромисс состоит в том, чтобы каждый раз при обновлении параметров случайным образом брать частичную выборку. В крайнем случае выборка случайным образом выбирается на каждой итерации, и одна выборка используется для обновления параметров этой итерации.Этот алгоритм называется стохастическим градиентным спуском (SGD) следующим образом:

\begin{aligned} &\quad Repeat \ k \ iterations \quad \{ \\ & \quad\quad randomly \ choose \ i \\ & \quad\quad for \ every \ \ w_j \ in \ \boldsymbol{w} \ w_j := w_j - \alpha (f(\boldsymbol{x_{i}})-y_i)x_{i,j}\quad \\ &\quad\} \end{aligned}

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

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

不同梯度下降法的收敛速度示意图

NumPy реализация градиентного спуска

Столько всего было выведено ранее, Обсуждайте дешево, Покажите какой-нибудь код. Затем мы реализуем модель линейной регрессии с использованием NumPy, используя пакетный градиентный спуск и стохастический градиентный спуск соответственно. В процессе реализации мы обнаружим, что некоторые проблемы являются техническими проблемами, которые не будут упомянуты при выводе формулы Например, данные в процессе расчета слишком велики и превышаютfloat64представимый диапазон. Инженерная реализация отражает разницу между теорией и практикой, и на самом деле часто именно эти инженерные детали определяют простоту использования фреймворков машинного обучения.

import numpy as np

class LinearRegression:

    def __init__(self):
        # the weight vector
        self.W = None

    def train(self, X, y, method='bgd', learning_rate=1e-2, num_iters=100, verbose=False):
        """
        Train linear regression using batch gradient descent or stochastic gradient descent

        Parameters
        ----------
        X: training data, shape (num_of_samples x num_of_features), num_of_samples rows of training sample, each training sample has num_of_features-dimension features.
        y: target, shape (num_of_samples, 1). 
        method: (string) 'bgd' for Batch Gradient Descent or 'sgd' for Stochastic Gradient Descent
        learning_rate: (float) learning rate or alpha
        num_iters: (integer) number of steps to iterate for optimization
        verbose: (boolean) if True, print out the progress

        Returns
        -------
        losses_history: (list) of losses at each training iteration
        """
        num_of_samples, num_of_features = X.shape

        if self.W is None:
            # initilize weights with values
            # shape (num_of_features, 1)
            self.W = np.random.randn(num_of_features, 1) * 0.001
        losses_history = []

        for i in range(num_iters):

            if method == 'sgd':
                # randomly choose a sample
                idx = np.random.choice(num_of_samples)
                loss, grad = self.loss_and_gradient(X[idx, np.newaxis], y[idx, np.newaxis])
            else:
                loss, grad = self.loss_and_gradient(X, y)
            losses_history.append(loss)

            # Update weights using matrix computing (vectorized)
            self.W -= learning_rate * grad

            if verbose and i % (num_iters / 10) == 0:
                print('iteration %d / %d : loss %f' %(i, num_iters, loss))
        return losses_history


    def predict(self, X):
        """
        Predict value of y using trained weights

        Parameters
        ----------
        X: predict data, shape (num_of_samples x num_of_features), each row is a sample with num_of_features-dimension features.

        Returns
        -------
        pred_ys: predicted data, shape (num_of_samples, 1)
        """
        pred_ys = X.dot(self.W)
        return pred_ys


    def loss_and_gradient(self, X, y, vectorized=True):
        """
        Compute the loss and gradients

        Parameters
        ----------
        The same as self.train function

        Returns
        -------
        tuple of two items (loss, gradient)
        loss: (float)
        gradient: (array) with respect to self.W 
        """
        if vectorized:
            return linear_loss_grad_vectorized(self.W, X, y)
        else:
            return linear_loss_grad_for_loop(self.W, X, y)


def linear_loss_grad_vectorized(W, X, y):
    """
    Compute the loss and gradients with weights, vectorized version
    """
    # vectorized implementation 
    num_of_samples = X.shape[0]
    # (num_of_samples, num_of_features) * (num_of_features, 1)
    f_mat = X.dot(W)

    # (num_of_samples, 1) - (num_of_samples, 1)
    diff = f_mat - y 
    loss = 1.0 / 2 * np.sum(diff * diff)
    
    # {(num_of_samples, 1).T dot (num_of_samples, num_of_features)}.T
    gradient = ((diff.T).dot(X)).T

    return (loss, gradient)


def linear_loss_grad_for_loop(W, X, y):
    """
    Compute the loss and gradients with weights, for loop version
    """
    
    # num_of_samples rows of training data
    num_of_samples = X.shape[0]
    
    # num_of_samples columns of features
    num_of_features = X.shape[1]
    
    loss = 0
    
    # shape (num_of_features, 1) same with W
    gradient = np.zeros_like(W) 
    
    for i in range(num_of_samples):
        X_i = X[i, :] # i-th sample from training data
        f = 0
        for j in range(num_of_features):
            f += X_i[j] * W[j, 0]
        diff = f - y[i, 0]
        loss += np.power(diff, 2)
        for j in range(num_of_features):
            gradient[j, 0] += diff * X_i[j]
            
    loss = 1.0 / 2 * loss

    return (loss, gradient)

использованная литература

  1. Andrew Ng: CS229 Lecture Notes
  2. Чжоу Чжихуа: «Машинное обучение»
  3. Post-Line Sequence.GitHub.IO/2015/03/31/…