如何向量化相關For-Loop

我正在研究一個函數，它需要一個1xn向量x作爲輸入並返回一個nxn矩陣L。
我想通過向量化循環來加快速度，但有一個令我困惑的問題：循環索引b取決於循環索引a。任何幫助，將不勝感激。如何向量化相關For-Loop

x = x(:); 
n = length(x); 
L = zeros(n, n); 
for a = 1 : n, 
    for b = 1 : a-1, 
     c = b+1 : a-1; 
     if all(x(c)' < x(b) + (x(a) - x(b)) * ((b - c)/(b-a))), 
      L(a,b) = 1; 
     end 
    end 
end

來源

2013-08-27 bluebox

旁註用於在Matlab中的變量：[http://stackoverflow.com/questions/14790740/using-i-and-j-as-variables-in-matlab](http://stackoverflow.com/questions/14790740/using-i-和-j-as-variables-in-matlab） –

感謝您的提示。我已將「i」，「j」，「k」更名爲「a」，「b」，「c」。 – bluebox

我不確定你是否意識到這一點，但所有（[]）評估爲真。所以，例如，當一個== 2，b將是1，並且c將是[]。這將最終使L（a，b）1。你可能需要考慮你的b和c範圍（我可能是錯的;也許這就是你打算如何工作）。 – ioums

從一個快速測試，它看起來像你正在做的事情只與下三角。您可能能夠使用醜陋的技巧，比如ind2sub和arrayfun類似的向量化這個

tril_lin_idx = find(tril(ones(n), -1)); 
[A, B] = ind2sub([n,n], tril_lin_idx); 
C = arrayfun(@(a,b) b+1 : a-1, A, B, 'uniformoutput', false); %cell array 
f = @(a,b,c) all(x(c{:})' < x(b) + (x(a) - x(b)) * ((b - c{:})/(b-a))); 
L = zeros(n, n); 
L(tril_lin_idx) = arrayfun(f, A, B, C);

我無法測試它，因爲我沒有x，我不知道預期的結果做。我通常喜歡矢量化的解決方案，但這可能會推動它太多:)。我會堅持你的明確的for循環，這可能會更清晰，哪些Matlab的JIT應該能夠輕鬆加速。你可以用L(a,b) = all(...)代替if。

EDIT1

更新版本，以防止C浪費~ n^3空間：

tril_lin_idx = find(tril(ones(n), -1)); 
[A, B] = ind2sub([n,n], tril_lin_idx); 
c = @(a,b) b+1 : a-1; 
f = @(a,b) all(x(c(a, b))' < x(b) + (x(a) - x(b)) * ((b - c(a, b))/(b-a))); 
L = zeros(n, n); 
L(tril_lin_idx) = arrayfun(f, A, B);

EDIT2

輕微變形，不使用ind2sub並且應該更容易在b的情況下將以更復雜的方式取決於a。我爲速度內聯c，似乎特別是調用函數句柄很貴。

[A,B] = ndgrid(1:n); 
v = B<A; % which elements to evaluate 
f = @(a,b) all(x(b+1:a-1)' < x(b) + (x(a) - x(b)) * ((b - (b+1:a-1))/(b-a))); 
L = false(n); 
L(v) = arrayfun(f, A(v), B(v));

來源

2013-08-27 21:45:47

謝謝，你的向量化代碼返回預期的結果 - 然而，（令我驚訝的）它並沒有像for循環版本那麼快。 Btw：你說得對，我只在一個矩陣三角形上循環 - 另一個三角形將像這樣創建：''L = tril（L，-1）+ tril（L，-1）''' – bluebox

Matlab的JIT非常善於優化'fortran風格'for循環，這是很難打敗的。一個常見的抱怨是，arrayfun有點慢... –

矢量化版本的另一個問題是，你必須預先計算'A'，'B'和'C'。特別是後者可能會佔用大量無用的內存，如果'x'很大。它的大小可能與'n^3'一致，哎喲。你的for循環不會遇到這個問題。爲了解決這個問題，你可以嘗試在每次迭代時計算'c'，我將編輯答案... –

如果我正確理解你的問題，L(a, b) == 1如果有<ç< B中所述C，（C，X（C））是「下面」的線連接（A，X（a））和（b，x（b）），對嗎？

這不是一個矢量化，但我找到了另一種方法。（a，b）中，我保存了從a到c的最大斜率，並將其用於（a，b + 1），而不是將所有c與每個新b的所有c進行比較。（我想只有一個方向，但我覺得用兩個方向也是可能的。）

x = x(:); 
n = length(x); 
L = zeros(n); 

for a = 1:(n - 1) 
    L(a, a + 1) = 1; 
    maxSlope = x(a + 1) - x(a); 
    for b = (a + 2):n 
    currSlope = (x(b) - x(a))/(b - a); 
    if currSlope > maxSlope 
     maxSlope = currSlope; 
     L(a, b) = 1; 
    end 
    end 
end

我不知道你的數據，但也有一些隨機數據，結果與原來的代碼相同（與轉置）。

來源

2013-08-28 10:23:07 dlimpid

不錯的工作，在我的筆記本電腦上（使用Matlab R2012b），它大約是原始代碼的兩倍。 – bluebox

真的嗎？當我檢查運行時時，可能是錯誤的'tic; toc'。 – dlimpid

對於''n = 100''我的代碼通常需要約0.04秒，而你的約0.02秒完成。 – bluebox

一個深奧的答案：你可以從1：n對每個a，b，c進行計算，排除不關心的事情，然後在c維上進行。

[a, b, c] = ndgrid(1:n, 1:n, 1:n); 

La = x(c)' < x(b) + (x(a) - x(b)) .* ((b - c)./(b-a)); 
La(b >= a | c <= b | c >= a) = true; 

L = all(La, 3);

儘管jit對於for循環可能會做得很好，因爲它們做的很少。

編輯：仍然使用所有的存儲器，但具有較少的數學

[A, B, C] = ndgrid(1:n, 1:n, 1:n); 

valid = B < A & C > B & C < A; 
a = A(valid); b = B(valid); c = C(valid); 

La = true(size(A)); 
La(valid) = x(c)' < x(b) + (x(a) - x(b)) .* ((b - c)./(b-a)); 
L = all(La, 3);

EDIT2：備用最後一行添加子句C無元素是真實的

L = all(La,3) | ~any(valid,3);

來源

2013-08-29 03:42:29 RyanD

幹得好！這需要大約3倍於我的機器上的循環版本...我想Mathworks的人真的做了很多工作來提高Matlab中循環的速度。否則，我無法向自己解釋爲什麼所有偉大的矢量化版本似乎都比較慢...... – bluebox

這個版本可能很慢的一個原因是'a'，'b'，'c'和'La'都採用' O（n^3）'空間，而for循環只使用'O（n^2）'空間。 –

如何向量化相關For-Loop

回答

相關問題