Метод наименьших квадратов06.07.2024 22:00

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

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

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

База

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

k1*X1 + k2*X2...+ km*Xm + C = y

и нахождение k и C по известным y и X.

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

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

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

Что важно:

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

Решением для n=m и n=m+1 будет нахождение идеально подходящих значений коэффициентов. Подходящих для этих 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 хотят выбрать функцию ошибки и минимизировать её пользуясь двумя фактами:
1) Производная любой функции это быстрота изменения значения этой функции при наращивании независимой переменной.
2 Я заранее знаю, чему функция ошибки должна быть равна — нулю. Я ищу коэффициенты там, где функция или ноль, или минимальна. Это то место, где производная по k_i слева отрицательна (увеличение k_i приводит к снижению функции ошибки), а справа положительна — приводит к росту значения ошибки. И равна производная чему-то между отрицательными и положительными значениями — »0».
Часто предстоит брать производные по стандартным функциям ошибок и как-то совать в них X_train, y_train и как-то усредняться.
В ML любят квадратичную ошибку.

За что?

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

Формулы часто выглядят как (Фигня — y)^2 и с непривычки кажутся одинаковыми. Надо смотреть в первую очередь на контекст, на то в какой момент Фигня подсчитана и зачем. Придётся различать.

Я тоже возьму квадратичную ошибку (мог бы и среднюю, так как производная у них одинаковая, но МНК этого не подразумевает) и минимизирую её квадрат. 2 признака + С.

$\sum_{i=1}^{} (y_i - k_1 x_{1i} + k_2 x_{2i} + C)^2 = 0$

X-ы у меня есть в исходных данных, k — нет. Я буду дифференцировать по k, а X у меня будет пока что неким коэффициентом:

$\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})$

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

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

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

$\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*}$

Суммарно это новая СЛАУ. Именно СЛАУ, так как x1_i хоть в первой хоть в сотой степени это просто число. Немного отредачу СЛАУ, а именно «развалю» 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}$

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

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

Подставив, перемножив и сложив, я получил:
15.2 = 5*C + 15*k1 + 15*k2
54.6 = 15*C + 55*k1 + 51*k2
50.0 = 15*C + 51*k1 + 55*k2

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

MSE = 0.0971

Посмотреть велосипед в 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]

Посмотреть матричный велосипед (Тут тоже немного сложно, но можно и пропустить)

Нашёл такую формулу:

МНК (ols) это число, соответственно все слагаемые должны быть числами, проверю хотя бы это:

y.T очевидно можно перемножить на y, это изначально столбец, так что в итоге число будет.
b.T dot X.T это строка (L=n) умноженная на таблицу n столбцов и m строк. Строка умножаемая на таблицу даёт строку.
b.T dot X.T dot y это строка n на столбец n, в итоге число
b.T dot X.T разбиралось на 2 строки выше. Строка n
b.T dot X.T dot X это строка n на таблицу высоты n. Строка n
b.T dot X.T dot X dot b это строка n на столбец n, число

Матричное дифференцирование похоже на обычное. У b снимаем степень и добавляем двойку перед слагаемым если надо. Слагаемые без b выкидываем

$X^\top \cdot y = X^\top \cdot X \cdot b$

Осталось выразить b = X.T. dot (X) / X.T. dot (y), но, с матрицами так нельзя. Вместо этого матрица, обратная делимому матрично умножается на частное
X_T_X_inv = np.linalg.inv(X.T.dot(X)) matrix_b = X_T_X_inv.dot(X.T.dot(y))Если матрица содержит строки, получаемые из других строк масштабированием, обратная матрица не найдется, для этого есть псведообратные матрицы и метод pinv.

Почему так?

Матрица это всего-лишь табличная форма записи СЛАУ (Почти всегда). СЛАУ с 3 неизвестными и с 3 уравнениями, где одно уравнение получается масштабированием другого тоже не решаема, а A.dot(inv(A)) = ЕД. МАТРИЦА это тоже по сути уравнение и ища обратную матрицу, я его именно решаю. И строки и столбцы равнозначимы в этом решении.

Теперь проверю через python:

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   
# Тут надо уже столбцами поставить
# И добавить столбец единиц. Добавляю в начало типо k0, но без разницы
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)  # вообще в ЛР обычно пишут b или w, математическое k не в почёте
#[[ 0.5275 ], [ 0.99375], [-0.15625]]

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

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

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