多变量函数极值-牛顿法

简介

之前我们已经采用过牛顿法进行一维搜索，这次我们将牛顿法推广到多维。利用牛顿法来求解多变量极值问题。跟一维问题类似，在局部采用二次函数$\phi (\mathbf{x})$来近似替代目标函数$f(\mathbf{x})$，然后用$\phi (\mathbf{x})$的极小点作为$f(\mathbf{x})$的近似极小点。设${\mathbf{x}}^k$为$f(x)$的一个近似极小点，将$f(\mathbf{x})$在${\mathbf{x}}^k$点利用Taylor公式展开，并略去高于二次的项，则可得:
$$f(\mathbf{x})\approx \phi (\mathbf{x})=f({\mathbf{x}}^k)+(▽f({\mathbf{x}}^k){)}^T(\mathbf{x}-{\mathbf{x}}^k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^k{)}^T{▽}^{2}f({\mathbf{x}}^k)(\mathbf{x}-{\mathbf{x}}^k)$$
故$f(\mathbf{x})$的极小点迭代公式如下：
$${\mathbf{x}}^{k+1}={\mathbf{x}}^k-(▽f({\mathbf{x}}^k){)}^{-1}{▽}^{2}f({\mathbf{x}}^k)$$
牛顿法优点是收敛速度很快，对于二次函数只需要一次迭代就能找到极小值点，缺点是$f(\mathbf{x})$必须是二阶连续可微的，并且初始点不能离极小点太远否则迭代可能不收敛。

优化问题

求解$f(x_1,x_2)=x_1^2+2x_2^2-4x_1-2x_1{x}_2$的极小值点，设置初始点$x_0=[0,0]$，收敛精度$\epsilon =1e^{-3}$，计算程序如下：

import numpy as np


def cal_fun_f(x1, x2):
    """
    函数值计算函数
    """
    return x1**2+2*x2**2-4*x1-2*x1*x2

def partial_derivative_x1(x1, x2):
    """
    求解函数对x1的偏导数
    :param x1:
    :param x2:
    """
    return 2*x1-4-2*x2

def partial_derivative_x2(x1, x2):
    """
    求解函数对x2的偏导数
    :param x1:
    :param x2:
    """
    return 4*x2-2*x1

def partial_derivative_x1x2(x1, x2):
    """
    求解对x1和x2的二阶偏导数
    :param x1:
    :param x2:
    """
    return -2

def partial_derivative_x1x1(x1, x2):
    """
    对x1的二阶偏导数
    :param x1:
    :param x2:
    """
    return 2

def partial_derivative_x2x2(x1, x2):
    """
    对x2的二阶偏导数
    :param x1:
    :param x2:
    """
    return 4

def get_derivative_matrix1(x1, x2):
    """
    计算一阶偏导数矩阵
    :param x1:
    :param x2:
    """
    derivative_matrix1 = np.full((2, 1), np.nan, dtype=float)
    derivative_matrix1[0, 0] = partial_derivative_x1(x1, x2)
    derivative_matrix1[1, 0] = partial_derivative_x2(x1, x2)
    return derivative_matrix1

def get_derivative_matrix2(x1, x2):
    """
    计算二阶偏导数矩阵
    :param x1:
    :param x2:
    """
    derivative_matrix2 = np.full((2, 2), np.nan, dtype=float)
    derivative_matrix2[0, 0] = partial_derivative_x1x1(x1, x2)
    derivative_matrix2[0, 1] = partial_derivative_x1x2(x1, x2)
    derivative_matrix2[1, 0] = partial_derivative_x1x2(x1, x2)
    derivative_matrix2[1, 1] = partial_derivative_x2x2(x1, x2)
    return derivative_matrix2

def newton_method(x, gap):
    i = 0
    while True:
        derivative_matrix1 = get_derivative_matrix1(x[0, 0], x[1, 0])
        derivative_matrix2 = get_derivative_matrix2(x[0, 0], x[1, 0])
        x = x-np.linalg.inv(derivative_matrix2)@derivative_matrix1
        i += 1
        print(f"第{i}次迭代极小值点", (x[0, 0], x[1, 0]))
        if np.sum(np.abs(get_derivative_matrix1(x[0, 0], x[1, 0])), axis=0) < gap:
            break

if __name__ == "__main__":
    newton_method(np.array([0, 0]).reshape(-1, 1), 1e-3)

计算结果如下：

第1次迭代极小值点 (4.0, 2.0)

可以看到仅仅经过一次迭代就找到了极小值点。