作者celestialgod (天)
看板MATLAB
標題Re: [問題]三維矩陣對二維矩陣的擬合
時間Tue Sep 8 19:52:51 2015
第三個是你的想法直接用兩層迴圈做的
前兩個可以加速不少，提供參考
第一個方法比較tricky的方式把你的想法用成矩陣方式去解會快不少
第二個就只是改成cellfun版本，不用preallocate D
% data generation
A = randi(255, 128, 128, 20);
B = rand(20 ,6);
% permute + for-loop
tic
D = zeros(128, 128, 6);
A2 = permute(A, [3,2,1]);
for i = 1:128
D(i, :, :) = (B \ squeeze(A2(:, :, i)))';
% 可以不要先permute 直接這樣做：D(i, :, :) = (B \ squeeze(A(i, :, :))')';
% 速度上沒差多少，看個人寫法而定
end
toc % Elapsed time is 0.035446 seconds.
% cellfun
tic
tmp_cell = mat2cell(A, ones(size(A, 1), 1), size(A,2), size(A,3));
tmp_cell = cellfun(@squeeze, tmp_cell, 'UniformOutput', false);
D2 = cellfun(@(x) B \ x', tmp_cell, 'UniformOutput', false);
D2 = permute(cat(3, D2{:}), [3,2,1]);
toc % Elapsed time is 0.047776 seconds.
% double for-loop
tic
D3 = zeros(128, 128, 6);
for i = 1:128
for j = 1:128
D3(i, j, :) = B \ squeeze(A(i, j, :));
end
end
toc % Elapsed time is 0.793275 seconds.
all(all(all(abs(D - D2) < 1e-4))) % 1
all(all(all(abs(D - D3) < 1e-4))) % 1 %% 精度會不同，allequal會給錯
※ 引述《victor6954 (維克)》之銘言：
: 大家好，我是初入matlab世界的新手
: 有個問題想要請教各位
: 我現在有一個 128 X 128 X 20 的三維 A 矩陣，
: 以及一個 6 X 20 的二維 B 矩陣，
: 想要由這兩個矩陣求出三維 D 矩陣，
: 其式為 A = B X D。
: 目前嘗試取出 A 矩陣的一條 Z 方向轉換為 A1 = 20 X 1 的矩陣後，
: 用D=B\A，
: 可以求出 6 X 1 的矩陣
: http://imgur.com/VIVjA27
: 然而，卻無法直接計算三維矩陣 (error : Input arguments must be 2-D.)
: 與二維的差別是，A與D矩陣，每個單一值變成128X128的矩陣
: 想請問要如何解出三維 D 矩陣
: http://imgur.com/NTAb6Ln
: A矩陣應該要如何進行轉換，或是是否有直接計算的語法
: 謝謝各位了
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 111.248.8.123
※ 文章網址: https://webptt.com/m.aspx?n=bbs/MATLAB/M.1441713174.A.4F7.html
1F：推 sunev: B的(pseudo) inverse可以先求。 09/08 23:16
2F：推 profyang: 雖然我同意樓上 但pinv(B)*A和B\A未必會有一樣的結果 09/08 23:35
3F：→ profyang: 但以這問題先求好pinv(B)然後去掃20*6這兩個維度應該比 09/08 23:36
4F：→ profyang: 較快 反正pinv和\都是求最小方差 09/08 23:36
5F：推 sunev: 都是求最小方差的話，答案應該一樣啊？ 09/08 23:51
6F：→ sunev: t=pinv(B)*reshape(permute(A,[3 1 2]),20,128*128); 09/08 23:52
7F：→ sunev: t=reshape(t,20,128,128); 09/08 23:53
8F：→ sunev: 上一行錯了，應為D=reshape(t,6,128,128); 09/08 23:54
10F：→ profyang: 你B本身不是full rank的話就不會一樣...雖然它20*6要不 09/09 00:01
11F：→ profyang: 是full rank難度也是頗大啦... 09/09 00:02
12F：→ profyang: 而且樓上都那樣打了 跟打t=B\reshape...不是差不多快? 09/09 00:03
13F：推 sunev: 也是 XD 09/09 00:12
14F：→ sunev: 所以 pinv 是最小方差解，"\" 是最多非零元素解，大概是 09/09 00:13
15F：→ sunev: gauss elimination演算法的結果吧 09/09 00:13
16F：推 profyang: 不是耶 是如果B不是full rank的時候(如網頁例子) 它的 09/09 00:15
17F：→ profyang: 最小方差解就不只一組 那這兩個方法結果會不一樣 但都是 09/09 00:15
18F：→ profyang: 最小方差 但是B\A的非0項會最少 至於它此時B\A怎麼算的 09/09 00:16
19F：→ profyang: 我就也不確定了... 09/09 00:16
20F：→ profyang: 上面忘了打 pinv是讓x的norm最小(不是方"差"最小喔~) 09/09 00:18
21F：推 sunev: 我把underdetermined 和overdetermined 搞混了 09/09 00:29
22F：→ sunev: overdetermined，二者都會給出唯一最小方差解。 09/09 00:29
23F：→ sunev: underdetermined，解有無限多組，pinv給出norm最小的解 09/09 00:30
24F：→ sunev: "\"給出非零元素最多的解 09/09 00:30
想不到這個問題 可以討論那麼熱烈@@...
辛苦p大跟s大了
※ 編輯: celestialgod (111.248.8.123), 09/09/2015 00:48:24
25F：推 victor6954: 解出來了!!非常感謝您!!! 09/09 17:36