是否可以在 scipy.optimize.minimize 中给出一个具有 2 个输出的函数作为目标函数梯度和约束梯度？

Question

我有一个非常基本的优化问题，它是一个 6 维 Rosenbrock 函数。我对目标函数和约束都使用单纯形梯度，但我想在单个函数内计算梯度，以使其计算成本更便宜。这是因为目标函数的梯度计算还包含约束梯度计算所需的数据。尽管对于这个问题在单个函数中计算它们没有意义，但它对我来说只是一个玩具问题。我无法在这里解释我真正的、复杂的问题，但如果你帮助我解决这个问题，我会将相同的逻辑应用于我的问题。

这是我之前尝试过的，但它不起作用。我还尝试将其转储到 pickle 文件中，然后尝试将其回调，但在这种情况下优化无法正常工作。

import numpy as np
from scipy.optimize import minimize, NonlinearConstraint


def rosenbrock_6d(x):
    print("obj is calculating!")
    result = 0
    for i in range(5):
        result += 100 * (x[i + 1] - x[i] ** 2) ** 2 + (1 - x[i]) ** 2
    
    return result

def nonlinear_constraint(x):
    print("Nonlinear constraint is calculating!")
    return 10*x[0]**2 + 11*x[1]**2 + 12*x[2]**2 + 13*x[3]**2 + 14*x[4]**2 + 15*x[5]**2


def StoSAG(x):
    print("OBj StoSAG is calculating!")
    CU = np.eye(6)
    n_pert=10
    u_norm_perturbed = np.random.multivariate_normal(x, CU, n_pert)
    u_norm_perturbed = np.transpose(u_norm_perturbed)
    J_ = np.zeros(n_pert)
    NL_array_pert=np.zeros((n_pert,1))
    for j in range(n_pert):
        J_[j] = rosenbrock_6d(u_norm_perturbed[:, j])
        NL_array_pert[j] = nonlinear_constraint(u_norm_perturbed[:, j])
    J_=np.reshape(J_, (n_pert, 1))
    j=J_- rosenbrock_6d(x)
    U= u_norm_perturbed-np.reshape(x, (6, 1))
    du=np.linalg.pinv(np.transpose(U))
    grad_f=np.matmul(du, j).ravel()
    c_nl=(NL_array_pert- nonlinear_constraint(x))
    g_nl=np.matmul(du, c_nl)
    grad_g=np.transpose(g_nl)
    return grad_f, grad_g

nlc = NonlinearConstraint(nonlinear_constraint, -np.inf,2,jac=StoSAG()[1])

# Example usage:
# Define a sample 6-dimensional point
x_values = (1, 2, 3, 4, 5, 6)
result = minimize(rosenbrock_6d, x_values, method='SLSQP', jac=StoSAG()[0],constraints=nlc,options={'maxiter':200,'disp':True,'iprint':2,'ftol':1e-2})

# Print the result
print("Optimal solution: x =", result.x)
print("Optimal objective function value:", result.fun)

Answer 1

有一长串 Numpy 使用问题，我不会在这里讨论，因为我认为它们超出了范围，但总的来说 - 向量化更多的东西。

不要使用泡菜。您可以使用 SQLite，但我认为在不了解更多有关您的用例的情况下这有点过分了。我认为完成此操作的唯一实用方法是 LRU 缓存：

import numpy as np
from scipy.optimize import minimize, NonlinearConstraint
from functools import lru_cache, wraps


# https://stackoverflow.com/a/52332109/313768
def np_cache(function):
    @lru_cache()
    def cached_wrapper(hashable_array):
        array = np.array(hashable_array)
        return function(array)

    @wraps(function)
    def wrapper(array):
        return cached_wrapper(tuple(array))

    # copy lru_cache attributes over too
    wrapper.cache_info = cached_wrapper.cache_info
    wrapper.cache_clear = cached_wrapper.cache_clear

    return wrapper


def rosenbrock_6d(x: np.ndarray) -> float:
    result = 0
    for i in range(5):
        result += 100*(x[i + 1] - x[i]**2)**2 + (1 - x[i])**2

    return result


def nonlinear_constraint(x: np.ndarray) -> float:
    return 10*x[0]**2 + 11*x[1]**2 + 12*x[2]**2 + 13*x[3]**2 + 14*x[4]**2 + 15*x[5]**2


@np_cache
def StoSAG(x: np.ndarray) -> tuple[
    np.ndarray,  # objective gradient
    np.ndarray,  # constraint gradient
]:
    print("OBj StoSAG is calculating!", flush=True)

    CU = np.eye(6)
    n_pert = 10
    u_norm_perturbed = np.random.multivariate_normal(x, CU, n_pert)
    u_norm_perturbed = np.transpose(u_norm_perturbed)
    J_ = np.zeros(n_pert)
    NL_array_pert = np.zeros((n_pert, 1))

    for j in range(n_pert):
        J_[j] = rosenbrock_6d(u_norm_perturbed[:, j])
        NL_array_pert[j] = nonlinear_constraint(u_norm_perturbed[:, j])

    J_ = np.reshape(J_, (n_pert, 1))
    j = J_ - rosenbrock_6d(x)
    U = u_norm_perturbed - np.reshape(x, (6, 1))
    du = np.linalg.pinv(np.transpose(U))
    grad_f = np.matmul(du, j).ravel()
    c_nl = (NL_array_pert - nonlinear_constraint(x))
    g_nl = np.matmul(du, c_nl)
    grad_g = np.transpose(g_nl)

    return grad_f, grad_g


def StoSAG_objective_grad(x: np.ndarray) -> np.ndarray:
    print('Objective gradient is calculating!', flush=True)
    return StoSAG(x)[0]


def StoSAG_constraint_grad(x: np.ndarray) -> np.ndarray:
    print('Constraint gradient is calculating!', flush=True)
    return StoSAG(x)[1]


def main() -> None:
    nlc = NonlinearConstraint(
        fun=nonlinear_constraint, lb=-np.inf, ub=2,
        jac=StoSAG_constraint_grad,
    )

    # Example usage:
    # Define a sample 6-dimensional point
    x_values = (1, 2, 3, 4, 5, 6)
    result = minimize(
        fun=rosenbrock_6d,
        x0=x_values,
        method='SLSQP',
        jac=StoSAG_objective_grad,
        constraints=nlc,
        options={
            'maxiter': 200,
            'ftol': 1e-2,
        },
    )

    print("Optimal solution: x =", result.x)
    print("Optimal objective function value:", result.fun)


if __name__ == '__main__':
    main()

Constraint gradient is calculating!
OBj StoSAG is calculating!
Objective gradient is calculating!
...