看到一堆点后试图绘制某种趋势的曲线的人。每个人都有这种想法。当只有几个点并且我绘制的曲线只是一条直线时，这很容易。但是每次我加更多的点，或者当我要找的曲线与直线不同时，它就会变得越来越难。在这种情况下，曲线拟合过程可以解决我所有的问题。输入一堆点并找到“完全”匹配趋势的曲线是令人兴奋的。但这如何工作？为什么拟合直线与拟合奇怪形状的曲线并不相同。每个人都熟悉线性最小二乘法，但是，当我们尝试匹配的表达式不是线性时，会发生什么？这使我开始了一段数学文章之旅，stack overflow发布了[1]一些深奥的数学表达式(至少对我来说是这样的!)，以及一个关于发现算法的有趣故事。这是我试图用最简单而有效的方式来解释这一切。
提出问题
在某些情况下，线性回归是不够的。有时需要将一系列数据调整为非线性表达式。在这些情况下，普通最小二乘对我们不起作用，我们需要求助于不同的方法。我第一次遇到这种情况是在我尝试将2D数据拟合到如下函数时:
幸运的是，我可以通过许多方法自动找到Beta的最佳值。任何熟悉MATLAB中的nlinfit或SciPy的curve_fit函数的人都知道，一旦您有了模型的数学表达式，这个非线性回归过程是简单的。所有这些库的工作方式都类似，它们使用迭代算法，试图找到参数或参数组合，使观测数据和模型响应之间的差异最小化。我们用一些方程来表示它。
假设我们有一个函数f它由一个自变量x和一组参数a决定，这是y= f(x,a)这个函数正在对我们已经知道输出ŷ的流程进行建模。目标是找到一组参数a，使y尽可能接近ŷ。衡量我们离ŷ有多近的一种方法是计算差的平方和。残差定义为y和ŷ在每一点上的差。这可以表示为:
在本例中，下标i指的是我们正在分析的数据点。如果我们试图用100个数据点调整一条曲线，那么我们需要计算每一个数据点的差。最后，我们会得到一个r1 r2 r3，等等，直到我们在这个例子中达到r100。差平方和对应于:
找到生成S可能的最低值的参数a组合，意味着参数a是从我们的模型计算出的y与ŷ值之间可能的最佳匹配。
用图形的方式来显示这个问题
下图用红色表示一些数据点，用紫色表示模型响应。如果我们想测量这个模型如何适应数据点，我们可以计算数据点(ŷ)和模型响应(y)之间的差异，然后将这些差异的平方和(残差)。这种思想可以外推到包含多个自变量(x1,x2，…，xn)的函数上。
解决方案
求函数最小值的一种常用方法是计算函数对特定变量的导数。在这种情况下，我们想找到使函数s最小的a值。可以写成:
下标j表示a可能有多个值，因为函数f依赖于自变量x和一个或多个参数a1, a2，…，aM。在这种情况下，我们需要根据每个参数部分推导函数。当函数的导数值为零时，函数的最小值才会出现。所以，我们之前的方程会是这样的:
注意我是如何展开ri的，只是为了提醒你这个差就是计算值和实际值之间的差。在这一点上，重要的是要有一个关于导数的图解解释，以及为什么当它们等于0时，我们可以说我们找到了一个最小值(或最大值)。
用导数使函数最小化的图解说明
一个导数可以被定义为一个函数相对于它的参数如何变化的度量。我们能找到的一个最简单的例子是y=mx类型的函数。这个函数关于x的导数(dy/dx)是m，这意味着x每改变一点，输出y就改变m次。所以这个函数的导数表示了x变化后y的变化量，直观上，这可以看作是函数中某一点上切线的斜率。
下图展示了一个与我们之前提到的直线完全不同的函数。函数的类型y = mx,变化量的比值对x总是不管x的值是相同的。在这种情况下,这一比率变化根据x。你可以看到每个点所示图有不同斜率切线(m)。这个斜率表示函数在某一点的导数。求函数的最小值和最大值的一种方法是寻找斜率为零的地方。在这种情况下，一个24.5的x将给我们一个最小值，而一个10的x将给我们一个最大值。
现在可能感觉事情开始变得复杂了，但是请耐心听我说。这并不像看起来那么难!如果你熟悉微积分和导数，你就会知道推导最后一个方程的结果是:
项df（xi，aj）/ daj对应于函数f关于每个参数a的导数。在这一点上要记住的是，我们的模型可以包含多个参数a，并且我们需要找到函数f相对于每个参数的导数。
请注意，此计算是针对数据中存在的每个点执行的。这就是为什么我们的函数f取决于xi和aj的原因：我们有x的i值和a的j值。我们可以将所有这些导数汇编成一个称为Jacobian的术语。雅可比行列式是一个矩阵，其中包含一个函数相对于每个参数的所有一阶偏导数。
记住，下标i代表一个特定的数据点。如果数据包含100个点那么雅可比矩阵就有100行3列因为我们有3个参数。
如果我们使用雅可比行列式的概念来重写最后找到的dS / da方程。我们将有：
注意我是如何用矩阵来表示这个方程的。我去掉了现在雅可比矩阵的和，剩余都用矩阵来写。记住，所有这些方程都是针对所有数据点同时求解的，所以使用矩阵是非常方便的。在这一点上，我将向您展示两种方法，我们可以解决这个方程，并找到参数更好地调整初始方程f。
梯度下降
你可能听过这个名字。梯度下降法是一种优化算法，用于寻找函数的局部最小值。它背后的理念并不难理解。因为我们要求最小值的函数是可微的我们知道任意点处的梯度。这意味着我们知道要继续往下走，我们需要走哪个方向。在每次迭代中，我们都会向函数的最小值移动一点。梯度下降法的两个重要方面是初始猜测和我们在每次迭代时采取的步骤的大小。这种方法的效率在这两个方面是非常可靠的。
这和非线性回归有什么关系?好的，我们可以使用梯度下降法来求函数s的最小值。在这种情况下，我们向最小值点所采取的每一步都可以表示为:
将此高阶差分添加到参数的初始估计中，并重复此过程，直到我们找到一个最小迭代次数或我们超过最大迭代次数为止。在最后一个方程中出现的α是用来增加或减少我们所采取的步骤的大小。正如我前面提到的，梯度下降法的性能与步骤的大小以及初始猜测有很大关系。
高斯牛顿法
梯度下降法是众所周知和广泛使用的，但它可能是相当缓慢并取决于参数的数量。另一种方法是高斯-牛顿法，它类似于梯度下降法，是一种迭代过程，我们采取多个步骤，直到我们接近正确的解。在本例中，我们通过以下方式得到一个新的参数组合:
hGN代表我们采用高斯-牛顿法的步骤。我们如何知道每次迭代的hGN值？
在高斯-牛顿法中，函数f是使用一阶泰勒展开式近似的，这意味着
还记得我们说过的术语dfi（a）/ daj也称为雅可比行列式，因此前面的等式也可以写成：
如果我们使用此表达式用f（an + 1）代替f（an），我们将得出：
可以重新组织为：
并使用以下公式计算步长：
下表适用于两种方法。在这两种情况下，都必须指定参数的初始猜测以及停止条件。在这种情况下，停止标准由最大迭代次数或平方误差的最小值组成。
Levenberg-Marquardt 或 damped least squares
请注意，hGD和hGN方程非常相似，这与Levenberg-Marquardt方法有很大关系。该方法根据我们与解的接近程度在梯度下降和高斯牛顿之间切换。Levenberg-Marquardt方法表示为：
在前面的等式中，I表示单位矩阵，并且λ被称为阻尼因子。此参数是允许在高斯牛顿或梯度下降更新之间进行更改的参数。当λ小时，该方法采用高斯-牛顿步长；当λ大时，该方法遵循梯度下降法。通常，λ的第一个值较大，因此第一步位于梯度下降方向[2]。其背后的逻辑是，高斯-牛顿法在最终迭代中更有效，而梯度下降法在过程开始时很有用，因为该过程仍距离理想解决方案还很远。
如您所见，Levenberg-Marquardt算法是梯度下降算法与高斯-牛顿算法的结合。因此，Levenberg-Marquardt算法的效率也高度依赖于初始猜测的选择以及阻尼系数[3]。另外，阻尼因子的增加和减少也影响算法的性能。在每次迭代中，阻尼系数将乘以或除以一个系数，具体取决于前一次迭代的质量。通常，lambda增加2倍，减少3倍。
作为一个有趣的历史记录，Levenberg在1944年提出了这种方法，但是直到19年后，Marquardt在一篇论文中提到后他的方法才被人们所知晓。Levenberg在陆军弹药厂Frankford Arsenal工作时发布了该算法。从某种意义上讲，可以说Marquardt重新发现了damped least squares，这就是为什么今天两个名称都被用作参考。Pujol [4]完整描述了Levenberg和Marquardt所做的工作，以及它们各自如何分别为我们今天所知的算法做出了贡献。
代码
Levenberg-Marquardt算法可以以多种形式实现[5]。在这种情况下，我将介绍一种ython实现此算法的非常简单的方法。我还在将我的结果与Scipy的curve_fit函数的结果进行比较。此函数对算法的实现更可靠，将比我向您展示的算法更好。但是，我认为这段代码对于任何更复杂的事情以及了解“幕后”正在发生的事情都是一个很好的起点。尽管此笔记本中显示的示例涉及到二维问题，但是该算法背后的逻辑可以应用于多种情况。
本笔记本中包含的示例称为DCA，这是石油工程界常用的方法。笔记本包含DCA的简要说明和一些示例。您可以在我的GitHub存储库中访问此代码。https://github.com/manfrezord/MediumArticles/blob/main/LM-Algorithm.ipynb
总结
科技发展使得执行几年前可能要花费很长时间的计算变得更加容易。但是，了解所有这些计算的来源始终很重要。进行线性和非线性回归是可以在数据分析和机器学习中完成的许多其他事情的基础。如今，当每个人都在注视着这些领域试图寻找答案或更有效地执行流程时，重要的是要了解基本原理。
