Как лучше всего подогнать квадратичный многочлен к p-мерным данным и вычислить его градиент и матрицу Гессе? ⇐ Python

[Anonymous]

Я пытался использовать библиотеку scikit-learn для решения этой проблемы. Примерно:
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# Make or load an n x p data matrix X and n x 1 array y of the corresponding
# function values.

poly = PolynomialFeatures(degree=2)
Xp = poly.fit_transform(X)
model = LinearRegression()
model.fit(Xp, y)

# Approximate the derivatives of the gradient and Hessian using the relevant
# finite-difference equations and model.predict.

Как показано выше, sklearn делает выбор в пользу разделения полиномиальной регрессии на PolynomialFeatures и LinearReгрессия, а не объединения их в одну функцию. У такого разделения есть концептуальные преимущества, но есть и серьезный недостаток: оно фактически не позволяет model предлагать методы градиента и гессиана, а модель была бы значительно полезнее, если бы это было так.
Мой текущий обходной путь использует уравнения конечных разностей и model.predict для аппроксимации элементов градиента и гессиана (как описано здесь). Но мне не нравится этот подход — он чувствителен к ошибкам с плавающей запятой, а «точная» информация, необходимая для построения градиента и гессиана, уже содержится в model.coef_.
Есть ли более элегантный и точный метод подбора p-мерного полинома и нахождения его градиента и гессиана в Python? Меня бы устроил тот, который использует другую библиотеку.

Обновить. На мой взгляд, самый простой метод (самый простой на программист), похоже, использует simpy, примерно так, как описано в ответе Джона. Вот мой окончательный код, который был отправлен как ожидающее редактирования ответа Джона.
from numpy.random import rand
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sympy import symbols, prod, diff
from numpy import sum, linspace

n = 200
p = 4
X = rand(n, p)
y = rand(n)

poly = PolynomialFeatures(degree=2)
Xp = poly.fit_transform(X)
model = LinearRegression()
model.fit(Xp, y)

powers = poly.powers_
coeffs = model.coef_
syms = symbols(f'x0:{p}')

f = sum([coeffs*prod(syms**powers) for i in range(len(coeffs))])

Df = [diff(f, xj) for xj in syms]
H = [[diff(djf, xk) for xk in syms] for djf in Df]
print('Gradient:', Df)
print('Hessian:', H)

x0 = linspace(0.2, 0.7, p) # test point
Df_x0 = [Df[k].evalf(subs={syms[j]:x0[j] for j in range(p)}) for k in range(p)]
print('x0:', x0)
print('Gradient at x0:', Df_x0)

Но использовать simpy не требуется, как показывает отличный ответ bb1. Вот альтернативный способ написать ответ bb1, который мне кажется немного легче читать.
from numpy import *

x0 = [0, 1]

# Define p, X, y, poly, model, powers, and coeffs.

grad = array(p * [float(0)])
for i in range(p):
gc = powers.T * coeffs # gradient coefficients
gp = powers.copy() # gradient powers
gp[:,i] = maximum(0, gp[:,i] - 1)
grad = sum([gc[j]*prod(x0**gp[j]) for j in range(len(coeffs))])
print('grad = ', grad)

hess = p * [None]
for i in range(p):
gc = powers.T * coeffs
gp = powers.copy()
gp[:,i] = maximum(0, gp[:,i] - 1)
hess = p * [float(0)]
for j in range(p):
hc = gp.T[j] * gc # Hessian coefficients
hp = gp.copy() # Hessian powers
hp[:,j] = maximum(0, hp[:,j] - 1)
hess[j] = sum([hc[k]*prod(x0**hp[k]) for k in range(len(coeffs))])
hess = array(hess)
print('hess = ', hess)


Подробнее здесь: https://stackoverflow.com/questions/790 ... ta-and-com