2017-01-27 65 views
1

矢量化我想“矢量化”这个循环在Matlab的计算效率矩阵二次方程式在MATLAB

for t=1:T 
    j=1; 
    for m=1:M 
     for n=1:N 
     y(t,j) = v{m,n} + data(t,:)*b{m,n} + data(t,:)*f{m,n}*data(t,:)'; 
     j=j+1; 
     end 
    end 
end 

哪里v是标量的(M X N)细胞。 b是(Kx1)个载体的(M xN)单元。 f是(K x K)矩阵的(M x N)单元格。 data是(T x K)数组。

为了让我的意思我曾经向量化同一回路中无需二次项的代码是一个示例:

B = [reshape(cell2mat(v)',1,N*M);cell2mat(reshape(b'),1,M*N)]; 
X = [ones(T,1),data]; 
y = X*B; 

谢谢!

回答

2

对于那些有兴趣在这里被我发现

f = f'; 
tMat = blkdiag(f{:})+(blkdiag(f{:}))'; 
y2BB = [reshape(cell2mat(v)',1,N*M);... 
     cell2mat(reshape(b',1,M*N));... 
     reshape(diag(blkdiag(f{:})),K,N*M);... 
     reshape(tMat((tril(tMat,-1)~=0)),sum(1:K-1),M*N)]; 
y2YBar = [ones(T,1),data,data.^2]; 

jj=1; 
kk=1; 
ll=1; 
for k=1:sum(1:K-1) 
    y2YBar = [y2YBar,data(:,jj).*data(:,kk+jj)]; 
    if kk<(K-ll) 
     kk=kk+1; 
    else 
     kk=1; 
     jj=jj+1; 
     ll=ll+1; 
    end 
end 
y = y2YBar*y2BB; 
+0

最后一个循环是草率的,但我的大脑伤害。如果您想到更有效的方式来添加这些元素,请告诉我。 – hipHopMetropolisHastings

0

解决方案下面是针对性能的最量化的形式 -

% Extract as multi-dim arrays 
vA = reshape([v{:}],M,N); 
bA = reshape([b{:}],K,M,N); 
fA = reshape([f{:}],K,K,M,N); 

% Perform : data(t,:)*f{m,n} for all iterations 
data_f_mult = reshape(data*reshape(fA,K,[]),T,K,M,N); 

% Now there are three parts : 
% v{m,n} 
% data(t,:)*b{m,n} 
% data(t,:)*f{m,n}*data(t,:)'; 

% Compute those parts one by one 
parte1 = vA(:).'; 
parte2 = data*reshape(bA,[],M*N); 

parte3 = zeros(T,M*N); 
for t = 1:T 
    parte3(t,:) = data(t,:)*reshape(data_f_mult(t,:,:),K,[]); 
end 

% Finally sum those up and to present in desired format permute dims 
sums = bsxfun(@plus, parte1, parte2 + parte3); 
out = reshape(permute(reshape(sums,T,M,N),[1,3,2]),[],M*N); 
+0

不错的工作!我会测试哪一种效率最高,并选择答案。考虑到我的应用程序T = 152,N = K = 3和M = 8,我的循环次数要少得多。 – hipHopMetropolisHastings

+0

@hipHopMetropolisHastings有没有更新? – Divakar

+0

是的,你的代码运行在0.025和我的0.008左右。再次感谢您的回答。 – hipHopMetropolisHastings