分段多项式求值的并行化

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 循环之外（以额外的内存为代价），并摆脱了ppval 内部的for 循环：

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%）。

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

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

Answer 2

使用parfor 命令进行并行循环。参见here，还将z向量预计算为z(j) = x(j:j+M-1)和parfor中的hcurrent以加快速度。

Answer 3

样条参数估计可以写成矩阵形式。

一旦您以矩阵形式编写并求解它，您就可以使用模型矩阵来使用矩阵乘法计算所有数据点上的样条曲线，这可能是 MATLAB 中最优化的操作。