使用scipy.optimize.minimize优化数组输入函数遇错求助

我希望用scipy.optimize.minimize优化含数组类型变量输入的函数，具体场景如下：

先生成含噪声的余弦叠加信号：

import numpy as np
time = np.linspace(0, 1, 501)
data = np.cos(2 * np.pi * 4 * time) + np.cos(2 * np.pi * 9 * time) + np.cos(2 * np.pi * 20 * time)
noise = np.sqrt(1 / 25) * np.random.randn(501)
signal = data + noise 

定义待优化的余弦叠加函数cos_sum：

def cos_sum(x, P):
    assert isinstance(P, np.ndarray)
    assert P.shape[0] == P.shape[1]

    sums = []
    for param in P:
        a, b, c = param
        sums.append(a * np.cos(b * (x - c)))
    sums = np.array(sums)

    return np.sum(sums, axis=0) 

创建残差函数resid用于优化：

def resid(params, x):
    assert isinstance(params, np.ndarray)
    fit = cos_sum(x, params)
    residual = np.sqrt(np.mean(np.abs(fit - signal)) ** 2)
    return residual

构造初始猜测参数：

guess_A = np.random.normal(1, .2, size=3)
guess_B = 2 * np.pi * np.array([4, 9, 20], dtype=float)
guess_C = np.random.normal(0, .2, size=3)
guess = np.array([guess_A, guess_B, guess_C]).T

运行优化代码时出错：

from scipy.optimize import minimize
optimization = minimize(resid, guess, args=(time))

报错信息：

Traceback (most recent call last):
  File "/Users/nickeisenberg/GitRepos/Python_Misc/Misc/minimize_curvefit_vector_variables.py", line 70, in <module>
    optimization = minimize(resid, guess, args=(time))
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_minimize.py", line 676, in minimize
    res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 1296, in _minimize_bfgs
    sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_optimize.py", line 263, in _prepare_scalar_function
    sf = ScalarFunction(fun, x0, args, grad, hess,
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 158, in __init__
    self._update_fun()
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 251, in _update_fun
    self._update_fun_impl()
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 155, in update_fun
    self.f = fun_wrapped(self.x)
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py", line 137, in fun_wrapped
    fx = fun(np.copy(x), *args)
  File "/Users/nickeisenberg/GitRepos/Python_Misc/Misc/minimize_curvefit_vector_variables.py", line 53, in resid
    fit = cos_sum(x, params)
  File "/Users/nickeisenberg/GitRepos/Python_Misc/Misc/minimize_curvefit_vector_variables.py", line 30, in cos_sum
    assert P.shape[0] == P.shape[1]
IndexError: tuple index out of range

请问这种优化需求是否可行？该如何解决当前错误？

可行性说明 #

这种优化需求完全可行，本质是用非线性最小二乘法拟合多分量余弦信号，scipy.optimize.minimize可以胜任这类任务，只要调整参数传递的格式符合函数要求即可。

错误原因及解决方法 #

1. 核心错误原因 #

scipy.optimize.minimize要求待优化的参数必须是一维数组，但你传入的guess是二维数组（形状(3,3)）。优化过程中，minimize会自动将二维数组展平为一维，但你的resid和cos_sum函数依然期望接收二维参数矩阵，导致访问P.shape[1]时出错（一维数组没有第二个维度）。

2. 具体修正步骤 #

步骤1：修改cos_sum的参数处理逻辑

允许输入一维参数数组，将其重塑为需要的二维形状：

def cos_sum(x, P):
    assert isinstance(P, np.ndarray)
    # 将一维参数数组重塑为N行3列的矩阵（N为余弦分量数量）
    P = P.reshape(-1, 3)
    sums = []
    for param in P:
        a, b, c = param
        sums.append(a * np.cos(b * (x - c)))
    sums = np.array(sums)
    return np.sum(sums, axis=0) 

步骤2：将初始猜测参数展平为一维

guess_flat = guess.flatten()

步骤3：修正minimize的args传递格式

args需要传入元组，(time)不是元组（仅为优先级括号），需添加逗号：

optimization = minimize(resid, guess_flat, args=(time,))

步骤4：优化残差计算（可选但推荐）

当前残差计算等价于L1损失，若使用更常用的L2损失（最小二乘），可修改为：

def resid(params, x):
    assert isinstance(params, np.ndarray)
    fit = cos_sum(x, params)
    # 均方根误差（L2损失）
    residual = np.sqrt(np.mean((fit - signal)**2))
    return residual

3. 完整修正代码示例 #

import numpy as np
from scipy.optimize import minimize

# 生成含噪声的信号
time = np.linspace(0, 1, 501)
data = np.cos(2 * np.pi * 4 * time) + np.cos(2 * np.pi * 9 * time) + np.cos(2 * np.pi * 20 * time)
noise = np.sqrt(1 / 25) * np.random.randn(501)
signal = data + noise 

# 修改后的余弦叠加函数
def cos_sum(x, P):
    assert isinstance(P, np.ndarray)
    P = P.reshape(-1, 3)
    sums = []
    for param in P:
        a, b, c = param
        sums.append(a * np.cos(b * (x - c)))
    sums = np.array(sums)
    return np.sum(sums, axis=0) 

# 修改后的残差函数
def resid(params, x):
    assert isinstance(params, np.ndarray)
    fit = cos_sum(x, params)
    residual = np.sqrt(np.mean((fit - signal)**2))
    return residual

# 构造初始猜测并展平
guess_A = np.random.normal(1, .2, size=3)
guess_B = 2 * np.pi * np.array([4, 9, 20], dtype=float)
guess_C = np.random.normal(0, .2, size=3)
guess = np.array([guess_A, guess_B, guess_C]).T
guess_flat = guess.flatten()

# 执行优化
optimization = minimize(resid, guess_flat, args=(time,))

# 提取优化后的参数并重塑为二维矩阵
optimized_params = optimization.x.reshape(-1, 3)
print("优化后的参数：")
print(optimized_params)

额外优化建议 #

• 添加参数边界约束：比如限制振幅a非负、角频率b为正，避免无意义参数，使用minimize的bounds参数实现。
• 向量化运算提升效率：若分量数量较多，替代循环用向量化计算：

def cos_sum(x, P):
    assert isinstance(P, np.ndarray)
    P = P.reshape(-1, 3)
    a, b, c = P[:,0], P[:,1], P[:,2]
    # 向量化生成所有分量的余弦值
    components = a[:, np.newaxis] * np.cos(b[:, np.newaxis] * (x - c[:, np.newaxis]))
    return np.sum(components, axis=0)

内容的提问来源于stack exchange，提问作者alpastor