Matlab：3维矩阵中的矢量化技术答案

【问题标题】：Matlab : Vectorize technics in 3 dimensions matrixMatlab：3维矩阵中的矢量化技术
【发布时间】：2017-03-06 17:55:59
【问题描述】：

我实际上对我的一个代码进行了矢量化处理，但我遇到了一些问题。

这是我的初始代码：

CoordVorBd = random(N+1,3)
CoordCP = random(N,3)
v = random(1,3)
for i = 1 : N
    for j = 1 : N
            ri1j = (-CoordVorBd (i,:) + CoordCP(j,:));

            vij(i,j,:) = cross(v,ri1j))/(norm(ri1j)
     end     

end

我已经开始向量化创建一些包含 3*1 向量的矩阵。我的矩阵大小是 N*N*3。

CoordVorBd1(1:N,:) = CoordVorBd(2:N+1,:);
CoordCP_x= CoordCP(:,1);
CoordCP_y= CoordCP(:,2);
CoordCP_z= CoordCP(:,3);
CoordVorBd_x = CoordVorBd([1:N],1);
CoordVorBd_y = CoordVorBd([1:N],2);
CoordVorBd_z = CoordVorBd([1:N],3);
CoordVorBd1_x = CoordVorBd1(:,1);
CoordVorBd1_y = CoordVorBd1(:,2);
CoordVorBd1_z = CoordVorBd1(:,3);
[X,Y] = meshgrid (1:N);

ri1j_x = (-CoordVorBd_x(X) + CoordCP_x(Y));
ri1j_y = (-CoordVorBd_y(X) + CoordCP_y(Y));
ri1j_z = (-CoordVorBd_z(X) + CoordCP_z(Y));

ri1jmat(:,:,1) = ri1j_x(:,:);
ri1jmat(:,:,2) = ri1j_y(:,:);
ri1jmat(:,:,3) = ri1j_z(:,:);
vmat(:,:,1) = ones(N)*v(1);
vmat(:,:,2) = ones(N)*v(2);
vmat(:,:,3) = ones(N)*v(3);

此代码有效，但在变量创建方面很繁重。我没有实现一次将矢量化应用于所有矩阵。

这样的公式

ri1jmat(X,Y,1:3) = (-CoordVorBd (X,:) + CoordCP(Y,:));

没有用... 如果有人有一些想法让一些更清洁的东西。

此时我有一个包含所有向量的 N*N*3 矩阵 ri1jmat。

我想计算作为向量范数的 N*N rij1norm 矩阵

rij1norm(i,j) = norm(ri1jmat(i,j,1:3))

能够向量化 vij 矩阵。

vij(:,:,1:3) = (cross(vmat(:,:,1:3),ri1jmat(:,:,1:3))/(ri1jmatnorm(:,:));

叉积有效。

我尝试了多种方法，但没有实现这个 rij1norm 矩阵，而无需执行双循环。

如果有人有一些技巧，在此先感谢。

【问题讨论】：

标签： matlab matrix vector vectorization

【解决方案1】：

这是一个矢量化版本。请注意，您的原始循环不包括CoordVorBd 的最后一列，因此如果这是故意的，您还需要将其从下面的代码中删除。我认为这是一个错误。

CoordVorBd = rand(N+1,3);
CoordCP = rand(N,3);
v = rand(1,3);

repCoordVor=kron(CoordVorBd', ones(1,size(CoordCP,1)))'; %based on http://stackoverflow.com/questions/16266804/matlab-repeat-every-column-sequentially-n-times
repCoordCP=repmat(CoordCP, size(CoordVorBd,1),1); %repeat matrix
V2=-repCoordVor + repCoordCP; %your ri1j 
nrm123=sqrt(sum(V2.^2,2)); %vectorized norm for each row
vij_unformatted=cat(3,(v(:,2).*V2(:,3) - V2(:,2).*v(:,3))./nrm123,(v(:,3).*V2(:,1) - V2(:,3).*v(:,1))./nrm123,(v(:,1).*V2(:,2) - V2(:,1).*v(:,2))./nrm123); % cross product, expanded, and each term divided by norm, could use bsxfun(@rdivide,cr123,nrm123) instead, if cr123 is same without divisions
vij=permute(reshape( vij_unformatted,N,N+1,3),[2,1,3]); %reformat to match your vij

【讨论】：

【解决方案2】：

这是使用 arrayfun 的另一种方法

% Define a meshgrid of indices to run over
[I, J] = meshgrid(1:N, 1:(N+1));
% Calculate ril for each index
rilj = arrayfun(@(x, y) -CoordVorBd (y,:) + CoordCP(x,:), I, J, 'UniformOutput', false);
%Calculate vij for each point
temp_vij1 = arrayfun(@(x, y) cross(v, rilj{x, y}) / norm(rilj{x, y}), J, I, 'UniformOutput', false);
%Reshape the matrix into desired format
temp_vij2 = cell2mat(temp_vij1);
vij = cat(3, temp_vij2(:, 1:3:end), temp_vij2(:, 2:3:end), temp_vij2(:, 3:3:end));

【讨论】：