用python 编写一个Newton-CG算法案例

在Python中，我们可以编写一个简单的牛顿- conjugate gradient (CG)算法用于求解二次函数的最小值问题。牛顿-CG算法是一种迭代优化方法，特别适用于大型稀疏矩阵的问题。这里我们先假设有一个函数f(x)，它的梯度和Hessian矩阵可以方便地计算。

import numpy as np

def f_and_gradient(x):
    # 假设f(x) = 0.5 * x.T @ A @ x + b.T @ x
    A = np.random.rand(10, 10)  # 随机生成一个对称正定矩阵
    b = np.random.randn(10)     # 随机生成一个向量
    return 0.5 * np.dot(np.dot(x, A), x) + np.dot(b, x), np.dot(A, x) + b

def conjugate_gradients(A, b, x0=np.zeros(10), maxiter=100, tol=1e-6):
    r = b - np.dot(A, x0)
    p = r
    rsold = np.linalg.norm(r)
    for _ in range(maxiter):
        Ap = np.dot(A, p)
        alpha = rsold / np.dot(p, Ap)
        x_new = x0 + alpha * p
        r_new = r - alpha * Ap
        rsnew = np.linalg.norm(r_new)
        if rsnew < tol:
            break
        beta = rsnew / rsold
        p = r_new + beta * p
        r = r_new
        x0 = x_new
        rsold = rsnew
    return x_new

A, b = f_and_gradient(np.zeros(10))
solution = conjugate_gradients(A, b)
print(f"Solution found after {len(solution)} iterations: {solution}")