找到这个函数的 Hessian 矩阵

问题描述

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

我需要一种方法来计算关于 (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，但我对如何进行感到困惑。

解决方法

PyTorch 具有 GPU 优化的 hessian 操作：

import torch

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

您可以使用常规的 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])（其中 { {1}} 和 x 切换位置以适合所需尺寸）。

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

干杯

convex-optimization numpy python sympy