我正在將一些代碼從 Python 轉換為 Matlab。我的代碼可以產生相同的結果,但我想知道是否有辦法在 Matlab 代碼中向量化我的一些 for 回圈,因為它需要很長時間才能運行。X在 Nxd 矩陣中,diff是 NxNxd 張量,kxy是 NxN 矩陣,gradK是 NxNx2 張量,sumkxy,dxkxy, 和obj都是 Nxd 矩陣。
這是原始的 Python 代碼:
diff = x[:, None, :] - x[None, :, :] # D_{ij, s}
kxy = np.exp(-np.sum(diff ** 2, axis=-1) / (2 * h ** 2)) / np.power(np.pi * 2.0 * h * h, d / 2) # -1 last dimension K_{ij]
gradK = -diff * kxy[:, :, None] / h ** 2 # N * N * 2
sumkxy = np.sum(kxy, axis=1)
dxkxy = np.sum(gradK, axis=1) # N * 2 sum_{i} d_i K_{ij, s}
obj = np.sum(gradK / sumkxy[None, :, None], axis=1) # N * 2
這是我最初的帶有所有 for 回圈的 Matlab 代碼:
diff = zeros([n,n,d]);
for i = 1:n
for j = 1:n
for k = 1:d
diff(i,j,k) = x(i,k) - x(j,k);
end
end
end
kxy = exp(-sum(dif.^2, 3)/(2*h^2))/((2*pi*h^2)^(d/2));
sumkxy = sum(kxy,2);
gradK = zeros([n,n,d]);
for i = 1:n
for j = 1:n
for k = 1:d
gradK(i,j,k) = -diff(i,j,k)*kxy(i, j)/h^2;
end
end
end
dxkxy = squeeze(sum(gradK,2));
a = zeros([n,n,d]);
for i =1:n
for j = 1:n
for k = 1:d
a(i,j,k) = gradK(i,j,k)/sumkxy(i);
end
end
end
obj = squeeze(sum(a, 2));
我知道一種更快的計算kxy術語的方法是使用以下代碼:
XY = x*x';
x2= sum(x.^2, 2);
X2e = repmat(x2, 1, n);
H = (X2e X2e' - 2*XY); % calculate pairwise distance
Kxy = exp(-H/(2*h^2))/((2*pi*h*h)^(d/2));
但是后來我努力尋找一種方法,然后在gradK沒有diff. 任何幫助或建議將不勝感激!
uj5u.com熱心網友回復:
如果您的目標是計算obj您甚至不需要計算gradK和a:
sx = sum(x.^2, 2);
H = sx - 2*x*x.' sx.';
kxy = exp(-H/(2*h^2))/((2*pi*h^2)^(d/2));
kh = kxy / h^2;
sumkxy = sum(kxy, 2);
khs = kh ./ sumkxy;
obj = khs * x - sum(khs, 2) .* x;
gradK并且dif可以這樣計算:
dif = reshape(x, n, 1, d) - reshape(x, 1, n, d);
gradK = -dif .* (kxy / h^2);.
uj5u.com熱心網友回復:
我喜歡嘗試通過將其分解為“子組件”來嘗試解決此類問題,其中包含一些可以快速執行且可用于測驗代碼功能的虛假資料。您可能開始的第一個子組件是您的第一個計算差異的嵌套回圈:
n = 100;
d = 50;
x = round(100*rand(n,d));
tic
diff = zeros([n,n,d]);
for i = 1:n
for j = 1:n
for k = 1:d
diff(i,j,k) = x(i,k) - x(j,k);
end
end
end
toc
首先,單獨考慮最里面的回圈:
...
for k = 1:d
diff(i,j,k) = x(i,k) - x(j,k);
end
...
看看這個回圈(至少對我來說!)大大簡化了事情。為了僅矢量化這個“子組件”,我們可以撰寫如下內容:
diff(i,j,:) = x(i,:) - x(j,:);
既然懸而未決的果實已經解決,讓我們考慮下一層回圈。做和以前一樣的技巧有用嗎?
diff(i,:,:) = x(i,:) - x; % where x(:,:) can just be written as x.
如果您不確定,您可以通過使用相同(強調相同)虛假資料運行嵌套回圈版本和上面的版本并使用 isequal() 檢查它們是否相等來檢查這一點。為了切入正題,它應該是一樣的,現在你的原始回圈歸結為:
tic
diff = zeros([n,n,d]);
for i = 1:n
diff(i,:,:) = x(i,:) - x;
end
toc
對于最后一點,您可以利用 matlab 的矩陣/陣列整形/置換函式。查看 reshape() 或 permute() 的檔案以獲取更多詳細資訊。簡而言之,如果您將 x 的一個副本的維度順序從 Nxd 重塑或更改為 1xNxd,則從另一個規則大小的矩陣中減去 x 將在 matlab 中按元素執行操作。例如:
diff = permute(x,[1,3,2]) - permute(x,[3,1,2]); % this is Nx1xd - 1xNxd
應該有效地計算您在第一個回圈中尋找的張量差!
如果您愿意,我可以擴展此答案以顯示其他回圈如何解決,但請先使用相同的邏輯嘗試其他回圈。希望您可以保持差異,然后更快地計算 kxy。在不知道你的原始矩陣有多大的情況下,我不能說你應該期待多少加速。
更新:
我應該補充一點,為確保您正在執行元素乘法、除法和轉置運算,請確保添加一個“.”。在每個命令之前。例如
gradK(i,j,:) = -diff(i,j,:).*kxy(i, j)/h^2;
有關更多資訊,請在 Matlab 中查找元素操作
轉載請註明出處,本文鏈接:https://www.uj5u.com/gongcheng/361460.html
