我有一个3维矩阵Y(i,j,w)。 我想得到一个行列式向量d(w),其中每个数字都是矩阵行列式的行列式(:,:,w)。返回行列式向量 - Matlab
它是否有一个优雅的语法,或者我只需要使用循环?
感谢
我有一个3维矩阵Y(i,j,w)。 我想得到一个行列式向量d(w),其中每个数字都是矩阵行列式的行列式(:,:,w)。返回行列式向量 - Matlab
它是否有一个优雅的语法,或者我只需要使用循环?
感谢
嗯,首先,你几乎从来没有真正要计算一个决定因素,你只是觉得你这样做。事实上,它几乎从来都不是一件好事,因为行列式的缩放比例很小。通常他们被用来推断矩阵的奇点状态,这在数值分析方面是一件可怕的事情。
在阐述对一般决定我的小咆哮......
选项1:
转换您的3-d阵列成方阵的单元阵列,与阵列为一体的每个平面细胞。 mat2cell将轻松高效地完成这项技巧。
接下来,在cell阵列上使用cellfun。 cellfun可以为每个单元格应用一个函数(@det),然后返回一个行列式向量。这是非常有效的吗?在循环中应用det可能不是一个巨大的收益,只要您在编写循环时提前预先分配vector。
选项2:
如果基质是小的,因此说2×2 3×3或基质,然后展开出来的行列式作为明确的矢量相乘的乘法。我觉得这是我写的,所以对于一个2×2的情况下,其中Y是2x2xn不明确:
d = Y(1,1,:).*Y(2,2,:) - Y(1,2,:).*Y(2,1,:);
当然,你看,这形成2x2的决定因素的向量矩阵Y的每一个面3x3的情况很简单,可以写成6个3路产品。我没有仔细检查下面的3x3情况,但它应该接近。
d = Y(1,1,:).*Y(2,2,:).*Y(3,3,:) + ...
Y(2,1,:).*Y(3,2,:).*Y(1,3,:) + ...
Y(3,1,:).*Y(1,2,:).*Y(2,3,:) - ...
Y(3,1,:).*Y(2,2,:).*Y(1,3,:) - ...
Y(2,1,:).*Y(1,2,:).*Y(3,3,:) - ...
Y(1,1,:).*Y(3,2,:).*Y(2,3,:);
正如你所看到的,选项2会很快,而且是矢量化的。
编辑:作为对克里斯的回应,所需时间存在显着差异。考虑一组1e5矩阵所需的时间。
p = 2;
n = 1e5;
Y = rand(p,p,n);
tic,
d0 = squeeze(Y(1,1,:).*Y(2,2,:) - Y(2,1,:).*Y(1,2,:));
toc
Elapsed time is 0.002141 seconds.
tic,
X = squeeze(mat2cell(Y,p,p,ones(1,n)));
d1= cellfun(@det,X);
toc
Elapsed time is 12.041883 seconds.
这两个调用将相同的值返回到浮点垃圾内。
std(d0-d1)
ans =
3.8312e-17
循环不会更好,事实上肯定会更糟。所以我要写一段代码来面对在数组中产生许多这样的矩阵的决定因素的任务,我会特别说明2x2和3x3矩阵的代码。我甚至可以写出4x4矩阵。是的,写出来很麻烦,但所需时间有很大差异。
一个原因是MATLAB的det使用了一个调用LU的因子分解矩阵。理论上这比中等大矩阵的乘法更好,但对于2x2或3x3,额外的开销是一个杀手。 (我不会猜测盈亏平衡点落在哪里,但可以轻松测试。)
我会使用arrayfun:
d = arrayfun(@(w) det(Y(:, :, w)), 1 : size(Y, 3));
编辑:速度测试:
p = 10;
n = 1e4;
Y = rand(p,p,n);
测试1:
>> tic, d1 = arrayfun(@(w) det(Y(:, :, w)), 1 : size(Y, 3)); toc
Elapsed time is 0.139030 seconds.
试验2(由木屑):
>> tic, X = squeeze(mat2cell(Y,p,p,ones(1,n))); d2= cellfun(@det,X); toc
Elapsed time is 1.318396 seconds.
试验3(天真的方法):
>> p = 10;
>> n = 1e4;
>> Y = rand(p,p,n);
>> tic; d = nan(n, 1); for w = 1 : length(d), d(w) = det(Y(:, :, w)); end; toc
Elapsed time is 0.069279 seconds.
结论:幼稚的做法是最快的。
我更喜欢选项1或者甚至选项2的循环。如果循环上解释器的速度损失是重要的,我宁愿使用更灵活的代码。想象一下,用这种方式做一个5x5的行列式 –
@ChrisA - 是12.041883/0.002141 = 5624.4一个令人震惊的区别? – 2012-06-02 17:08:16
+1是!好的回应。也许有一个更灵活的方法来做到这一点...因为它使得术语的数量变得更大,所以'p!' –