找到这个函数的 Hessian 矩阵

Question

您好，我有以下功能：

sum from 1 to 5000 -log(1−(xi)^2) -log(1-(a_i)^t*x)，其中 a_i 是一个随机向量，我们试图通过 Netwon 最小化这个函数的值方法。

我需要一种方法来计算关于 (x1, x2, x3, ...) 的 Hessian 矩阵。我尝试了自动渐变，但花了太多时间。这是我现在的时间。

from autograd import elementwise_grad as egrad
from autograd import jacobian
import autograd.numpy as np

x=np.zeros(5000);
a = np.random.rand(5000,5000)
def f (x):
  sum = 0;
  for i in range(5000):
      sum += -np.log(1 - x[i]*x[i]) - np.log(1-np.dot(x,a[i]))
  
  return sum;

df = egrad(f)
d2f = jacobian(egrad(df)); 
print(d2f(x));

我已尝试调查 sympy，但我对如何继续感到困惑。

Answer 1

PyTorch 具有 GPU 优化的 hessian 操作：

import torch

torch.autograd.functional.hessian(func, inputs)

Answer 2

您可以使用常规的NumPy 向量化数组操作，这将显着加快程序的执行速度：

from autograd import elementwise_grad as egrad
from autograd import jacobian
import autograd.numpy as np
from time import time
import datetime


n = 5000
x = np.zeros(n)
a = np.random.rand(n, n)

f = lambda x: -1 * np.sum(np.log(1-x**2) + np.log(1-np.dot(a, x)))

t_start = time()
df = egrad(f)
d2f = jacobian(egrad(df))
t_end = time() - t_start
print('Execution time: ', datetime.datetime.fromtimestamp(t_end).strftime('%H:%M:%S'))

输出

Execution time: 02:02:27 

一般来说，每次处理数字数据时，都应该尽量避免使用loops进行计算，因为它们的头部和计数器变量的维护通常会成为程序的瓶颈。

另一方面，NumPy 为每个 array 使用一个非常短的标头，并且如您所料，针对数字计算进行了高度优化。

注意x**2 将x 的每一项平方而不是x[i]*x[i]，而np.dot(a, x) 仅在一个命令中执行np.dot(x, a[i])（其中x 和a 切换位置以适应所需的尺寸）。

您可以参考this 很棒的电子书，它会更详细地解释这一点。

干杯