共轭梯度法——正定函数和一般函数python代码

方法不赘述，代码如下：

正定函数问题

import numpy as np


def grad(A, x):
    return A @ x


def func_PR(k, k1):
    if k.T @ k1 != 0:
        return k.T @ k / k1.T @ k1
    else:
        return 0


def conjugate_gradient(A, x0):
    """
    parameter A: 二次函数矩阵
    parameter x0: 初值
    """
    count = 0  # 计数器
    print("起始点为：", x0)
    eps = 1e-100  # 允许误差
    gs = []  # 储存梯度值
    while True:
        g = grad(A, x0)  # 梯度
        p = - g  # 搜索方向
        gs.append(g)
        if np.linalg.norm(g) < eps:
            return False
        if count == 0:
            alpha = 0
        else:
            alpha = func_PR(gs[count - 1], gs[count])
            # alpha = func_PR(g, g)
            # alpha = g.T @ A @ p / p.T@ A@p
        print("alpha", alpha)
        p = p + alpha * p  # 更新搜索方向
        beta = -(g.T @ p) / (p.T @ A @ p)
        x0 = x0 + beta * p  # 更新初值
        print('%d次的梯度：' % count, g)
        print('gs为：', gs)

        print('x0的值为： {}'.format(x0))
        count += 1
        print("迭代次数为：", count)


# 测试函数为：f(x) = x1**2 + 4* x2**2   x0=[1,1] 最优值【0,0】
A = np.array([[2, 0], [0, 8]])
x0 = np.array([[1], [1]])

conjugate_gradient(A, x0)

一般函数

from sympy import *


x, y = symbols('x, y')

# 对一般函数求梯度


def get_grad(f, X):
    # 一阶导数
    fx = diff(f, x)
    fy = diff(f, y)
    grad = Matrix([[fx], [fy]])
    return Matrix([[fx.subs([(x, X[0]), (y, X[1])])],
                   [fy.subs([(x, X[0]), (y, X[1])])]])


# PR公式
def func_PR(k, k1):
    if k.dot(k1) != 0:
        return k.dot(k) / k1.dot(k1)
    else:
        return 0


def conjugate_gradient(f, X0):
    count = 0  # 计数器
    print("起始点为：", X0)
    eps = 1e-10  # 允许误差
    
    while True:
        g = get_grad(f, X0)
        p = - g
        # print(g.norm())
        if g.norm() < eps:
            return False
        if count == 0:
            alpha = 0
        else:
            alpha = func_PR(g, g)
        p = p + alpha * p  # 更新搜索方向

        beta = 0.01
        X0 = X0 + beta * p  # 更新初值
     
        print('%d次的梯度：' % count, g)
        print('x0的值为： {}'.format(X0))
        count += 1
        print("迭代次数为：", count)


f = (x+2*y - 7)**2+(2*x+y-5)**2  # 最优解f(1,3)=0
# f = 2*x**2-1.05*x**4+(1/6)*x**6+x*y+y**2 # f(0,0) = 0


X0 = Matrix([[-1], [1]])
conjugate_gradient(f, X0)

note

这里迭代点的参数beta应该用 线搜索 的方法确定，这里直接给出了，不好，有空改进。