输入 numpy 数组无法解决 optimize.curve_fit 中的值错误

Question

我正在尝试使用scipy.optimize.curve_fit 将一些数据拟合到函数中，但不断收到值错误，类似于this post 中描述的错误，但是按照 cmets 中的建议将所有输入放入 numpy 数组中并没有似乎解决了我的问题。

这是我的基本代码：

def RamanHamiltonian(k, omega, delta, epsilon):
    H = np.array([[(k+2.0)**2.0 - delta, omega/2.0, 0.0],
                  [omega/2.0, k**2.0-epsilon,omega/2.0],
                  [0.0,omega/2.0,(k-2.0)**2.0 + delta]])
    return H

def propagateHamiltonianTest(t, omega, delta, epsilon):  
    k = 0.0
    psi0 = np.array([0+1j*0.0, 1.0+1j*0.0, 0.0+1j*0.0])
    H = RamanHamiltonian(k, omega, delta ,epsilon)
    Energy, V = LA.eig(H)

    V = V + 1j*0.0
    Vinv = np.conjugate(np.transpose(V))

    U = np.diag(np.exp(-1j*np.array(Energy)*t))
    a  =np.dot(Vinv,psi0)
    b = np.dot(U,a)
    psi = np.dot(V,b)
    pop0 = np.absolute(psi[0])**2.0

    return pop0

popt, pcov = optimize.curve_fit(propagateHamiltonianTest,
                               np.array(tRecoils), 
                               np.array(frac0), 
                               p0=(3.0,0.05,0.03))

下面是tRecoils 和frac0 的值，它们的长度都是24 个数组：

tRecoils = array([ 2.88597836,  1.15439135,  1.73158702,  2.19334356,  0.23087827,
    2.07790442,  0.11543913,  1.50070875,  2.77053923,  1.61614788,
    1.03895221,  2.42422183,  0.92351308,  0.80807394,  0.3463174 ,
    0.57719567,  2.6551001 ,  0.46175654,  1.84702615,  1.38526961,
    2.53966096,  0.69263481,  2.30878269,  1.96246529])
frac0 = array([ 0.15761062,  0.02044625,  0.17275937,  0.02236243,  0.07388558,
    0.00967176,  0.01886309,  0.20412516,  0.21667489,  0.21783697,
    0.00173812,  0.14038657,  0.03145599,  0.08644404,  0.13153078,
    0.18794377,  0.2139092 ,  0.17141201,  0.13021916,  0.12671806,
    0.21090369,  0.1611094 ,  0.08732627,  0.05764911])

我收到此错误：

U = np.diag(np.exp(-1j*np.array(能量)*t))

ValueError: 操作数不能与形状一起广播 (3) (24)

因此，由于某种原因，curve_fit 读取整个 tRecoil 数组，而不是按元素进行读取，我似乎无法通过更改输入格式来解决这个问题。

Answer 1

我认为问题在于如何评估您的函数 propagateHamiltonianTest 存在歧义。它既可以应用于t 的单个值（这就是您的想法），也可以应用于整个数组，在这种情况下，numpy 将整个数组传递给参数t 并尝试向量化。

事实上，例如，当我在数组tRecoil[0] 的单个元素上测试你的函数时，它可以工作，但是当你在整个数组上尝试它时，它会失败并出现同样的错误。所以问题与curve_fit无关。如果你可以让你的函数接收整个数组并返回一个数组，那么它应该与curve_fit一起工作。

为此，我不得不使用一些 hacky 技巧，使用 scipy.linalg.block_diag。这使我们能够创建一个更大的块对角矩阵，它基本上具有沿对角线的每个时间步产生的每个 3x3 矩阵。

下面是带有 cmets 的 propagateHamiltonianTest 函数的修改版本。我已经测试过了，它应该可以工作。

from scipy.linalg import block_diag

def propagateHamiltonianTest(t, omega, delta, epsilon):  
    k = 0.0
    psi0 = np.array([0+1j*0.0, 1.0+1j*0.0, 0.0+1j*0.0])
    H = RamanHamiltonian(k, omega, delta ,epsilon)
    Energy, V = LA.eig(H)

    V = V + 1j*0.0
    Vinv = np.conjugate(np.transpose(V))

    # np.outer(t, Energy).flatten() creates a matrix for all t
    U = np.diag(np.exp(-1j*np.outer(t, Energy).flatten()))  
    a = np.dot(Vinv, psi0)
    # This repeats a so that the shape is consitent with U
    aa = block_diag(*([a]*t.size))                          
    # Have to add the transpose to make shapes match 
    b = np.dot(U, aa.T)                                     
    # Same block diagonal trick for eigenvector matrix
    VV = block_diag(*([V]*t.size))                          
    psi = np.dot(VV, b)
    pop0 = np.absolute(psi)**2.0                       
    # Since you want the first value, need to take every 3rd row 
    # and extract the values you want from the diagonal
    return np.diag(pop0[::3])

输出：

print('popt: ', popt, '\n\npcov: ', pcov)

popt:  [ 2.59126808 -0.13125704  0.52597681] 

pcov:  [[ 0.00114962  0.00072806 -0.0005113 ]
        [ 0.00072806  0.00076575 -0.00020831]
        [-0.0005113  -0.00020831  0.00063511]]