Линейная регрессия. Часть 3

Материалы > Обучение модели

Продолжим изучать линейную регрессию.

Линейная регрессия как проекция

Представим данные как систему линейных уравнений.

На рисунке выше представлена обычная задача простой линейной регрессии с немного измененной нотацией. У нас есть пять наблюдений, каждое наблюдение описано признаками $a_0$ и $a_1,$ при этом $a_0$ заполнен единицами. Подставив некоторые веса $\theta_0$ и $\theta_1,$ мы получим для каждой точки значение целевой переменной $b,$ ранее мы ее обозначали через $y.$

Единственного решения для такой системы не существует, поскольку мы не сможем подобрать такие $\theta_0$ и $\theta_1,$ которые бы удовлетворяли всем $a_0$ и $a_1$ в каждой строке (для каждого наблюдения).

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

Перепишем эту систему с помощью векторов и матриц. В общем случае у нас конечно может быть больше признаков $\mathbf a_0, \mathbf a_1, \mathbf a_2,.., \mathbf a_d,$ однако систему по-прежнему будем считать переопределенной, то есть $n > d.$ Также заменим $\boldsymbol \theta$ на $\mathbf x,$

$$ \begin{bmatrix} \vdots & \vdots & \vdots & \vdots \\ \mathbf a_1 & \mathbf a_2 & \dots & \mathbf a_d \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_d \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix} $$

$$ \underset{n \times d}{ \mathbf A} \mathbf x = \mathbf b. $$

При этом в данном случае $\mathbf A, \mathbf x \in \mathbb{R}^d,$ а $\mathbf b \in \mathbb{R}^n.$ То есть вектор $\mathbf b$ находится в пространстве большей размерности, чем $\mathbf {x}.$ Например, вектор $\mathbf x$ может находиться на плоскости (двумерный вектор), а вектор $\mathbf b$ в трехмерном пространстве.

Можно также сказать, что не существует линейной комбинации значений вектора $\mathbf x$ (веса модели) и векторов $\mathbf a_0, \mathbf a_1, \mathbf a_2,.., \mathbf a_d,$ которые преобразовывались бы в вектор $\mathbf {b},$

$$ x_1 \mathbf a_1 + x_2 \mathbf a_2 + x_3 \mathbf a_3 + … + x_d \mathbf a_d \not= \mathbf b. $$

Наконец, справедливо, что $ \mathbf b $ не находится в пространстве столбцов $\mathbf A,$ $\mathbf b \not\in col(\mathbf A)$ (то есть не находится на плоскости), а вектор $ \mathbf x $ как раз лежит на этой плоскости, и это значит, что мы никак не можем перевести с помощью матрицы $\mathbf A$ вектор $ \mathbf x $ из двумерного пространства в трехмерное.

С другой стороны, мы можем попробовать получить наилучшее возможное решение, найдя такой вектор $ \mathbf x^*{,} $ который будет максимально приближен к вектору $ \mathbf b $ (т.е. будет иметь минимальное расстояние до него). Расстояние же между векторами можно определить как разницу двух векторов.

Для того чтобы положительные и отрицательные значения не взаимоудалялись, возведем значения в квадрат,

$$ min || \mathbf b-\mathbf A \mathbf x^* ||^2. $$

Очевидно, что вектор (назовем его $ \mathbf p $), получившийся в результате $\mathbf A \mathbf x^*,$ будет в пространстве $\mathbb{R}^d,$ что то же самое, что $ \mathbf p \in col(\mathbf A) $ (то есть в данном случае на плоскости).

Далее, наименьшее расстояние от вектора $ \mathbf b $ до $ \mathbf p $ можно определить как ортогональную проекцию $ \mathbf b $ на пространство столбцов $\mathbf A,$

$$ \mathbf A \mathbf x^* = \mathbf p = proj_{col(A)} \mathbf b. $$

Попробуем найти решение относительно $\mathbf b,$

