在 Matlab 中优化 fsolve 的大量调用

Question

我正在使用 MATLAB 2016b 中的fsolve() 为约十亿体素的数据集中的每个体素求解一对非线性方程。

我已经完成了我所知道的所有“简单”优化。内存本地化没问题，我使用的是parfor，方程的数值相当简单。积分的所有不连续性都馈送到integral()。我正在使用具有良好起始值和合适起始阻尼常数的 Levenberg-Marquardt 算法，它平均收敛 6 次迭代。

我现在每个体素约为 6 毫秒，这很好，但还不够好。我需要减少一个数量级才能使该技术可行。在开始敲定准确性之前，我只能想到几件事要改进：

方程中的样条用于对复杂方程进行快速采样。每个方程有两个，一个在“复杂的非线性方程”内。它们代表两个方程，一个具有大量项但平滑且没有不连续性，另一个近似于从频谱中绘制的直方图。正如编辑建议的那样，我正在使用griddedInterpolant()。

有没有更快的方法从预先计算的分布中采样点？

parfor i=1:numel(I1)
    sols = fsolve(@(x) equationPair(x, input1, input2, ... 
         6 static inputs, fsolve options)
    output1(i) = sols(1); output2(i) = sols(2)
end

调用fsolve 时，我使用 Mathworks 建议的“参数化”来输入变量。我有一种挥之不去的感觉，在这一点上为每个体素定义一个匿名函数会占用大量时间。这是真的吗，一次又一次地定义匿名函数有比较大的开销吗？我有什么方法可以矢量化对fsolve 的调用吗？

有两个输入变量不断变化，所有其他输入变量都保持不变。我需要为每个输入对求解一个方程对，因此我无法将其变成一个庞大的系统并立即求解。除了fsolve 之外，我还有其他选择来求解非线性方程组吗？

如果不是，一些静态输入相当大。有没有办法使用 MATLAB 的persistent 将输入保持为持久变量，这会提高性能吗？我只看到了如何加载持久变量的示例，我怎样才能使它们只被输入一次，并且将来的函数调用可以避免大量输入的假定较大开销？

编辑：

完整形式的原始方程如下所示：

地点：

和：

其他一切都是已知的，求解 x_1 和 x_2。 f_KN 由样条曲线近似。 S_low (E) 和 S_high(E) 是样条曲线，它们的直方图如下所示：

Answer 1

罗迪的答案是正确的。提供雅可比矩阵是最大的一个因素。尤其是矢量化版本，提供雅可比行列式的速度差异有 3 个数量级，而没有提供。

我在网上找不到有关此主题的信息时遇到了麻烦，因此我将在此处将其拼写出来以供将来参考：可以使用 fsolve() 对独立的平行方程进行矢量化，并获得很大的收益。

我还做了一些内联 fsolve() 的工作。在提供雅可比行列式并对方程更聪明之后，我的代码的串行版本的开销主要是每个体素 ~1*10^-3 秒。那时，函数内部的大部分时间都花在传递 options -struct 并创建永远不会发送的错误消息 + 很多未使用的东西，假设用于优化函数中的其他优化函数（对我来说是 levenberg-marquardt）。我成功地处理了函数 fsolve 和它调用的一些函数，将我机器上每个体素的时间减少到 ~1*10^-4s。因此，如果您坚持使用串行实现，例如由于必须依赖以前的结果，因此很有可能内联 fsolve() 并获得良好的结果。

在我的例子中，矢量化版本提供了最好的结果，每个体素约为 5*10^-5 秒。

Answer 2

所以，我想到了一些事情：

查找表

因为你的函数中的积分不依赖于x以外的任何参数，你可以用它们制作一个简单的二维查找表：

% assuming simple (square) range here, adjust as needed
[x1,x2]  = meshgrid( linspace(0, xmax, N) );

LUT_high = zeros(size(x1));
LUT_low  = zeros(size(x1));

for ii = 1:N        

    LUT_high(:,ii) = integral(@(E) Fhi(E, x1(1,ii), x2(:,ii)), ...
                              0, E_high, ...
                              'ArrayValued', true);

    LUT_low(:,ii) = integral(@(E) Flo(E, x1(1,ii), x2(:,ii)), ...
                             0, E_low, ...
                             'ArrayValued', true);

end 

其中Fhi 和Flo 是计算这些积分的辅助函数，在此示例中使用标量x1 和向量x2 进行矢量化。将N 设置为内存允许的最高值。

然后将这些查找表作为参数传递给equationPair()（它允许parfor 分发数据）。然后在equationPair()中使用interp2：

F(1) = I_high - interp2(x1,x2,LUT_high, x(1), x(2));
F(2) = I_low  - interp2(x1,x2,LUT_low , x(1), x(2));

因此，您不必每次重新计算整个积分，而是在x 的预期范围内对其进行一次评估，然后重复使用结果。

您可以指定使用的插值方法，默认为linear。如果您真的关心准确性，请指定cubic。

粗/细

如果由于某种原因无法使用查找表方法（内存限制，以防x 的可能范围太大），您可以做另一件事：将整个过程分成两部分，我'将调用 coarse 和 fine。

coarse 方法的目的是真正快速地改进您的初始估计，但可能不是那么准确。 到目前为止，近似积分的最快方法是通过rectangle method：

• 不用spline近似S，只需使用原始列表数据（所以S_high/low = [S_high/low@E0, S_high/low@E1, ..., S_high/low@E_high/low]

• 在 S 数据（E0、E1、...）使用的 E 的相同值下，计算 x 处的指数：

  Elo = linspace(0, E_low, numel(S_low)).';
  integrand_exp_low = exp(x(1)./Elo.^3 + x(2)*fKN(Elo));
  
  Ehi = linspace(0, E_high, numel(S_high)).';
  integrand_exp_high = exp(x(1)./Ehi.^3 + x(2)*fKN(Ehi));
  

  然后使用矩形方法：

  F(1) = I_low  - (S_low  * Elo) * (Elo(2) - Elo(1));
  F(2) = I_high - (S_high * Ehi) * (Ehi(2) - Ehi(1));
  

像这样对所有I_low 和I_high 运行fsolve 将提高您的初始估计x0 可能接近“实际”收敛的程度。

或者，您可以使用trapz (trapezoidal method) 而不是矩形方法。有点慢，但可能更准确。

请注意，如果(Elo(2) - Elo(1)) == (Ehi(2) - Ehi(1))（步长相等），您可以进一步减少计算次数。在这种情况下，两个被积函数的第一个 N_low 元素是相同的，因此指数的值只会在 N_low + 1 : N_high 元素中不同。那么只需计算integrand_exp_high，并将integrand_exp_low 设置为等于integrand_exp_high 的第一个N_low 元素。

fine 方法然后使用您的原始实现（使用实际的integral()s），然后从粗略步骤的更新初始估计开始。

这里的整个目标是尝试将所需的迭代总数从大约 6 次减少到少于 2 次。也许您甚至会发现 trapz 方法已经提供了足够的准确性，从而呈现整个 很好步骤是不必要的。

矢量化

上述粗略步骤中的矩形方法很容易矢量化：

% (uses R2016b implicit expansion rules)

Elo = linspace(0, E_low, numel(S_low));
integrand_exp_low = exp(x(:,1)./Elo.^3 + x(:,2).*fKN(Elo));

Ehi = linspace(0, E_high, numel(S_high));
integrand_exp_high = exp(x(:,1)./Ehi.^3 + x(:,2).*fKN(Ehi));

F = [I_high_vector - (S_high * integrand_exp_high) * (Ehi(2) - Ehi(1))
     I_low_vector  - (S_low  * integrand_exp_low ) * (Elo(2) - Elo(1))];

trapz 也适用于矩阵；它将整合矩阵中的每一列。

您可以调用equationPair()，然后使用x0 = [x01; x02; ...; x0N]，然后fsolve 将收敛到[x1; x2; ...; xN]，其中N 是体素的数量，每个x0 是1×2（@987654376 @)，所以x0 是N×2。

