牛顿迭代法求解多元函数的最小值--以二元函数为例

【python】牛顿迭代法求解多元函数的最小值–以二元函数为例

本文转载链接：[https://blog.csdn.net/qq_45726331/article/details/115804812]

一元函数到多元函数的牛顿迭代法
python代码实现过程
一元函数到多元函数的牛顿迭代法
参考 多元函数的牛顿迭代和高斯牛顿法怎么推导？

python代码实现过程
计算梯度函数
使用了sympy库的diff函数计算导数，再代入具体数值进行计算，返回X处的梯度grad

求解梯度值

def get_grad(f, X):
    # 计算一阶导数
    f1 = diff(f, x1)
    f2 = diff(f, x2)
    # 代入具体数值计算
    grad = np.array([[f1.subs([(x1, X[0]), (x2, X[1])])],
                   [f2.subs([(x1, X[0]), (x2, X[1])])]])
    return grad

计算Hession矩阵函数
同样使用diff函数计算二次偏导，组成Hession矩阵（里面的括号、中括号要盯对清楚）

求解Hession矩阵

def get_hess(f, X):
    # 计算二次偏导
    f1 = diff(f, x1)
    f2 = diff(f, x2)
    f11 = diff(f,x1,2)
    f22 = diff(f,x2,2)
    f12 = diff(f1,x2)
    f21 = diff(f2,x1)
    # 计算具体数值计算
    hess = np.array([[f11.subs([(x1,X[0]), (x2,X[1])]),
                        f12.subs([(x1,X[0]), (x2,X[1])])],

                        [f21.subs([(x1,X[0]), (x2,X[1])]),
                        f22.subs([(x1,X[0]), (x2,X[1])])]])
    # 转换数值类型为了后续求逆矩阵
    hess = np.array(hess, dtype = 'float')
    return hess

牛顿迭代过程
迭代公式里使用的是Hession矩阵的逆，所以还要求个逆

牛顿迭代

def newton_iter(X0, err):
    # 记录迭代次数
    count = 0
    X1 = np.array([[0],[0]])
    while True:
        # 得到两次迭代结果差值
        X2 = X0 - X1
        if sqrt(X2[0]**2 + X2[1]**2) <= err:
            break
        else:
            hess = get_hess(f, X0)
            # 得到Hession矩阵的逆
            hess_inv = np.linalg.inv(hess)
            grad = get_grad(f, X0)
            X1 = X0 #保存上次迭代结果
            # 迭代公式
            X0 = X1 - np.dot(hess_inv, grad)
            count += 1
            print('第',count,'次迭代:','x1=',X0[0],'x2=',X0[1])
    print('迭代次数为:',count)
    print('迭代结果为',X0)

全部代码如下

from sympy import *
import numpy as np

定义符号

x1,x2 = symbols('x1, x2')

定义所求函数

f = 60 - 10*x1 - 4*x2 + x1**2 + x2**2 - x1*x2

求解梯度值

def get_grad(f, X):
    # 计算一阶导数
    f1 = diff(f, x1)
    f2 = diff(f, x2)
    # 代入具体数值计算
    grad = np.array([[f1.subs([(x1, X[0]), (x2, X[1])])],
                   [f2.subs([(x1, X[0]), (x2, X[1])])]])
    return grad

求解Hession矩阵

def get_hess(f, X):
    # 计算二次偏导
    f1 = diff(f, x1)
    f2 = diff(f, x2)
    f11 = diff(f,x1,2)
    f22 = diff(f,x2,2)
    f12 = diff(f1,x2)
    f21 = diff(f2,x1)
    # 计算具体数值计算
    hess = np.array([[f11.subs([(x1,X[0]), (x2,X[1])]),
                        f12.subs([(x1,X[0]), (x2,X[1])])],

                        [f21.subs([(x1,X[0]), (x2,X[1])]),
                        f22.subs([(x1,X[0]), (x2,X[1])])]])
    # 转换数值类型为了后续求逆矩阵
    hess = np.array(hess, dtype = 'float')
    return hess

牛顿迭代

def newton_iter(X0, err):
    # 记录迭代次数
    count = 0
    X1 = np.array([[0],[0]])
    while True:
        # 得到两次迭代结果差值
        X2 = X0 - X1
        if sqrt(X2[0]**2 + X2[1]**2) <= err:
            break
        else:
            hess = get_hess(f, X0)
            # 得到Hession矩阵的逆
            hess_inv = np.linalg.inv(hess)
            grad = get_grad(f, X0)
            X1 = X0 #保存上次迭代结果
            # 迭代公式
            X0 = X1 - np.dot(hess_inv, grad)
            count += 1
            print('第',count,'次迭代:','x1=',X0[0],'x2=',X0[1])
    print('迭代次数为:',count)
    print('迭代结果为',X0)

设置初始点

X0 = np.array([[1],[1]])
err = 1e-10
newton_iter(X0, err)

原文链接：https://blog.csdn.net/qq_45726331/article/details/115804812