理查德外推法计算偏导数近似值-python实现

一阶偏导数

思路

自变量： 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