w0lkolak 6 июл в 21:47

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

Простой

7 мин

13K

Python*Математика*Машинное обучение*

Туториал

Я прохожу онлайн курс по ML, а здесь я пишу заметки, в которых, как мне кажется, я нуждался неделю назад.

class GradientDescentMse:
    """
    Базовый класс для реализации градиентного спуска в задаче линейной МНК регрессии 
    """

Прочтя такое задание, я почувствовал как мой мир рушится. Всё, что я знал о ML, а это два с половиной урока кажутся пройденными зря - ~~МНК это не метод решения для линейной регрессии, а что то другое..~~. И я начал разбираться заново. Какую статью не открою - везде простыня с формулами. В итоге самые короткие и простые я прочёл. Но можно ещё проще. Перейти сразу к МНК

База

Линейная регрессия это сведение зависимости y(X) к линейному уравнению вида

$y = k_1 X_1 + k_2 X_2 + \cdots + k_m X_m + C$

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

Во-первых какова геометрия, описываемая уравнением?

Если признак X в количестве одной штуки (таблица имеет один столбец входных данных), то это прямая, делящая двумерную плоскость на то, что выше и то, что ниже.
Если столбцов в исходных данных два, то это плоскость, делящая 3х-мерное пространство на то, что выше и то что ниже.
Короче, это m мерная геометрическая штука, рассекающая m+1 мерное пространство, ибо это пространство ещё имеет ось искомого таргета y.

Сразу определюсь, что m, как и в формуле - это количество признаков (так же это количество столбцов в таблице с данными), n - количество строк или количество точек в геометрическом представлении. Почему речь о прямой, плоскости и тд? Потому что линейные уравнения не описывают парабол, гипербол, сфер или цилиндров в пространстве. Они описывают m-мерные гиперплоскости(это общее название, подходящее даже для прямой). Можно попытаться представить 3х мерное пространство, секущее 4х мерное как временной слепок местоположения объектов, однако представить 4х-мерное пространство секущее 5-мерное я не могу, но мне и не нужно, так как принципиального отличия от работы с 2d плоскостью нет.

Что важно для понимания описанной геометрии:

Если C=0, то секущая проходит через начало координат. Если C $\neq$ 0, то гиперплоскость смещена по оси таргета от начала системы координат. C - это величина, на которую гиперплоскость скользит от нуля по оси таргета.
Если n=m, то решение СЛАУ не даст свободного коэффициента C.
Если n=m+1, то для решения СЛАУ мне придётся добавить свободный коэффициент C. Почему?

Пример линейной зависимости таргета y от признаков ,

Решением для n=m и n=m+1 будет нахождение идеально подходящих значений коэффициентов. Идеальных в том смысле, что они опишут изначальные данные с нулевой ошибкой. Однако такие коэффициенты вряд-ли останутся идеальными или хотя бы приемлемыми для любых других точек: ожидать, что объекты реального мира подчинены точной линейной зависимости было бы наивно. Способность коэффициентов предсказывать результат по данным, не использованным для их нахождения, называется обобщающей способностью. В случае, с n $\approx$ m, она не будет низкой.

Почему n=m+1 недостаточно, если СЛАУ решается?

Даже если случайности не существует, существует признак, который я не учёл. А те признаки, что учёл, все равно с ошибкой. Это результат или неверно работающего измерительного прибора или пьяные/уставшие/тапающие_хомяка человеки вбивали данные в БД и время от времени косячили. Это если данные снимались и вбивались специально обученной обезьяной в рамках рабочего процесса. Если данные - результат опроса, то никто ничего не помнит, а то, что помнят - все равно могут исказить, так как стремно признаться, что ты полтора метра и полтора центнера одновременно.

Побороть влияние искажений можно увеличив количество данных. Увеличивать m можно только пока есть какая-то корреляция между новым признаком и таргетом (И пока сложность модели соизмерима со сложностью задачи). Увеличивать n можно почти безгранично и чем больше, тем лучше. Больше данных ~~Богу Данных~~ Data Scientistу. Он разберётся кого в X_train, кого в корзину.

Итак стандартная ситуация это n>>m, но можно же решить все СЛАУ, и найти среднее.

Можно, конечно, решить СЛАУ для m первых уравнений, потом сместиться на step>=1 решать снова и снова, а потом найти средне арифметические значения всех коэффициентов. Если данных много, то это даже может сработать. Я провёл тест с теми данными что рассматриваются в течении всей статьи: и получил
[0.7, 0.8, -0.336] что случайными данными не назовешь, но и от минимизации MSE это сильно отличается. Лучшая MAE тоже ближе к идеально возможной MSE - [0.46, .0033, -0.147]. Придуманный мной на ходу велосипед едет, но не очень. Тем более на очень скромном наборе данных, что я рассматриваю.

Важно для ML в целом:

