如何使用numpy函数计算以下Hessian以加快计算速度？

Question

我需要实现一个功能，来计算逻辑损失的海森矩阵（hessian），这个海森矩阵是通过对一些指数项的对数求和得到的。我在Python中实现了以下这个函数：

def hessian(self,w,hess_trick=0):
        hess = 0
        for x_i,y_i in zip(self.data, self.labels):
            hess += np.exp(y_i * np.dot(w.T, x_i))/((1 + np.exp(y_i * np.dot(w.T,x_i)))**2) *       np.outer(x_i,x_i.T)
        return hess + lambda_reg * np.identity(w.shape[0]) + hess_trick * 10**(-12) *    np.identity(w.shape[0])

我的问题是，如何能写一个等效的、但更快的函数，而不使用慢吞吞的Python循环？

因为我对numpy不太自信，所以我尝试写了以下这个函数：

    def new_hessian(self, w, hess_trick=0):
        exp_term = np.exp(self.labels * np.dot(self.data, w))
        sigmoid_term = 1 + exp_term
        inv_sigmoid_sq = 1 / sigmoid_term ** 2

        diag_elements = np.sum((exp_term * inv_sigmoid_sq)[:, np.newaxis] * self.data ** 2, axis=0)
        off_diag_elements = np.dot((exp_term * inv_sigmoid_sq) * self.data.T, self.data)
        hess = np.diag(diag_elements) + off_diag_elements
        regularization = lambda_reg * np.identity(w.shape[0])

        hess += hess_trick * 1e-12 * np.identity(w.shape[0])

        return hess + regularization

在调试这个函数时，我发现了一个根本性的问题。当特征数量很少（比如少于200）时，这两个海森矩阵的实现结果并不相等。但是当特征数量增加时，这两个函数的结果似乎就相等了。问题在于，当我用牛顿法来优化逻辑损失时，那个更快的实现需要更多的迭代才能收敛，而第一个（但运行较慢的）实现则收敛得更快。

Answer 1

这里是同样内容的向量化版本：

a = np.exp(y * (X @ w))

## Method 1
b = np.sqrt(a)*X.T/(1 + a)
b @ b.T 

## method 2
X.T @ np.diag(a/(1 + a)**2) @ X 

## method 3
np.einsum('ij,i,ik',X, a/(1 + a)**2, X)

编辑

你提供的链接中的公式是错误的。请注意，计算对数似然函数的二阶导数时，和权重有关的部分并不需要标签。所以你不需要标签。这里有一个关于逻辑回归的正确海森矩阵的链接。

现在有了这个正确的公式，我们可以这样计算海森矩阵：

 a = np.exp(X @ w)
 - X.T @ np.diag(a/(1 + a)**2) @ X

真实示例：

import pandas as pd, numpy as np
dat = pd.read_csv("https://raw.githubusercontent.com/oonyambu/ECON430/main/admit.csv", index_col = 0)

X = np.c_[np.ones(dat.shape[0]), dat.iloc[:,1:3].to_numpy()]
y = dat.admit.to_numpy()

from statsmodels.formula.api import logit
mod = logit("admit~gre+gpa", dat).fit()

# Hessian = -inverse(covariance)
np.linalg.inv(-mod.normalized_cov_params)
array([[-8.25030780e+01, -4.98883263e+04, -2.84169634e+02],
       [-4.98883263e+04, -3.11768300e+07, -1.72965203e+05],
       [-2.84169634e+02, -1.72965203e+05, -9.89840314e+02]])

mod.model.hessian(mod.params)
array([[-8.25030780e+01, -4.98883263e+04, -2.84169634e+02],
       [-4.98883263e+04, -3.11768300e+07, -1.72965203e+05],
       [-2.84169634e+02, -1.72965203e+05, -9.89840314e+02]]) 

## Direct Computation. (No y)
a = np.exp(X @ mod.params)
b = np.sqrt(a)*X.T/(1 + a)
-b @ b.T
array([[-8.25030780e+01, -4.98883263e+04, -2.84169634e+02],
       [-4.98883263e+04, -3.11768300e+07, -1.72965203e+05],
       [-2.84169634e+02, -1.72965203e+05, -9.89840314e+02]])

或者甚至可以这样：

k = 1/(1 + np.exp(-X @ mod.params))
-k*(1-k)*X.T @ X
array([[-8.25030780e+01, -4.98883263e+04, -2.84169634e+02],
       [-4.98883263e+04, -3.11768300e+07, -1.72965203e+05],
       [-2.84169634e+02, -1.72965203e+05, -9.89840314e+02]])