use*_*865 6 matlab vectorization multidimensional-array
假设我有两个矩阵A和B.
A = rand(4,5,3);
B = rand(4,5,6)
Run Code Online (Sandbox Code Playgroud)
我想应用函数'corr2'来计算相关系数.
corr2(A(:,:,1),B(:,:,1))
corr2(A(:,:,1),B(:,:,2))
corr2(A(:,:,1),B(:,:,3))
...
corr2(A(:,:,1),B(:,:,6))
...
corr2(A(:,:,2),B(:,:,1))
corr2(A(:,:,2),B(:,:,2))
...
corr2(A(:,:,3),B(:,:,6))
Run Code Online (Sandbox Code Playgroud)
如何避免使用循环来创建这样的矢量化?
入侵m文件corr2以创建用于处理3D阵列的自定义矢量化版本.这里提出两种方法bsxfun(当然!)
方法#1
szA = size(A);
szB = size(B);
a1 = bsxfun(@minus,A,mean(mean(A)));
b1 = bsxfun(@minus,B,mean(mean(B)));
sa1 = sum(sum(a1.*a1));
sb1 = sum(sum(b1.*b1));
v1 = reshape(b1,[],szB(3)).'*reshape(a1,[],szA(3));
v2 = sqrt(sb1(:)*sa1(:).');
corr3_out = v1./v2; %// desired output
Run Code Online (Sandbox Code Playgroud)
corr3_out商店corr2的所有3D片之间的结果A和B.
因此,对于A = rand(4,5,3), B = rand(4,5,6),我们将corr3_out作为一个6x3数组.
方法#2
保存少量呼叫sum和mean使用reshape相反的方法略有不同-
szA = size(A);
szB = size(B);
dim12 = szA(1)*szA(2);
a1 = bsxfun(@minus,A,mean(reshape(A,dim12,1,[])));
b1 = bsxfun(@minus,B,mean(reshape(B,dim12,1,[])));
v1 = reshape(b1,[],szB(3)).'*reshape(a1,[],szA(3));
v2 = sqrt(sum(reshape(b1.*b1,dim12,[])).'*sum(reshape(a1.*a1,dim12,[])));
corr3_out = v1./v2; %// desired output
Run Code Online (Sandbox Code Playgroud)
基准代码 -
%// Create random input arrays
N = 55; %// datasize scaling factor
A = rand(4*N,5*N,3*N);
B = rand(4*N,5*N,6*N);
%// Warm up tic/toc
for k = 1:50000
tic(); elapsed = toc();
end
%// Run vectorized and loopy approach codes on the input arrays
%// 1. Vectorized approach
%//... solution code (Approach #2) posted earlier
%// clear variables used
%// 2. Loopy approach
tic
s_A=size(A,3);
s_B=size(B,3);
out1 = zeros(s_B,s_A);
for ii=1:s_A
for jj=1:s_B
out1(jj,ii)=corr2(A(:,:,ii),B(:,:,jj));
end
end
toc
Run Code Online (Sandbox Code Playgroud)
结果 -
-------------------------- With BSXFUN vectorized solution
Elapsed time is 1.231230 seconds.
-------------------------- With loopy approach
Elapsed time is 139.934719 seconds.
Run Code Online (Sandbox Code Playgroud)
MATLAB-JIT爱好者在这里展示了一些爱!:)