Метод нормальных уравнений и метод Ньютона

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

Помимо градиентного спуска, для наших целей можно применить и другие подходы: метод нормальных уравнений и метод Ньютона. Этот пост – о них.

Метод нормальных уравнений

Когда я был в Массачусетском технологическом институте, я часто любил подшучивать над людьми. Однажды в кабинете черчения какой-то шутник поднял лекало (кусок пластмассы для рисования гладких кривых — забавно выглядящая штука в завитушках) и спросил: «Имеют ли кривые на этих штуках какую-либо формулу?»

Я немного подумал и ответил: «Несомненно. Это такие специальные кривые. Дай-ка я покажу тебе. — Я взял свое лекало и начал его медленно поворачивать. — Лекало сделано так, что, независимо от того, как ты его повернешь, в наинизшей точке каждой кривой касательная горизонтальна».

Все парни в кабинете начали крутить свои лекала под различными углами, подставляя карандаш к нижней точке и по-всякому прилаживая его. Несомненно, они обнаружили, что касательная горизонтальна. Все были крайне возбуждены от этого открытия, хотя уже много прошли по математике и даже «выучили», что производная (касательная) в минимуме (нижней точке) для любой кривой равна нулю (горизонтальна). Они не совмещали эти факты. Они не знали даже того, что они уже «знали».

(с) Вы, конечно, шутите, мистер Фейнман

Да, мы знаем, что производная функции в точке минимума равна нулю. Этот минимум нам-то и нужен. Функция стоимости имеет вид некой «чаши», и минимум у нее всегда один:

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

Но перед этим давайте посмотрим, как можно представить функцию стоимости в матричном виде – это сильно упростит вычисления, а кроме того, позволит реализовать метод программно более просто и ускорить выполнение кода, используя библиотеки для работы с матрицами вроде Python Numpy.

Итак, наша функция стоимости:

J ( \theta_{0},\theta_{1} )=\frac{1}{2}\sum_{i=1}^{m}\left ( h_{\theta }(x^{(i)})-y^{(i)} \right )^{2}

Эту же формулу можно записать в векторной форме. Сумма от 1 до m – это разность между предсказанной стоимостью дома и реальной ценой для всех m домов, которые мы продали в прошлом. Предсказанную стоимость мы вычисляем по формуле h(x)=\theta_{0}+\theta_{1}x_{1}, где x_{1} – площадь дома. Давайте добавим в список параметров, по которым мы оцениваем дом, кроме его площади, дополнительный параметр x_{0}, который всегда будет равен единице. Тогда нашу формулу можно будет переписать так:

h(x)=\theta_{0}x_{0}+\theta_{1}x_{1}

Функция стоимости же будет выглядеть следующим образом:

J ( \theta_{0},\theta_{1} )=\frac{1}{2}((\theta_{0}x_{0}^{(2)}-y^{(1)})^{1}+(\theta_{0}x_{0}^{(2)}-y^{(2)})^{2}+...+(\theta_{0}x_{0}^{(m)}-y^{(m)})^{2})

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

A^{(n\times n)}*B^{(n\times 1)}=\begin{bmatrix}    \sum_{i=1}^{n}A_{1,i}*B_{i}\\    \sum_{i=1}^{n}A_{2,i}*B_{i}\\    ...\\    \sum_{i=1}^{n}A_{n,i}*B_{i}    \end{bmatrix}

Разница двух векторов выглядит как:

A^{(n\times 1)}-B^{(n\times 1)}=\begin{bmatrix}    A_{1}-B_{1}\\    A_{2}-B_{2}\\    ...\\    A_{n}-B_{n}    \end{bmatrix}

Удобно, не правда ли? Идеально применимо для наших целей. Запишем параметры наших домов как матрицу X. Не забудьте, что мы добавили дополнительный параметр 1, и теперь у нас каждый дом характеризуется двумя параметрами: единицей и площадью. Матрица будет иметь размерность m колонок на 2 столбца, где m – количество проданных ранее домов.

X=\begin{bmatrix}    1 & 2104\\    1 & 1416\\    ... & ...\\    1 & size \; of\; house\; m    \end{bmatrix}

Параметры \theta зададим как вектор:

\theta =\begin{bmatrix}    \theta_{0}\\    \theta_{1}    \end{bmatrix}

Количество проданных ранее домов (исторических примеров) равно m. Если умножить матрицу X размером m x 2 на вектор \theta размером 2 x 1, мы получим вектор размером m x 1, каждый элемент которого будет предсказанным значением стоимости дома с использованием параметров модели \theta.

Если от этого вектора отнять вектор y, размером m x 1, получим вектор размера m x 1, в котором каждый элемент будет разницей между предсказанной стоимостью и реальной ценой.

