共軛梯度法（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()}')

運行結果