我编写了这段代码来执行二维矩阵值函数的一维卷积（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
问问题
3460 次
3 回答
4
新解决方案：
这是一个基于旧解决方案的新解决方案，它解决了先前给出的公式。问题中的代码实际上是对该公式的修改，其中第三维中两个矩阵之间的重叠被重复移动（类似于沿数据第三维的卷积）。我之前给出的解决方案只计算了问题中代码的最后一次迭代的结果（即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
即使使用 just kend = 100，这个解决方案对我来说比问题中的解决方案快大约 30 倍，比纯基于 for 循环的解决方案（这将涉及5 个循环！）快大约 4 倍。请注意，下面关于浮点精度的讨论仍然适用，因此您会看到解决方案之间在相对浮点精度上的细微差别是正常的，并且可以预期。
旧解决方案：
根据您在评论中链接到的公式：
看来您实际上想要做的事情与您在问题中提供的代码不同。假设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适当地重塑和重组来执行此操作而无需任何 for 循环B：
%# 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，我会得到这些结果（在几秒钟内执行给定解决方案 MATLAB R2010b 的 10,000 次迭代所需的时间）：
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 左右时，向量化矩阵乘法解决方案开始胜出，与kfor 循环解决方案相比，随着增加的扩展更好。由于我不确定的原因，混合的 for-loop/ matrix -multiply 解决方案的规模非常大。
所以，如果你希望k很大，我会选择矢量化矩阵乘法解决方案。要记住的一件事是，您从每个解决方案（矩阵C）获得的结果将略有不同（在浮点精度级别上），因为每个解决方案执行的加法和乘法的顺序不同，因此导致舍入误差累积的差异。简而言之，这些解决方案的结果之间的差异应该可以忽略不计，但您应该意识到这一点。
于 2011-04-25T15:15:09.443 回答
2
你有没有研究过 Matlab 的conv方法？
我无法将它与您提供的代码进行比较，因为您提供的内容让我在尝试访问 A 的第零个元素时遇到了问题。（何时k=1，k-1=0。）
于 2011-04-23T16:38:13.117 回答
1
您是否考虑过使用 FFT 进行卷积？卷积运算只是频域中的逐点乘法。你必须对有限序列采取一些预防措施，因为如果你不小心，你最终会得到循环卷积（但这很容易处理）。
这是一维案例的简单示例。
>> 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（当我对其进行点 FFT 时，这会自动处理）。b4+3-1=6
编辑
尽管您显示的方程类似于循环卷积，但并非完全如此。所以你可以放弃 FFT 方法和内置conv*函数。要回答您的问题，以下是在没有显式循环的情况下完成的相同操作：
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 的循环已经走过了漫长的道路。当然，矢量化操作正在飞速发展。但并非所有内容都可以向量化，有时，循环比这种复杂的匿名函数更有效。通过优化我的构造（或者可能采取不同的方法），我可能会在这里和那里再削减十分之一，但我不会这样做。
请记住，好的代码应该是可读的，以及以可读性为代价的高效和微小的优化对任何人都没有好处。虽然我写了上面的代码，但如果我在一个月后重新访问它，我肯定无法破译它的作用。您的循环代码清晰、易读且快速，我建议您坚持使用它。
于 2011-04-23T19:42:48.940 回答