$$ \mathbf A \mathbf x^* = proj_{col(A)} \mathbf b. $$

Вычтем вектор $\mathbf b$ из обеих частей,

$$ \mathbf A \mathbf x^*-\mathbf b = proj_{col(A)} \mathbf b-\mathbf b. $$

Заметим, что вектор $ proj_{col(A)} \mathbf b-\mathbf b $ (на рисунке представлен красным вектором) ортогонален к плоскости $col(\mathbf A)$. Как следствие, $\mathbf A \mathbf x^*-\mathbf b$ ортогонально $col(\mathbf A).$

Можно также сказать, что $\mathbf A \mathbf x^*-\mathbf b$ является ортогональным дополнением пространства $col(\mathbf A).$ Запишем это как $\mathbf A \mathbf x^*-\mathbf b = col(\mathbf A)^{\perp}.$

Одновременно ортогональное дополнение пространства столбцов матрицы $\mathbf A$ равно ядру $\mathbf A^\top,$ то есть $col(\mathbf A)^{\perp} = null(\mathbf A^\top)$. Тогда,

$$ \mathbf A \mathbf x^*-\mathbf b = col(\mathbf A)^{\perp} = null(\mathbf A^\top). $$

Если умножить матрицу $\mathbf A^\top$ на ее ядро $ \mathbf A \mathbf x^*-\mathbf {b}, $ то мы получим нулевой вектор,

$$ \begin{split} \mathbf A^\top (\mathbf A \mathbf x^*-\mathbf b) &= \mathbf 0 \\[1ex] \mathbf A^\top \mathbf A \mathbf x^*-\mathbf A^\top \mathbf b &= \mathbf 0 \\[1ex] \mathbf A^\top \mathbf A \mathbf x^* &= \mathbf A^\top \mathbf b. \end{split} $$

Таким образом можно найти $\mathbf x^*,$ которое минимизирует квадрат расстояния между вектором $ \mathbf b $ и вектором проекции $\mathbf A \mathbf x^*.$

---

Уверен, вы узнали в этом выражении нормальные уравнения,

$$ \mathbf X^\top \mathbf X \boldsymbol \theta = \mathbf X^\top \mathbf y. $$

Нормальными эти уравнения называются, потому что $ \mathbf A \mathbf x^*-\mathbf b \perp col(\mathbf A), $ а нормалью⧉ в геометрии как раз считается обобщенное понятие перпендикуляра к поверхности.

---

Более того, можно сказать, что минимизация расстояний от точек до прямой вдоль оси $y$ одновременно приводит к минимизации длины перпендикуляров к проекциям точек.

Пример

Возьмем систему уравнений,

$$ \begin{cases} 2x &-& y &=& 2 \\ x &+& 2y &=& 1 \\ x &+& y &=& 4. \end{cases} $$

Такая система опять же могла бы описывать данные, где x и y представляют собой некоторые коэффициенты без сдвига, например, $\theta_1$ и $\theta_2,$ а сами уравнения — объекты (наблюдения).

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

Продолжим работу в том же блокноте⧉

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

Опять же решить систему методом Гаусса у нас не получится. То есть $ \mathbf A \mathbf x = \mathbf b $ решения не имеет. Однако мы можем найти приближенное решение через $ \mathbf A^\top \mathbf A \mathbf x^* = \mathbf A^\top \mathbf b. $ Найдем $\mathbf A^\top \mathbf A$ и $\mathbf A^\top \mathbf b.$

Таким образом, мы получили

$$ \begin{bmatrix} 6 & 1 \\ 1 & 6 \end{bmatrix} \mathbf x^* = \begin{bmatrix} 9 \\ 4 \end{bmatrix}. $$

Решим систему (а такая система уже решается) методом Гаусса. Применим следующие элементарные преобразования,

• $R_1 = \frac{R_1}{6};$
• $R_2=R_2-R_1;$
• $R_2=\frac{6R_2}{35};$
• $R_1=R_1-\frac{R_2}{6}.$

