最小二乘（模型拟合）算法

本页对应的英文页面已更新，但尚未翻译。若要查看最新内容，请点击此处访问英文页面。

最小二乘定义

一般情况下，最小二乘问题求的是使某一函数局部最小的向量 x，函数具有平方和的形式，求解可能需要满足一定的约束：

$\min_{x} {‖ F (x) ‖}_{2}^{2} = \min_{x} \sum_{i} F_{i}^{2} (x)$

满足 A·x ≤ b, Aeq·x = beq, lb ≤ x ≤ ub。

Optimization Toolbox™ 提供几种求解器，适用于各种类型的 F(x) 和各种类型的约束：

求解器	F(x)	约束
`mldivide`	C·x – d	无
`lsqnonneg`	C·x – d	x ≥ 0
`lsqlin`	C·x – d	边界、线性
`lsqnonlin`	一般形式的 F(x)	边界
`lsqcurvefit`	F(x, xdata) – ydata	边界

除了 mldivide 中使用的算法外，Optimization Toolbox 求解器中还有五种最小二乘算法：

lsqlin 内点
lsqlin 活动集
信赖域反射（非线性或线性最小二乘）
Levenberg-Marquardt（非线性最小二乘）
lsqnonneg 使用的算法

除 lsqlin 活动集外，所有算法均为大规模算法；请参阅大规模算法与中等规模算法。有关非线性最小二乘方法的一般概况，请参阅 Dennis 的论著·[8]。Levenberg-Marquardt 方法的具体细节可见于 Moré 的论著 [28]。

线性最小二乘：内点或活动集

lsqlin 'interior-point' 算法使用 interior-point-convex quadprog 算法，lsqlin 'active-set' 算法使用活动集 quadprog 算法。quadprog 问题定义用于最小化以下二次函数

$\min_{x} \frac{1}{2} x^{T} H x + c^{T} x$

且需要满足线性约束和边界约束。lsqlin 函数在满足线性约束和边界约束的前提下最小化向量 Cx – d 的 2-范数平方。换句话说，lsqlin 最小化下式

$\begin{matrix} {‖ C x - d ‖}_{2}^{2} = {(C x - d)}^{T} (C x - d) \\ = (x^{T} C^{T} - d^{T}) (C x - d) \\ = (x^{T} C^{T} C x) - (x^{T} C^{T} d - d^{T} C x) + d^{T} d \\ = \frac{1}{2} x^{T} (2 C^{T} C) x + {(- 2 C^{T} d)}^{T} x + d^{T} d . \end{matrix}$

将上式的 2C^TC 替换为 H 矩阵，并将 (–2C^Td) 替换为 c 向量，即符合 quadprog 框架。（加法项 d^Td 对最小值的位置没有影响。）在重新表示 lsqlin 问题后，quadprog 会计算解。

注意

quadprog 'interior-point-convex' 算法有两种代码路径。当 Hessian 矩阵 H 是由双精度值组成的普通（满）矩阵时，它采用一种代码路径；当 H 是稀疏矩阵时，它采用另一种代码路径。有关稀疏数据类型的详细信息，请参阅稀疏矩阵。通常，对于非零项相对较少的大型问题，该算法在 H 指定为 sparse 时执行更快。同样，对于较小或相对稠密的问题，该算法在 H 指定为 full 时执行更快。

信赖域反射最小二乘

信赖域反射最小二乘算法

Optimization Toolbox 求解器中使用的许多方法都基于信赖域，这是一个简单而功能强大的优化概念。

要理解信赖域优化方法，请考虑无约束最小化问题，最小化 f(x)，该函数接受向量参数并返回标量。假设您现在位于 n 维空间中的点 x 处，并且您要寻求改进，即移至函数值较低的点。基本思路是用较简单的函数 q 来逼近 f，该函数需能充分反映函数 f 在点 x 的邻域 N 中的行为。此邻域是信赖域。试探步 s 是通过在 N 上进行最小化（或近似最小化）来计算的。以下是信赖域子问题，

