在 Matlab 中优化大量的 fsolve 调用
Optimizing huge amounts of calls of fsolve in Matlab
我正在使用 MATLAB 2016b 中的 fsolve() 为约十亿体素的数据集中的每个体素求解一对非线性方程。
我已经完成了我所知道的所有 'easy' 优化。内存本地化没问题，我使用的是 parfor，方程式的数字形式相当简单。积分的所有不连续性都被馈送到 integral()。我使用的 Levenberg-Marquardt 算法具有良好的起始值和合适的起始阻尼常数，它平均收敛 6 次迭代。
我现在每个体素大约 6 毫秒，这很好，但还不够好。我需要一个数量级的减少才能使该技术可行。在开始提高准确性之前，我只能考虑改进几件事：
等式中的样条用于复杂等式的快速采样。每个方程有两个，一个在'complicated nonlinear equation'里面。它们表示两个等式，一个具有大量项但平滑且没有间断点，另一个近似于从频谱中绘制的直方图。我按照编辑的建议使用 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 建议的 'parametrization' 输入变量。我有一种挥之不去的感觉，为每个体素定义一个匿名函数在这一点上占用了大量时间。是这样吗，反复定义匿名函数的开销是不是比较大？我有什么方法可以向量化对 fsolve 的调用吗？
有两个输入变量不断变化，所有其他输入变量保持不变。我需要为每个输入对求解一个方程对，所以我不能把它变成一个庞大的系统并立即求解。除了 fsolve 之外，我还有其他选择来求解非线性方程对吗？
如果不是，一些静态输入相当大。有没有办法使用 MATLAB 的 persistent 将输入保持为持久变量，这会提高性能吗？我只看到了如何加载持久变量的示例，我怎样才能使它们只被输入一次，并且将来的函数调用可以避免大输入的假定大开销？
编辑：
完整形式的原始方程式如下：
其中：
和：
其他一切都已知，求解x_1和x_2。 f_KN 由样条近似。 S_low (E) 和 S_high(E) 是样条曲线，它们来自的直方图如下所示：
所以，我想到了一些事情：
查找table
因为函数中的积分不依赖于 x 以外的任何参数，您可以根据它们进行简单的二维查找 table：
% 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 设置为内存允许的最高值。
然后您将这些查找 table 作为参数传递给 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。
Coarse/Fine
如果查找 table 方法由于某种原因（内存限制，以防 x 的可能范围太大）而无法实现，您可以做另一件事：拆分整个过程分为两部分，我称之为 coarse 和 fine。
粗略 方法的目的是真正 快速改进您的初始估计，但可能不是那么准确。到目前为止近似积分的最快方法是通过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));
运行 fsolve 像这样对所有 I_low 和 I_high 将改进你的初始估计 x0 可能接近 "actual"收敛。
或者，您可以使用 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），但随后从粗略步骤开始更新初始估计。
这里的全部objective就是尝试将所需的迭代总数从大约6次减少到2次以下。也许你甚至会发现trapz方法已经提供了足够的准确性，使整个 fine 步骤变得不必要。
向量化
上面概述的粗略步骤中的矩形方法很容易向量化：
% (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 应该能够相当轻松地将所有这些切分给池中的所有工作人员。
同样，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 is the derivative w.r.t. x_1, and 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 and pass this on to fsolve via the 'SpecifyObjectiveGradient' option (via optimoptions)。 'CheckGradients' 选项在那里会派上用场。
因为 fsolve 通常花费大部分时间通过有限差分计算雅可比行列式，手动计算它的值通常会大大加快算法速度 .
会更快，因为
fsolve 不必进行额外的函数评估即可进行有限差分- 由于雅可比矩阵的精度提高，收敛速度会提高
特别是如果您像上面那样使用矩形方法或 trapz，您可以重用您已经为函数值本身完成的许多计算，也就是说，甚至更多 speed-up。
罗迪的回答是正确的。提供雅可比行列式是最大的单一因素。尤其是矢量化版本，提供和不提供 Jacobian 的速度有 3 个数量级的差异。
我在网上找不到关于这个主题的信息，所以我会在这里拼写以供将来参考：可以用 fsolve() 向量化独立的平行方程并获得很大的收益。
我还做了一些内联 fsolve() 的工作。在提供了雅可比行列式并对方程式进行了更智能的处理之后，我的代码的串行版本主要是开销为每个体素 ~1*10^-3 秒。在那一点上，函数内部的大部分时间都花在传递一个 options -struct 和创建 error-messages 上，这些 error-messages 永远不会发送 + 许多未使用的东西假定用于优化函数内的其他优化函数（levenberg-marquardt为了我）。我成功地扼杀了函数 fsolve 和它调用的一些函数，将我机器上每个体素的时间减少到 ~1*10^-4s。因此，如果您坚持使用串行实现，例如因为必须依赖以前的结果，所以内联 fsolve() 很可能会得到很好的结果。
矢量化版本在我的案例中提供了最好的结果，每个体素大约 5*10^-5 秒。
我正在使用 MATLAB 2016b 中的 fsolve() 为约十亿体素的数据集中的每个体素求解一对非线性方程。
我已经完成了我所知道的所有 'easy' 优化。内存本地化没问题，我使用的是 parfor，方程式的数字形式相当简单。积分的所有不连续性都被馈送到 integral()。我使用的 Levenberg-Marquardt 算法具有良好的起始值和合适的起始阻尼常数，它平均收敛 6 次迭代。
我现在每个体素大约 6 毫秒，这很好，但还不够好。我需要一个数量级的减少才能使该技术可行。在开始提高准确性之前，我只能考虑改进几件事：
等式中的样条用于复杂等式的快速采样。每个方程有两个，一个在'complicated nonlinear equation'里面。它们表示两个等式，一个具有大量项但平滑且没有间断点，另一个近似于从频谱中绘制的直方图。我按照编辑的建议使用 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 建议的 'parametrization' 输入变量。我有一种挥之不去的感觉，为每个体素定义一个匿名函数在这一点上占用了大量时间。是这样吗，反复定义匿名函数的开销是不是比较大？我有什么方法可以向量化对 fsolve 的调用吗？
有两个输入变量不断变化，所有其他输入变量保持不变。我需要为每个输入对求解一个方程对，所以我不能把它变成一个庞大的系统并立即求解。除了 fsolve 之外，我还有其他选择来求解非线性方程对吗？
如果不是，一些静态输入相当大。有没有办法使用 MATLAB 的 persistent 将输入保持为持久变量，这会提高性能吗？我只看到了如何加载持久变量的示例，我怎样才能使它们只被输入一次，并且将来的函数调用可以避免大输入的假定大开销？
编辑：
完整形式的原始方程式如下：
其中：
和：
其他一切都已知，求解x_1和x_2。 f_KN 由样条近似。 S_low (E) 和 S_high(E) 是样条曲线，它们来自的直方图如下所示：
所以，我想到了一些事情：
查找table
因为函数中的积分不依赖于 x 以外的任何参数，您可以根据它们进行简单的二维查找 table：
% 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 设置为内存允许的最高值。
然后您将这些查找 table 作为参数传递给 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。
Coarse/Fine
如果查找 table 方法由于某种原因（内存限制，以防 x 的可能范围太大）而无法实现，您可以做另一件事：拆分整个过程分为两部分，我称之为 coarse 和 fine。
粗略 方法的目的是真正 快速改进您的初始估计，但可能不是那么准确。到目前为止近似积分的最快方法是通过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));
运行 fsolve 像这样对所有 I_low 和 I_high 将改进你的初始估计 x0 可能接近 "actual"收敛。
或者，您可以使用 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），但随后从粗略步骤开始更新初始估计。
这里的全部objective就是尝试将所需的迭代总数从大约6次减少到2次以下。也许你甚至会发现trapz方法已经提供了足够的准确性，使整个 fine 步骤变得不必要。
向量化
上面概述的粗略步骤中的矩形方法很容易向量化：
% (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 应该能够相当轻松地将所有这些切分给池中的所有工作人员。
同样，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 is the derivative w.r.t. x_1, and 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 and pass this on to fsolve via the 'SpecifyObjectiveGradient' option (via optimoptions)。 'CheckGradients' 选项在那里会派上用场。
因为 fsolve 通常花费大部分时间通过有限差分计算雅可比行列式，手动计算它的值通常会大大加快算法速度 .
会更快，因为
fsolve不必进行额外的函数评估即可进行有限差分- 由于雅可比矩阵的精度提高，收敛速度会提高
特别是如果您像上面那样使用矩形方法或 trapz，您可以重用您已经为函数值本身完成的许多计算，也就是说，甚至更多 speed-up。
罗迪的回答是正确的。提供雅可比行列式是最大的单一因素。尤其是矢量化版本，提供和不提供 Jacobian 的速度有 3 个数量级的差异。
我在网上找不到关于这个主题的信息，所以我会在这里拼写以供将来参考：可以用 fsolve() 向量化独立的平行方程并获得很大的收益。
我还做了一些内联 fsolve() 的工作。在提供了雅可比行列式并对方程式进行了更智能的处理之后，我的代码的串行版本主要是开销为每个体素 ~1*10^-3 秒。在那一点上，函数内部的大部分时间都花在传递一个 options -struct 和创建 error-messages 上，这些 error-messages 永远不会发送 + 许多未使用的东西假定用于优化函数内的其他优化函数（levenberg-marquardt为了我）。我成功地扼杀了函数 fsolve 和它调用的一些函数，将我机器上每个体素的时间减少到 ~1*10^-4s。因此，如果您坚持使用串行实现，例如因为必须依赖以前的结果，所以内联 fsolve() 很可能会得到很好的结果。
矢量化版本在我的案例中提供了最好的结果，每个体素大约 5*10^-5 秒。