MATLAB：矩阵值函数的卷积

Question

我编写了这段代码来执行二维矩阵值函数的一维卷积（k 是我的时间索引，kend 大约是 10e3）。有没有更快或更简洁的方法来做到这一点，也许使用内置函数？

for k=1:kend
  C(:,:,k)=zeros(3);
  for l=0:k-1
      C(:,:,k)=C(:,:,k)+A(:,:,k-l)*B(:,:,l+1);
  end
end

Answer 1

新解决方案：

这是基于旧解决方案的新解决方案，它解决了先前给出的公式。问题中的代码实际上是对该公式的修改，其中第三维中两个矩阵之间的重叠被重复移动（类似于沿数据第三维的卷积）。我之前给出的解决方案只计算了问题中代码的 last 迭代的结果（即k = kend）。所以，这里有一个完整的解决方案，它应该比问题中kend 大约 1000 的代码更高效：

kend = size(A,3);                         %# Get the value for kend
C = zeros(3,3,kend);                      %# Preallocate the output
Anew = reshape(flipdim(A,3),3,[]);        %# Reshape A into a 3-by-3*kend matrix
Bnew = reshape(permute(B,[1 3 2]),[],3);  %# Reshape B into a 3*kend-by-3 matrix
for k = 1:kend
  C(:,:,k) = Anew(:,3*(kend-k)+1:end)*Bnew(1:3*k,:);  %# Index Anew and Bnew so
end                                                   %#   they overlap in steps
                                                      %#   of three

即使只使用kend = 100，这个解决方案对我来说比问题中的解决方案快大约 30 倍，比纯基于 for 循环的解决方案快大约 4 倍（这将涉及 5 个循环！）。请注意，下面关于浮点精度的讨论仍然适用，因此这是正常的，并且预计您会看到在 relative floating-point accuracy 顺序上的解决方案之间存在细微差异。

---

旧解决方案：

根据您在评论中链接到的公式：

看来您实际上想要做的事情与您在问题中提供的代码不同。假设 A 和 B 是 3×3×k 矩阵，则结果 C 应该是一个 3×3 矩阵，并且链接中的公式写成一组嵌套的 for 循环将看起来像这样：

%# Solution #1: for loops
k = size(A,3);
C = zeros(3);
for i = 1:3
  for j = 1:3
    for r = 1:3
      for l = 0:k-1
        C(i,j) = C(i,j) + A(i,r,k-l)*B(r,j,l+1);
      end
    end
  end
end

现在，可以通过适当地重塑和重组 A 和 B 来执行此操作而无需任何 for 循环：

%# Solution #2: matrix multiply
Anew = reshape(flipdim(A,3),3,[]);        %# Create a 3-by-3*k matrix
Bnew = reshape(permute(B,[1 3 2]),[],3);  %# Create a 3*k-by-3 matrix
C = Anew*Bnew;                            %# Perform a single matrix multiply

您甚至可以修改问题中的代码，以创建一个具有单个循环的解决方案，该循环执行您的 3×3 子矩阵的矩阵乘法：

%# Solution #3: mixed (loop and matrix multiplication)
k = size(A,3);
C = zeros(3);
for l = 0:k-1
  C = C + A(:,:,k-l)*B(:,:,l+1);
end

那么现在的问题是：这些方法中哪一种更快/更清洁？

嗯，“更干净”是非常主观的，老实说，我无法告诉您以上哪些代码更容易理解操作正在做什么。第一个解决方案中的所有循环和变量使得跟踪正在发生的事情有点困难，但它清楚地反映了公式。第二种解决方案将其全部分解为一个简单的矩阵运算，但很难看出它与原始公式有何关系。第三种解决方案似乎介于两者之间。

所以，让我们把速度作为决胜局。如果我为k 的多个值对上述解决方案计时，我会得到这些结果（以秒为单位执行给定解决方案的 10,000 次迭代，MATLAB R2010b）：

 k   |  loop  | matrix multiply |  mixed
-----+--------+-----------------+--------
 5   | 0.0915 |      0.3242     | 0.1657
 10  | 0.1094 |      0.3093     | 0.2981
 20  | 0.1674 |      0.3301     | 0.5838
 50  | 0.3181 |      0.3737     | 1.3585
 100 | 0.5800 |      0.4131     | 2.7311   * The matrix multiply is now fastest
 200 | 1.2859 |      0.5538     | 5.9280

