共轭梯度法的二元函数形式，python实现

共轭梯度法，找一个二元函数的最小值，在整个二维平面中。需要一个起始点。目前只能做到二元函数。

import sympy as sy
import numpy as np


def differ(f):
    fx = sy.diff(f, x)
    fy = sy.diff(f, y)
    fxx = sy.diff(fx, x)
    fyx = sy.diff(fx, y)
    fxy = sy.diff(fy, x)
    fyy = sy.diff(fy, y)

    return fx, fy, fxx, fyx, fxy, fyy


def FV(x0, f):
    result = float(f.evalf(subs={x: x0[0], y: x0[1]}))

    return result


def Con_Gra(x, f, eps=10e-7):
    times = 0
    gx, gy, fxx, fyx, fxy, fyy = differ(f)
    gx_v = float(FV(x, gx))
    gy_v = float(FV(x, gy))
    g = np.array([gx_v, gy_v]).T
    d = -g
    while True:
        Q = np.array([[FV(x, fxx), FV(x, fxy)], [FV(x, fyx), FV(x, fyy)]])
        alpha = -np.matmul(g.T, d)/np.matmul(d.T, np.matmul(Q, d))
        x = x + alpha * d
        print(times, g, d, alpha, x)
        times += 1
        gap = np.matmul(g.T, g)
        gx_v = float(FV(x, gx))
        gy_v = float(FV(x, gy))
        g = np.array([gx_v, gy_v]).T
        beta = np.matmul(g.T, np.matmul(Q, d))/np.matmul(d.T, np.matmul(Q, d))
        d = -g+beta * d

        if gap < eps:  # 精度满足要求，结束
            break

        elif times > 1000:  # 如果迭代超过10000次，结束
            print("10000次迭代仍不收敛")
            break

    print(x)
    print(times)


if __name__ == '__main__':  # 当模块被直接运行时，以下代码块将被运行，当模块是被导入时，代码块不被运行。

    x, y = sy.symbols("x y")
    f = x*x+2*y*y-x*y+3*x-2*y  # 公式
    x0 = np.array([10, 10]).T  # x=0 y=0 为初值
    Con_Gra(x0, f)