В результате получим следующую матрицу,

$$ \begin{bmatrix} 1 & 0 & \frac{10}{7} \\ 0 & 1 & \frac{3}{7} \end{bmatrix} \approx \begin{bmatrix} 1 & 0 & 1{,}43 \\ 0 & 1 & 0{,}43 \end{bmatrix}. $$

Таким образом, вектор $ \begin{bmatrix} 1{,}43 \\ 0{,}43 \end{bmatrix} $ будет решением этой системы уравнений.

Обозначим его на графике.

Решим эту задачу на Питоне. Вначале с помощью функции np.linalg.lstsq().

Теперь с помощью класса LinearRegression библиотеки sklearn (напомню, что этот класс также использует метод наименьших квадратов).

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

Обратимость матрицы $\mathbf A^\top \mathbf A$

Общее решение требует нахождения обратной матрицы $(\mathbf A^\top \mathbf A)^{-1},$

$$ \mathbf x = (\mathbf A^\top \mathbf A)^{-1} \mathbf A^\top \mathbf b. $$

Посмотрим, в каких случаях матрица $\mathbf A^\top \mathbf A$ обратима.

Такая матрица будет обратима, если столбцы $\mathbf A$ линейно независимы (то есть матрица имеет полный ранг).

Приведем простой пример двух матриц $\mathbf A_1$ и $\mathbf A_2$ размерностью $3 \times 2$ с линейно независимыми (полный ранг) и линейно зависимыми (ранг 1) столбцами,

$$ \mathbf A_1 = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 5 \end{bmatrix}, \hspace{5pt} \mathbf A_2 = \begin{bmatrix} 1 & 3 \\ 1 & 3 \\ 1 & 3 \end{bmatrix}. $$

Найдем $\mathbf A_1^\top \mathbf A_1$ и $\mathbf A_2^\top \mathbf A_2,$

$$ \mathbf A_1^\top \mathbf A_1 = \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 5 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 5 \end{bmatrix} = \begin{bmatrix} 3 & 8 \\ 8 & 30 \end{bmatrix} $$

$$ \mathbf A_2^\top \mathbf A_2 = \begin{bmatrix} 1 & 1 & 1 \\ 3 & 3 & 3 \end{bmatrix} \begin{bmatrix} 1 & 3 \\ 1 & 3 \\ 1 & 3 \end{bmatrix} = \begin{bmatrix} 3 & 9 \\ 9 & 27 \end{bmatrix}. $$

Очевидно, в первом случае матрица обратима, во втором — нет. Приведем и более формальное доказательство.

Предположим, что $\mathbf A^\top \mathbf A \mathbf x = \mathbf {0}.$ Тогда, чтобы столбцы $\mathbf A^\top \mathbf A$ были линейно независимы, должно выполняться условие, при котором ядро содержит только нулевой вектор $null(\mathbf A) = \{ \mathbf 0 \}$ и $ \mathbf x = \mathbf {0} .$ Умножим обе части на $ \mathbf x^\top,$

$$ \mathbf x^\top \mathbf A^\top \mathbf A \mathbf x = \mathbf x^\top \mathbf 0. $$

Получим умножение двух векторов (их результатом будет скаляр, а не нулевой вектор),

$$ \begin{split} (\mathbf A \mathbf x)^\top (\mathbf A \mathbf x) &= 0 \\[1ex] || \mathbf A \mathbf x ||^2 &= 0. \end{split} $$

Значит,

$$ \mathbf A \mathbf x = \mathbf 0. $$

Так как из $\mathbf A^\top \mathbf A \mathbf x = \mathbf 0$ следует, что $\mathbf A \mathbf x = \mathbf {0},$ то можно сказать, что если $\mathbf x \in null(\mathbf A^\top \mathbf A),$ то $\mathbf x \in null(\mathbf A).$ Отсюда,

$$ null(\mathbf A^\top \mathbf A) = null(\mathbf A). $$

