非线性模型拟合：从基础到高级方法

# 非线性模型拟合：从基础到高级方法 ## 1. 非线性模型拟合概述 在处理模型与一组 M 个未知参数 $a_k$（$k = 1, 2, \cdots, M$）呈非线性依赖关系的拟合问题时，我们采用与之前类似的方法，即定义一个 $\chi^2$ 优度函数，并通过最小化该函数来确定最佳拟合参数。由于存在非线性依赖，最小化过程必须迭代进行。给定参数的试验值后，我们会开发一个改进试验解的程序，然后重复该程序，直到 $\chi^2$ 停止（或实际上停止）减小。 与一般的非线性函数最小化问题相比，这里有一个关键区别。在足够接近最小值的情况下，$\chi^2$ 函数可以用二次形式很好地近似： $\chi^2(a) \approx \gamma - d \cdot a + \frac{1}{2}a \cdot D \cdot a$ 其中 $d$ 是一个 M 维向量，$D$ 是一个 $M \times M$ 矩阵。如果这个近似效果良好，我们可以从当前试验参数 $a_{cur}$ 一步跳到最小化参数 $a_{min}$： $a_{min} = a_{cur} + D^{-1} \cdot [-\nabla\chi^2(a_{cur})]$ 然而，如果上述二次近似在 $a_{cur}$ 处效果不佳，我们只能像最速下降法那样沿梯度方向迈出一步： $a_{next} = a_{cur} - constant \times \nabla\chi^2(a_{cur})$ 其中常数要足够小，以免耗尽下坡方向。 要使用上述两个公式，我们必须能够计算任意参数集 $a$ 下 $\chi^2$ 函数的梯度。使用第一个公式时，还需要 $\chi^2$ 优度函数在任意 $a$ 处的二阶导数矩阵（Hessian 矩阵）$D$。与之前处理的问题不同的是，这里我们确切知道 $\chi^2$ 的形式，因为它基于我们自己指定的模型函数，所以 Hessian 矩阵是已知的。因此，只要愿意，我们可以自由使用第一个公式。只有当第一个公式无法改善拟合效果，即二次近似效果不佳时，才会使用第二个公式。 ## 2. 梯度和 Hessian 矩阵的计算 ### 2.1 模型和 $\chi^2$ 优度函数 要拟合的模型为： $y = y(x; a)$ $\chi^2$ 优度函数为： $\chi^2(a) = \sum_{i = 1}^{N} \left[\frac{y_i - y(x_i; a)}{\sigma_i}\right]^2$ ### 2.2 梯度计算 $\chi^2$ 关于参数 $a$ 的梯度在 $\chi^2$ 最小值处为零，其分量为： $\frac{\partial\chi^2}{\partial a_k} = -2 \sum_{i = 1}^{N} \frac{[y_i - y(x_i; a)]}{\sigma_i^2} \frac{\partial y(x_i; a)}{\partial a_k}$，$k = 1, 2, \cdots, M$ ### 2.3 Hessian 矩阵计算 再求一次偏导数可得： $\frac{\partial^2\chi^2}{\partial a_k\partial a_l} = 2 \sum_{i = 1}^{N} \frac{1}{\sigma_i^2} \left[\frac{\partial y(x_i; a)}{\partial a_k} \frac{\partial y(x_i; a)}{\partial a_l} - [y_i - y(x_i; a)]\frac{\partial^2 y(x_i; a)}{\partial a_l\partial a_k}\right]$ 为了简化，通常定义： $\beta_k \equiv -\frac{1}{2} \frac{\partial\chi^2}{\partial a_k}$ $\alpha_{kl} \equiv \frac{1}{2} \frac{\partial^2\chi^2}{\partial a_k\partial a_l}$ 这样，在最小二乘的背景下，矩阵 $[\alpha]$（等于 Hessian 矩阵的一半）通常被称为曲率矩阵。上述公式可以重写为一组线性方程： $\sum_{l = 1}^{M} \alpha_{kl} \delta a_l = \beta_k$ 求解这组方程得到增量 $\delta a_l$，将其加到当前近似值上，即可得到下一个近似值。最速下降公式则转化为： $\delta a_l = constant \times \beta_l$ ### 2.4 二阶导数的处理 Hessian 矩阵的分量 $\alpha_{kl}$ 既依赖于基函数关于其参数的一阶导数，也依赖于二阶导数。在某些情况下，可以忽略二阶导数项。二阶导数项出现是因为梯度已经依赖于 $\frac{\partial y}{\partial a_k}$，所以下一次求导必然包含涉及 $\frac{\partial^2 y}{\partial a_l\partial a_k}$ 的项。当二阶导数项为零（如在线性情况下），或者与一阶导数项相比足够小可以忽略时，或者在实际中该项乘以的 $[y_i - y(x_i; a)]$ 项在求和时相互抵消时，可以忽略二阶导数项。如果模型拟合不佳或存在异常点，包含二阶导数项可能会导致不稳定。因此，从这里开始，我们将 $\alpha_{kl}$ 定义为： $\alpha_{kl} = \sum_{i = 1}^{N} \frac{1}{\sigma_i^2} \left[\frac{\partial y(x_i; a)}{\partial a_k} \frac{\partial y(x_i; a)}{\partial a_l}\right]$ 需要注意的是，对 $[\alpha]$ 进行微小（甚至重大）的调整不会影响最终达到的参数集 $a$，只会影响迭代的路径。$\chi^2$ 最小值处 $\beta_k = 0$ 的条件与 $[\alpha]$ 的定义无关。 ## 3. Levenberg - Marquardt 方法 ### 3.1 方法原理 Levenberg - Marquardt 方法是一种在逆 Hessian 方法和最速下降方法之间平滑过渡的优雅方法。在远离最小值时使用最速下降方法，随着接近最小值，逐渐过渡到逆 Hessian 方法。该方法基于两个重要的见解： - **尺度确定**：考虑最速下降公式中的“常数”，梯度本身无法提供该常数的信息，它只告诉我们斜率，而不告诉我们斜率延伸的距离。Hessian 矩阵的分量虽然不能精确使用，但可以提供问题量级尺度的一些信息。$\chi^2$ 是无量纲的，而 $\beta_k$ 的维度是 $1/a_k$。$\beta_k$ 和 $\delta a_k$ 之间的比例常数必须具有 $a_k^2$ 的维度，而 $[\alpha]$ 中唯一具有此维度的明显量是对角元素的倒数 $1/\alpha_{kk}$。为了避免尺度过大，我们将常数除以一个无量纲的调节因子 $\lambda$，即： $\delta a_l = \frac{1}{\lambda\alpha_{ll}} \beta_l$ 或 $\lambda \alpha_{ll} \delta a_l = \beta_l$ - **方法融合**：通过定义一个新矩阵 $\alpha'$： $\alpha'_{jj} \equiv \alpha_{jj}(1 + \lambda)$ $\alpha'_{jk} \equiv \alpha_{jk}$（$j \neq k$） 将最速下降方法和逆 Hessian 方法的公式统一为： $\sum_{l = 1}^{M} \alpha'_{kl} \delta a_l = \beta_k$ 当 $\la