最优化方法Python计算：n元函数梯度与Hesse阵的数值计算

对于连续可微的$n$元函数$f(\mathbf{x})$，$\mathbf{x}=\left(\begin{array}{c}{x}_1\\ ⋮\\ x_i\\ ⋮\\ x_n\end{array}\right)\in {\text{R}}^n$，其梯度是由$f(\mathbf{x})$对$\mathbf{x}$的各分量$x_i$的偏导数$\partial f/\partial x_i$（$i=1,2\cdots ,n$）构成。事实上，偏导数$\partial f/\partial x_i$是将$f(\mathbf{x})$中的$x_i$视为唯一变量，其他均为常量的一元函数的导数。所以可利用中心差分（详见博文《最优化方法Python计算：一元函数导数的数值计算》）计算偏导数
$$\frac{\partial f({\mathbf{x}}_0)}{\partial x_i}\approx \frac{f({\mathbf{x}}_1+{\mathbf{e}}_i\Delta x/2)-f({\mathbf{x}}_1-{\mathbf{e}}_i\Delta x/2)}{\Delta x},i=1,2,\cdots ,n.$$
其中，
$$\begin{array}{cc}{\mathbf{e}}_i=\left(\begin{array}{c}0\\ ⋮\\ 1\\ ⋮\\ 0\end{array}\right)& i\end{array}$$
即第$i$个分量为1其他分量均为0的向量，称为第$i$个基本单位向量，$i=1,2,\cdots ,n$。利用上式，可定义如下计算$n$元可微函数$f(\mathbf{x})$在点${\mathbf{x}}_1$处梯度$\nabla f({\mathbf{x}}_1)$的Python函数。

import numpy as np						#导入numpy
def grad(f,x1):							#f为函数f(x)，x1为自变量值
   dx=1e-7								#设置Δx=0.0000001
   n=x1.size							#x的维度
   e=np.eye(n)							#n个基本单位向量
   g0=np.zeros(n)						#梯度初始化为0向量
   for i in range(n):					#逐一计算偏导数
      g0[i]=(f(x1+e[i]*dx/2)-f(x1-e[i]*dx/2))/dx
return g0

程序的第2～9行定义计算n元可微函数梯度的Python函数grad，参数f表示函数$f(\mathbf{x})$，x1表示自变量（$n$维向量）。第3行设置自变量增量$\Delta x_i=\frac{1}{10^7}$为dx。第4行读取自变量所含分量个数n。第5行调用numpy的eye函数设置n个单位向量${\mathbf{e}}_1,{\mathbf{e}}_2,\cdots ,{\mathbf{e}}_n$。第6行调用numpy的zeros函数将梯度g0初始化为0向量。第7～8行的for语句按前式逐一计算偏导数$\frac{\partial f({\mathbf{x}}_1)}{\partial x_i}$，$i=1,2,\cdots ,n$。
相仿地，利用对广义导数的中心差分（详见博文《最优化方法Python计算：一元函数导数的数值计算》），可得二阶可微函数$f(\mathbf{x})$的二阶偏导数的计算公式
$$\begin{gathered} \frac{{\partial }^{2}f({\mathbf{x}}_1)}{\partial x_i\partial x_j}\approx [f\left({\mathbf{x}}_1+({\mathbf{e}}_i+{\mathbf{e}}_j)\frac{\Delta x}{2}\right)+f\left({\mathbf{x}}_1-({\mathbf{e}}_i+{\mathbf{e}}_j)\frac{\Delta x}{2}\right) \\ -f\left({\mathbf{x}}_1+({\mathbf{e}}_i-{\mathbf{e}}_j)\frac{\Delta x}{2}\right)-f\left({\mathbf{x}}_1-({\mathbf{e}}_i-{\mathbf{e}}_j)\frac{\Delta x}{2}\right)]/\Delta x^2,1\le i,j\le n \end{gathered}$$
利用上式，定义计算$n$元二阶可微函数$f(\mathbf{x})$在点${\mathbf{x}}_1$处Hesse阵的Python函数如下

import numpy as np							#导入numpy
def hess(f,x1):								#f表示函数，x1表示自变量
   dx=7.5e-8								#设置自变量增量Δx
   n=x1.size								#向量维数
   e=np.eye(n)								#基本单位向量
   H0=np.zeros((n,n))						#Hesse阵初始化为0阵
   for i in range(n):						#逐一计算二阶偏导
      for j in range(n):
         H0[i,j]=(f(x1+(e[i]+e[j])*dx/2)+f(x1-(e[i]+e[j])*dx/2)-
	          f(x1+(e[i]-e[j])*dx/2)-f(x1-(e[i]-e[j])*dx/2))/dx**2
   return H0

程序的第2～1行定义$n$元二阶可微函数$f(\mathbf{x})$在点${\mathbf{x}}_1$处Hesse矩阵$\mathbf{H}({\mathbf{x}}_1)$的计算函数hess。第3行设置自变量增量$\Delta x=7.5\cdot 10^{-8}$为dx。第4行读取作为自变量的${\mathbf{x}}_1$的维数n。第5行调用numpy的eye函数设置n个单位向量${\mathbf{e}}_1,{\mathbf{e}}_2,\cdots ,{\mathbf{e}}_n$。第6行调用numpy的zeros函数将Hesse阵H0初始化为0阵。第7～10行的两重for循环按上式逐一计算$f(\mathbf{x})$在${\mathbf{x}}_1$的二阶偏导数$\frac{{\partial }^{2}f({\mathbf{x}}_1)}{\partial x_i\partial x_j}$并置于H[i,j]。
例1:计算函数$f(x,y)=x_1^4+x_2^4-4x_1^2{x}_2^2,(x_1,x_2)\in {\text{R}}^2$在${\mathbf{x}}_1=\left(\begin{array}{c}\frac{1}{2}\\ \frac{1}{2}\end{array}\right)$的梯度和Hesse矩阵。
解：下列代码调用上列程序定义的函数grad和hess计算$\nabla f({\mathbf{x}}_1)$和$\mathbf{H}({\mathbf{x}}_1)$。

import numpy as np											#导入numpy
from utility import grad,hess								#导入grad,hess
f=lambda x:(x[0]**4)+(x[1]**4)-4*(x[0]**2)*(x[1]**2)		#设置函数
x1=np.array([0.5,0.5])										#设置自变量x1
print(grad(f,x1))											#计算梯度
print(hess(f,x1))											#计算Hesse阵

利用代码内部的注释信息不难理解上述程序。运行程序，输出

[-0.5 -0.5]
[[ 1.00166788 -4.00667154]
[-4.00667154  1.00166788]]

输出的第1行，表示$f(\mathbf{x})$在${\mathbf{x}}_1=\left(\begin{array}{c}\frac{1}{2}\\ \frac{1}{2}\end{array}\right)$的梯度$\nabla f({\mathbf{x}}_1)=\left(\begin{array}{c}-\frac{1}{2}\\ -\frac{1}{2}\end{array}\right)$的近似，输出的第2～3行表示Hesse矩阵$\mathbf{H}({\mathbf{x}}_1)=\left(\begin{array}{cc}1& -4\\ -4& 1\end{array}\right)$的近似。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！