Из этого следует, что

$$ dim(null(\mathbf A^\top \mathbf A)) = dim(null(\mathbf A)) \implies rank(\mathbf A^\top \mathbf A) = rank(\mathbf A). $$

В целом, как известно, матрица обратима, если она представляет собой квадратную матрицу, и ее столбцы независимы (она имеет полный ранг). Мы знаем, что $\mathbf A^\top \mathbf A$ квадратная и имеет тот же ранг, что и $\mathbf A.$ Соответственно, если $\mathbf A$ имеет полный ранг, то $\mathbf A^\top \mathbf A$ обратима.

Добавлю, что если матрицу $\mathbf A^\top \mathbf A$ возможно привести к упрощенному ступенчатому виду (единичной матрице), т.е. $ rref(\mathbf A^\top \mathbf A) = \mathbf I, $ то она обратима.

Про матрицу проекции

Выше мы сказали, что $ \mathbf p = \mathbf A \mathbf x^* .$ Одновременно, решив нормальные уравнения, получим $ \mathbf x^* = (\mathbf A^\top \mathbf A)^{-1} \mathbf A^\top \mathbf {b}. $ Как следствие,

$$ \mathbf p = \mathbf A \mathbf x^* = \mathbf A (\mathbf A^\top \mathbf A)^{-1} \mathbf A^\top. $$

Свойства матрицы проекции

Так мы получили матрицу проекции (projection matrix). Обозначим ее как $ \mathbf P .$ Посмотрим на несколько интересных свойств. Во-первых, такая матрица симметрична,

$$ \mathbf P^\top = \mathbf P. $$

Во-вторых, квадрат матрицы проекции равен самой матрице проекции $ \mathbf P^2 = \mathbf P ,$

$$ \mathbf A (\mathbf A^\top \mathbf A)^{-1} \cancel{\mathbf A^\top \mathbf A} \cancel{(\mathbf A^\top \mathbf A)^{-1}} \mathbf A^\top = \mathbf A (\mathbf A^\top \mathbf A)^{-1} \mathbf A^\top. $$

Это же справедливо для любого количества степеней.

$\mathbf b \in col(\mathbf A)$

Что интересно, если $\mathbf A \mathbf x = \mathbf b $ имеет решение, то есть $\mathbf b \in col(\mathbf A),$ то $\mathbf A$ — будет квадратной обратимой матрицей, а значит, по свойству $(\mathbf A\mathbf B)^{-1} = \mathbf A^{-1}\mathbf B^{-1},$

$$ \mathbf A\mathbf A^{-1} (\mathbf A^\top)^{-1} \mathbf A^\top = \mathbf I. $$

Другими словами, если $\mathbf A \mathbf x = \mathbf b $ проецирует $\mathbf b$ на $\mathbf A$ и $\mathbf b \in col(\mathbf A),$ то проекцией будет сам вектор $\mathbf {b},$ а скалярной проекцией единица (единичная матрица),

$$ \mathbf P \mathbf b = \mathbf I \mathbf b = \mathbf b. $$

Можно также сказать, что $ \mathbf A \mathbf x = \mathbf b, $

$$ \mathbf p = \mathbf P \mathbf b = (\mathbf A (\mathbf A^\top \mathbf A)^{-1} \mathbf A^\top) \mathbf A \mathbf x = \mathbf A \cancel{(\mathbf A^\top \mathbf A)^{-1}} \cancel{(\mathbf A^\top \mathbf A)} \mathbf x = \mathbf A \mathbf x. $$

$\mathbf b \perp col(\mathbf A)$

Если вектор $\mathbf b$ ортогонален пространству столбцов $col(\mathbf A),$ то он находится в ядре $\mathbf A^\top.$ Логично, что проекция $\mathbf b$ на $col(\mathbf A)$ представляет собой нулевой вектор,

$$ \mathbf P \mathbf b = \mathbf 0. $$

Метод максимального правдоподобия