В ML ВСЕГДА или почти ВСЕГДА не хотят решать исходную СЛАУ: где
В ML хотят выбрать функцию ошибки и минимизировать её пользуясь двумя фактами:
1) Производная любой функции это быстрота изменения значения этой функции при наращивании независимой переменной.
2) Я заранее знаю, чему функция ошибки должна быть равна - в идеале нулю. Я ищу коэффициенты там, где функция или ноль, или минимальна. Это то место, где производная по слева отрицательна (увеличение приводит к снижению функционала ошибки), а справа положительна - увеличение приводит к росту ошибки. И равна производная чему-то между отрицательными и положительными значениями - "0".
Часто предстоит брать производные по стандартным функциям ошибок и как-то совать в них X_train, y_train и как-то усредняться.
В ML любят квадратичную ошибку.

За что?

Она сильнее штрафует за абсолютную ошибку q, чем за 2 ошибки величиной q/2. А ещё у неё точно есть производная в нуле. А у модуля от ошибки производную в ноле не посчитаешь - 0/|0| стремится к -1 с одной стороны и 1 с другой. MAE может вызывать как минимум внутренний дискомфорт. А ошибку без квадратов и модулей вообще рассматривать нет смысла, так как ошибись ты на одних данных в среднем на +9000, на других на -9000 получишь нулевую ошибку и нулевую предсказательную силу.

МНК (метод наименьших квадратов) aka OLS (ordinary least squares)

Если я возьму квадратичную ошибку (мог бы и среднюю, так не изменит расчетов, но МНК этого всё таки не подразумевает) и минимизирую её, решив новую СЛАУ из производных по всем и , то это и будет применение МНК. Приступаю, выдумав пример из 2 признаков и свободного коэффициента .

$\sum_{i=1}^{n} (y_i - k_1 x_{1i} + k_2 x_{2i} + C)^2 \rightarrow min.$

*Верное уточнение от Dimaush : а $\rightarrow min.$

-ы у меня есть в исходных данных, а - нет. Я буду дифференцировать по и , а у меня будет пока что неким коэффициентом, МНК я запишу английской аббревиатурой OLS:

$\frac{\partial \text{OLS}}{\partial C} = -2 \sum_{i=1}^{n} (y_i - C - k_1 x_{1i} - k_2 x_{2i})$ $\frac{\partial \text{OLS}}{\partial k_1} = -2 \sum_{i=1}^{n} x_{1i} (y_i - C - k_1 x_{1i} - k_2 x_{2i})$ $\frac{\partial \text{OLS}}{\partial k_2} = -2 \sum_{i=1}^{n} x_{2i} (y_i - C - k_1 x_{1i} - k_2 x_{2i})$

Поскольку я ищу экстремум, точку покоя (не путать с точкой перелома), как было сказано в базе в искомой точке производная по равна нулю:

$\begin{align*} \frac{\partial \text{OLS}}{\partial C} &= 0, & \frac{\partial \text{OLS}}{\partial k_1} &= 0, & \frac{\partial \text{OLS}}{\partial k_2} &= 0 \end{align*}$

Получается, есть 3 выражения с тремя неизвестными, приравненные к нолю. Это искомая СЛАУ. Именно СЛАУ, так как хоть в первой хоть в сотой степени - просто число. Немного отредачу СЛАУ, а именно "развалю" $\sum$ и выкину -2 за ненадобностью. Тут нет ничего сложного, я просто раскрыл скобки и поставил SUM каждому слагаемому по-отдельности. В аналитической части ничего сложнее уже не будет.

$\sum_{i=1}^{n} x_{1i}y_{i} = C \sum_{i=1}^{n} x_{1i} - k_{1} \sum_{i=1}^{n} x_{1i}^{2} - k_{2} \sum_{i=1}^{n} x_{1i}x_{2i}$

Аналогично для k2 и C

$\sum_{i=1}^{n} x_{2i}y_{i} = C \sum_{i=1}^{n} x_{2i} - k_{1} \sum_{i=1}^{n} x_{1i}x_{2i} - k_{2} \sum_{i=1}^{n} x_{2i}^{2}$ $\sum_{i=1}^{n} y_{i} = nC + k_{1} \sum_{i=1}^{n} x_{1i} + k_{2} \sum_{i=1}^{n} x_{2i}$

В новой СЛАУ всё - числа, кроме коэффициентов. Представляю, что это неизвестные и просто решаю. Точнее там будут числа, если у меня появятся какие-то исходные данные:

X1	X2	y
1	2	1.5
2	3	1.8
3	1	3.2
4	5	3.6
5	4	5.1

Подставив, перемножив и сложив, я получил:

Посмотреть решение в онлайн калькуляторе СЛАУ

MSE = 0.0971

Матричная форма

Ранее я использовал обозначения и для коэффициентов линейной регрессии, но если обозначить коэффициенты как , а также обозначить как то можно будет увидеть смысл в такой формуле OLS, с вектором-столбцом:

