一阶偏导数

一阶偏导数

思路

自变量： x0,...,xi+h,...,xn

带入理查德外推法公式，求对xi的偏导数

计算gradi: grad[grad0,...,gradi,...,null]

实现

import numpy as np

def computeGrad(pointFun, x, h):

temp = np.array(x)

grad = np.array(x)

for i in range(len(x)):

temp[i] += 0.5 * h

grad[i] = 4 * pointFun(temp) / (3 * h)

temp[i] -= h

grad[i] -= 4 * pointFun(temp) / (3 * h)

temp[i] += 3 * h / 2

grad[i] -= pointFun(temp) / (6 * h)

temp[i] -= 2 * h

grad[i] += pointFun(temp) / (6 * h)

temp[i] = x[i]

return grad

# def fun(x):

# return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2

#

#

# grad = computeGrad(fun, np.array([-1.2, 1.0]), 1e-3)

#

# print "grad:", grad

调用

def fun(x):

return 100 * (x[1] - x[0] ** 2) ** 2 + (1 - x[0]) ** 2

# 手动计算偏导数：

def gradient(x):

return np.array([-400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0]), 200 * (x[1] - x[0] ** 2)])

# 理查德外推法计算偏导数

def gradient(x):

grad = computeGrad(fun, x, 1e-3)

return grad

二阶偏导数

二阶偏导

思路

计算主对角线hessian矩阵元素fxixi''

计算下三角阵hessian矩阵元素fxixj''

实现

import numpy as np

def computeHes(pointFun, x, epsilon):

temp = np.array(x)

Hess = np.array(x)

for index in range(len(x) - 1):

Hess = np.vstack((Hess, x))

for i in range(len(x)):

for j in range(len(x)):

Hess[i, j] = 0

# compute fxx''

for i in range(len(x)):

temp[i] += epsilon / 2

Hess[i, i] = 16 * pointFun(temp)

temp[i] -= epsilon

Hess[i, i] += 16 * pointFun(temp)

temp[i] += 3 * epsilon / 2

Hess[i, i] -= pointFun(temp)

temp[i] -= 2 * epsilon

Hess[i, i] -= pointFun(temp)

Hess[i, i] -= 30 * pointFun(x)

Hess[i, i] /= 3 * epsilon * epsilon

temp[i] = x[i]

# compute Zij

for i in range(len(x)):

for j in range(len(x)):

if j > i:

temp[i] += epsilon / 2

temp[j] += epsilon / 2

Hess[i, j] = 16 * pointFun(temp)

temp[i] -= epsilon

temp[j] -= epsilon

Hess[i, j] += 16 * pointFun(temp)

temp[i] += 3 * epsilon / 2

temp[j] += 3 * epsilon / 2

Hess[i, j] -= pointFun(temp)

temp[i] -= 2 * epsilon

temp[j] -= 2 * epsilon

Hess[i, j] -= pointFun(temp)

Hess[i, j] -= 30 * pointFun(x)

Hess[i, j] /= (3 * epsilon * epsilon)

Hess[i, j] -= Hess[i, i]

Hess[i, j] -= Hess[j, j]

Hess[i, j] /= 2

temp[i] = x[i]

temp[j] = x[j]

# compute fij'' use fii'' & Zij

for i in range(len(x)):

for j in range(len(x)):

if j > i:

Hess[j, i] -= Hess[i, j]

return Hess

测试

def fun(x):

return x[0] * np.exp(x[0] - x[1] * x[1])

point = np.array([1.0, 1.0])

Hessian = computeHes(fun, point, 1e-3)

print "Hessian:", Hessian

测试结果

调用

def hessianMatrix(x):

# Hes = np.array([[-400 * (x[1] - 3 * x[0] ** 2) + 2, -400 * x[0]], [-400 * x[0], 200]])

Hes = computeHes(fun, x, 1e-3)

return Hes