使用“scipy.optimize.least_squares()”用可变数量的参数拟合两条曲线

Question

我正在使用

scipy.optimize.least_squares()

函数编写代码来适应系列中相关函数对的参数，但测试参数似乎没有正确传递给函数。这是问题的简化说明。

import numpy as np
import numpy.linalg as la
import scipy.optimize
import math

Ha = {'H': lambda x, a, b : a + b*x, 'nParams': 2}
Hb = {'H': lambda x, a, b, c : a + b*x + c*x*x, 'nParams': 3}
Hc = {'H': lambda x, a, b, c, d : a*math.cos(b*(x-c)) + d, 'nParams': 4}

points1 = [(1,0), (2,4), (3,5), (7,8)]
points2 = [(1,1), (2,6), (4,1), (6,6)]

def coupled_1 (x, g, *params, **kwargs):
    H1, H2 = kwargs.values() # H1 and H2 are dictionaries containing the two curves
                             # to be fit, 'H', and their respective numbers of
                             # parameters, 'nParams'
    matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
            + np.diag([g],k=1) + np.diag([g],k=-1)
    return la.eigh(matrix)['eigenvalues'][0]

def coupled_2 (x, g, *params, **kwargs):
    H1, H2 = kwargs.values()
    matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
            + np.diag([g],k=1) + np.diag([g],k=-1)
    return la.eigh(matrix)['eigenvalues'][1]

def getParameters(H1, H2, pts1, pts2):
    nParams = H1['nParams'] + H2['nParams']
    # These points sometimes come in pairs and sometimes don't so I've found no better way to zip them
    pts     = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
    res     = lambda *params, **kwargs: \
                     np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
                     + (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
    result  = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
    return result['x']

params = getParameters(Ha, Hc, points1, points2)

我从

Ha()

收到以下错误

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[5], line 37
     34     result  = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
     35     return result['x']
---> 37 params = getParameters(Ha, Hc, points1, points2)
     38 print(params)

Cell In[5], line 34, in getParameters(H1, H2, pts1, pts2)
     30 pts     = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
     31 res     = lambda *params, **kwargs: \
     32                  np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
     33                  + (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
---> 34 result  = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
     35 return result['x']

File ~\anaconda3\lib\site-packages\scipy\optimize\_lsq\least_squares.py:830, in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
    827 if method == 'trf':
    828     x0 = make_strictly_feasible(x0, lb, ub)
--> 830 f0 = fun_wrapped(x0)
    832 if f0.ndim != 1:
    833     raise ValueError("`fun` must return at most 1-d array_like. "
    834                      "f0.shape: {0}".format(f0.shape))

File ~\anaconda3\lib\site-packages\scipy\optimize\_lsq\least_squares.py:825, in least_squares.<locals>.fun_wrapped(x)
    824 def fun_wrapped(x):
--> 825     return np.atleast_1d(fun(x, *args, **kwargs))

Cell In[5], line 32, in getParameters.<locals>.<lambda>(*params, **kwargs)
     29 # These points sometimes come in pairs and sometimes don't so I've found no better way to zip them
     30 pts     = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
     31 res     = lambda *params, **kwargs: \
---> 32                  np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
     33                  + (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
     34 result  = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
     35 return result['x']

Cell In[5], line 32, in <listcomp>(.0)
     29 # These points sometimes come in pairs and sometimes don't so I've found no better way to zip them
     30 pts     = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
     31 res     = lambda *params, **kwargs: \
---> 32                  np.sqrt(sum( [( p[0][1]-coupled_1(p[0][0], *params, **kwargs) )**2 \
     33                  + (p[1][1]-coupled_2(p[0][0], *params, **kwargs) )**2 for p in pts] ))
     34 result  = scipy.optimize.least_squares(res,[1] + [0]*nParams, kwargs={'H1':H1,'H2':H2})
     35 return result['x']

Cell In[5], line 17, in coupled_1(x, g, *params, **kwargs)
     14 H1, H2 = kwargs.values() # H1 and H2 are dictionaries containing the two curves
     15                          # to be fit, 'H', and their respective numbers of
     16                          # parameters, 'nParams'
---> 17 matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
     18         + np.diag([g],k=1) + np.diag([g],k=-1)
     19 return la.eigh(matrix)['eigenvalues'][0]

TypeError: <lambda>() missing 2 required positional arguments: 'a' and 'b'

似乎

*params

列表没有进入

Ha()

，但我不确定为什么。

Answer 1

这是 res 函数的另一种方法，这对我有用。

import numpy as np
import numpy.linalg as la
import scipy.optimize
import math

# Define the dictionaries for the functions
Ha = {'H': lambda x, a, b: a + b*x, 'nParams': 2}
Hb = {'H': lambda x, a, b, c: a + b*x + c*x*x, 'nParams': 3}
Hc = {'H': lambda x, a, b, c, d: a*math.cos(b*(x-c)) + d, 'nParams': 4}

points1 = [(1,0), (2,4), (3,5), (7,8)]
points2 = [(1,1), (2,6), (4,1), (6,6)]

def coupled_1(x, g, *params, **kwargs):
    H1, H2 = kwargs['H1'], kwargs['H2']
    matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
            + np.diag([g], k=1) + np.diag([g], k=-1)
    return la.eigh(matrix)[0][0]  # Get the smallest eigenvalue

def coupled_2(x, g, *params, **kwargs):
    H1, H2 = kwargs['H1'], kwargs['H2']
    matrix = np.diag([H1['H'](x, *params[:H1['nParams']]), H2['H'](x, *params[H1['nParams']:])]) \
            + np.diag([g], k=1) + np.diag([g], k=-1)
    return la.eigh(matrix)[0][1]  # Get the second smallest eigenvalue

def getParameters(H1, H2, pts1, pts2):
    nParams = H1['nParams'] + H2['nParams']
    # Filter points with matching x-values
    pts = [(p1, p2) for p1 in pts1 for p2 in pts2 if p1[0] == p2[0]]
    
    def res(params, **kwargs):
        error_sum = 0
        for p1, p2 in pts:
            x_val = p1[0]
            y1 = p1[1]
            y2 = p2[1]
            error_sum += (y1 - coupled_1(x_val, *params, **kwargs))**2
            error_sum += (y2 - coupled_2(x_val, *params, **kwargs))**2
        return np.sqrt(error_sum)
    
    # Initial guess for parameters
    initial_guess = [1] + [0]*nParams
    
    # Optimize parameters using least squares method
    result = scipy.optimize.least_squares(res, initial_guess, kwargs={'H1': H1, 'H2': H2})
    
    return result['x']

# Obtain parameters
params = getParameters(Ha, Hc, points1, points2)
print(params)