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

最小二乘定义

一般情况下，最小二乘问题求的是使某一函数局部最小的向量 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)

其中 t_i 是等间距的。在此问题中，向量 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)‖ 可能很小，因为通常的做法是设置现实可行的目标轨迹。虽然可以使用一般的无约束最小化方法来最小化公式 7中的函数（如无约束优化的基础知识中所述），但问题的某些特征经常可以被用来提高求解过程的迭代效率。公式 7 的梯度和 Hessian 矩阵具有特殊结构。

将 F(x) 的 m×n Jacobian 矩阵表示为 J(x)，f(x) 的梯度向量表示为 G(x)，f(x) 的 Hessian 矩阵表示为 H(x)，并将每个 F_i(x) 的 Hessian 矩阵表示为 D_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)‖ 很小时，一个高效方法是使用高斯-牛顿方向作为优化过程的基础。

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

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

(11)

当 Q(x) = 0 时，从该方法推导出的方向等效于牛顿方向。算法可将搜索方向 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 的幅值和方向，而 diag(A) 表示 A 中对角项的矩阵。将选项 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 方法使用的搜索方向综合了高斯-牛顿方向和最陡下降方向。

Levenberg-Marquardt 方法的另一项优势体现在当 Jacobian 矩阵 J 秩亏时。在这种情况下，高斯-牛顿法可能会有数值问题，因为公式 11 中的最小化问题不适定。相反，Levenberg-Marquardt 方法在每次迭代中都有满秩，因此可避免这些问题。

下图显示 Levenberg-Marquardt 方法在最小化 Rosenbrock 函数时的迭代，这是众所周知的以最小二乘法形式出现的最小化难题。

对 Rosenbrock 函数使用 Levenberg-Marquardt 方法

有关此图的更完整说明，包括生成迭代点的脚本，请参阅香蕉函数的最小化。

Levenberg-Marquardt 方法中的边界约束

当问题包含边界约束时，lsqcurvefit 和 lsqnonlin 会修改 Levenberg-Marquardt 迭代。如果建议的迭代点 x 位于边界之外，算法会将该步投影到最近的可行点上。换句话说，在 P 定义为将不可行点投影到可行域上的投影算子的情况下，算法会将建议的点 x 修改为 P(x)。根据定义，投影算子 P 根据以下条件独立地对每个分量 x_i 进行运算

$P (x_{i}) = {\begin{array}{l} {lb}_{i} & if x_{i} < {lb}_{i} \\ {ub}_{i} & if x_{i} > {ub}_{i} \\ x_{i} & otherwise \end{array}$

或者采用以下等效命令，

$P (x_{i}) = \min (\max (x_{i}, {lb}_{i}), {ub}_{i}) .$

该算法修改一阶最优性度量的停止条件。假设 x 是建议的迭代点。在无约束情况下，停止条件为

{‖ \nabla f (x) ‖}_{\infty} \leq tol,

(14)

其中 tol 是最优性容差值，为 1e-4*FunctionTolerance。在有界情况下，停止条件是

{‖ x - P (x - \nabla f (x)) ‖}_{\infty}^{2} \leq tol {‖ \nabla f (x) ‖}_{\infty} .

(15)

要理解此条件，首先要注意，如果 x 在可行域的内部，则算子 P 不起作用。因此，终止条件变为

${‖ x - P (x - \nabla f (x)) ‖}_{\infty}^{2} = {‖ \nabla f (x) ‖}_{\infty}^{2} \leq tol {‖ \nabla f (x) ‖}_{\infty},$

这与原始的无约束停止条件 ${‖ \nabla f (x) ‖}_{\infty} \leq tol$ 相同。如果边界约束处于活动状态，这意味着 x – ∇f(x) 不可行，则在算法应停止的点上，边界上某点的梯度垂直于边界。因此，点 x 等于最陡下降步长的投影 P(x – ∇f(x))，如下图所示。