matlab - 在Matlab中优化大量的fsolve调用

Question

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



我已经完成了我所知道的所有“轻松”优化。内存本地化还可以，我使用的是parfor，方程式的数值形式非常简单。积分的所有不连续性均馈入integral()。我使用的Levenberg-Marquardt算法具有良好的起始值和适当的起始阻尼常数，平均经过6次迭代即可收敛。

我现在每个体素大约为6ms，这很好，但是还不够好。我需要降低一个数量级才能使这项技术可行。在开始精打细算之前，我只有几件事需要改进：

方程中的样条曲线用于快速采样复杂方程。每个方程有两个，一个在“复杂非线性方程”中。它们代表两个方程，一个方程具有大量项，但平滑且没有间断，而另一个方程则近似于从频谱绘制的直方图。我正在使用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建议的“参数化”来输入变量。我有点a的感觉是，在此时，为每个体素定义一个匿名函数将花费大量时间。是真的吗，一次又一次地定义匿名函数是否有相对较大的开销？我有什么方法可以对 fsolve的呼叫进行矢量化处理？

有两个输入变量保持不变，所有其他输入变量保持不变。我需要为每个输入对求解一个方程对，所以不能使其成为一个庞大的系统并立即求解。除了 fsolve之外，我还有其他选择用于求解非线性方程对吗？

如果不是，则某些静态输入会相当大。有没有一种方法可以使用MATLAB的 persistent将输入保留为持久变量，这会提高性能吗？我只看到了有关如何加载持久变量的示例，如何使它们只能输入一次，而将来的函数调用可以从大输入的假定庞大开销中省下来呢？

编辑：

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





哪里：



和：



解决x_1和x_2的其他所有知识都是已知的。 f_KN由样条曲线近似。 S_low（E）和S_high（E）是样条曲线，它们的直方图如下所示：

Answer 1

因此，我想到了几件事：

查找表

由于函数中的积分不依赖于x以外的任何参数，因此您可以根据它们创建一个简单的2D查找表：

% 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分发数据）。然后只需在 interp2中使用 equationPair()：

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的可能范围太大）无法使用查找表方法，那么您可以做另一件事：将整个过程分为两部分，我称之为粗而细。

粗略方法的目的是真正快速地提高您的初始估计，但可能没有那么精确。到目前为止，最接近该积分的最快方法是通过 rectangle method：


不要使用 S近似 spline，只需使用原始列表数据即可（因此 S_high/low = [S_high/low@E0, S_high/low@E1, ..., S_high/low@E_high/low]
在与 E数据使用的 S值相同的值（ E0， E1，...）下，计算 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));

对所有 fsolve和 I_low像这样运行 I_high将会改善您的初始估计 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设置为等于 N_low的前 integrand_exp_high个元素。

好的方法然后使用您的原始实现（带有实际的 integral()），但随后从粗略步骤的更新后的初始估算开始。

此处的总体目标是尝试使所需的迭代总数从大约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（ [x(1) x(2)]），因此 x0是 N×2。

parfor应该能够轻松地将所有这些内容分割到池中的所有工作人员上。

同样，精细方法的矢量化也应该是可能的。只需对 '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并将其通过 fsolve选项（ via 'SpecifyObjectiveGradient')传递给 optimoptions。 'CheckGradients'选项将派上用场。

因为 fsolve通常花费大量时间通过有限差分来计算Jacobian，所以手动为其计算值通常会极大地加快算法的速度。

它将更快，因为


fsolve不必进行额外的函数求值来进行有限差分
由于雅可比行列式的精度提高，收敛速度将提高


特别是如果您使用上述矩形方法或 trapz，则可以重用已经为函数值本身完成的许多计算，这意味着甚至可以提高速度。