矢量化插值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

我在工作，没时间熟悉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。

如果你不能让它工作，请告诉我，我会在回家后尝试整理一些更明确的东西。另外，如果你确实让它工作了，但它并没有加快速度，我很想知道，所以请在此处发布结果。

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 和 10000 个元素的 TI，以下版本在我的计算机上快了大约 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

如果没有数据并运行分析器，很难说任何事情，但如果您的数据已排序，您可以使用 interp1q 而不是 interp ，它不会对数据进行任何检查。

取自 Matlab 帮助：

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