共轭梯度法法（Python实现）

使用共轭梯度法，分别使用Armijo准则和Wolfe准则来求步长

求解方程

\(f(x_1,x_2)=(x_1^2-2)^4+(x_1-2x_2)^2\)的极小值

import numpy as np


# import tensorflow as tf


def gfun(x):  # 梯度
    # x = tf.Variable(x, dtype=tf.float32)
    # with tf.GradientTape() as tape:
    #     tape.watch(x)
    #     z = fun(x)
    # return tape.gradient(z, x).numpy()  # 这里使用TensorFlow来求梯度，直接手算梯度返回也行
    return np.array([4 * (x[0] - 2) ** 3 + 2 * (x[0] - 2 * x[1]), -4 * (x[0] - 2 * x[1])]).reshape(len(x), 1)


def fun(x):  # 函数
    return (x[0] - 2) ** 4 + (x[0] - 2 * x[1]) ** 2


def frcg_armijo(x0):
    maxk = 5000  # 最大迭代次数
    rho = .6  # Armijo准测参数
    sigma = .4
    k = 0
    epsilon = 1e-4
    n = len(x0)  # 输入的维度
    while k < maxk:  # 最大迭代次数
        g = gfun(x0)  # 计算梯度
        itern = k % n
        if itern == 0:  # 迭代n（维度）次后，重新选取负梯度方向作为搜索方向
            d = -g
        else:
            beta = (g.T @ g) / (g0.T @ g0)  # 计算beta
            d = -g + beta * d0
            gd = g.T @ d
            if gd >= 0:  # 必要条件，要小于0，取负梯度方向
                d = -g
        if np.linalg.norm(g) < epsilon:  # 满足精度则结束循环
            break
        m = 0
        mk = 0
        while m < 20:  # 使用Armijo搜索（非精确线搜索）
            if fun(x0 + rho ** m * d) < fun(x0) + sigma * rho ** m * g.T @ d:
                mk = m
                break
            m += 1
        x0 = x0 + rho ** mk * d
        g0 = g
        d0 = d
        k += 1
    val = fun(x0)
    return x0, val, k


def frcg_wolfe(x0):
    maxk = 5000  # 最大迭代次数
    k = 0
    epsilon = 1e-4
    n = len(x0)  # 输入的维度
    while k < maxk:  # 最大迭代次数
        g = gfun(x0)  # 计算梯度
        itern = k % n
        if itern == 0:  # 迭代n（维度）次后，重新选取负梯度方向作为搜索方向
            d = -g
        else:
            beta = (g.T @ g) / (g0.T @ g0)  # 计算beta
            d = -g + beta * d0
            gd = g.T @ d
            if gd >= 0:  # 必要条件，要小于0，取负梯度方向
                d = -g
        if np.linalg.norm(g) < epsilon:  # 满足精度则结束循环
            break
        rho = 0.4
        sigma = 0.5
        a = 0
        b = np.inf
        alpha = 1
        while True:
            if not ((fun(x0) - fun(x0 + alpha * d)) >= (-rho * alpha * gfun(x0).T @ d)):
                b = alpha
                alpha = (a + alpha) / 2
                continue
            if not (gfun(x0 + alpha * d).T @ d >= sigma * gfun(x0).T @ d):
                a = alpha
                alpha = np.min([2 * alpha, (alpha + b) / 2])
                continue
            break

        x0 = x0 + alpha * d
        g0 = g
        d0 = d
        k += 1
    x = x0
    val = fun(x)
    return x, val, k


if __name__ == '__main__':
    x0 = np.array([[0], [3]])
    x0, val, k = frcg_armijo(x0)  # 使用armijo准则
    print(f'近似最优点：\n{x0}\n迭代次数：{k}\n目标函数值：{val.item()}')
    x0 = np.array([[-1.2], [-1]])
    x0, val, k = frcg_wolfe(x0)  # 使用wolfe准则
    print(f'近似最优点：\n{x0}\n迭代次数：{k}\n目标函数值：{val.item()}')

运行结果