方程求解算法

方程求解定义

给定一组（n 个）非线性函数 F_i(x)，其中 n 是向量 x 中的分量数，方程求解的目标是找到使所有 F_i(x) = 0 的向量 x。

fsolve 尝试通过最小化分量的平方和来求解方程组。如果平方和为零，则方程组得到求解。fsolve 有三种算法：

信赖域
信赖域 dogleg
莱文贝格-马夸特

所有算法均为大规模的；请参阅大规模算法与中等规模算法。

fzero 函数求解单个一维方程。

mldivide 函数对一个线性方程组求解。

信赖域算法

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

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

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

如果 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 Δ},

(1)

其中，g 是 f 在当前点 x 处的梯度，H 是黑塞矩阵（二阶导数的对称矩阵），D 是对角缩放矩阵，Δ 是正标量，‖ . ‖ 是 2-范数。为了求解公式 1，一种算法（请参阅[48]）会计算 H 的所有特征值，然后将牛顿法应用于以下久期方程

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

这种算法提供公式 1 的精确解。但是，这需要耗费与 H 的几个分解成比例的时间。因此，处理信赖域问题需要另一种方法。文献（[42] 和 [50]）提出了几种基于公式 1 的逼近和启发式策略。Optimization Toolbox 求解器采用一种逼近方法，将信赖域子问题限制在二维子空间 S 内（[39] 和 [42]）。在求解器计算出子空间 S 后，求解公式 1 的工作量是微不足道的，因为在子空间中，问题只是二维的。现在的主要工作转移到子空间的确定上。

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

$H \cdot s_{2} = - g,$

或是负曲率的方向，

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

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

现在，我们可以很容易地指定基于信赖域方法的无约束最小化过程：

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

求解器重复这四个步骤，直到收敛，从而根据标准规则调整信赖域维度 Δ。具体来说，当 f(x + s) ≥ f(x) 时，求解器不接受试探步，并减小信赖域大小。有关这方面的讨论，请参阅[46] 和 [49]。

Optimization Toolbox 求解器用专用函数处理 f 的重要实例：非线性最小二乘、二次函数和线性最小二乘。然而，其底层算法思路与一般情况相同。

预条件共轭梯度法

求解大型对称正定线性方程组 Hp = –g 的一种常用方法是预条件共轭梯度法 (PCG)。这种迭代方法需要能够计算 H·v 形式的矩阵-向量积，其中 v 是任意向量。对称正定矩阵 M 是 H 的预条件子。也就是说，M = C²，其中 C^–1HC^–1 是良态矩阵或具有聚类特征值的矩阵。

在最小化上下文中，您可以假设黑塞矩阵 H 是对称矩阵。然而，只有在强最小解的邻域内才能保证 H 是正定矩阵。算法 PCG 在遇到负（或零）曲率方向（即 d^THd ≤ 0）时退出。PCG 输出方向 p 要么是负曲率方向，要么是牛顿方程组 Hp = –g 的近似解。无论哪种情况，p 都有助于定义在非线性最小化信赖域方法中讨论的信赖域方法中使用的二维子空间。

信赖域 dogleg 算法

另一种方法是求解线性方程组以找到搜索方向。牛顿法指定求解搜索方向 d_k，满足

J(x_k)d_k = –F(x_k)
x_{k + 1} = x_k + d_k，

其中，J(x_k) 是 n×n 雅可比矩阵

$J (x_{k}) = [\begin{matrix} \nabla F_{1} {(x_{k})}^{T} \\ \nabla F_{2} {(x_{k})}^{T} \\ ⋮ \\ \nabla F_{n} {(x_{k})}^{T} \end{matrix}] .$

牛顿法可能存在问题。J(x_k) 可能是奇异的，在这种情况下，牛顿步 d_k 甚至没有定义。精确的牛顿步 d_k 的计算成本也可能很高昂。此外，如果起点与解的距离较远，牛顿法可能不会收敛。