parfor 应该能够相当轻松地对池中的所有工作人员进行切片。

同样，fine 方法的矢量化也应该是可能的；只需使用'ArrayValued' 选项到integral()，如上所示：

F = [I_high_vector - integral(@(E) S_high(E) .* exp(x(:,1)./E.^3 + x(:,2).*fKN(E)),...
                              0, E_high,...
                              'ArrayValued', true);
    I_low_vector   - integral(@(E) S_low(E) .* exp(x(:,1)./E.^3 + x(:,2).*fKN(E)),...
                              0, E_low,...
                              'ArrayValued', true);
    ];

雅可比

对你的函数求导非常容易。 Here 是 w.r.t 的导数。 x_1 和 here w.r.t. x_2。然后，您的雅可比矩阵必须是 2×2 矩阵

J = [dF(1)/dx(1)  dF(1)/dx(2)
     dF(2)/dx(1)  dF(2)/dx(2)]; 

不要忘记前面的减号 (F = I_hi/lo - g(x) → dF/dx = -dg/dx)

使用上面列出的一种或两种方法，您可以实现一个函数来计算Jacobian matrix，并通过'SpecifyObjectiveGradient' 选项（via optimoptions)）将其传递给fsolve。'CheckGradients' 选项将会出现在那里派上用场。

由于fsolve 通常将其大部分时间用于通过有限差分计算雅可比行列式，因此手动为其手动计算值通常会极大地加速算法。

会更快，因为

1. fsolve 不需要做额外的函数评估来做有限差分
2. 由于雅可比行列精度的提高，收敛速度会提高

特别是如果您使用矩形方法或trapz 像上面那样，您可以重用您已经为函数值本身完成的许多计算，这意味着，甚至可以更快地加速。