好吧，事实证明，对于较小的 k 值（大约 50 或更小），for 循环解决方案实际上胜出，再次表明 for 循环并不像过去认为的那样“邪恶”旧版本的 MATLAB。在某些情况下，它们可能比聪明的矢量化更有效。然而，当k 的值大于 100 左右时，向量化矩阵乘法解决方案开始胜出，随着k 的增加，其缩放比 for 循环解决方案要好得多。由于我不确定的原因，混合的 for-loop/matrix-multiply 解决方案严重扩展。

所以，如果您希望 k 很大，我会选择矢量化矩阵乘法解决方案。要记住的一件事是，您从每个解决方案（矩阵C）中获得的结果将略有不同（在floating-point precision 的级别上），因为每个解决方案执行的加法和乘法的顺序是不同的，从而导致rounding errors的累积差异。简而言之，这些解决方案的结果之间的差异应该可以忽略不计，但您应该意识到这一点。

Answer 2

你研究过 Matlab 的 conv 方法吗？

我无法将它与您提供的代码进行比较，因为您提供的内容让我无法尝试访问 A 的第零个元素。（当 k=1、k-1=0 时。）

Answer 3

您是否考虑过使用 FFT 进行卷积？卷积操作是simply a point-wise multiplication in the frequency domain。你必须对有限序列采取一些预防措施，因为如果你不小心，你最终会得到循环卷积（但这很容易处理）。

这是一个一维案例的简单示例。

>> a=rand(4,1);
>> b=rand(3,1);
>> c=conv(a,b)

c =

    0.1167
    0.3133
    0.4024
    0.5023
    0.6454
    0.3511

同样使用 FFT

>> A=fft(a,6);
>> B=fft(b,6);
>> C=real(ifft(A.*B))

C =

    0.1167
    0.3133
    0.4024
    0.5023
    0.6454
    0.3511

M 点向量和N 点向量的卷积得到M+N-1 点向量。因此，在进行 FFT 之前，我已经用零填充了每个向量 a 和 b（这在我对其进行 4+3-1=6 点 FFT 时会自动处理）。

编辑

虽然您展示的方程类似于循环卷积，但并非完全如此。因此，您可以放弃 FFT 方法和内置的 conv* 函数。为了回答您的问题，以下是在没有 explicit 循环的情况下完成的相同操作：

dim1=3;dim2=dim1;
dim3=10;
a=rand(dim1,dim2,dim3);
b=rand(dim1,dim2,dim3);


mIndx=cellfun(@(x)(1:x),num2cell(1:dim3),'UniformOutput',false);
fun=@(x)sum(reshape(cell2mat(cellfun(@(y,z)a(:,:,y)*b(:,:,z),num2cell(x),num2cell(fliplr(x)),'UniformOutput',false)),[dim1,dim2,max(x)]),3);
c=reshape(cell2mat(cellfun(@(x)fun(x),mIndx,'UniformOutput',false)),[dim1,dim2,dim3]);

• mIndx 这里是一个单元格，其中ith 单元格包含一个向量1:i。这是您的l 索引（正如其他人所指出的，请不要使用l 作为变量名）。
• 下一行是执行卷积运算的匿名函数，利用了k 索引只是l 索引翻转的事实。这些操作是在单个单元上进行的，然后进行组装。
• 最后一行实际上是对矩阵进行运算。

答案与循环获得的答案相同。但是，您会发现循环解决方案实际上是一个数量级快（我的代码平均为 0.007 秒，循环平均为 0.0006 秒）。这是因为循环非常简单，而在这种嵌套结构中，有大量的函数调用开销和重复的整形会减慢它的速度。

自从早期人们对循环感到恐惧以来，MATLAB 的循环已经走过了漫长的道路。当然，矢量化操作正在飞速发展。但并非所有内容都可以向量化，有时，循环比这种复杂的匿名函数更有效。通过优化我的结构（或者可能采取不同的方法），我可能会在这里和那里再削减十分之一，但我不会这样做。

请记住，好的代码应该是可读的，以及以牺牲可读性为代价的高效和微小的优化对任何人都没有好处。虽然我写了上面的代码，但如果我在一个月后重新访问它，我肯定无法破译它的作用。您的循环代码清晰、可读且速度很快，我建议您坚持下去。