使用 scipy.optimize 最小化两个变量的函数

Question

我用Python编写了以下代码：

def f_multi_para(q_sq, alpha_par):
    return f_0 / (q_sq * alpha_par)

def dB_dqsq_model2_para(q_sq, V_ub_par, alpha_par):
    sec1 = G_F**2 * V_ub_par
    sec2 = (1 - m_e**2 / q_sq)**2
    sec3 = (q_sq * (1 + m_e**2 / q_sq) * f_multi_para(q_sq, alpha_par)**2)
    return sec1 * sec2 * sec3

def chi_sq(params):
    V_ub_par, alpha_par = params
    return np.sum(((dB_dqsq_arr - np.array([dB_dqsq_model2_para(v, V_ub_par, alpha_par) for v in q_sq_arr])) / dBdqsq_err_arr)**2)

initial_guess = [0.0037, 0.54]
result = optimize.minimize(chi_sq, initial_guess)
if result.success:
    fitted_params = result.x
    print(fitted_params)
else:
    raise ValueError(result.message)

其中

q_sq

、

V_ub_par

、

alpha_par

为函数参数，其余为存储变量。我试图最小化关于参数

chi_sq(params)

和

V_ub_par

的

alpha_par

，同时

q_sq

通过 for 循环获取多个值。我找到了使用

scipy.optimize

here 最小化多元函数的解决方案，并且遵循了相同的方法；但是，当我运行代码时，会弹出以下错误：

ValueError: The user-provided objective function must return a scalar value.

我不确定这个错误意味着什么，以及如何修复它（对于上下文，我对 python 和一般编程相对较新）。

我正在使用以下库：

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy.integrate import quad
import scipy.optimize as optimize

上面代码中的变量要么是整数，要么是数组，给出如下：

f_0 = 0.261
G_F = 1.166 * 10**(-5)
m_e = 0.511 * 10**(-3)
dB_dqsq_arr = np.array([7.2,7.14,6.7,7.56,6.44,7.17,6.67,6.33,6.2,4.32,4.25,3.4,1.17])
dBdqsq_err_arr = np.array([0.70,0.45,0.39,0.43,0.43,0.45,0.47,0.48,0.44,0.43,0.41,0.42,0.26])
q_sq_arr = np.array([1,3,5,7,9,11,13,15,17,19,21,23,25])

Answer 1

您实际上是在要求 scipy 一次优化多个问题，因为您正在评估

q_sq

的多个值的结果。您收到错误是因为

chi_sq

返回一个列表而不是单个值，即标量。相反，您需要对每个

q_sq

值运行优化。您没有提供足够的信息来自己运行代码，但它可能看起来像这样：

# above code is unchanged

def chi_sq(params, q_sq):
    return ((dB_dqsq_arr - dB_dqsq_model2_para(q_sq, *params)) / dBdqsq_err_arr)**2
    
results = []
for q_sq in q_sq_arr:
    initial_guess = [0.0037, 0.54]
    result = optimize.minimize(chi_sq, initial_guess, args=(q_sq,))
    if result.success:
        fitted_params = result.x
        print(fitted_params)
        results.append(fitted_params)
    else:
        raise ValueError(result.message)

Answer 2

为了进一步调试您的代码，我构建了这个分析函数和标量通用类型函数，用它来进一步评估您的值。

from typing import Union, TypeVar
import numpy as np
from scipy.optimize import OptimizeResult


# Define a generic scalar type Scalar
scalar_types = (int, float, complex)
Scalar = TypeVar('Scalar', *scalar_types)


def scalar(value: Scalar) -> Union[Scalar, Exception]:
    """ construct a scalar using generic typing """

    if isinstance(value, scalar_types):
        return value
    
    elif isinstance(value, np.ndarray):
        raise Exception("Expected a scalar (int, complex, float) value but received np.ndarray")

    elif isinstance(value, OptimizeResult):
        raise Exception("Expected a scalar but got an scipy.optimize.OptimizeResult")

    try:
        print(dir(value))
    except:
        
        try:
            print(vars(value))
        except:
            print("dir and vars can not be used on ", type(value))
        
    return Exception(f"Expected a scalar (int, complex, float) value but received {type(value)}")