\min_{s} {q (s), s \in N} .

(1)

如果 f(x + s) < f(x)，当前点更新为 x + s；否则，当前点保持不变，信赖域 N 缩小，算法再次计算试探步。

在定义特定信赖域方法以最小化 f(x) 的过程中，关键问题是如何选择和计算逼近 q（在当前点 x 上定义）、如何选择和修改信赖域 N，以及如何准确求解信赖域子问题。本节重点讨论无约束问题。后面的章节讨论由于变量约束的存在而带来的额外复杂性。

在标准信赖域方法 ([48]) 中，二次逼近 q 由 F 在 x 处的泰勒逼近的前两项定义；邻域 N 通常是球形或椭圆形。以数学语言表述，信赖域子问题通常写作

\min {\frac{1}{2} s^{T} H s + s^{T} g such that ‖ D s ‖ \leq Δ},

(2)

其中，g 是 f 在当前点 x 处的梯度，H 是 Hessian 矩阵（二阶导数的对称矩阵），D 是对角缩放矩阵，Δ 是正标量，∥ . ∥ 是 2-范数。存在求解公式 2 的好算法（请参阅[48]）；此类算法通常涉及计算 H 的所有特征值，并将牛顿法应用于以下久期方程

$\frac{1}{Δ} - \frac{1}{‖ s ‖} = 0.$

此类算法提供公式 2 的精确解。但是，它们要耗费与 H 的几个分解成比例的时间。因此，对于信赖域问题，需要采取另一种方法。文献（[42] 和 [50]）中提出了几种基于公式 2 的逼近和启发式方法建议。Optimization Toolbox 求解器采用的逼近方法是将信赖域子问题限制在二维子空间 S 内（[39] 和 [42]）。一旦计算出子空间 S，即使需要完整的特征值/特征向量信息，求解公式 2 的工作量也不大（因为在子空间中，问题只是二维的）。现在的主要工作已转移到子空间的确定上。

二维子空间 S 是借助下述预条件共轭梯度法确定的。求解器将 S 定义为由 s₁ 和 s₂ 确定的线性空间，其中 s₁ 是梯度 g 的方向，s₂ 是近似牛顿方向，即下式的解

H \cdot s_{2} = - g,

(3)

或是负曲率的方向，

s_{2}^{T} \cdot H \cdot s_{2} < 0.

(4)

以此种方式选择 S 背后的理念是强制全局收敛（通过最陡下降方向或负曲率方向）并实现快速局部收敛（通过牛顿步，如果它存在）。

现在，我们可以很容易地给出基于信赖域的无约束最小化的大致框架：

构造二维信赖域子问题。
求解公式 2 以确定试探步 s。
如果 f(x + s) < f(x)，则 x = x + s。
调整 Δ。

重复这四个步骤，直到收敛。信赖域维度 Δ 根据标准规则进行调整。具体来说，它会在试探步不被接受（即 f(x + s) ≥ f(x)）时减小。有关这方面的讨论，请参阅[46] 和 [49]。

Optimization Toolbox 求解器用专用函数处理 f 的一些重要特例：非线性最小二乘、二次函数和线性最小二乘。然而，其底层算法思路与一般情况相同。这些特例将在后面的章节中讨论。

大规模非线性最小二乘

f(x) 的一个重要特例是非线性最小二乘问题

\min_{x} \sum_{i} f_{i}^{2} (x) = \min_{x} {‖ F (x) ‖}_{2}^{2},

(5)

其中，F(x) 是向量值函数，F(x) 的分量 i 等于 f_i(x)。用于求解此问题的基本方法与非线性最小化信赖域方法中所述的一般情况相同。不过，我们可以充分利用非线性最小二乘问题的结构来提高效率。具体来说，就是使用近似高斯-牛顿方向，即下式的解 s