Рассмотрим модель простой линейной регрессии с независимыми одинаково распределенными наблюдениями,

$$ y_i = \theta_0 + \theta_1 x_i + \varepsilon_i, \quad \varepsilon_i \sim \mathcal{N}(0, \sigma^2), \; i = 1,\ldots,n. $$

Проиллюстрируем эту модель для трех наблюдений $x_1, x_2, x_3.$

Применим метод максимального правдоподобия (maximum likelihood estimation) для оценки параметров $\theta_0, \theta_1$ и $\sigma^2,$

$$ \begin{split} P(y \mid \theta_0, \theta_1,\sigma^2) &= \prod_{i=1}^n P(y_i \mid \theta_0,\theta_1,\sigma^2) \\ &= \prod_{i=1}^n \mathcal{N}(y_i; \theta_0 + \theta_1 x_i, \sigma^2) \\ &= \prod_{i=1}^n \frac{1}{\sqrt{2 \pi \sigma}} \cdot \exp \left[ -\frac{(y_i-\theta_0-\theta_1 x_i)^2}{2 \sigma^2} \right] \\ &= \frac{1}{\sqrt{(2 \pi \sigma^2)^n}} \cdot \exp\left[ -\frac{1}{2 \sigma^2} \sum_{i=1}^n (y_i-\theta_0-\theta_1 x_i)^2 \right]. \end{split} $$

Возьмем логарифм от функции правдоподобия,

$$ \begin{split} \ell(\theta_0,\theta_1,\sigma^2) &= \log P(y \mid \theta_0,\theta_1,\sigma^2) \\ &= -\frac{n}{2} \log(2\pi)-\frac{n}{2} \log (\sigma^2) -\frac{1}{2 \sigma^2} \sum_{i=1}^n (y_i-\theta_0 — \theta_1 x_i)^2. \end{split} $$

Найдем производную $\ell$ относительно $\theta_0,$

$$ \frac{\partial \ell (\theta_0,\theta_1,\sigma^2)}{\partial \theta_0} = \frac{1}{\sigma^2} \sum_{i=1}^n (y_i-\theta_0-\theta_1 x_i). $$

Приравняем к нулю и выразим $\theta_0,$

$$ \begin{split} \frac{\partial \ell (\hat{\theta}_0,\hat{\theta}_1,\hat{\sigma}^2)}{\partial \theta_0} &= 0 \\ 0 &= \frac{1}{\hat{\sigma}^2} \sum_{i=1}^n (y_i-\hat{\theta}_0-\hat{\theta}_1 x_i) \\ 0 &= \sum_{i=1}^n y_i-n \hat{\theta}_0-\hat{\theta}_1 \sum_{i=1}^n x_i \ \hat{\theta}_0 &= \frac{1}{n} \sum_{i=1}^n y_i-\hat{\theta}_1 \frac{1}{n} \sum_{i=1}^n x_i \\ \hat{\theta}_0 &= \bar{y}-\hat{\theta}_1 \bar{x}. \end{split} $$

Найдем производную $\ell$ со значением $\hat{\theta}_0$ относительно $\theta_1,$

$$ \frac{\partial \ell (\hat{\theta}_0,\theta_1,\sigma^2)}{\partial \theta_1} = \frac{1}{\sigma^2} \sum_{i=1}^n (x_i y_i-\hat{\theta}_0 x_i-\theta_1 x_i^2). $$

Приравняем к нулю и выразим $\theta_1,$

