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

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

import numpy as np						#导入numpy
def grad(f,x1):							#f为函数f(x)，x1为自变量值dx=1e-7								#设置Δx=0.0000001n=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行设置自变量增量 Δ x i = 1 1 0 7 \\Delta x_i=\\frac{1}{10^7} Δxi​=1071​为dx。第4行读取自变量所含分量个数n。第5行调用numpy的eye函数设置n个单位向量 e 1 , e 2 , ⋯ , e n \\boldsymbol{e}_1,\\boldsymbol{e}_2,\\cdots,\\boldsymbol{e}_n e1​,e2​,⋯,en​。第6行调用numpy的zeros函数将梯度g0初始化为0向量。第7～8行的for语句按前式逐一计算偏导数 ∂ f ( x 1 ) ∂ x i \\frac{\\partial f(\\boldsymbol{x}_1)}{\\partial x_i} ∂xi​∂f(x1​)​，$i=1,2,\cdots ,n$。
相仿地，利用对广义导数的中心差分（详见博文《最优化方法Python计算：一元函数导数的数值计算》），可得二阶可微函数$f(\mathbf{x})$的二阶偏导数的计算公式
∂ 2 f ( x 1 ) ∂ x i ∂ x j ≈ [ f ( x 1 + ( e i + e j ) Δ x 2 ) + f ( x 1 − ( e i + e j ) Δ x 2 ) − f ( x 1 + ( e i − e j ) Δ x 2 ) − f ( x 1 − ( e i − e j ) Δ x 2 ) ] / Δ x 2 , 1 ≤ i , j ≤ n \\frac{\\partial^2f(\\boldsymbol{x}_1)}{\\partial x_i\\partial x_j}\\approx\\left[ f\\left(\\boldsymbol{x}_1+(\\boldsymbol{e}_i+\\boldsymbol{e}_j)\\frac{\\Delta x}{2}\\right)+f\\left(\\boldsymbol{x}_1-(\\boldsymbol{e}_i+\\boldsymbol{e}_j)\\frac{\\Delta x}{2}\\right)\\right.\\\\ \\quad-\\left.f\\left(\\boldsymbol{x}_1+(\\boldsymbol{e}_i-\\boldsymbol{e}_j)\\frac{\\Delta x}{2}\\right)-f\\left(\\boldsymbol{x}_1-(\\boldsymbol{e}_i-\\boldsymbol{e}_j)\\frac{\\Delta x}{2}\\right)\\right]\\big/\\Delta x^2,1\\leq i,j\\leq n ∂xi​∂xj​∂2f(x1​)​≈[f(x1​+(ei​+ej​)2Δx​)+f(x1​−(ei​+ej​)2Δx​)−f(x1​+(ei​−ej​)2Δx​)−f(x1​−(ei​−ej​)2Δx​)]/Δx2,1≤i,j≤n
利用上式，定义计算$n$元二阶可微函数$f(\mathbf{x})$在点 x 1 \\boldsymbol{x}_1 x1​处Hesse阵的Python函数如下

import numpy as np							#导入numpy
def hess(f,x1):								#f表示函数，x1表示自变量dx=7.5e-8								#设置自变量增量Δxn=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**2return H0

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

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})$在 x 1 = ( 1 2 1 2 ) \\boldsymbol{x}_1=\\begin{pmatrix} \\frac{1}{2}\\\\\\frac{1}{2} \\end{pmatrix} x1​=(21​21​​)的梯度 ∇ f ( x 1 ) = ( − 1 2 − 1 2 ) \\nabla f(\\boldsymbol{x}_1)=\\begin{pmatrix} -\\frac{1}{2}\\\\-\\frac{1}{2} \\end{pmatrix} ∇f(x1​)=(−21​−21​​)的近似，输出的第2～3行表示Hesse矩阵 H ( x 1 ) = ( 1 − 4 − 4 1 ) \\boldsymbol{H}(\\boldsymbol{x}_1)=\\begin{pmatrix} 1&-4\\\\ -4&1 \\end{pmatrix} H(x1​)=(1−4​−41​)的近似。

岩板材料