Линейная регрессия

Матрица объекты-признаки

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

В общем случае, линейная регрессия - это просто поиск наилучшей аппроксимации среди моделей вида

\[ \check y (x) = g_\mathbf{w}(x) = \sum^m_{k=1} \text{w}_k \cdot f_k (x) \]

где \(\{f_k(x)\}^m_{k=1}\) - какой-то набор из \(m\) функций (крайне желательно независимых).

В идеале (для минимизации любого loss-a) , мы должны получить полное совпадение на лэйблах (парах \((x_i, y_i))\):

\[ \forall x_i \; \rightarrow \; \check y (x_i) = \sum^m_{k=1} \text{w}_k \cdot f_k (x_i) \stackrel{\text{want}}{=} y_i \]

Переписывая в матричном виде:

\[\begin{split} \left(\begin{array}{c} \check y_1 \\ \check y_2 \\ \ldots \\ \check y_N \end{array}\right)= \left(\begin{array}{cccc} f_1\left(x_1\right) & f_2\left(x_1\right) & \ldots & f_m\left(x_1\right) \\ f_1\left(x_2\right) & f_2\left(x_2\right) & \ldots & f_m\left(x_2\right) \\ \ldots & \ldots & \ldots & \ldots \\ f_1\left(x_N\right) & f_2\left(x_N\right) & \ldots & f_m\left(x_N\right) \end{array}\right)\left(\begin{array}{c} \text{w}_1 \\ \text{w}_2 \\ \ldots \\ \text{w}_N \end{array}\right)\stackrel{\text{want}}{=} \left(\begin{array}{c} y_1 \\ y_2 \\ \ldots \\ y_N \end{array}\right) \end{split}\]

Или короче:

\[ \check Y = X \cdot \mathbf{w} \stackrel{\text{want}}{=} Y \]

Определение. Назовём строки матрицы \(X\) - объектами, а столбцы - признаками. Легко принять, что тогда матрицу \(X\) назовём матрицей объекты-признаки.

То есть множество всех признаков для \(x_i\) - \(f_1(x_i), \ldots, f_m(x_i)\) - это объект.

Далее будем часто работать с компонентами этой матрицы (тензорные обозначения) \(X = (x_{ij}) = (f_j(x_i))\), где \(i\) - номер строки, \(j\) - номер столбца.

Примечание. Символ \(\stackrel{\text{want}}{=}\) здесь не означает строго равенства (его может не быть какие бы ни брали \(\mathbf{w}\)). Понимать его стоит с точки срезния минимизации конкретного loss-а (в частности, далее получим аналитический результат для MSE).

Точное решение линейной регрессии для MSE

Начнём с анекдота.

Внимание! Рассуждения ниже не надо воспринимать серьезно.

Хотим найти такой \(\mathbf{w}^*\), чтобы \(X \cdot \mathbf{w}^* = Y\). Домножим слева на \(X^\top\), чтобы сделать матрицу в левой части уравнения квадратной:

\[ X^\top \cdot X \cdot \mathbf{w}^* = X^\top Y \]

Надеемся\верим\предполагаем, что \((X^\top \cdot X)\) обратима, а тогда можем домножить на обратную слева

\[ \mathbf{w}^* = \left( X^\top \cdot X \right)^{-1} \cdot X^\top \cdot Y \]

Вуаля, нашли искомый столбец параметров и рЕшИлИ задачу регрессии в общем случае для лЮбОгО loss-а и лЮбОгО набора функций \(\{f_k(x)\}^m_{k=1}\), т.к. равенство достигнуто 😏😏😏

Конец внимания! Рассуждения ниже нужно воспринимать серьезно.

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

Теорема. В случае линейной регрессии c MSE функции ошибок, для заданного набора \(\{f_k(x)\}^m_{k=1}\), наилучшим набором параметров будет вектор:

\[ \mathbf{w}^* = \left( X^\top \cdot X \right)^{-1} \cdot X^\top \cdot Y \]

Доказательство:

Распишем MSE (под конец приведём всё к тензорным обозначениям для удобства дифференциирования):