解决方案 您的函数 chi_sq 不返回标量。它返回一个 np.ndarray。重新设计 chi_sq 以确保它返回 int、complex 或 float 类型。

def chi_sq(params):
    V_ub_par, alpha_par = params
    vals = (dB_dqsq_arr - np.array([dB_dqsq_model2_para(v, V_ub_par, alpha_par) for v in q_sq_arr]))
    chi = [x / dBdqsq_err_arr for x in vals]
    return np.sum(chi)

结果代码

from math import *
import numpy as np
from scipy.integrate import quad
import scipy.optimize as optimize


# DATA
f_0 = 0.261 # Double
G_F = 1.166 * 10**(-5) # Double
m_e = 0.511 * 10**(-3) # Double

# TYPE array[float]
dB_dqsq_arr = np.array([7.2,7.14,6.7,7.56,6.44,7.17,6.67,6.33,6.2,4.32,4.25,3.4,1.17])
# TYPE array[float]
dBdqsq_err_arr = np.array([0.70,0.45,0.39,0.43,0.43,0.45,0.47,0.48,0.44,0.43,0.41,0.42,0.26])
# TYPE array[int]
q_sq_arr = np.array([1,3,5,7,9,11,13,15,17,19,21,23,25])



def f_multi_para(q_sq, alpha_par):
    return f_0 / (q_sq * alpha_par)

def dB_dqsq_model2_para(q_sq, V_ub_par, alpha_par):
    sec1 = G_F**2 * V_ub_par
    sec2 = (1 - m_e**2 / q_sq)**2
    sec3 = (q_sq * (1 + m_e**2 / q_sq) * f_multi_para(q_sq, alpha_par)**2)
    return sec1 * sec2 * sec3

def chi_sq(params):
    V_ub_par, alpha_par = params
    vals = (dB_dqsq_arr - np.array([dB_dqsq_model2_para(v, V_ub_par, alpha_par) for v in q_sq_arr]))
    chi = [x / dBdqsq_err_arr for x in vals]
    return np.sum(chi)


initial_guess = [0.0037, 0.54]
print(type(chi_sq(initial_guess)))
result = optimize.minimize(chi_sq, initial_guess)
if result.success:
    fitted_params = result.x
    print(fitted_params)
else:
    raise ValueError(result.message)

Answer 3

以下内容修复了矢量化间隙并成功运行。发现的最小值等于最初的猜测是值得怀疑的，但我认为这是一个可分离的问题。

import numpy as np
from scipy import optimize

f_0 = 0.261
G_F = 1.166e-5
m_e = 0.511e-3
dB_dqsq_arr = np.array((7.2,7.14,6.7,7.56,6.44,7.17,6.67,6.33,6.2,4.32,4.25,3.4,1.17))
dBdqsq_err_arr = np.array((0.70,0.45,0.39,0.43,0.43,0.45,0.47,0.48,0.44,0.43,0.41,0.42,0.26))
q_sq_arr = np.arange(1, 26, 2)


def f_multi_para(q_sq: int | np.ndarray, alpha_par: float) -> float:
    return f_0 / (q_sq * alpha_par)


def dB_dqsq_model2_para(q_sq: int | np.ndarray, V_ub_par: float, alpha_par: float) -> float:
    sec1 = G_F**2 * V_ub_par
    sqdiff = 1 - m_e**2 / q_sq
    sec2 = sqdiff**2
    sec3 = q_sq * sqdiff * f_multi_para(q_sq, alpha_par)**2
    return sec1 * sec2 * sec3


def chi_sq(params: np.ndarray) -> float:
    V_ub_par, alpha_par = params
    chi = (
        dB_dqsq_arr -
        dB_dqsq_model2_para(q_sq_arr, V_ub_par, alpha_par)
    ) / dBdqsq_err_arr
    return chi.dot(chi)


def regression_test() -> None:
    actual = chi_sq((0.004, 0.54))
    assert np.isclose(actual, 2307.989692506141, rtol=0, atol=1e-12)


def solve() -> None:
    initial_guess = (0.0037, 0.54)
    result = optimize.minimize(fun=chi_sq, x0=initial_guess)
    if not result.success:
        raise ValueError(result.message)
    fitted_params = result.x
    print(fitted_params)


if __name__ == '__main__':
    regression_test()
    solve()