\min {‖ J s + F ‖}_{2}^{2},

(6)

（其中 J 是 F(x) 的 Jacobian 矩阵）来定义二维子空间 S。此处不使用分量函数 f_i(x) 的二阶导数。

每次迭代中，算法都使用预条件共轭梯度法来近似求解标准方程，即，

$J^{T} J s = - J^{T} F,$

（虽然并不显式构造标准方程）。

大规模线性最小二乘

在这种情况下，要求解的函数 f(x) 是

$f (x) = {‖ C x + d ‖}_{2}^{2},$

可能需要满足线性约束。算法生成严格可行的迭代，在极限情况下收敛于局部解。每次迭代都涉及大型线性方程组的近似解（方程组阶数为 n，其中 n 是 x 的长度）。迭代矩阵具有矩阵 C 的结构。具体来说，算法使用预条件共轭梯度法来近似求解标准方程，即，

$C^{T} C x = - C^{T} d,$

（虽然并不显式构造标准方程）。

此处使用子空间信赖域方法来确定搜索方向。然而，它不采用非线性最小化情形中的做法，将步限制为（可能）一个反射步，而是采用二次情形中的做法，在每次迭代中进行分段反射线搜索。请参阅[45] 了解线搜索的详细信息。最终，这些线性方程组表示一种牛顿法，该方法在解处捕获一阶最优性条件，实现快速局部收敛。

Jacobian 矩阵乘法函数. lsqlin 无需显式构造矩阵 C 即可求解线性约束最小二乘问题。它使用 Jacobian 矩阵乘法函数 jmfun，

W = jmfun(Jinfo,Y,flag)

（该函数由您提供。）该函数必须计算矩阵 Y 的下列乘积：

如果 flag == 0，则 W = C'*(C*Y)。
如果 flag > 0，则 W = C*Y。
如果 flag < 0，则 W = C'*Y。

如果 C 很大，但其包含的结构足以让您在不显式构造 C 的情况下编写 jmfun，这种函数将非常有用。有关示例，请参阅Jacobian Multiply Function with Linear Least Squares。

Levenberg-Marquardt 方法

最小二乘问题对平方和函数 f(x) 进行最小化。

\min_{x} f (x) = {‖ F (x) ‖}_{2}^{2} = \sum_{i} F_{i}^{2} (x) .

(7)

这种类型的问题出现在大量的实际应用中，特别是在对数据进行模型函数拟合（即非线性参数估计）的过程中。它们在下面这类控制系统中也很普遍：对于向量 x 和标量 t，您希望输出 y(x,t) 遵循一定的连续模型轨迹 φ(t)。此问题可以表达为

\min_{x \in ℜ^{n}} \int_{t_{1}}^{t_{2}} {(y (x, t) - φ (t))}^{2} d t,

(8)

其中 y(x,t) 和 φ(t) 是标量函数。

当使用合适的求积公式对积分进行离散化时，上述问题可表示为最小二乘问题：

\min_{x \in ℜ^{n}} f (x) = \sum_{i = 1}^{m} {(y (x, t_{i}) - φ (t_{i}))}^{2},

(9)

其中 $\bar{y}$ 和 $\bar{φ}$ 包含求积方案的权重信息。请注意，在此问题中，向量 F(x) 为

$F (x) = [\begin{matrix} y (x, t_{1}) - φ (t_{1}) \\ y (x, t_{2}) - φ (t_{2}) \\ ... \\ y (x, t_{m}) - φ (t_{m}) \end{matrix}] .$

在这类问题中，最优情形下的残差 ∥F(x)∥ 可能很小，因为通常的做法是设置现实可行的目标轨迹。虽然最小二乘中的函数可以使用一般的无约束最小化方法来最小化（如无约束优化的基础知识中所述），但通常可以利用问题的某些特性来提高求解过程的迭代效率。最小二乘的梯度和 Hessian 矩阵具有特殊结构。