使用信赖域方法（在非线性最小化信赖域方法中介绍）处理 J(x_k) 为奇异时的情况，并提高在起点与解的距离较远时的稳健性。要使用信赖域策略，您需要评价函数来决定 x_{k + 1} 比 x_k 更好还是更差。一种可能的选择是

$\min_{d} f (d) = \frac{1}{2} F {(x_{k} + d)}^{T} F (x_{k} + d) .$

但 f(d) 的最小值不一定是 F(x) 的根。

牛顿步 d_k 是下式的根

M(x_k + d) = F(x_k) + J(x_k)d，

因此它也是 m(d) 的最小值，其中

\begin{matrix} \min_{d} m (d) = \frac{1}{2} {‖ M (x_{k} + d) ‖}_{2}^{2} = \frac{1}{2} {‖ F (x_{k}) + J (x_{k}) d ‖}_{2}^{2} \\ = \frac{1}{2} F {(x_{k})}^{T} F (x_{k}) + d^{T} J {(x_{k})}^{T} F (x_{k}) + \frac{1}{2} d^{T} J {(x_{k})}^{T} J (x_{k}) d . \end{matrix}

(2)

m(d) 是优于 f(d) 的评价函数选择，因此信赖域子问题是

\min_{d} [\frac{1}{2} F {(x_{k})}^{T} F (x_{k}) + d^{T} J {(x_{k})}^{T} F (x_{k}) + \frac{1}{2} d^{T} J {(x_{k})}^{T} J (x_{k}) d],

(3)

满足 ‖D·d‖ ≤ Δ。您可以使用 dogleg 策略高效求解此子问题。

有关信赖域方法的概述，请参阅 Conn 的论著 [4] 和 Nocedal 的论著 [31]。

信赖域 dogleg 实现

信赖域 dogleg 算法的关键特征是使用鲍威尔 dogleg 过程来计算步 d，从而最小化公式 3。有关详细说明，请参阅鲍威尔的论著·[34]。

该算法基于柯西步（沿最陡下降方向的步）和高斯-牛顿步的凸组合来为 f(x) 构造步 d。柯西步的计算如下

d_C = –αJ(x_k)^TF(x_k)，

其中 α 最小化公式 2。

高斯-牛顿步通过求解下式来计算

J(x_k)·d_GN = –F(x_k)，

在计算中使用 MATLAB^® mldivide（矩阵左除）运算符。

该算法选择步 d，满足

d = d_C + λ(d_GN – d_C)，

其中 λ 是区间 [0,1] 中的最大值，使得 ‖d‖ ≤ Δ。如果 J_k 是（或几乎是）奇异的，则 d 就是柯西方向。

信赖域 dogleg 算法很高效，因为它每次迭代只需进行一次线性求解（用于计算高斯-牛顿步）。此外，该算法比基于高斯-牛顿法的线搜索更稳健。

莱文贝格-马夸特方法

莱文贝格-马夸特算法（[25] 和 [27]）使用的搜索方向是以下线性方程组的解

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

(4)

也可以是以下方程的解

(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}),

(5)

其中标量 λ_k 控制 d_k 的模和方向。将 fsolve 选项 ScaleProblem 设置为 'none' 以使用公式 4，或将此选项设置为 'jacobian' 以使用公式 5。

当 λ_k 为零时，d_k 方向是高斯-牛顿方向。随着 λ_k 趋向于无穷，d_k 趋向于最陡下降方向，模趋向于零。这意味着，对于一些足够大的 λ_k，项 F(x_k + d_k) < F(x_k) 成立。因此，即便有二阶项对高斯-牛顿法的效率造成限制，该算法仍能控制项 λ_k 以确保下降。所以，莱文贝格-马夸特算法使用的搜索方向综合了高斯-牛顿方向和最陡下降方向。有关详细信息，请参阅最小二乘法文档中的莱文贝格-马夸特方法。