求此函数的Hessian矩阵

Question

嗨，我有以下功能：

从1到5000的总和-对数（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

您可以使用常规的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一本很棒的电子书，它将更详细地解释这一点

干杯

Answer 2

PyTorch具有GPU优化的^{}操作：

import torch

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