\[\begin{split} \begin{aligned} \mathcal{L} & =\frac{1}{N} \sum_{i=1}^N\left(\check{y}_i-y_i\right)^2=\left\{\check{y}_i=x_{i j} \text{w}_i\right\}=\frac{1}{N}(X \cdot \mathbf{w}-Y)^{\top}(X \cdot \mathbf{w}-Y)= \\ & =\frac{1}{N}\left(\mathbf{w}^{\top} X^{\top} X \mathbf{w}-Y^{\top} X \cdot \mathbf{w}-\mathbf{w}^{\top} \cdot X^{\top} \cdot Y+Y^{\top} Y\right)= \\ & =\frac{1}{N}\left(\text{w}_j x_{i j} x_{i s} \text{w}_s-y_i \cdot x_{i j} \cdot \text{w}_j-x_{i j} \text{w}_j \cdot y_i+y_i y_i\right) \end{aligned} \end{split}\]

Для нахождения минимума функции многих переменных надо найти нули производной:

\[\begin{split} \begin{aligned} \frac{\partial \mathcal{L}}{\partial \text{w}_k} & =\left\{\begin{array}{l} \frac{\partial \text{w}_i}{\partial \text{w}_k}=\delta_{i k} \\ A_{k s} \delta_{s m}=A_{k m} \end{array}\right\}=\frac{1}{N}\left(\delta_{j k} x_{i j} x_{i s} \text{w}_s+\text{w}_j x_{i j} x_{i s} \delta_{s k}-y_i x_{i j} \delta_{j k}-x_{i j} \delta_{j k} \cdot y_i+0\right)= \\ & =\frac{2}{N}\left(x_{i k} x_{i s} \text{w}_s-x_{i k} y_i\right)=0 \end{aligned} \end{split}\]

А тогда получаем уже знакомое выражение (в этот раз переход к нему абсолютно «законный»):

\[ X^{\top} X \cdot \mathbf w-X^{\top} Y=0 \]

Опять-таки вынуждены предположить обратимость \((X^\top \cdot X)\), итого,

\[ \mathbf w ^ *=\left(X^{\top} X\right)^{-1} X^{\top} Y \]

Определение. Для произвольной матрицы \(X\) (даже прямоугольной) назовём псевдообратной матрицу \(\left(X^{\top} X\right)^{-1} X^{\top}\)

• Пример с независимыми признаками

Пусть имеем набор \(\left\{\left(\mathbf{x}_i, y_i\right)\right\}_{i=1}^N\), где \(\mathbf{x}_i \in \mathbb{R}^3\).

Пусть этот набор определён функцией:

\[ y \left(x^1, x^2, x^3 \right) = \mathbf{3} \cdot x^1 - \mathbf{2} \cdot x^2 + \mathbf{5} \cdot x^3 - \mathbf{4} \]

(Сделаем вид, что мы не знаем эти коэффициенты и как раз-таки их надо восстановить) Здесь обозначение \(x^k\) - не степень, а номер оси.

То есть хотим зафитить какую-то функцию трёх переменных. Вполне себе законно выбрать функции \(f_1, f_2, f_3\) как проекции вектора \(\mathbf{x}\) на три оси из \(\mathbb{R}^3\) - это будут наши признаки. Тем не менее даже видим, что их не хватит, и надо добавить \(f_4(x) = 1\) как четвёртый признак (intercept или bias).

Итого, матрица объекты-признаки будет выглядеть как

\[\begin{split} X=\left(\begin{array}{cccc} x^1_1 & x^2_1 & x^3_1 & 1 \\ x^1_2 & x^2_2 & x^3_2 & 1 \\ \ldots & \ldots & \ldots & \ldots \\ x^1_N & x^2_N & x^3_N & 1 \end{array}\right) \end{split}\]

Ожидаем, что сможем восстановить все коээфициенты из полученной ранее формулы:

\[ \mathbf w ^ *=\left(X^{\top} X\right)^{-1} X^{\top} Y \]

import numpy as np

# Генерируем выборку для исследования

x1 = np.random.random(100).reshape((-1, 1))
x2 = np.random.random(100).reshape((-1, 1))
x3 = np.random.random(100).reshape((-1, 1))

y_label = 3*x1 - 2*x2 + 5*x3 - 4 # Наша цель дальше - восстановить эти коэффициенты

# Создаём матрицу объекты-признаки
X = np.concatenate([x1, x2, x3, np.ones(x1.shape)], axis = 1)

# Считаем псевдообратную и сразу же умножаем на y_label для получения коэффициентов
w = np.linalg.pinv(X) @ y_label # Лучше использовать pinv, чем самому писать выражение
w                               # для почти вырожденных матриц сработает лучше

