Minimizing a function of two variables using scipy.optimize

Question

I have written the following code in 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)

Where q_sq, V_ub_par, and alpha_par are the function parameters, and the rest are stored variables. I am trying to minimize chi_sq(params) with respect to the parameters V_ub_par and alpha_par whilst q_sq takes several values through a for loop. I found a solution to minimizing a multivariate function using scipy.optimize here and I followed the same methodology; however, when I run the code, the following error pops up:

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

I am not sure what this error means, and how to fix it (for context, I am relatively new to python, and programming in general).

I am using the following libraries:

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

And the variables in the above code are either integers or arrays, given as follows:

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

To further debug your code i constructed this analyzing function and scalar generic typing functions, use it to evalute your values to a further extent.

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)}")

SOLUTION Your function chi_sq does not return a scalar. it returns a np.ndarray. Rework the chi_sq to ensure that it returns a type of int, complex or 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)

Resulting Code

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 2

The following fixes vectorisation gaps and runs successfully. It's suspicious that the minimum found is equal to the initial guess, but I consider that to be a separable problem.

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()