$$ \begin{split} \frac{\partial \ell (\hat{\theta}_0,\hat{\theta}_1,\hat{\sigma}^2)}{\partial \theta_1} &= 0 \\ 0 &= \frac{1}{\hat{\sigma}^2} \sum_{i=1}^n (x_i y_i-\hat{\theta}_0 x_i-\hat{\theta}_1 x_i^2) \\ 0 &= \sum_{i=1}^n x_i y_i-\hat{\theta}_0 \sum_{i=1}^n x_i-\hat{\theta}_1 \sum_{i=1}^n x_i^2) \\ 0 &= \sum_{i=1}^n x_i y_i-(\bar{y}-\hat{\theta}_1 \bar{x}) \sum_{i=1}^n x_i-\hat{\theta}_1 \sum_{i=1}^n x_i^2 \\ 0 &= \sum_{i=1}^n x_i y_i-\bar{y} \sum_{i=1}^n x_i + \hat{\theta}_1 \bar{x} \sum_{i=1}^n x_i-\hat{\theta}_1 \sum_{i=1}^n x_i^2 \\ 0 &= \sum_{i=1}^n x_i y_i-n \bar{x} \bar{y} + \hat{\theta}_1 n \bar{x}^2-\hat{\theta}_1 \sum_{i=1}^n x_i^2 \\ \hat{\theta}_1 &= \frac{\sum_{i=1}^n x_i y_i-\sum_{i=1}^n \bar{x} \bar{y}}{\sum_{i=1}^n x_i^2-\sum_{i=1}^n \bar{x}^2} \\ \hat{\theta}_1 &= \frac{\sum_{i=1}^n (x_i-\bar{x}) (y_i-\bar{y})}{\sum_{i=1}^n (x_i-\bar{x})^2}. \end{split} $$

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

Осталось получить MLE-оценку $\sigma^2.$ Найдем производную $\ell$ со значением $(\hat{\theta}_0, \hat{\theta}_1)$ относительно $\theta_1,$

$$ \frac{\partial \ell (\hat{\theta}_0,\hat{\theta}_1,\sigma^2)}{\partial \sigma^2} = -\frac{n}{2\sigma^2} + \frac{1}{2(\sigma^2)^2} \sum_{i=1}^n (y_i-\hat{\theta}_0-\hat{\theta}_1 x_i)^2. $$

Приравняем производную к нулю и выразим $\sigma^2,$

$$ \begin{split} \frac{\partial \ell (\hat{\theta}_0,\hat{\theta}_1,\hat{\sigma}^2)}{\partial \sigma^2} &= 0 \\ 0 &= -\frac{n}{2\hat{\sigma}^2} + \frac{1}{2(\hat{\sigma}^2)^2} \sum_{i=1}^n (y_i-\hat{\theta}_0-\hat{\theta}_1 x_i)^2 \\ \frac{n}{2\hat{\sigma}^2} &= \frac{1}{2(\hat{\sigma}^2)^2} \sum_{i=1}^n (y_i-\hat{\theta}_0-\hat{\theta}_1 x_i)^2 \\ \hat{\sigma}^2 &= \frac{1}{n} \sum_{i=1}^n (y_i-\hat{\theta}_0-\hat{\theta}_1 x_i)^2. \end{split} $$

Теперь рассмотрим модель множественной линейной регрессии с независимыми одинаково распределенными наблюдениями,

$$ \mathbf y = \mathbf X \boldsymbol \theta + \varepsilon, \varepsilon \sim \mathcal{N}(0, \sigma^2 \mathbf I), $$

где остатки следуют многомерному нормальному распределению с нулевым средним и дисперсией, равной $\sigma^2 \mathbf I,$

$$ \begin{split} P(\mathbf y \mid \boldsymbol \theta, \sigma^2) &= \mathcal{N}(\mathbf y; \mathbf X \boldsymbol \theta, \sigma^2 \mathbf I) \\ &= \sqrt{\frac{1}{(2\pi)^n |\sigma^2 \mathbf I|}} \cdot \exp\left[ -\frac{1}{2} (\mathbf y-\mathbf X \boldsymbol \theta)^\top (\sigma^2 \mathbf I)^{-1} (\mathbf y-\mathbf X \boldsymbol \theta) \right]. \end{split} $$

Используя свойство умножения определителя матрицы на скаляр $| \sigma^2 \mathbf I | = (\sigma^2)^n |\mathbf I|, $ зададим логарифмическую функцию правдоподобия,

