n维函数最优化算法代码实现

由于n维函数最优化算法有很多种，这里只列举其中几种常见的算法的代码实现。

1. 梯度下降法

def gradient_descent(f, df, x0, learning_rate, max_iters):
    x = x0
    for i in range(max_iters):
        grad = df(x)
        x = x - learning_rate * grad
        if abs(grad).max() < 1e-6:
            break
    return x, f(x)

其中，f表示目标函数，df表示目标函数的梯度函数，x0表示初始点，learning_rate表示学习率，max_iters表示最大迭代次数。

2. 牛顿法

def newton_method(f, df, ddf, x0, max_iters):
    x = x0
    for i in range(max_iters):
        grad = df(x)
        hess = ddf(x)
        x = x - np.linalg.inv(hess) @ grad
        if abs(grad).max() < 1e-6:
            break
    return x, f(x)

其中，f表示目标函数，df表示目标函数的梯度函数，ddf表示目标函数的Hessian矩阵函数，x0表示初始点，max_iters表示最大迭代次数。

3. 共轭梯度法

def conjugate_gradient(f, df, x0, max_iters):
    x = x0
    r = df(x)
    p = -r
    for i in range(max_iters):
        alpha = np.sum(r ** 2) / np.sum(p @ df(x + p))
        x = x + alpha * p
        new_r = r + alpha * df(x + p)
        beta = np.sum(new_r ** 2) / np.sum(r ** 2)
        p = -new_r + beta * p
        r = new_r
        if abs(r).max() < 1e-6:
            break
    return x, f(x)

其中，f表示目标函数，df表示目标函数的梯度函数，x0表示初始点，max_iters表示最大迭代次数。

4. BFGS算法

def bfgs(f, df, x0, max_iters):
    x = x0
    H = np.eye(len(x))
    for i in range(max_iters):
        grad = df(x)
        d = -H @ grad
        alpha = line_search(f, df, x, d)
        s = alpha * d
        new_x = x + s
        new_grad = df(new_x)
        y = new_grad - grad
        rho = 1 / (y @ s)
        H = (np.eye(len(x)) - rho * np.outer(s, y)) @ H @ (np.eye(len(x)) - rho * np.outer(y, s)) + rho * np.outer(s, s)
        x = new_x
        if abs(grad).max() < 1e-6:
            break
    return x, f(x)

其中，f表示目标函数，df表示目标函数的梯度函数，x0表示初始点，max_iters表示最大迭代次数。

5. L-BFGS算法

def lbfgs(f, df, x0, max_iters, m=5):
    x = x0
    s_list = []
    y_list = []
    rho_list = []
    for i in range(max_iters):
        grad = df(x)
        if i == 0:
            d = -grad
        else:
            q = grad
            alpha_list = []
            for j in range(len(s_list) - 1, -1, -1):
                alpha = rho_list[j] * s_list[j] @ q
                alpha_list.append(alpha)
                q = q - alpha * y_list[j]
            r = H0 @ q
            for j in range(len(s_list)):
                beta = rho_list[j] * y_list[j] @ r
                r = r + s_list[j] * (alpha_list[-j - 1] - beta)
            d = -r
        alpha = line_search(f, df, x, d)
        s = alpha * d
        new_x = x + s
        new_grad = df(new_x)
        y = new_grad - grad
        rho = 1 / (y @ s)
        if len(s_list) == m:
            s_list.pop(0)
            y_list.pop(0)
            rho_list.pop(0)
        s_list.append(s)
        y_list.append(y)
        rho_list.append(rho)
        H0 = np.eye(len(x))
        q = new_grad
        alpha_list = []
        for j in range(len(s_list) - 1, -1, -1):
            alpha = rho_list[j] * y_list[j] @ q
            alpha_list.append(alpha)
            q = q - alpha * s_list[j]
        r = H0 @ q
        for j in range(len(s_list)):
            beta = rho_list[j] * s_list[j] @ r
            r = r + y_list[j] * (alpha_list[-j - 1] - beta)
        H0 = (s @ y) / (y @ y) * np.eye(len(x)) + (r @ s.T + s @ r.T) / (s @ y) - (y @ r.T + r @ y.T) / (y @ y)
        x = new_x
        if abs(grad).max() < 1e-6:
            break
    return x, f(x)

其中，f表示目标函数，df表示目标函数的梯度函数，x0表示初始点，max_iters表示最大迭代次数，m表示存储的历史信息的个数。