共轭梯度法python实现

Cherry_csc

已于 2023-05-11 20:00:02 修改

阅读量636

点赞数 1

文章标签：算法 python

于 2023-05-11 19:57:35 首次发布

本文链接：https://blog.csdn.net/Cherry_csc/article/details/130624196

版权

两种共轭梯度法的算法框架：

Fletcher-Reeves

1. 任取 $X_0$ $\in R^n$ ,令 $k=0$

2. 如果 $||\bigtriangledown f(X_k)|| \leq \varepsilon$ ,停止计算

3. 如果k/n等于0或整数，令 $D_k=-\bigtriangledown f(X_k)$ ,否则令 $D_k=-\bigtriangledown f(X_k)+\alpha D_{k-1}$ ，其中

$\alpha _{k-1}=\frac{||\bigtriangledown f(X_k)||^2}{||\bigtriangledown f(X_{k-1})||^2}$

4. 进行精确搜索获得 $X_{k+1}=X_k+t_kD_k$

5. 用k+1替代k，回到2，继续迭代

Polak-Ribiere

1. 任取 $X_0$ $\in R^n$ ,令 $k=0$

2. 如果 $||\bigtriangledown f(X_k)|| \leq \varepsilon$ ,停止计算

3. 如果k/n等于0或整数，令 $D_k=-\bigtriangledown f(X_k)$ ,否则令 $D_k=-\bigtriangledown f(X_k)+\alpha D_{k-1}$ ，其中

$\alpha _{k-1}=\frac{\bigtriangledown ^Tf(X_k)(\bigtriangledown f(X_k)-\bigtriangledown f(X_{k-1}))}{||\bigtriangledown f(X_{k-1})||^2}$

4. 进行精确搜索获得 $X_{k+1}=X_k+t_kD_k$

5. 用k+1替代k，回到2，继续迭代

实例

考虑无约束优化问题

$minimize (1-x-y)^2+2(y-x^2)^2$

取初始点（0,0），终止条件 $||\bigtriangledown f(x)||_2\leq 10^{-4}$

画出目标函数等值线，给出每种算法求得最优解和最优值，画出不同算法函数值随迭代次数增加变化曲线，以及迭代过程中决策变量在等值线上的变化曲线

Fletcher-Reeves代码：

import numpy as np


def grad_f(X):
    x = X[0]
    y = X[1]
    d = np.array([
        -2 * (1 - x - y) - 8 * x * (y - x * x),
        -2 * (1 - x - y) + 4 * (y - x * x)
    ])
    return d


def one_d_search(D, X):  # 使用Newton法精确搜索
    x = X[0]
    y = X[1]
    dx = D[0]
    dy = D[1]
    t = 0

    def first_deriv(t):
        return -2 * (dx + dy) * (1 - x - y - t * dx - t * dy) + 4 * (
            y + t * dy - (x + t * dx) * (x + t * dx)) * (dy - 2 *
                                                         (x + t * dx) * dx)

    def sec_deriv(t):
        return 2 * (dx + dy) * (dx + dy) + 4 * (dy - 2 * (x + t * dx) * dx) * (
            dy - 2 * (x + t * dx) * dx) - 8 * (y + t * dy - (x + t * dx) *
                                               (x + t * dx)) * dx * dx

    while abs(first_deriv(t)) >= 1e-10:
        t = t - first_deriv(t) / sec_deriv(t)
    return t


def FR_conjugate_gradient(x0, tol=1e-4, max_iter=1000):
    x = x0.copy()
    g = grad_f(x)
    d = -g
    alpha = 0
    it = 0
    while np.linalg.norm(g) > tol and it < max_iter:
        step = one_d_search(d, x)
        x = x + step * d
        print("inter", it, ":", "x", x)
        g_new = grad_f(x)
        alpha = np.dot(g_new, g_new) / np.dot(g, g)
        d = -g_new + alpha * d
        g = g_new
        it += 1
    return x


x0 = np.zeros(2)
print(FR_conjugate_gradient(x0))

运行结果：

Polak-Ribiere代码：

def PR_conjugate_gradient(x0, tol=1e-4, max_iter=1000):
    x = x0.copy()
    g = grad_f(x)
    d = -g
    alpha = 0
    it = 0
    while np.linalg.norm(g) > tol and it < max_iter:
        step = one_d_search(d, x)
        x = x + step * d
        print("inter", it, ":", "x", x)
        g_new = grad_f(x)
        alpha = np.dot(g_new, g_new - g) / np.dot(g, g)
        d = -g_new + alpha * d
        g = g_new
        it += 1
    return x

运行结果:

绘图

import matplotlib.pyplot as plt

def f(x, y):
    return (1 - x - y) * (1 - x - y) + 2 * (y - x * x) * (y - x * x)

# 绘图
dx = 0.01
dy = 0.01
x = np.arange(0.3, 0.8, dx)
y = np.arange(0.1, 0.6, dy)
X, Y = np.meshgrid(x, y)
C = plt.contour(X, Y, f(X, Y), levels=10, colors='black')  # 生成等值线图
plt.contourf(X, Y, f(X, Y), 8)
plt.clabel(C, inline=1, fontsize=10)
plt.colorbar()
# 在图上标注点
points_x = PR[1]
points_y = PR[2]
plt.scatter(points_x, points_y, c='r', marker='o', s=5)
print(PR[0][0], PR[0][1])
plt.scatter(PR[0][0], PR[0][1], c='r', marker='*', s=150)
plt.plot(points_x, points_y, color='blue')
plt.show()

等值线