$$ \begin{split} \ell (\boldsymbol \theta,\sigma^2) =& \; \log P(\mathbf y \mid \boldsymbol \theta, \sigma^2) \\ =& -\frac{n}{2} \log(2\pi)-\frac{n}{2} \log (\sigma^2)-\frac{1}{2} \log (|\mathbf I|) \\ & -\frac{1}{2 \sigma^2} (\mathbf y-\mathbf X \boldsymbol \theta)^\top \mathbf I^{-1} (\mathbf y-\mathbf X \boldsymbol \theta) \\ =& -\frac{n}{2} \log(2\pi)-\frac{n}{2} \log(\sigma^2)-\frac{1}{2} \log(|\mathbf I|) \\ & -\frac{1}{2 \sigma^2} \left(\mathbf y^\top \mathbf I^{-1} \mathbf y-2 \boldsymbol \theta^\top \mathbf X^\top \mathbf I^{-1} \mathbf y + \boldsymbol \theta^\top \mathbf X^\top \mathbf I^{-1} \mathbf X \boldsymbol \theta \right). \end{split} $$

Найдем производную $\ell$ относительно $\boldsymbol \theta,$

$$ \begin{split} \frac{ \partial \ell (\boldsymbol \theta, \sigma^2) }{ \partial \boldsymbol \theta } &= \frac{\partial}{\partial \boldsymbol \theta} \left( -\frac{1}{2 \sigma^2} \left( \mathbf y^\top \mathbf I^{-1} \mathbf y-2 \boldsymbol \theta^\top \mathbf X^\top \mathbf I^{-1} \mathbf y + \boldsymbol \theta^\top \mathbf X^\top \mathbf I^{-1} \mathbf X \boldsymbol \theta\right) \right) \\ &= \frac{1}{2 \sigma^2} \frac{\partial}{\partial \boldsymbol \theta} \left( 2 \boldsymbol \theta^\top \mathbf X^\top \mathbf I^{-1} \mathbf y-\boldsymbol \theta^\top \mathbf X^\top \mathbf I^{-1} \mathbf X \boldsymbol \theta \right) \\ &= \frac{1}{2 \sigma^2} \left( 2 \mathbf X^\top \mathbf I^{-1} \mathbf y-2 \mathbf X^\top \mathbf I^{-1} \mathbf X \boldsymbol \theta \right) \\ &= \frac{1}{\sigma^2} \left( \mathbf X^\top \mathbf I^{-1} \mathbf y-\mathbf X^\top \mathbf I^{-1} \mathbf X \boldsymbol \theta \right). \end{split} $$

Приравняем производную к нулю и выразим $\boldsymbol \theta,$

$$ \begin{split} \frac{\partial \ell (\hat{\boldsymbol \theta},\sigma^2)}{\partial \boldsymbol \theta} &= 0 \\ 0 &= \frac{1}{\sigma^2} \left( \mathbf X^\top \mathbf I^{-1} \mathbf y-\mathbf X^\top \mathbf I^{-1} \mathbf X \hat{\boldsymbol \theta} \right) \\ 0 &= \mathbf X^\top \mathbf I^{-1} \mathbf y-\mathbf X^\top \mathbf I^{-1} \mathbf X \hat{\boldsymbol \theta} \\ \mathbf X^\top \mathbf I^{-1} \mathbf X \hat{\boldsymbol \theta} &= \mathbf X^\top \mathbf I^{-1} \mathbf y \\ \hat{\boldsymbol \theta} &= \left( \mathbf X^\top \mathbf I^{-1} \mathbf X \right)^{-1} \mathbf X^\top \mathbf I^{-1} \mathbf y \\ \hat{\boldsymbol \theta} &= \left( \mathbf X^\top \mathbf X \right)^{-1} \mathbf X^\top \mathbf y. \end{split} $$

