如何從這些循環中只獲得矩陣的對角線？

Question

我對8維經驗copula有以下代碼，它創建了一個8d矩陣，但我只需要在此代碼中名爲EC的矩陣的對角線。由於此代碼非常慢，是否有無需計算所有「ecop」就可以獲得「EC」？

function EC = ecopula8d(x) 

[m n] = size(x); 

y = sort(x); 

for r=1:m 
     for q=1:m 
      for p=1:m 
       for o=1:m 
        for l=1:m 
         for k=1:m 
          for j=1:m 
           for i=1:m 
      ecop(i,j,k,l,o,p,q,r) = sum((x(:,1)<=y(i,1)).*(x(:,2)<=y(j,2)).*(x(:,3)<=y(k,3)).*(x(:,4)<=y(l,4))... 
      .*(x(:,5)<=y(o,5)).*(x(:,6)<=y(p,6)).*(x(:,7)<=y(q,7)).*(x(:,8)<=y(r,8)))/(m+1); 
           end 
          end 
         end 
        end 
       end 
      end 
     end 
end 

for i=1:m 
EC(i,1)=ecop(i,i,i,i,i,i,i,i); 
end 

Answer 1

您可以採用完全與bsxfun量化的解決方案 -

EC = squeeze(sum(all(bsxfun(@le,x,permute(y,[3 2 1])),2),1))/(m+1) 

神奇與使用permute使我們能夠繼續vectorization全油門在這裏發生。

這裏的一個基準測試來比較這種方法和運行時效率other partially vectorized approach with bsxfun - 由此獲得

x = rand(2000,2000); 
y = sort(x); 
m = size(x,1); 

%// Warm up tic/toc. 
for k = 1:100000 
    tic(); elapsed = toc(); 
end 

disp('----------- With completely vectorized solution') 
tic 
EC1 = squeeze(sum(all(bsxfun(@le,x,permute(y,[3 2 1])),2),1))/(m+1); 
toc, clear EC1 

disp('----------- With partial vectorized solution') 
tic 
for i = 1:m 
    EC2(i,1) = sum(all(bsxfun(@le, x, y(i,:)), 2), 1)/(m+1); 
end 
toc 

的運行時間分別爲 -

----------- With completely vectorized solution 
Elapsed time is 2.883594 seconds. 
----------- With partial vectorized solution 
Elapsed time is 4.752508 seconds. 

人們可以預先分配的其他部分向量化方法 -

EC2 = zeros(m,1); 
for i = 1:m 
    EC2(i,1) = sum(all(bsxfun(@le, x, y(i,:)), 2), 1)/(m+1); 
end 

運行時間因此得到的沒有那麼不同 -

----------- With completely vectorized solution 
Elapsed time is 2.963835 seconds. 
----------- With partial vectorized solution 
Elapsed time is 4.620455 seconds. 

Answer 2

我會使用方法。一旦是ND數組轉換成方形2- d（如果可能的話），然後簡單地提取對角項，因爲它們應該在兩種情況下等於：

EC=diag(reshape(ecop,size1,size2)); 

我建議嘗試使用Python，因爲numpy有非常好的高效線性代數包來處理ND數組。 Matlab在添加和更新庫時相當緩慢。

Answer 3

如果您最初的計算是正確的我沒有檢查（如相比，與維基百科文章的公式您的實現），但你的代碼應該是等同於：

[m n] = size(x); 
y = sort(x); 
for i = 1:m 
    EC(i,1) = sum(all(bsxfun(@le, x, y(i,:)), 2), 1)/(m+1); 
end