решать проблему в лоб сильно затратно по памяти: лучше гляньте как задачу регрессии оптимизируют mashine learning библиотеки. но в лоб
на numpy это выглядит так:
XtX = X.T * X
result = inv(XtX).dot(X.T).dot(Y)
X = матрица MxN*K где строкам соответствуют точкам а столбцам координаты возведённые в соответствующие степени
Y = вектор столбец из M нулей
Matrix.T - транспонирование матрицы
inv(Matrix) -- функция inv вычисляет обратную матрицу
MatrixA.dot(MatrixB) -- вычисляет произведение матриц
На выходе получаете вектор result -- коэффициэнты полинома соответствующие колонкам матрицы X (Соответственно, если хотите, чтобы свободный коэфициент был, столбец единиц в матрице X должен быть)
ещё раз обращу внимание, что матрица XtX имеет размерность MxM что при попытке обработать любые реальные данные с миллионами измерений просто сожжёт вам то, на чём вы там считаете.
За нотацию извиняюсь -- классическую математическую подзабыл, и честно говоря вспоминать не очень хочется.
UPD:
После некоторых посиделок оформил всё в такую конструкцию:
import numpy as np
K = 2 #степень целевой плоскости
points = np.array([[a, 5 + 2*a + 6 * a**2] for a in range(0, 60)]) # точки, по которым строим
class Symbolic:
def __init__(self, letter=None):
self.letters = dict() if letter is None else {letter : 1}
def addLetter(self, letter, power=1):
if letter not in self.letters:
self.letters[letter] = 0
self.letters[letter] += power
def __pow__(self, pow):
res = Symbolic()
for letter in self.letters:
res.letters[letter] = self.letters[letter] * pow
return res
def __mul__(self, right):
res = Symbolic()
for letter in self.letters:
res.addLetter(letter, self.letters[letter])
for letter in right.letters:
res.addLetter(letter, right.letters[letter])
return res
def __repr__(self):
return '*'.join(map(lambda letter: f'{letter}**{self.letters[letter]}', self.letters.keys()))
def pointsToPower(points, power):
result = points ** power
if points.shape[1] > 1 and power > 1:
for startpower in range(power-1, 0, -1):
for index in range(points.shape[1] - 1) :
base = np.array(points[:,index]) ** startpower
submatrix = pointsToPower(
points[:, index+1:],
power - startpower
)
result = np.hstack([result, np.array([submatrix[lineIndex, :] * base[lineIndex] for lineIndex in range(base.shape[0])])])
return result
def pointsToX(points, K=1):
return np.hstack([
# np.array([[1]] * points.shape[0]),
*[pointsToPower(points, power=i+1) for i in range(K)]
])
def krammer(X, sigma=1):
XtX = X.T.dot(X)
baseDet = np.linalg.det(XtX)
if baseDet == 0:
return [[None]]
Y = np.array([[sigma] * X.shape[0]]).T
return np.linalg.inv(XtX).dot(X.T).dot(Y)
X = pointsToX(points,K)
values = krammer(X).T[0]
if values[0] is None:
print('''Обратной матрицы не существует, и решить систему матричным методом невозможно.
В этом случае система решается методом исключения неизвестных (методом Гаусса).
Но его реализовывать мне влом''')
symbols = pointsToX(np.array([[Symbolic(chr(ord('a') + i)) for i in range(points.shape[1])]]),K)[0]print('Уравнение целевой плоскости:')
print(*[f'{value} * {symbol}' for symbol,value in zip(symbols, values)], sep=' + \n', end=' = 1')
