如何在Python中拟合p维二次多项式并精准计算梯度与海森矩阵？

你完全可以利用model.coef_中的精确系数直接推导梯度和海森矩阵，无需依赖易受浮点误差影响的有限差分。以下是几种可行方案：

方案一：基于scikit-learn手动解析系数计算 #

scikit-learn的PolynomialFeatures支持获取特征对应的多项式项名称，结合model.coef_就能精准计算梯度和海森矩阵：

代码示例 #

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
import numpy as np

# 生成示例数据
np.random.seed(42)
X = np.random.rand(100, 2)  # 2维输入
y = 2 + 3*X[:,0] + 5*X[:,1] + 1.2*X[:,0]**2 + 0.8*X[:,0]*X[:,1] + 2.1*X[:,1]**2 + np.random.randn(100)*0.1

# 拟合多项式回归
poly = PolynomialFeatures(degree=2, include_bias=False)
Xp = poly.fit_transform(X)
model = LinearRegression()
model.fit(Xp, y)

# 获取多项式特征项与对应系数
feature_names = poly.get_feature_names_out(['x0', 'x1'])
coef = model.coef_

def compute_gradient(x):
    """计算给定点x的梯度"""
    grad = np.zeros_like(x)
    for name, c in zip(feature_names, coef):
        # 解析多项式项的变量与指数
        var_counts = {}
        for term in name.split(' '):
            if '^' in term:
                var, exp = term.split('^')
                var_counts[var] = int(exp)
            else:
                var_counts[term] = var_counts.get(term, 0) + 1
        # 计算每个变量的偏导贡献
        for var, exp in var_counts.items():
            idx = int(var[-1])
            if exp == 0:
                continue
            # 偏导系数：c * 指数
            deriv_c = c * exp
            # 偏导项的乘积部分
            deriv_term = 1.0
            for v, e in var_counts.items():
                if v == var:
                    deriv_term *= x[int(v[-1])] ** (e-1) if e > 1 else 1
                else:
                    deriv_term *= x[int(v[-1])] ** e
            grad[idx] += deriv_c * deriv_term
    return grad

def compute_hessian(x):
    """计算给定点x的海森矩阵"""
    p = len(x)
    hess = np.zeros((p, p))
    for name, c in zip(feature_names, coef):
        var_counts = {}
        for term in name.split(' '):
            if '^' in term:
                var, exp = term.split('^')
                var_counts[var] = int(exp)
            else:
                var_counts[term] = var_counts.get(term, 0) + 1
        # 1. 同一变量的二阶偏导
        for var, exp in var_counts.items():
            if exp >= 2:
                idx = int(var[-1])
                deriv_c = c * exp * (exp-1)
                deriv_term = 1.0
                for v, e in var_counts.items():
                    if v == var:
                        deriv_term *= x[int(v[-1])] ** (e-2) if e > 2 else 1
                    else:
                        deriv_term *= x[int(v[-1])] ** e
                hess[idx, idx] += deriv_c * deriv_term
        # 2. 不同变量的混合偏导
        vars_list = list(var_counts.keys())
        if len(vars_list) >= 2:
            for i in range(len(vars_list)):
                for j in range(i+1, len(vars_list)):
                    var1, var2 = vars_list[i], vars_list[j]
                    exp1, exp2 = var_counts[var1], var_counts[var2]
                    idx1, idx2 = int(var1[-1]), int(var2[-1])
                    deriv_c = c * exp1 * exp2
                    deriv_term = 1.0
                    for v, e in var_counts.items():
                        if v == var1:
                            deriv_term *= x[int(v[-1])] ** (e-1) if e > 1 else 1
                        elif v == var2:
                            deriv_term *= x[int(v[-1])] ** (e-1) if e > 1 else 1
                        else:
                            deriv_term *= x[int(v[-1])] ** e
                    hess[idx1, idx2] += deriv_c * deriv_term
                    hess[idx2, idx1] += deriv_c * deriv_term
    return hess

# 测试计算
x_test = np.array([0.5, 0.5])
print("梯度：", compute_gradient(x_test))
print("海森矩阵：\n", compute_hessian(x_test))

方案二：使用PyTorch自动求导（更简洁） #

PyTorch的自动求导机制可以直接对拟合的多项式模型计算梯度和海森矩阵，无需手动推导偏导公式：

代码示例 #

import torch
import torch.nn as nn
import numpy as np
from sklearn.preprocessing import PolynomialFeatures

# 生成示例数据
np.random.seed(42)
X_np = np.random.rand(100, 2).astype(np.float32)
y_np = 2 + 3*X_np[:,0] + 5*X_np[:,1] + 1.2*X_np[:,0]**2 + 0.8*X_np[:,0]*X_np[:,1] + 2.1*X_np[:,1]**2 + np.random.randn(100)*0.1
y_np = y_np.reshape(-1, 1).astype(np.float32)

# 转换为PyTorch张量
X = torch.tensor(X_np)
y = torch.tensor(y_np)

# 定义多项式模型
class PolynomialModel(nn.Module):
    def __init__(self, degree=2, input_dim=2):
        super().__init__()
        self.poly = PolynomialFeatures(degree=degree, include_bias=False)
        self.linear = nn.Linear(self.poly.fit_transform(np.zeros((1, input_dim))).shape[1], 1)
    
    def forward(self, x):
        x_poly = torch.tensor(self.poly.transform(x.detach().numpy()), dtype=torch.float32)
        return self.linear(x_poly)

# 拟合模型
model = PolynomialModel(degree=2, input_dim=2)
criterion = nn.MSELoss()
optimizer = torch.optim.SGD(model.parameters(), lr=0.1)

for epoch in range(1000):
    optimizer.zero_grad()
    outputs = model(X)
    loss = criterion(outputs, y)
    loss.backward()
    optimizer.step()

# 计算梯度与海森矩阵
x_test = torch.tensor([[0.5, 0.5]], requires_grad=True)
y_pred = model(x_test)

# 梯度
y_pred.backward(create_graph=True)
gradient = x_test.grad.clone()
print("梯度：", gradient.numpy())

# 海森矩阵
hessian = []
for i in range(x_test.shape[1]):
    grad_i = gradient[0, i]
    x_test.grad.zero_()
    grad_i.backward(retain_graph=True)
    hessian_row = x_test.grad.clone().numpy()[0]
    hessian.append(hessian_row)
hessian = np.array(hessian)
print("海森矩阵：\n", hessian)

内容的提问来源于stack exchange，提问作者SapereAude