2010-05-19 82 views
9

我想向量化以下MATLAB代碼。我認爲這一定很簡單,但我覺得它很混亂。向量化矩陣中不同對角線的總和

r = some constant less than m or n 
[m,n] = size(C); 
S = zeros(m-r,n-r); 
for i=1:m-r+1 
    for j=1:n-r+1 
     S(i,j) = sum(diag(C(i:i+r-1,j:j+r-1))); 
    end 
end 

的代碼計算得分的表,小號,對於動態規劃算法,從另一個得分表,Ç
對角線求和是針對所有可能的部分(大小爲r)生成用於生成的各個數據片段的分數,所述數據用於生成C

在此先感謝您的任何答案!很抱歉,如果這一塊應該是顯而易見的......

注意
內置CONV2竟然是比convnfft更快,因爲我的眼睛(R)是相當小的(5 < = R < = 20) 。 convnfft.m指出r應該是> 20以表示任何好處。

回答

10

如果我理解正確,你試圖計算C的每個子數組的對角線總和,你已經刪除C的最後一行和一列(如果你不應該刪除行/列,你需要循環到m-r + 1,並且您需要將整個數組C傳遞給我的解決方案中的函數)。

您可以通過convolution做這個手術,就像這樣:

S = conv2(C(1:end-1,1:end-1),eye(r),'valid'); 

如果C和R都很大,你可能想看看CONVNFFT從MATLAB文件交換來加快計算。

+0

太棒了,謝謝。 明天我會看看CONVNFFT。在我用於測試的數據子集(比實際數據大約小500倍)上,與循環相比,內置的conv2函數實現了69,652次調用次數減少和34.56次執行時間減少(23.5 vs 0.68秒)。 – 2010-05-19 12:33:12

2

我認爲你可能需要將C重新排列成3D矩陣,然後沿着其中一個維度求和它。我會盡快發佈答案。

編輯

我沒能找到一個方法來清潔它vectorise,但我沒有找到功能accumarray,這可能有一定的幫助。當我回家時,我會更詳細地看它。

EDIT#2

實測值利用線性索引更簡單的解決方案,但是這可能是存儲器密集型的。在C(1,1)處,我們想要求和的索引是1+ [0,m + 1,2 * m + 2,3 * m + 3,4 * m + 4,...] ,或(0:R-1)+(0:M:(R-1)* M)

sum_ind = (0:r-1)+(0:m:(r-1)*m); 

創建S_offset,一個(MR)通過(NR)中的R矩陣,使得S_offset(: ,:,1)= 0,S_offset(:,:2)= m + 1,S_offset(:,:3)= 2 * m + 2,依此類推。

S_offset = permute(repmat(sum_ind, [m-r, 1, n-r]), [1, 3, 2]); 

創建S_base,基地陣列地址的矩陣,從中計算偏移量。

S_base = reshape(1:m*n,[m n]); 
S_base = repmat(S_base(1:m-r,1:n-r), [1, 1, r]); 

最後,使用S_base+S_offset解決C.

S = sum(C(S_base+S_offset), 3); 

你可以,當然,使用bsxfun等方法,使之更加高效;在這裏我選擇了爲了清晰起見。我還沒有對此進行基準測試,以瞭解它如何與雙循環方法進行比較;我需要先回家吃晚飯!

+0

嗯,與從二維索引對角線成爲該索引的第三維? – 2010-05-19 09:36:27

+1

@JS:好主意。 'im2col'可能可以爲你節省幾行代碼。 – Jonas 2010-05-20 01:30:58

+0

@Jonas:我添加了一個完全相同的解決方案! – Amro 2010-05-21 21:02:12

3

基於對JS的想法,並且如Jonas在評論所指出的,這可以在使用IM2COL一些數組操作兩行來完成:

B = im2col(C, [r r], 'sliding'); 
S = reshape(sum(B(1:r+1:end,:)), size(C)-r+1); 

基本上B包含所有滑動塊的元素的矩陣r-by-r在矩陣C上。然後我們將這些塊的每個塊的對角線上的元素B(1:r+1:end,:),計算它們的總和,並將結果重塑爲預期的大小。


比較這由喬納斯基於卷積的解決方案,這並不執行任何矩陣乘法,只有索引...

+0

如果你能負擔得起內存,im2col通常是一個很好的選擇。 – Jonas 2010-05-21 22:34:30

+0

im2col返回值而不是索引,所以第二行應該有'reshape(sum(B(1:r + 1:end,:))),size(C)-r + 1);'。 – 2010-11-26 00:27:21

+0

@reve_etrange:你說的很對,謝謝指出 – Amro 2010-11-26 05:20:44

1

這是你在找什麼?該函數添加對角線並將它們放入一個類似於函數'sum'將矩陣中的所有列相加並將其放入向量中的方式的向量。

function [diagSum] = diagSumCalc(squareMatrix, LLUR0_ULLR1) 
% 
% Input: squareMatrix: A square matrix. 
%  LLUR0_ULLR1: LowerLeft to UpperRight addition = 0  
%      UpperLeft to LowerRight addition = 1 
% 
% Output: diagSum: A vector of the sum of the diagnols of the matrix. 
% 
% Example: 
% 
% >> squareMatrix = [1 2 3; 
%     4 5 6; 
%     7 8 9]; 
% 
% >> diagSum = diagSumCalc(squareMatrix, 0); 
% 
% diagSum = 
% 
%  1 6 15 14 9 
% 
% >> diagSum = diagSumCalc(squareMatrix, 1); 
% 
% diagSum = 
% 
%  7 12 15 8 3 
% 
% Written by M. Phillips 
% Oct. 16th, 2013 
% MIT Open Source Copywrite 
% Contact [email protected] fmi. 
% 

if (nargin < 2) 
    disp('Error on input. Needs two inputs.'); 
    return; 
end 

if (LLUR0_ULLR1 ~= 0 && LLUR0_ULLR1~= 1) 
    disp('Error on input. Only accepts 0 or 1 as input for second condition.'); 
    return; 
end 

[M, N] = size(squareMatrix); 

if (M ~= N) 
    disp('Error on input. Only accepts a square matrix as input.'); 
    return; 
end 

diagSum = zeros(1, M+N-1); 

if LLUR0_ULLR1 == 1 
    squareMatrix = rot90(squareMatrix, -1); 
end 

for i = 1:length(diagSum) 
    if i <= M 
     countUp = 1; 
     countDown = i; 
     while countDown ~= 0 
      diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i); 
      countUp = countUp+1; 
      countDown = countDown-1; 
     end 
    end 
    if i > M 
     countUp = i-M+1; 
     countDown = M; 
     while countUp ~= M+1 
      diagSum(i) = squareMatrix(countUp, countDown) + diagSum(i); 
      countUp = countUp+1; 
      countDown = countDown-1; 
     end 
    end 
end 

乾杯