И теперь, если мы возведем этот вектор в квадрат и умножим на \frac{1}{2}, мы получим вектор, снова-таки, размером m x 1, каждый элемент которого будет функцией стоимости для конкретного проданного дома.

Матрицы возводятся в квадрат по формуле: A^{T}A. Если не верите, посчитайте сами 🙂

Итого, наша функция стоимости может быть записана в виде:

J ( \theta_{0},\theta_{1} )=\frac{1}{2}\sum_{i=1}^{m}\left ( h_{\theta }(x^{(i)})-y^{(i)} \right )^{2}=(X\theta-y)^{T}(X\theta-y)

Дальше, нам надо взять от нее градиент.

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

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

Для начала, нам понадобится оператор trace. Применительно к матрице, trace – это сумма всех ее диагональных элементов. Результат ее – скалярное значение.

У этого оператора есть ряд интересных свойств:

trABC=trCAB=trBCA\; (1) (мы можем переместить последний элемент на место первого)

trA=trA^{T}\; (2)

tr(A+B)=trA+trB\; (3)

tra*A=a*trA\; (4)

И от этого же оператора очень удобно считать градиенты:

\bigtriangledown _{A^{T}}f(A)=(\bigtriangledown _{A}f(A))^{T}\; (5)

\bigtriangledown _{A}ABA^{T}C=CAB+C^{T}AB^{T}\; (6)

\bigtriangledown _{A}trAB=B^{T}\; (7)

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

\bigtriangledown _{\theta} J(\theta )=\bigtriangledown _{\theta}(X\theta-y)^{T}(X\theta-y)=\bigtriangledown _{\theta}((X\theta )^{T}-y^{T})(X\theta-y)=(\theta ^{T}X^{T}-y^{T})(X\theta-y)=\bigtriangledown _{\theta}(\theta ^{T}X^{T}X\theta-\theta ^{T}X^{T}y-y^{T}X\theta +y^{T}y)

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

\bigtriangledown _{\theta}(\theta ^{T}X^{T}X\theta-\theta ^{T}X^{T}y-y^{T}X\theta +y^{T}y)=\bigtriangledown _{\theta}tr(\theta ^{T}X^{T}X\theta-\theta ^{T}X^{T}y-y^{T}X\theta +y^{T}y)

Мы берем градиент от параметров \theta. Последний элемент полинома y^{T}y не зависит от \theta, так что его производная от любого из \theta будет равна нулю, и мы можем исключить его из уравнения.

Дальше, используя формулу (3):

\bigtriangledown _{\theta}tr(\theta ^{T}X^{T}X\theta-\theta ^{T}X^{T}y-y^{T}X\theta)=\bigtriangledown _{\theta}(tr\theta ^{T}X^{T}X\theta-tr\theta ^{T}X^{T}y-try^{T}X\theta)

Дальше, используя формулу (2) и правила работы с транспонированными матрицами:

\bigtriangledown _{\theta}(tr\theta ^{T}X^{T}X\theta-tr\theta ^{T}X^{T}y-try^{T}X\theta)=\bigtriangledown _{\theta}(tr\theta ^{T}X^{T}X\theta-tr\theta ^{T}X^{T}y-tr(y^{T}X\theta)^{T})=\bigtriangledown _{\theta}(tr\theta ^{T}X^{T}X\theta-tr\theta ^{T}X^{T}y-tr\theta ^{T}X^{T}y)=\bigtriangledown _{\theta}(tr\theta ^{T}X^{T}X\theta-2tr\theta ^{T}X^{T}y)

Используя формулу (5) и считая A=\theta, B=X^{T}X, C=I — identity matrix, и формулу (7):

\bigtriangledown _{\theta}(tr\theta ^{T}X^{T}X\theta-2tr\theta ^{T}X^{T}y)=\frac{1}{2}(X^{T}X\theta + X^{T}X\theta - 2X^{T}y)=X^{T}X\theta-X^{T}y

Все, с этого места можете снова начинать читать 🙂

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

X^{T}X\theta-X^{T}y=0

И выразить из этого уравнения \theta:

X^{T}X\theta=X^{T}y

\theta=(X^{T}X)^{-1}X^{T}y

Вот по этой формуле, зная X и y, вы можете в один ход аналитически посчитать нужные вам значения \theta.

В чем плюсы и минусы этого метода, сравнительно с градиентным спуском? Плюс: он считает минимум в один ход и не требует итераций. Минус: он требует вычисления обратной матрицы. Сложность этого вычисления – порядка O(n^{3}) , и на каком-то этапе считать обратную матрицу становится сложнее, чем использовать градиентный спуск. Специалисты в МЛ говорят, что использовать метод линейных уравнений стоит, если у вас имеется до 10 тыс. исторических примеров. Если же больше – имеет смысл подумать об итерационных методах.

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

