Проблема линейной регрессии
кого заинтересовал мой заголовок, сейчас я попробую как можно подробнее и интереснее описать проблему с которой я столкнулся в процессе изучения основ машинного обучения. Решая учебную задачу с kaggle, я захотел сравнить различные модели в решении одной и той же задачи, чтобы сравнить их эффективность в различных условиях, но использовать реализацию из sklearn не интересно, поэтому я решил написать собственные аналоги моделей, и решено было начать с линейной регрессии.
Класс "model" должен используя функцию Lineal_Regression(x, y), принимающую на вход:
x - numpy массив из массивов записей нашего датасета
y - numpy массив из результатов, которые мы хотим предсказывать для тестовых данных
import numpy as np
class model():
def Lineal_Regression(self, x, y):...
def factorial(self, n):...
def RepComb(self, m, n):...
def RSS(self, x, y):...
Начнём с вспомогательных функций, как не трудно догадаться factorial(n) ищет факториал, а RepComb(m, n) перестановки с повторениями.
Как и в LinearRegression() из sklearn, было решено взять за функцию ошибки RSS т.е. так называемые наименьшие квадраты. Мы создаём функцию, которая считает сумму квадратов отклонений и пытаемся её минимизировать и тем самым найти заветные коэффициенты
В моей реализация функция принимает на вход вышеописанные x и y и выдаёт массив коэффициентов вида: [ω0^2, ω1^2, ..., ωn^2, ω0ω1, ω0ω2, ..., ω0ωn, ω1ω2, ..., ω1*ωn, ..., ω(n-1)ωn, ω0, ω1, ..., ωn, a]
Где ω0^2 - коэффициент стоящий перед β0^2, ω1 - коэффициент перед β1 и так далее. А самый последний элемент это константное значение функции RSS. Проще говоря он раскрывает скобки и находит коэффициенты при неизвестных.
Например для x = [1, 2], y = [2, 3]
SSR = (1 * ω1 + ω0)^2 + (2 * ω1 + ω0)^2
На что моя функция выдаст:
[ 2. 5. 6. -10. -16. 13.]
Что означает SSR = 2ω0^2 + 5ω1^2 + 6ω0ω1 - 10ω0 - 16ω1 + 13
Вот её полная реализация:
def RSS(self, x, y):
x = np.insert(x, 0, [1] * x.shape[0], axis=1)
res = np.array([0.0] * int(self.RepComb(2, x.shape[1]) + x.shape[1] + 1))
for k in range(x.shape[0]):
res[-1] += y[k]**2
for i in range(x.shape[1]):
res[i] += x[k][i]**2
i = x.shape[1]
for j in range(x.shape[1] - 1):
for l in range(j + 1, x.shape[1]):
res[i] += 2 * x[k][j] * x[k][l]
i += 1
for j in range(x.shape[1]):
res[i] += -2 * x[k][j] * y[k]
i += 1
return res
И теперь мы переходим к самому главному, к реализации самой Lineal_Regression(x, y) перед которой ставится задача отыскать ω0, ..., ωn такие, что функция RSS является минимальной. Первая моя идея заключалась в том, чтобы найти частные производные функции и прировнять их к нулю, мы получаем (n + 1) переменную и (n + 1) неизвестную, находим стационарную точку, которая исходя из здравого смысла будет минимальной, так как устремляя ω0, ..., ωn в разные стороны на бесконечность, значение RSS будет только увеличиваться.
Зная коэффициенты, задача показалась мне не сложной, её графическое значение я попытался изобразить от руки, несмотря на то, что не одарён художественными навыками:
Так как матрица производных будет симметричной, достаточно найти только верхний диагональный вид, после чего прибавить транспонированную и вычесть одну диагональку.
Так же необходимо найти вектор решений:
После чего решить обычную СЛАУ
И случилось чудо, моя линейная регрессия начала выдавать результаты идентичные результатам LinealRegression() из пакета sklearn, например:
x = np.array([[1.7], [1.62], [1.7], [1.85], [1.85], [1.63], [1.5]])
y = np.array([63, 56, 54, 75, 80, 59, 49])
model1 = model()
res1 = model1.Lineal_Regression(x, y)
m1 = linear_model.LinearRegression()
m1.fit(x, y)
print(res1)
print(m1.coef_)
Программа выдаёт:
[-78.06313282 82.90649196]
[82.90649196]
Так же хорошо работает и для многомерного случая:
x = np.array([[1, 2, 4, 6, 5], [2, 3, 3, 4, 5], [4, 2, 5, 4, 3], [0, 0, 0, 0, 0], [1, 4, 4, 5, 2], [6, 8, 1, 6, 5]])
y = np.array([2, 4, 2.5, 0, 3, 2])
model1 = model()
res1 = model1.Lineal_Regression(x, y)
m1 = linear_model.LinearRegression()
m1.fit(x, y)
print(res1)
print(m1.coef_)
Вывод:
[ 1.19371180e-14 -5.62072893e-01 1.02961276e+00 1.12243736e+00
-1.32744875e+00 7.95558087e-01]
[-0.56207289 1.02961276 1.12243736 -1.32744875 0.79555809]
Но самое больное началось, когда я начал применять модель в реальных боевых условиях, при большом количестве записей она начинает ломаться даже в случае одной переменной, вот пример для датасета:
import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
df = pd.read_csv("C://Users/79178/Desktop/House Prices/train.csv")
fig, ax = plt.subplots()
ax.scatter(df["LotArea"], df["SalePrice"])
x = np.array([[i] for i in df["LotArea"]])
y = np.array([i for i in df["SalePrice"]])
model3 = model()
res3 = model3.Lineal_Regression(x, y)
ax.plot([0, 200000], [0 * res3[1] + res3[0], 500 * res3[1] + res3[0]])
m3 = linear_model.LinearRegression()
m3.fit(x, y)
r3 = m3.coef_
ax.plot([0, 200000], [0 * r3[0], 500 * res3[1]])
print(res3)
print(r3)
Вывод:
[-4.28049172e+07 4.08733870e+03]
[2.09997195]
И выходит предупреждение:
C:\Users\79178\AppData\Local\Temp\ipykernel_14316\861510427.py:50: RuntimeWarning: overflow encountered in long_scalars
res[i] += -2 * x[k][j] * y[k]
По графику видно, что произошла ошибка как по переменной ω1 так и по ω0. Мной было замечено, что повышая количество записей, всё работает хорошо до какого-то момента, в моём случае до 54ой записи всё совпадает с sklearn решением и на графике выглядит адекватно, но потом, что-то идёт не по плану и решение ломается:
Вот полный код реализации
class model():
def Lineal_Regression(self, x, y):
if x.shape[0] == y.shape[0]:
coefs_of_Rss = self.RSS(x, y)
DifMatrix = np.zeros((x.shape[1] + 1, x.shape[1] + 1))
vector_sol = -1 * coefs_of_Rss[(x.shape[1]+2) * (x.shape[1] + 1) // 2: (x.shape[1]+2) * (x.shape[1] + 1) // 2 + x.shape[1] + 1]
sdvig = 0
for i in range(x.shape[1] + 1):
for j in range(i, x.shape[1] + 1):
if i == j:
DifMatrix[i][j] = 2 * coefs_of_Rss[i]
else:
DifMatrix[i][j] = coefs_of_Rss[x.shape[1] + j + i * x.shape[1] - sdvig]
sdvig += i + 1
DifMatrix = DifMatrix + DifMatrix.transpose() - (np.eye(x.shape[1] + 1) * DifMatrix)
return (np.linalg.solve(DifMatrix, vector_sol))
def factorial(self, n):
if n <= 1:
return 1
else:
return self.factorial(n - 1) * n
def RepComb(self, m, n):
return self.factorial(m + n - 1) / (self.factorial(m) * self.factorial(n - 1))
def RSS(self, x, y):
x = np.insert(x, 0, [1] * x.shape[0], axis=1)
res = np.array([0.0] * int(self.RepComb(2, x.shape[1]) + x.shape[1] + 1), dtype=np.float64)
for k in range(x.shape[0]):
res[-1] += y[k]**2
for i in range(x.shape[1]):
res[i] += x[k][i]**2
i = x.shape[1]
for j in range(x.shape[1] - 1):
for l in range(j + 1, x.shape[1]):
res[i] += 2 * x[k][j] * x[k][l]
i += 1
for j in range(x.shape[1]):
res[i] += -2 * x[k][j] * y[k]
i += 1
return res
Но на этом я не отчаялся и погуглив понял, что можно реализовать данное решение при помощи итерационного метода градиентного спуска, реализация которого совпадает с первой до нахождения коэффициентов матрицы частных производных, отличается только часть функции Lineal_Regression(x,y)
def Lineal_Regression(self, x, y):
if x.shape[0] == y.shape[0]:
coefs_of_Rss = self.RSS(x, y)
DifMatrix = np.zeros((x.shape[1] + 1, x.shape[1] + 1))
vector_sol = np.array([coefs_of_Rss[(x.shape[1]+2) * (x.shape[1] + 1) // 2: (x.shape[1]+2) * (x.shape[1] + 1) // 2 + x.shape[1] + 1]])
v = np.zeros((x.shape[1] + 1, 1))
vector_sol = vector_sol.T
sdvig = 0
for i in range(x.shape[1] + 1):
for j in range(i, x.shape[1] + 1):
if i == j:
DifMatrix[i][j] = 2 * coefs_of_Rss[i]
else:
DifMatrix[i][j] = coefs_of_Rss[x.shape[1] + j + i * x.shape[1] - sdvig]
sdvig += i + 1
DifMatrix = DifMatrix + DifMatrix.transpose() - (np.eye(x.shape[1] + 1) * DifMatrix)
for _ in range(10000):
v = v - 0.001*((np.dot(DifMatrix, v)) + vector_sol)
return v
Шаг в 0.001 и количество итераций 10000 были взяты из головы, я сначала хотел убедиться в эффективности метода, а потом его оптимизировать, но как и в прошлом случае, на маленьком наборе данных всё работает хорошо с заданной точностью, но при увеличении этих самых данных, происходит ошибка: v уходит на бесконечность по всем координатам, чего я уже никак не могу понять, даже если шаг подобран неудачно, v всё равно должен "прыгать" вокруг правильного решения, но никак не уходить на бесконечность вот пример для того же датасета:
import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
df = pd.read_csv("C://Users/79178/Desktop/House Prices/train.csv")
fig, ax = plt.subplots()
ax.scatter(df["LotArea"], df["SalePrice"])
x = np.array([[i] for i in df["LotArea"]], dtype=np.float64)
y = np.array([i for i in df["SalePrice"]], dtype=np.float64)
model3 = model()
res3 = model3.Lineal_Regression(x, y)
ax.plot([0, 200000], [0 * res3[1] + res3[0], 500 * res3[1] + res3[0]])
m3 = linear_model.LinearRegression()
m3.fit(x, y)
r3 = m3.coef_
ax.plot([0, 200000], [0 * r3[0], 500 * res3[1]])
print(res3)
print(r3)
Вывод:
[[nan]
[nan]]
[2.09997195]
Предупреждение:
C:\Users\79178\AppData\Local\Temp\ipykernel_14316\78557657.py:26: RuntimeWarning: invalid value encountered in subtract
v = v - 0.001*((np.dot(DifMatrix, v)) + vector_sol)
c:\users\79178\appdata\local\programs\python\python39\lib\site-packages\numpy\core\shape_base.py:65: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
ary = asanyarray(ary)
Если вы дочитали до этого момента, то мне было бы очень интересно услышать ваши предположения, с чем могут быть связанны мои ошибки, а так же, у меня есть глупый вопрос, на который я пока не нашёл ответ: В чём преимущество использования метода градиентного спуска перед обычным решением СЛАУ в данной задаче? Ведь реализация МГС встречалась чаще в процессе гуглёшки.