矢量化矩陣乘法

我有一個關於在Octave上運行SOR算法的學校項目，但是我的效率很低。所以我有這段代碼：矢量化矩陣乘法

for ii=1:n 
    r = 1/A(ii,ii); 
    for jj=1:n 
     if (ii!=jj) 
      A(ii,jj) = A(ii,jj)*r; 
     end; 
    end; 
    b(ii,1) = b(ii,1)*r; 
    x(ii,1) = b(ii,1); 
end;

我該如何實現這一點？我的第一次嘗試是這樣的：

for ii=1:n 
    r = 1/A(ii,ii); 
    A(find(eye(length(A))!=1)) = A(find(eye(length(A))!=1))*r; 
    b(ii,1) = b(ii,1)*r; 
    x(ii,1) = b(ii,1); 
end;

但我不確定它有多大的幫助。有沒有更好的和/或更有效的方法來做到這一點？

謝謝！

來源

2013-04-15 dccarmo

你能解釋代碼實際在做什麼嗎？（一些矩陣分解可能......） – tashuhka

它創建了SOR算法的迭代矩陣。 – dccarmo

發佈了一個基於bsxfun的解決方案[here]（http：// stackoverflow。com/a/26128876/3293881）如果你想看看！ – Divakar

你完全可以避免循環我相信。你必須看到你正在採取A的對角元素的倒數，那麼你爲什麼要使用循環。直接做。這是第一步。刪除Inf，現在您想將r乘以相應行的相應非對角元素，對嗎？因爲您將與（1,2），（1,3），...，（1，n）相同的r相乘，因此使用repmat來構造這樣的矩陣，其將沿着列具有複製的元素。但是R具有非零對角線元素。因此使它們爲零。現在你會得到你的A，但對角元素將爲零。因此，您只需將其從原始A中添加回來即可。這可以通過A=A.*R+A.*eye(size(A,1))完成。

矢量化來自經驗，最重要的是分析你的代碼。想想在每一步是否要使用循環，如果不是用等效命令替換該步驟，其他代碼將會跟隨（例如，我構建了一個矩陣R，而您正在構建單個元素r。所以我只考慮轉換r - >R然後剩下的代碼就會落到位）。

的代碼如下：

R=1./(A.*eye(size(A,1))); %assuming matrix A is square and it does not contain 0 on the main diagonal 
R=R(~isinf(R)); 
R=R(:); 
R1=R; 
R=repmat(R,[1,size(A,2)]); 
R=R.*(true(size(A,1))-eye(size(A,1))); 
A=A.*R+A.*eye(size(A,1)); %code verified till here since A comes out to be the same 

b = b.*R1; 
x=b;

來源

2013-04-16 00:35:44

太棒了，Parag！它確實有用！你有什麼資源可以學習更多關於矢量化的知識嗎？我很難說明你做了什麼，並且還想將我在代碼中的其他循環向量化！ – dccarmo

@dccarmo我剛查過。 'b'不出來是正確的。你能明白爲什麼，因爲你比我更瞭解代碼。或者我犯了一些錯誤。我會嘗試通過編輯我的答案來告訴邏輯 –

我剛剛在這裏檢查，而且b和原始代碼一樣，我檢查它是否錯誤？ – dccarmo

我假設有矩陣：

A (NxN) 
b (Nx1)

的代碼：

d = diag(A); 
A = diag(1 ./ d) * A + diag(d - 1); 
b = b ./ d; 
x = b;

來源

2013-04-16 17:58:27 Serg

隨機碰着到這個問題和第一眼看起來考慮到這個問題被標記爲vectorization問題，所以感興趣。

我能夠想出一個bsxfun的矢量化解決方案，也使用diagonal indexing。這個解決方案似乎給了我一個3-4x超過循環代碼與體面大小的輸入。

假設你對於加速這個問題的速度有所提高仍然有點感興趣，那麼我會非常想知道你將會獲得哪些加速。這裏的代碼 -

diag_ind = 1:size(A,1)+1:numel(A); 
diag_A = A(diag_ind(:)); 
A = bsxfun(@rdivide,A,diag_A); 
A(diag_ind) = diag_A; 
b(:,1) = b(:,1)./diag_A; 
x(:,1) = b(:,1);

讓我知道！

來源

2014-09-30 19:47:04 Divakar

矢量化矩陣乘法

回答

相關問題