Метод наименьших квадратов

Все курсы > Линейная алгебра > Занятие 7

Посмотрим на задачу регрессии с точки зрения линейной алгебры.

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

Общее решение

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

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

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

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

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

$$ \begin{bmatrix} \vdots & \vdots & \vdots & \vdots \\ a_1 & a_2 & \dots & a_k \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_k \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix} $$

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

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

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

$$ x_1 \mathbf a_1 + x_2 \mathbf a_2 + x_3 \mathbf a_3 + … + x_k \mathbf a_k \not= \mathbf b $$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$$ A \mathbf x^*-\mathbf b \in null(A^{T}) $$

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

$$ A^T (A \mathbf x^*-\mathbf b) = \mathbf 0 $$

$$ A^TA \mathbf x^*-A^T \mathbf b = \mathbf 0 $$

$$ A^TA \mathbf x^* = A^T \mathbf b $$

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

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

$$ X^TX\theta = X^Ty $$

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

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

Пример

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

$$ \left\{\begin{matrix} 2x-y=2 \\ x+2y=1 \\ x + y = 4 \end{matrix}\right. $$

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

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

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

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

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

Решим систему (а такая система уже решается) методом Гаусса.

Вектор $ \begin{bmatrix} 1.43 \\ 0.43 \end{bmatrix} $ и будет решением. Обозначим его на графике с помощью точки.

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

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

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

Обратимость матрицы $A^TA$

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

$$ \mathbf x = (A^TA)^{-1}A^T \mathbf b $$

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

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

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

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

Найдем $A_1^T A_1$ и $A_2^T A_2$.

$$ A_1^T 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} $$

$$ A_2^T 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} $$

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

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

$$ \mathbf x^T A^TA \mathbf x = \mathbf x^T \mathbf 0 $$

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

$$ (A \mathbf x)^T (A \mathbf x) = 0 $$

$$ || A \mathbf x ||^2 = 0 $$

Значит,

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

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

$$ null(A^TA) = null(A) $$

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

$$ dim(null(A^TA)) = dim(null(A)) \implies rank(A^TA) = rank(A) $$

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

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

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

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

$$ \mathbf p = A \mathbf x^* = A (A^TA)^{-1} A^T $$

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

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

$$ P^T = P $$

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

$$ A (A^TA)^{-1} \cancel{A^T A} \cancel{(A^TA)^{-1}} A^T = A (A^TA)^{-1} A^T $$

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

$\mathbf b \in col(A)$

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

$$ AA^{-1} (A^T)^{-1}A^T = I $$

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

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

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

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

$\mathbf b \perp col(A)$

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

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

Подведем итог

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

Вновь рассмотрим смену базиса.