array([[ 3.],
       [-2.],
       [ 5.],
       [-4.]])

• Пример с многочленами

Пусть имеем набор \(\left\{\left(x_i, y_i\right)\right\}_{i=1}^N\), где \(x_i \in \mathbb{R}\).

Пусть этот набор определён функцией:

\[ y \left(x \right) = - \mathbf{4} + \mathbf{3} \cdot x^1 - \mathbf{2} \cdot x^2 + \mathbf{5} \cdot x^3 \]

(Сделаем вид, что мы не знаем эти коэффициенты и как раз-таки их надо восстановить) А здесь обозначение \(x^k\) - степень, а не номер оси.

Мы ожидаем, что зафитим эту зависимость каким-то многочленом, но заранее его степень не знаем. Поэтому нашими признаками будут многочлены до степени \(m\) включительно. Опять-таки не забываем добавить intercept.

Матрица объекты-признаки, таким образом, будет:

\[\begin{split} X=\left(\begin{array}{} x^1_1 & x^2_1 & \ldots & x^m_1 & 1 \\ x^1_2 & x^2_2 & \ldots & x^m_2 & 1 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ x^1_N & x^2_N & \ldots & x^m_N & 1 \end{array}\right) \end{split}\]

Ожидаем, что сможем восстановить все коээфициенты из полученной ранее формулы:

\[ \mathbf w ^ *=\left(X^{\top} X\right)^{-1} X^{\top} Y \]

А параметр \(m\) (такие «внешние» парметры называются гиперпараметрами модели, т.к. непосредственно их мы не фитим) будем определять по величине ошибки MSE (то, где меньше, там и лучше, казалось бы 😏, - построим график MSE от \(m\)).

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

import numpy as np

# Генерируем выборку для исследования

x = 2*np.random.random(200).reshape((-1, 1)) - 1
noise = 0.1*np.random.normal(size = x.shape)
y_label = -4 + 3*x - 2*x**2 + 5*x**3 + noise

# Напишем функцию для создания матрицы объекты-признаки
def create_X(x, m):
    X = np.ones(x.shape) # Сразу же добавим intercept
    for i in range(1, m + 1):
        x.reshape(X.shape[0], 1)
        X = np.concatenate([X, x**i], axis = 1)
    return X


# Напишем функцию для определения весов при заданном m
def fit(x, y_label, m):
    X = create_X(x, m)
    w = np.linalg.pinv(X) @ y_label
    return w

# Функция вычисления y_predict для любого x модели с весами w
def predict(x, w):
    return create_X(x, np.size(w) - 1) @ w

# Функция ошибки MSE для модели (которая полностью определяется предсказаниями модели и y_label)
def MSE(y, y_label):
    return np.sum((y - y_label)**2) / np.size(y)


fit(x, y_label, 20)

array([[-3.99537082e+00],
       [ 2.52945011e+00],
       [-2.87013798e+00],
       [ 2.37031935e+01],
       [ 3.18064489e+01],
       [-2.92517087e+02],
       [-3.62608739e+02],
       [ 2.27117748e+03],
       [ 1.89589890e+03],
       [-9.75462693e+03],
       [-5.34239920e+03],
       [ 2.46783719e+04],
       [ 8.58441046e+03],
       [-3.76788116e+04],
       [-7.70537344e+03],
       [ 3.41192311e+04],
       [ 3.31380385e+03],
       [-1.68682118e+04],
       [-2.15665236e+02],
       [ 3.50712515e+03],
       [-1.99141226e+02]])

# Построим MSE от m
for m in range(100):
    w = fit(x, y_label, m)
    y_predict = predict(x, w)
    print(f"for {m} MSE is {MSE(y_predict, y_label)}")

