作者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/cn.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