Curve_fit 对于 returns numpy 数组的函数

Question

我知道 scipy 的库 curve_fit 及其拟合曲线的能力。我在这里和文档中阅读了很多示例，但我无法解决我的问题。例如，我有 10 个文件（化学结构但没关系）和十个实验能量值。我在 class 中有一个函数，它为每个结构计算某些参数的理论能量，它 return 是一个具有理论能量值的 numpy 数组。

我想找到最佳参数，使理论值最接近实验值。我将在这里提供我的代码的最小示例

这是 class 函数，它读取实验能量文件，提取正确的子字符串，并 return 将值作为 numpy 数组。 self.path 只是目录，self.nPoints = 10。不是很重要，但为了完整起见，我提供了

def experimentalValues(self):
        os.chdir(self.path)
        energy = np.zeros(self.nPoints)
        for i in range(1, self.nPoints):
            f = open("p_" + str(i + 1) + ".xyz", "r")
            energy[i] = float(f.readlines()[1].split()[1])
            f.close()
        os.chdir('..')
        return energy

我用这个 class 函数计算理论值，它以两个 numpy 数组作为参数，假设

sigma = np.full(nSubstrate, 2.)
epsilon = np.full(nSubstrate, 0.15)

其中 nSubstrate = 9

这里有class函数。它读取文件并执行两个嵌套循环来计算每个文件的理论值并将其 return 到一个 numpy 数组。

def theoreticalEnergy(self, epsilon, sigma):
        os.chdir(self.path)
        cE = np.zeros(self.nPoints)
        for n in range(0, self.nPoints):
            filenameXYZ = "p_" + str(n + 1) + "_extended.xyz"

            allCoordinates = np.loadtxt(filenameXYZ, skiprows = 0, usecols = (1, 2, 3))
            substrate = allCoordinates[0:self.nSubstrate]
            surface = allCoordinates[self.nSubstrate:]
            for i in range(0, substrate.shape[0]):
                positionAtomI = np.array(substrate[i][:])
                for j in range(0, surface.shape[0]):
                    positionAtomJ = np.array(surface[j][:])
                    distanceIJ = self.distance(positionAtomI, positionAtomJ)
                    cE[n] += self.LennardJones(distanceIJ, epsilon[i], sigma[i])
                
        os.chdir('..')
        return cE

同样，为了完整起见，Lennard Jones class 函数定义为

def LennardJones(self, distance, epsilon, sigma):
        repulsive = (sigma/distance) ** 12.
        attractive = (sigma/distance) ** 6.
        potential = 4. * epsilon* (repulsive - attractive)
        return potential

在这种情况下，所有参数都是作为 return 值的标量。总结问题介绍我有 3 个成分：

带有实验数据的 numpy 数组
两个带有对参数 sigma 和 epsilon 的猜测的 numpy 数组
一个接受最后一个参数的函数，return一个带有待拟合值的 numpy 向量。

如何像文档https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html中描述的方法那样解决这个问题？

Answer 1

曲线拟合

curve_fit 通过找到最小化 sum((f(w, x[i] - y[i])**2 for i in range(n)) 的 w 将函数 f(w, x[i]) 拟合到点 y[i]。正如您将在函数定义后的第一行中看到的那样

[It uses] non-linear least squares to fit a function, f, to data.

它指的是 least_squares 它指出的地方

Given the residuals f(x) (an m-D real function of n real variables) and the loss function rho(s) (a scalar function), least_squares finds a local minimum of the cost function F(x):

曲线拟合是一种凸成本多objective优化。由于每个单独的成本都是凸的，你可以将它们全部相加，这仍然是一个凸函数。请注意决策变量（要优化的参数）在每个点都是相同的。

你的问题

根据我的理解，每个能级都有一组不同的参数，如果将其写为曲线拟合问题，objective 函数可以表示为 sum((f(w[i], x[i]) - y[i])**2 ...), where y[i] is determined by the energy level. Since each of the terms in the sum is independent on the other terms, this is equivalent to finding each group of parametersw[i]separately minimizing(f(w[i], x[i]) - y[i])**2`.

凸性

凸性是一种非常方便的优化属性，因为它确保您在参数space中只有一个最小值。我没有做详细分析，但对你的能量函数的凸性有合理的怀疑。

Lennard Jones 函数具有排斥力和吸引力的区别，两者在距离上都具有负偶数指数，仅此一项就不太可能是凸的。
以不同位置为中心的多个局部函数之和没有定义的凸性
分子能量，或 crystal 能量，或蛋白质折叠众所周知是非凸的。

几天前（在骑自行车时）我在想这个问题，分子将如何以全局最小能量配置，我想知道它是否由于量子隧道效应而如此迅速地找到该配置。

非凸优化

非凸（全局）优化与（非线性）最小二乘法不同，在某种意义上，当找到局部最小值时，过程不会立即 return，它开始不同区域搜索的新尝试 spaces。如果函数是平滑的，你仍然可以利用基于梯度的局部优化方法，但复杂度仍然是 NP。

经典的全局优化方法是Simulated annenaling, if you have a chemical background I think you will have some insights reading about it. Once upon a time，scipy.optimize中提供了模拟退火。

你会发现一些global optimization methods in scipy.optimize. I would encourage you to try Basin hopping, since it was successfully applied to similar problems, as you can read in the references。

我希望这能让您找到解决问题的正确方法。但是，请注意，您可能需要花钱学习如何使用该功能，并且需要做出一些决定。您需要在准确性、简单性和效率之间找到平衡点。

如果您想要更好的解决方案，请花时间推导成本函数的梯度（您可以 return 两个值 f 和 df，其中 df 是 f 相对于决策变量的梯度） .

Curve_fit 对于 returns numpy 数组的函数

Curve_fit for a function that returns a numpy array

python

numpy

curve-fitting

scipy

曲线拟合

你的问题

凸性

非凸优化