将 F(x) 的 m×n Jacobian 矩阵表示为 J(x)，f(x) 的梯度向量表示为 G(x)，f(x) 的 Hessian 矩阵表示为 H(x)，并将每个 F_i(x) 的 Hessian 矩阵表示为 H_i(x)，您将得到

\begin{array}{l} G (x) = 2 J (^{x) T} F (x) \\ H (x) = 2 J (^{x) T} J (x) + 2 Q (x), \end{array}

(10)

其中

$Q (x) = \sum_{i = 1}^{m} F_{i} (x) \cdot D_{i} (x) .$

矩阵 Q(x) 的属性是：当残差 ∥F(x)∥ 随着 x_k 接近解而趋于零时，Q(x) 也趋于零。因此，当解处的 ∥F(x)∥ 很小时，一个非常有效的方法是使用高斯-牛顿方向作为优化过程的基础。

在高斯-牛顿法中，每次主迭代 d_k 都会得到一个搜索方向 k，该方向是线性最小二乘问题的解：

\min_{d_{k} \in ℜ^{n}} {‖ J (x_{k}) d_{k} + F (x_{k}) ‖}_{2}^{2} .

(11)

当 Q(x) 的项可以忽略时，从该方法推导出的方向等效于牛顿方向。搜索方向 d_k 可用作线搜索策略的一部分，以确保在每次迭代中，函数 f(x) 会减小。

当二阶项 Q(x) 很显著时，高斯-牛顿法经常遇到问题。克服此问题的一种方法是 Levenberg-Marquardt 方法。

Levenberg-Marquardt [25] 和 [27] 方法使用的搜索方向是线性方程组的解

(J {(x_{k})}^{T} J (x_{k}) + λ_{k} I) d_{k} = - J {(x_{k})}^{T} F (x_{k}),

(12)

也可以是以下方程的解

(J {(x_{k})}^{T} J (x_{k}) + λ_{k} diag (J {(x_{k})}^{T} J (x_{k}))) d_{k} = - J {(x_{k})}^{T} F (x_{k}),

(13)

其中标量 λ_k 控制 d_k 的模和方向。将选项 ScaleProblem 设置为 'none' 可选择公式 12，将 ScaleProblem 设置为 'Jacobian' 可选择公式 13。

您可以使用 InitDamping 选项设置参数 λ₀ 的初始值。个别情况下，此选项的默认值 0.01 可能不合适。如果您发现 Levenberg-Marquardt 算法的初始进展甚微，请尝试将 InitDamping 设置为不同于默认值的值（例如 1e2）。

当 λ_k 为零时，方向 d_k 与高斯-牛顿法的方向相同。随着 λ_k 趋向于无穷大，d_k 趋向于最陡下降方向，模趋向于零。这意味着，对于一些足够大的 λ_k，项 F(x_k + d_k) < F(x_k) 成立。因此，即使遇到会限制高斯-牛顿法效率的二阶项，也可以控制项 λ_k 以确保下降。当步成功（给出较低的函数值）时，该算法设置 λ_k+1 = λ_kλk/10。当步不成功时，该算法设置 λ_k+1 = λ_k*10。

Levenberg-Marquardt 算法在内部使用最优性容差（停止条件）1e-4 乘以函数容差。

因此，Levenberg-Marquardt 方法使用的搜索方向综合了高斯-牛顿方向和最陡下降方向。图 12-1：对 Rosenbrock 函数使用 Levenberg-Marquardt 方法就是一个例子。Rosenbrock 函数的解在 90 次函数计算后收敛，而高斯-牛顿法为 48 次。效率较低的部分原因是，当解处的残差为零时，高斯-牛顿法通常更高效。然而，此类信息并不能总能事先知晓，而 Levenberg-Marquardt 方法更高的稳健性弥补了它偶尔发生的相对低效。

图 12-1：对 Rosenbrock 函数使用 Levenberg-Marquardt 方法