Найдем производную функции $\ell$ со значением $\hat{\boldsymbol \theta}$ относительно $\sigma^2,$

$$ \begin{split} \frac{\partial \ell (\hat{\boldsymbol \theta},\sigma^2)}{\partial \sigma^2} &= \frac{\partial}{\partial \sigma^2} \left( -\frac{n}{2} \log (\sigma^2)-\frac{1}{2 \sigma^2} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top \mathbf I^{-1} (\mathbf y-\mathbf X\hat{\boldsymbol \theta}) \right) \\ &= -\frac{n}{2} \frac{1}{\sigma^2} + \frac{1}{2 (\sigma^2)^2} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top \mathbf I^{-1} (\mathbf y-\mathbf X\hat{\boldsymbol \theta}) \\ &= — \frac{n}{2 \sigma^2} + \frac{1}{2 (\sigma^2)^2} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top \mathbf I^{-1} (\mathbf y-\mathbf X\hat{\boldsymbol \theta}). \end{split} $$

Приравняем производную к нулю и выразим $\sigma^2,$

$$ \begin{split} \frac{\partial \ell (\hat{\boldsymbol \theta},\hat{\sigma}^2)}{\partial \sigma^2} &= 0 \\ 0 &= -\frac{n}{2 \hat{\sigma}^2} + \frac{1}{2 (\hat{\sigma}^2)^2} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top \mathbf I^{-1} (\mathbf y-\mathbf X\hat{\boldsymbol \theta}) \\ \frac{n}{2 \hat{\sigma}^2} &= \frac{1}{2 (\hat{\sigma}^2)^2} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top \mathbf I^{-1} (y-X\hat{\boldsymbol \theta}) \\ \frac{2 (\hat{\sigma}^2)^2}{n} \cdot \frac{n}{2 \hat{\sigma}^2} &= \frac{2 (\hat{\sigma}^2)^2}{n} \cdot \frac{1}{2 (\hat{\sigma}^2)^2} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top \mathbf I^{-1} (\mathbf y-\mathbf X\hat{\boldsymbol \theta}) \\ \hat{\sigma}^2 &= \frac{1}{n} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top \mathbf I^{-1} (\mathbf y-\mathbf X\hat{\boldsymbol \theta}) \\ \hat{\sigma}^2 &= \frac{1}{n} (\mathbf y-\mathbf X\hat{\boldsymbol \theta})^\top (\mathbf y-\mathbf X\hat{\boldsymbol \theta}). \end{split} $$

Масштабирование признаков

Рассмотрим линейную регрессию с двумя независимыми переменными с разным масштабом. Например, пусть один признак находится в диапазоне от 0 до 1000, а второй — от 0 до 1. Для того чтобы настроиться на такие признаки, модель выберет небольшой первый и большой второй коэффициенты.

Проблема подобной модели в том, что на этапе обучения скорость оптимизации таких коэффициентов (весов, параметров) не будет одинаковой.

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

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

Пропишем функцию потерь и найдем ее частные производные.

Попробуем найти оптимальные веса с помощью алгоритма градиентного спуска.

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

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

Посмотрим, что будет, если привести признаки к одному масштабу. Вначале найдем минимум методом градиентного спуска.

Выведем результат на графике.

Посмотрим на обновления весов.

Мы видим только оранжевую линию, потому что одна кривая перекрывает другую.

Полиномиальная регрессия

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

Попробуем «наложить» на данные модель линейной регрессии.

Посмотрим на получившиеся коэффициенты.

Выведем результат на графике и оценим качество через RMSE и $R^2$.

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

$$ \hat{y} = \theta_0 + \theta_1x + \theta_2x^2. $$

Обратите внимание, что в квадрат возводится не коэффициент $\theta$, а признак $x$. Поэтому все, что нам нужно, это создать еще один признак, равный квадрату исходного признака. Для этого в sklearn есть класс PolynomialFeatures.

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

Второй столбец — квадрат первого.

Построим модель.

Два важных замечания,

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