矢量化插值matlab

Question

我在S个时间点进行了一组N次测量(不同测量的时间点不同)。我有两个矩阵：

V-表示测量值的NxS矩阵

T-表示测量时间的NxS矩阵

我想生成一个矩阵VI，它表示时间TI处的线性插值测量。代码的非矢量化版本如下：

tic;
VI = zeros([size(V,1), size(TI,2)]);
for j = 1:size(V,1)
VI(j,:) = interp1(T(j,:),V(j,:),TI);
end
toc;

我想重写这段代码以消除for循环，从而使用矩阵运算和函数来实现它。它可以矢量化吗？

Answer 1

如果没有数据并运行探查器，很难说什么，但如果对数据进行排序，则可以使用interp1q而不是interp，后者不会对数据进行任何检查。

取自Matlab帮助：

interp1q要正常工作，x必须是单调递增的列向量。Y必须是长度为(x)的列向量或矩阵行。xi必须是列向量

Answer 2

Matlab将数据打包在列中，而不是行中，因此我怀疑您只需在循环中从遍历行更改为遍历列，就会看到速度的提高：

[N, S] = size(V);
VT = V';                             % Value series in columns
TT = T';                             % Time series in columns
VIT = zeros(length(TI), N);          % Interpolated value series in columns
for j = 1:N
VIT(:, j) = interp1(VT(:, j), TT(:, j), TI);
end
VI = VIT';                           % Interpolated value series in rows

不幸的是，我不认为可以做太多来进一步提高此代码的性能，因为不同的值序列没有相关的时间序列。如果它们有相同的时间，这样我们就可以将T折叠成具有length(T) = S的向量，那么这样做就很容易了：

VIT = interp1(VT, T, TI);            % (see VIT and VT from above)

Answer 3

如果您对所有测量都有不同的测量时间(T是矩阵而不是向量)，您可以通过调用arrayfun来完成您想要的操作，如下所示：

VI = arrayfun(@(x)(interp1(T(x,:),V(x,:),TI)), 1:size(V, 1), 'UniformOutput', false);
VI = cell2mat(VI');

arrayfun类似于一个循环，但由于它是一个内部matlab函数，它可能会更快。它返回一个向量单元格，所以第二行确保有一个矩阵作为输出。你可能不需要它——这取决于你以后对数据做什么。

另一方面，如果在同一时间对不同的N值进行测量(T是大小为S的向量，而不是矩阵，或者换句话说，T的所有行都相等)，则可以在对interp1的一次调用中进行插值。

VI = interp1(T(1,:), V', TI)

这里你必须转置V，因为interp1在列内插值。这是因为MATLAB按列存储矩阵(列在内存中是连续的)。如果将V作为SxN矩阵传递，则可能会使interp1的并行化更加高效，因为所有CPU都可以以更高效的方式访问内存。因此，我建议您在整个代码中转置矩阵，当然，除非出于性能原因，您在其他地方依赖这种精确的数据布局。

编辑由于矩阵的列布局，可以通过转换矩阵和按列工作来改进原始代码。对于N=1000、S=10000和TI为10000的元素，以下版本在我的计算机上的速度大约快20%。由于更高效的内存带宽利用率，它可能会随着系统大小而增长。

tic;
VI = zeros(size(TI,2), size(V,2));
for j = 1:size(V,2)
VI(:,j) = interp1(T(:,j),V(:,j),TI);
end
toc;

Answer 4

我正在上班，没有时间熟悉interp1(我以前从未使用过它)，所以如果下面的答案不合适，我提前道歉。

话虽如此，由于循环的迭代并不相互依赖，因此矢量化应该是可能的。在我看来，使用mat2cell和cellfun可以去掉显式循环。

我的意思的一个简单例子如下：

NumRow = 4;
NumCol = 3;
V = randn(NumRow, NumCol);
VCell = mat2cell(V, ones(NumRow, 1), NumCol);
A = cellfun(@sum, VCell);

当然，我所做的相当于sum(V, 2)。但我认为我使用的方法适用于你的情况。mat2cell函数将V转换为单元格列向量，其中每个单元格包含一行V。然后对cellfun的调用将sum函数应用于VCell中的每个单元格，并在A中返回结果。您可以简单地将@sum替换为@interp1，当然，还可以适当地将输入调整为cellfun。

如果你不能让它发挥作用，请告诉我，我回家后会尽量做一些更明确的事情。此外，如果你真的能让它发挥作用，但它不会加快速度，我很有兴趣知道这一点，所以请在这里发布结果。