如何使用 3 个生成矩阵的循环对代码进行矢量化答案

【问题标题】：How to vectorize code with 3 loops that generates a matrix如何使用 3 个生成矩阵的循环对代码进行矢量化
【发布时间】：2017-12-17 00:45:50
【问题描述】：

您能帮我向量化这个构造矩阵A、维度为MNx(2+N-1)xR 的Matlav 代码以加快速度吗？目前大约需要。 8 秒

正在初始化

R=200;
M=400;
N=20;
B=[kron((1:1:M)', ones(N,1)) repmat((1:1:N)', M,1)]; %(MN)x(2)

B 看起来像

B=[1 1;
   1 2;
   ...;
   1 20;
   2 1;
   2 2;
   ...
   2 20;
   ...
   400 20]

代码

A=[ repmat(B,1,1,R)  zeros(M*N,N-1,R)]; %Allocate the space
                                        %(MN)x(2+N-1)x(R)
                                        % I want A(:,1:2,r)=B for r=1,...,R


%Fill the remaining columns of A in the following way
for r=1:R
    utemp1=randn(M*N, N-1); %Generate a matrix of random 
                            %numbers of dimension (M*N)x(N-1)
    utemp2=randn(M, N);  %Generate a matrix of random 
                         %numbers of dimension (M)x(N)
    utemp3=zeros(M*N,N-1); %Generate a matrix of random 
                           %numbers of dimension(M)x(N-1)
    for m=1:M
        for j=1:N
            utemp3((m-1)*N+j,:)= utemp2(m,j)+[utemp2(m,1:j-1) utemp2(m,j+1:N)]; %(1)x(N-1)
              %e.g. if m=2, j=3, I want to fill the 23th row of utemp3
              %with utemp2(2,3)+[utemp2(2,1:2) utemp2(m,4:20)];

              %e.g. if m=4, j=1, I want to fill the 61st row of utemp3
              %with utemp2(4,1)+[utemp2(4,2:20)];
        end
    end
    A(:,3:end,r)=utemp1+utemp3;     %sum utemp1
end

一些解释

为r=1,...,R

A 是这样的

对于m=1,...,M 和对于j=1,...,N

A(:,:,r) 中前两列中以[m j] 开头的行被填充到剩余的(N-1) 列中，每个h~=j 用uj+uh+ujh 填充，其中uj, uh, ujh 是独立同分布的标准高斯数，可以是在utemp1 和utemp2 中找到。

【问题讨论】：

还请描述您要解决的问题，而不仅仅是转储代码并要求我们对其进行矢量化...如果我们知道您试图做什么，也许我们可以建议一个内置或 FEX 函数...
...然后你想用 A 做什么？
这是一段很长的代码（模拟矩估计方法）的一部分。简而言之，我用它来估计使用一些数据的模型的参数。我想提高效率，因为我需要重复该步骤 500 次。
如果你想对每一行的元素求和，我认为有一些快速的内置插件可用

标签： matlab performance for-loop multidimensional-array vectorization

【解决方案1】：

您可以预先计算索引并在 500 次迭代中使用它们：

idx = repmat(reshape(1:M*N,M,N).',1,N);
idx = reshape(idx(logical(kron(~eye(N),ones(1,M)))),N-1,[]).';
for k=1:500
    for r=1:R
        utemp2=randn(M*N,1);
        A(:,3:end,r)=randn(M*N, N-1)+bsxfun(@plus,utemp2,utemp2(idx) );
    end
end

然而，分配大型矩阵，特别是重复元素，并对其进行矢量化操作并不总是最有效的方法。有内置函数可以直接对原始数组进行操作，避免数组的重复元素。

【讨论】：