3

кого заинтересовал мой заголовок, сейчас я попробую как можно подробнее и интереснее описать проблему с которой я столкнулся в процессе изучения основ машинного обучения. Решая учебную задачу с 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)

Если вы дочитали до этого момента, то мне было бы очень интересно услышать ваши предположения, с чем могут быть связанны мои ошибки, а так же, у меня есть глупый вопрос, на который я пока не нашёл ответ: В чём преимущество использования метода градиентного спуска перед обычным решением СЛАУ в данной задаче? Ведь реализация МГС встречалась чаще в процессе гуглёшки.

2
  • 1
    Градиентный спуск нужен чтобы не решать уравнения. Не всегда функция такая хорошая, что можно просто решить линейное уравнение, поэтому и нужен градиентный спуск. А ещё лучше - стохастический градиентный спуск, за счёт стохастичности он лучше обходит локальные минимумы.
    – CrazyElf
    17 фев в 0:03
  • 1
    Вставлю свои пять копеек. Можно же найти коэффициенты при параметрах с помощью матричных операций: Составляется матрица А размерностью n*k, где n - число наблюдений, k - индекс последнего коэффициента с учётом того, что свободный член имеет нулевой индекс. Затем решается что-то типа такого ` np.linalg.inv(A.transpose() @ A) @ A.transpose() @ df.Y` На выходе получим вектор коэффициентов. Как то так)) 17 фев в 9:07

1 ответ 1

1

Перечитывая вопрос, я понял как решить проблему в методе решения обычной СЛАУ, достаточно было всего лишь указать 64-битный тип, dtype=numpy.float64

x = np.array([[i] for i in df["LotArea"]], dtype=np.float64)
y = np.array([i for i in df["SalePrice"]], dtype=np.float64)

И проблема решилась :) Но в методе решения через градиентный спуск, всё ещё осталась

3
  • 1
    xD это метод утки 16 фев в 22:57
  • 1
    @Senior Pomidor, ахах, забавно, но в процессе написания вопроса, на моём столе действительно была резиновая утка))
    – Gari
    16 фев в 23:07
  • 1
    А ведь есть ещё и np.float128 )
    – CrazyElf
    17 фев в 0:01

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.