$\sum_{i=1}^{n} \left( y_i - f(x_i, B) \right)^2\rightarrow min$

Однако, подробнее формула несколько сложнее:

$y^T y - 2B^T X^T y + B^T X^T X B\rightarrow min$

это число, равное $\sum_{i=1}^{n}y_i^2$ , с остальным не так очевидно, однако последовательное применение произведений и транспонирований позволяет по крайней мере понять, что два других слагаемых тоже являются числами. Предположу, что формула верна и попробую матрично продифференцировать по вектору ( исчезли, также исчезли слагаемые без ):

$X^\top y = X^\top X B$

Осталось выразить ~~B = X.T.dot(X) / X.T.dot(y)~~, но, с матрицами так нельзя. Вместо этого матрица, обратная делимому матрично умножается на частное: $B = \left( X^T X \right)^{-1} X^\top y$

Если матрица содержит строки, получаемые из других строк масштабированием, обратная матрица не найдется, для этого есть псведообратные матрицы и метод pinv в пакете np.linalg , например.

Почему обратная матрица может не найтись?

Матрица это всего-лишь табличная форма записи СЛАУ(Почти всегда). СЛАУ с 3 неизвестными и с 3 уравнениями, где одно уравнение получается масштабированием другого тоже не решаема, а $A \cdot A^{-1} = \mathbf{I}$ это тоже по сути уравнение и ища обратную матрицу, я его именно решаю. И строки и столбцы равнозначимы в этом решении.

Посмотреть как это делается в онлайн калькуляторе
(нужно указать foo в функцию и k_1 в аргумент)

Посмотреть велосипед в python с составлением и решением СЛАУ

import numpy as np
from scipy.linalg import solve


y = np.array([[1.5, 1.8, 3.2, 3.6, 5.1]]) # Всё лежит на боку, мне так удобно
X = np.array([[1,2,3,4,5], [2,3,1,5,4]])  # Если хочешь привычную форму - транспонируй X и y

SUM_x1 = X[0].sum()
SUM_x2 = X[1].sum()
SUM_y = y.sum()

SUM_x1_x1 = (X[0]**2).sum()
SUM_x2_x2 = (X[1]**2).sum()
SUM_x1_x2 = (X[0] * X[1]).sum()

SUM_x1_y = (X[0] * y).sum()
SUM_x2_y = (X[1] * y).sum()
#Дальше СЛАУ
"""
SUM_y = len(y[0])*C + k1*SUM_x1 + k2*SUM_x2
SUM_x1_y  = C*SUM_x1 - k1*SUM_x1_x1 - k2*SUM_x1_x2
SUM_x2_y = C*SUM_x2  - k1*SUM_x1_x2 - k2*SUM_x2_x2
"""

X_ = np.array([[len(y[0]), SUM_x1, SUM_x2],  # Тут уже нормально - одна строка это все признаки одной точки
                  [SUM_x1, SUM_x1_x1, SUM_x1_x2],
                  [SUM_x2, SUM_x1_x2, SUM_x2_x2]])
b = np.array([SUM_y, SUM_x1_y, SUM_x2_y])
k_all = solve(X_,b)
print(f'k_all = {k_all}\n')  # [ 0.5275   0.99375 -0.15625]

Посмотреть матричное решение

X = np.array([[1,1,1,1,1], [1,2,3,4,5], [2,3,1,5,4]]).T  #Сразу транспонирую
y = np.array([[1.5, 1.8, 3.2, 3.6, 5.1]]).T
X_T_X_inv = np.linalg.inv(X.T.dot(X))
matrix_b = X_T_X_inv.dot(X.T.dot(y))
print(f'b = {matrix_b}')  # b = [ 0.5275   0.99375 -0.15625]

Посмотреть заводское решение, python

import scipy
X = np.array([[1,1,1,1,1], [1,2,3,4,5], [2,3,1,5,4]]).T   
y = np.array([[1.5, 1.8, 3.2, 3.6, 5.1]]).T
b, squared_error_sum, matrix_rank, SVD_ = scipy.linalg.lstsq(X, y)
print(b)
#[[ 0.5275 ], [ 0.99375], [-0.15625]]

Нельзя не добавить:

solve и lstsq есть не только в scipy.linalg, но и в numpy.linalg
Это (наличие lstsq в numpy и scipy) логично, так как МНК и вообще ЛР появились задолго до ML. Учебники по статистике и эконометрике могут быть полезны.
Градиентный спуск, который я тоже запустил, чисто для сравнения с learning_rate = 0.005 до идеального значения, полученного МНК добрался только за 13к итераций. а с learning_rate = 0.01 его вовсе разболтало - если МНК считается и считается за уместное время, то лучше воспользоваться им.

И это всё про МНК, что мне нужно знать. Надеюсь.

Теги:

Хабы:

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

База

МНК (метод наименьших квадратов) aka OLS (ordinary least squares)

Матричная форма

Публикации

Истории

Работа

Ближайшие события