for 0 MSE is 15.071678459840461
for 1 MSE is 0.9811942013228904
for 2 MSE is 0.572496049338976
for 3 MSE is 0.010089320129317025
for 4 MSE is 0.009950532084079434
for 5 MSE is 0.009950438485090852
for 6 MSE is 0.009949288146964675
for 7 MSE is 0.009907487326393086
for 8 MSE is 0.009906623794384319
for 9 MSE is 0.00982684225541596
for 10 MSE is 0.00982628883908817
for 11 MSE is 0.00977639336886851
for 12 MSE is 0.009679743693416055
for 13 MSE is 0.009677648940186575
for 14 MSE is 0.009671455905351034
for 15 MSE is 0.009664807043356672
for 16 MSE is 0.00965398933027318
for 17 MSE is 0.009616458251452582
for 18 MSE is 0.00958586457141497
for 19 MSE is 0.00955119782228293
for 20 MSE is 0.009551174362673207
for 21 MSE is 0.00954717509357125
for 22 MSE is 0.00952359082311598
for 23 MSE is 0.009506641413195625
for 24 MSE is 0.009493684488390301
for 25 MSE is 0.009349161417467618
for 26 MSE is 0.009349065952548707
for 27 MSE is 0.009331775389885227
for 28 MSE is 0.009180371521815875
for 29 MSE is 0.009169115778730497
for 30 MSE is 0.009046648689574877
for 31 MSE is 0.009035476377882871
for 32 MSE is 0.008952961538979131
for 33 MSE is 0.008926551312570644
for 34 MSE is 0.008922990653483166
for 35 MSE is 0.00892064625801604
for 36 MSE is 0.008910131303252314
for 37 MSE is 0.008905600255871775
for 38 MSE is 0.008905742606756245
for 39 MSE is 0.008929104314532144
for 40 MSE is 0.00889978329197983
for 41 MSE is 0.008890893090345374
for 42 MSE is 0.008773839789685616
for 43 MSE is 0.008835703093760482
for 44 MSE is 0.008748268642354005
for 45 MSE is 0.008741537550430872
for 46 MSE is 0.008761499651783017
for 47 MSE is 0.008749519273825713
for 48 MSE is 0.008729568169513426
for 49 MSE is 0.00875494530855859
for 50 MSE is 0.008795783611986172
for 51 MSE is 0.008724639748509502
for 52 MSE is 0.008738708143867784
for 53 MSE is 0.008627729660765611
for 54 MSE is 0.008618597874757786
for 55 MSE is 0.008666925149220288
for 56 MSE is 0.008565977872597537
for 57 MSE is 0.00858472422223916
for 58 MSE is 0.008568268507908735
for 59 MSE is 0.008536388392472973
for 60 MSE is 0.008634546358496442
for 61 MSE is 0.00880435411742616
for 62 MSE is 0.008555952046914913
for 63 MSE is 0.008582576851050072
for 64 MSE is 0.008498726264527015
for 65 MSE is 0.008501797219667472
for 66 MSE is 0.00849994713611994
for 67 MSE is 0.008496975762834572
for 68 MSE is 0.008495357160327767
for 69 MSE is 0.008493619489219714
for 70 MSE is 0.008546805782612278
for 71 MSE is 0.008528966871771439
for 72 MSE is 0.008522890345794688
for 73 MSE is 0.008474639923292454
for 74 MSE is 0.008469993960783835
for 75 MSE is 0.008529326358066415
for 76 MSE is 0.008460488301038824
for 77 MSE is 0.008428647168531433
for 78 MSE is 0.008442609881198429
for 79 MSE is 0.008397549905914754
for 80 MSE is 0.008412822432991227
for 81 MSE is 0.008405035896233875
for 82 MSE is 0.008416672944564411
for 83 MSE is 0.008380465673201916
for 84 MSE is 0.008352994898924816
for 85 MSE is 0.008331062039667953
for 86 MSE is 0.008341615166082908
for 87 MSE is 0.008843693919905735
for 88 MSE is 0.008367421504380282
for 89 MSE is 0.008486968061458332
for 90 MSE is 0.008310989306543367
for 91 MSE is 0.00832032085742849
for 92 MSE is 0.008540373770987564
for 93 MSE is 0.008348175488294583
for 94 MSE is 0.008344238425524154
for 95 MSE is 0.008335237971722023
for 96 MSE is 0.008321390948480785
for 97 MSE is 0.008324486879171013
for 98 MSE is 0.008326665339188901
for 99 MSE is 0.008328174678272966

Ошибка уменьшается с ростом m, но как же так, у нас же не многочлен большой степени?

Да просто показатели коэффициенты при больших степенях будут маленькие:

fit(x, y_label, 5) # Значит ли это, что можем ставить m сколь угодно большим? Да и не то чтобы они совсем уже маленькие

array([[-3.9937026 ],
       [ 2.98968055],
       [-2.10265803],
       [ 5.01935966],
       [ 0.1536202 ],
       [-0.00795813]])

Переобучение

from bokeh.io import output_notebook
from bokeh.resources import INLINE
output_notebook(INLINE)

Loading BokehJS ...