multithreading - 分段多项式求值的并行化

Question

我正在尝试评估大型分段多项式中的点，该多项式是从三次样条曲线获得的。这需要很长时间才能完成，我想加快速度。

这样，我想用并行过程而不是顺序地评估分段多项式上的点。

代码:

z = zeros(1e6, 1) ;    % preallocate some memory for speed
Y = rand(11220,161) ;  %some data, rand for generating a working example
X = 0 : 0.0125 : 2 ;   % vector of data sites
pp = spline(X, Y) ;    % get the piecewise polynomial form of the cubic spline. 

所得的结构很大。

for t = 1 : 1e6  % big number
    hcurrent = ppval(pp,t); %evaluate the piecewise polynomial at t
    z(t) = sum(x(t:t+M-1).*hcurrent,1) ; % do some operation of the interpolated value. Most likely not relevant to this question.
end

不幸的是，使用矩阵形式并使用:

hcurrent = flipud(ppval(pp, 1: 1e6 ))

需要太多的内存来处理，因此无法完成。有没有一种方法可以批量处理此代码以加快处理速度？

Answer 1

对于标量第二个参数，如您的示例所示，您要处理两个问题。首先，有大量的函数调用开销和冗余计算(例如，每次循环迭代都会调用unmkpp(pp))。其次， ppval 被编写为通用的，因此它没有完全向量化，并且执行了许多您不需要的事情。

下面是矢量化的代码，它们利用了问题的某些结构(例如，t是大于0的整数)，避免了函数调用的开销，将一些计算移出了主for循环之外(付出了一些代价)。 )，并摆脱了for内部的ppval循环:

n = 1e6;
z = zeros(n,1);
X = 0:0.0125:2;
Y = rand(11220,numel(X));
pp = spline(X,Y);

[b,c,l,k,dd] = unmkpp(pp);
T = 1:n;
idx = discretize(T,[-Inf b(2:l) Inf]); % Or: [~,idx] = histc(T,[-Inf b(2:l) Inf]);
x = bsxfun(@power,T-b(idx),(k-1:-1:0).').';

idx = dd*idx;
d = 1-dd:0;
for t = T
    hcurrent = sum(bsxfun(@times,c(idx(t)+d,:),x(t,:)),2);
    z(t) = ...;
end

生成的代码花费 n=1e6示例的时间的〜34％。注意，由于矢量化，计算以不同的顺序执行。由于浮点数学的性质，这将导致 ppval和我的优化版本的输出之间略有差异。任何差异都应为 eps(hcurrent)的几倍。您仍然可以尝试使用 parfor进一步加快计算速度(在四个已经运行的工作程序中，我的系统仅花费了代码原始时间的12％)。

我认为以上是概念证明。如果您的示例与您的实际代码和数据不太吻合，则我可能过度优化了上面的代码。在这种情况下，建议您创建自己的优化版本。您可以通过在命令窗口中键入 ppval来查看 edit ppval的代码来开始。您可以通过查看问题的结构以及 z向量中最终想要的内容来实现进一步的优化。

在内部， ppval仍然使用 histc ，该名称已被弃用。上面的代码使用 discretize 执行相同的任务 as suggested by the documentation。