MATLAB中没有循环的矩阵计算

Question

我对执行某些数组操作的代码有疑问。它太慢了，因为我使用循环并且输入数据很大。这对我来说是最简单的方法，但现在我正在寻找比 for 循环更快的方法。我试图优化或重写代码，但没有成功。非常感谢您的帮助。

在我的代码中，我有三个数组x1、y1（网格中点的坐标）、g1（点中的值），例如它们的大小为 300 x 300。我将每个矩阵视为组合9，我计算中间的点。例如，我从g1(101,101) 开始，但我使用的是来自g1(1:201,1:201)=g2 的数据。我需要计算从g1(1:201,1:201) 的每个点到g1(101,101)（ll 矩阵）的距离，然后我计算代码中的 nn，接下来我从 nn 中找到 g1(101,101) 的值并将其放入 @987654330 @ 大批。然后我转到g1(101,102)，依此类推，直到g1(200,200)，在最后一种情况下g2=g1(99:300,99:300)。

正如我所说，这段代码效率不高，即使我必须使用比示例中给出的更大的数组，也需要花费太多时间。我希望我能足够清楚地解释我对代码的期望。我正在考虑使用arrayfun，但我从未使用过这个函数，所以我不知道应该如何使用它，但在我看来它无法处理。也许还有其他解决方案，但是我找不到任何合适的解决方案。

tic
x1=randn(300,300);
y1=randn(300,300);
g1=randn(300,300);
m=size(g1,1);
n=size(g1,2);
w=1/3*m;
k=1/3*n;
N=zeros(w,k);
for i=w+1:2*w 
    for j=k+1:2*k 
        x=x1(i,j);
        y=y1(i,j);
        x2=y1(i-k:i+k,j-w:j+w);
        y2=y1(i-k:i+k,j-w:j+w);
        g2=g1(i-k:i+k,j-w:j+w);
        ll=1./sqrt((x2-x).^2+(y2-y).^2);
        ll(isinf(ll))=0;
        nn=ifft2(fft2(g2).*fft2(ll));
        N(i-w,j-k)=nn(w+1,k+1);
      end
  end
  czas=toc;

Answer 1

不管怎样，arrayfun() 只是一个 for 循环的包装器，所以它不会带来任何性能改进。另外，您可能在x2 的定义中有错字，我假设它取决于x1。否则，这将是一个多余的变量。此外，您的i<->w/k、j<->k/w 配对似乎不一致，您也应该检查一下。 另外 同样，仅使用 tic/toc 进行计时很少准确。分析代码时，将其放入函数中并多次运行计时，并从计时中排除变量生成。更好的是：使用内置分析器。

免责声明：由于其巨大的内存需求，此解决方案可能无法解决您的实际问题。对于 300x300 矩阵的输入，这适用于大小为 300x300x100x100 的数组，这通常是不可行的。尽管如此，它还是在这里以较小的输入大小作为参考。我想添加一个基于 nlfilter() 的解决方案，但您的问题似乎太复杂而无法使用。

与向量化一样，如果您可以为它腾出内存，您可以更快地完成它。您正在尝试为每个 [i,j] 索引使用大小为 [2*k+1,2*w+1] 的矩阵。这需要形状为[2*k+1,2*w+1,w,k] 的4d 数组。对于每个元素[i,j]，您都有一个带有索引[:,:,i,j] 的矩阵与x1 和y1 的相应元素一起处理。 fft2 接受多维数组也有帮助。

这就是我的意思：

tic
x1 = randn(30,30);  %// smaller input for tractability
y1 = randn(30,30);
g1 = randn(30,30);
m = size(g1,1);
n = size(g1,2);
w = 1/3*m;
k = 1/3*n;

%// these will be indexed on the fly:    
%//x = x1(w+1:2*w,k+1:2*k);     %// size [w,k]
%//y = x1(w+1:2*w,k+1:2*k);     %// size [w,k]

x2 = zeros(2*k+1,2*w+1,w,k); %// size [2*k+1,2*w+1,w,k]
y2 = zeros(2*k+1,2*w+1,w,k); %// size [2*k+1,2*w+1,w,k]
g2 = zeros(2*k+1,2*w+1,w,k); %// size [2*k+1,2*w+1,w,k]

%// manual definition for now, maybe could be done smarter:
for ii=w+1:2*w       %// don't use i and j as variables
    for jj=k+1:2*k   %// don't use i and j as variables
        x2(:,:,ii-w,jj-k) = x1(ii-k:ii+k,jj-w:jj+w);  %// check w vs k here
        y2(:,:,ii-w,jj-k) = y1(ii-k:ii+k,jj-w:jj+w);  %// check w vs k here
        g2(:,:,ii-w,jj-k) = g1(ii-k:ii+k,jj-w:jj+w);  %// check w vs k here
    end
end

%// use bsxfun to operate on [2*k+1,2*w+1,w,k] vs [w,k]-sized arrays
%// need to introduce leading singletons with permute() in the latter
%// in order to have shape [1,1,w,k] compatible with the first array
ll = 1./sqrt(bsxfun(@minus,x2,permute(x1(w+1:2*w,k+1:2*k),[3,4,1,2])).^2 ...
           + bsxfun(@minus,y2,permute(y1(w+1:2*w,k+1:2*k),[3,4,1,2])).^2);
ll(isinf(ll)) = 0;

%// compute fft2, operating on [2*k+1,2*w+1,w,k]
%// will return fft2 for each index in the [w,k] subspace
nn = ifft2(fft2(g2).*fft2(ll));

%// we need nn(w+1,k+1,:,:) which is exactly of size [w,k] as needed
N = reshape(nn(w+1,k+1,:,:),[w,k]);  %// quicker than squeeze()
N = real(N);  %// this solution leaves an imaginary part of around 1e-12

czas=toc;