4

我在 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 循环，以便使用矩阵运算和函数实现它。可以矢量化吗？

4

4 回答 4

2

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

取自 Matlab 帮助：

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

于 2012-09-20T07:50:15.107 回答
1

如您所说，如果您对所有测量有不同的测量时间（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;
于 2012-09-20T07:39:07.367 回答
1

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)
于 2012-09-20T04:45:37.820 回答
0

我在工作，所以没有时间熟悉interp1（我以前从未使用过），所以如果以下答案不恰当，我提前道歉。

话虽如此，由于循环的迭代不相互依赖，向量化应该是可能的。mat2cell在我看来，您可以使用and摆脱显式循环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应用于 中的每个单元格，并返回 中的结果。您可以简单地替换为，当然，适当地调整输入。sumVCellA@sum@interp1cellfun

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

于 2012-09-20T03:02:56.043 回答