2

我有一个循环，在其中我使用 ppval 从分段多项式样条中评估一组值。插值很容易成为循环中最耗时的部分，我正在寻找一种提高函数效率的方法。

更具体地说，我使用有限差分方案来计算摩擦焊缝中的瞬态温度分布。为此，我需要在每个时间步重新计算材料属性（作为温度和位置的函数）。速率限制因素是这些值的插值。我可以使用另一种有限差分方案（在时域中限制较少），但如果可能的话，我宁愿坚持我所拥有的。

我在下面包含了一个 MWE：

x=0:.1:10;
y=sin(x);
pp=spline(x,y);
tic
for n=1:10000
    x_int=10*rand(1000,1);
    y_int=ppval(pp,x_int);
end
toc

plot(x,y,x_int,y_int,'*') % plot for sanity of data

Elapsed time is 1.265442 seconds.

编辑-我可能应该提到，我对值之间的简单线性插值非常满意，但 interp1 函数比 ppval 慢

x=0:.1:10;
y=sin(x);
tic
for n=1:10000
    x_int=10*rand(1000,1);
    y_int=interp1(x,y,x_int,'linear');
end
toc

plot(x,y,x_int,y_int,'*') % plot for sanity of data

Elapsed time is 1.957256 seconds.
4

2 回答 2

2

interp1比 慢有点令人惊讶ppval，但是快速查看它的源代码，它似乎必须检查许多特殊情况并且必须循环所有点，因为它无法确定步长是否是恒定的。

我没有检查时间，但我想如果你能保证你的表 x 中的步骤是恒定的，并且要插值的值严格在给定范围内，我想你可以大大加快线性插值，所以你不必做任何检查。在这种情况下，线性插值可以转换为一个简单的查找问题，如下所示：

%data to be interpolated, on grid with constant step
x = 0:0.5:10;
y = sin(x);

x_int = 0:0.1:9.9;

%make sure it is interpolation, not extrapolation
assert(all(x(1) <= x_int & x_int < x(end)));

% compute mapping, this can be precomputed for constant grid
slope = (length(x) - 1) / (x(end) - x(1));
offset = 1 - slope*x(1); 

%map x_int to interval 1..lenght(i)
xmapped = offset + slope * x_int;
ind = floor(xmapped);
frac = xmapped - ind;
%interpolate by taking weighted sum of neighbouring points
y_int = y(ind) .* (1 - frac) + y(ind+1) .* frac;

% make plot to check correctness
plot(x, y, 'o-', x_int, y_int, '.')
于 2013-09-05T19:23:52.593 回答
2

这很慢，因为您遇到了 JIT最烦人的一个限制。这是在 SO 上的 MATLAB 标记中出现许多许多问题的原因：

MATLAB 的 JIT 加速器无法加速调用非内置函数的循环。

两者ppval和interp1都不是内置的（用或type ppval检查edit interp1）。它们的实现并不是特别慢，只是放在循环中时并不快。

现在我的印象是它在最近版本的 MATLAB 中变得更好，但是“内联”和“非内联”循环之间仍然存在很大差异。为什么他们的 JIT 不能通过简单地递归到非内置函数来自动执行此任务，我真的不知道。

无论如何，要解决这个问题，您应该将发生的事情的本质复制粘贴ppval到循环体中：

% Example data
x = 0:.1:10;
y = sin(x);
pp = spline(x,y);


% Your original version
tic
for n = 1:10000
    x_int = 10*rand(1000,1);
    y_int = ppval(pp, x_int);
end
toc


% "inlined" version

tic

br = pp.breaks.';
cf = pp.coefs;

for n = 1:10000

    x_int = 10*rand(1000,1);

    [~, inds] = histc(x_int, [-inf; br(2:end-1); +inf]); 

    x_shf = x_int - br(inds);    
    zero  = ones(size(x_shf));
    one   = x_shf;
    two   = one .* x_shf;
    three = two .* x_shf;

    y_int = sum( [three two one zero] .* cf(inds,:), 2);
end
toc

探查器：

在此处输入图像描述

我蹩脚的机器上的结果：

Elapsed time is 2.764317 seconds.  % ppval
Elapsed time is 1.695324 seconds.  % "inlined" version

差异实际上比我预期的要小，但我认为这主要是由于sum()- 对于这种ppval情况，我通常只需要每次迭代评估一个站点，你可以不用histc（但使用简单的矢量化代码）和矩阵/vector 乘法x*y（BLAS）而不是sum(x.*y)（快速，但不是 BLAS-fast）。

哦，好吧，减少约 60% 还不错:)

于 2013-09-06T09:27:51.183 回答