Метод Ньютона

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

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

Итерационно? Итерационно. Удобно? Еще как! Минусы: касательная к производной – это производная второго порядка. Производная второго порядка от функции многих переменных называется гессиан.

Давайте посмотрим на него вблизи. У нас есть градиент. Градиент от функции – это вектор, имеющий столько же элементов, сколько у функции переменных:

\bigtriangledown _{\theta }J(\theta )=\begin{bmatrix}    \frac{\partial }{\partial \theta _{1}}J(\theta )\\    \frac{\partial }{\partial \theta _{2}}J(\theta )\\    ...\\    \frac{\partial }{\partial \theta _{n}}J(\theta )    \end{bmatrix}

Теперь нам надо посчитать Гессиан от этого градиента. Для этого нужно от каждой частной производной в градиенте взять частные производные, снова-таки, по каждому из параметров \theta. Получим матрицу размерностью n x n , где n – количество переменных в \theta:

\bigtriangledown_{\theta }^{2}J(\theta )=\begin{bmatrix}    \frac{\partial }{\partial \theta _{1}\partial \theta _{1}}J(\theta ) & \frac{\partial }{\partial \theta _{1}\partial \theta _{2}}J(\theta ) & ... & \frac{\partial }{\partial \theta _{1}\partial \theta _{n}}J(\theta )\\    \frac{\partial }{\partial \theta _{2}\partial \theta _{1}}J(\theta ) & \frac{\partial }{\partial \theta _{2}\partial \theta _{2}}J(\theta ) & ... & \frac{\partial }{\partial \theta _{2}\partial \theta _{n}}J(\theta )\\    ... & ... & ... & ...\\    \frac{\partial }{\partial \theta _{n}\partial \theta _{1}}J(\theta ) & \frac{\partial }{\partial \theta _{n}\partial \theta _{2}}J(\theta ) & ... & \frac{\partial }{\partial \theta _{n}\partial \theta _{n}}J(\theta )    \end{bmatrix}

Мы знаем, что градиент функции стоимости для линейной регрессии равен:

\sum_{j=1}^{m}(h_{\theta }(x^{(j)})-y^{(j)})x_{i}

\bigtriangledown_{\theta }J(\theta ) = \begin{bmatrix}    \sum_{j=1}^{m}(h_{\theta }(x^{(j)})-y^{(j)})x_{1}\\    \sum_{j=1}^{m}(h_{\theta }(x^{(j)})-y^{(j)})x_{2}\\    ...\\    \sum_{j=1}^{m}(h_{\theta }(x^{(j)})-y^{(j)})x_{n}    \end{bmatrix}

Если брать вторую частную производную и считать Гессиан, получим матрицу:

\bigtriangledown _{\theta}^{2} J(\theta ) = \begin{bmatrix}    x_{1}^{2} & x_{1}x_{2} & ... & x_{1}x_{n}\\    x_{2}x_{1} & x_{2}^{2} & ... & x_{2}x_{n}\\    ... & ... & ... & ...\\    x_{n}x_{1} & x_{n}x_{2} & ... & x_{n}^{2}    \end{bmatrix}

Итого – гессиан нашей функции стоимости – квадратная симметричная матрица размером n x n. В матричном виде ее можно представить как \bigtriangledown _{\theta}^{2} J(\theta ) = X^{T}X

Возвращаясь к реальности: все это нам нужно было, чтобы найти минимум функции J(\theta ) . Улучшенный метод Ньютона для функций от нескольких переменных называется методом Ньютона-Рафсона. После того, как вы высчитали градиент и гессиан функции, вы итерационно обновляете параметры \theta:

\theta := \theta - H^{-1}\bigtriangledown _{\theta }J(\theta )

Метод Ньютона обычно сходится быстрее, чем (batch) градиентный спуск, но каждый его «шаг» — более сложен вычислительно, так как нужно инвертировать гессиан. Тем не менее, пока количество переменных не очень велико, он показывает лучшую производительность в поиске минимума функции.

5 комментариев на “Метод нормальных уравнений и метод Ньютона

      1. tassadarius Автор записи

        Нет, вот этот:
        http://cs229.stanford.edu/

        Курс Курсеры сравнительно с CS229:
        «Compared to CS229 (Machine Learning), we cover fewer learning algorithms, and also spend less time on the math and theory of machine learning, but spend much more time on the pratical, hands-on skills (and «dirty tricks») for getting this stuff to work well on an application.»

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход /  Изменить )

Google photo

Для комментария используется ваша учётная запись Google. Выход /  Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход /  Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход /  Изменить )

Connecting to %s