6

我有一个具有 3 维 Y(i,j,w) 的矩阵。我想得到一个行列式向量 d(w)，其中每个数字都是矩阵 Y(:,:,w) 的行列式。

是否有一个优雅的语法，或者我只需要使用一个循环？

谢谢

4

2 回答 2

8

嗯，首先，你实际上从来没有真正想要计算一个行列式，你只是认为你这样做了。事实上，这几乎从来都不是一件好事，因为行列式的尺度太差了。它们经常被用来推断矩阵的奇点状态，这在数值分析方面是一件很糟糕的事情。

说了我对一般决定因素的迷你咆哮......

选项1：

将您的 3-d 数组转换为方阵元胞数组，阵列的每个平面为一个元胞。mat2cell 将轻松高效地完成任务。

接下来，在元胞数组上使用 cellfun。cellfun 可以将函数 (@det) 应用于每个单元格，然后将返回一个行列式向量。这是非常有效的吗？只要在编写循环时提前预先分配向量，它可能不会比在循环中应用 det 获得巨大收益。

选项 2：

如果矩阵很小，例如 2x2 或 3x3 矩阵，则将行列式的乘法展开为显式向量乘法。我认为这在我写它的时候还不清楚，所以对于一个 2x2 的情况，其中 Y 是 2x2xn：

d = Y(1,1,:).*Y(2,2,:) - Y(1,2,:).*Y(2,1,:);

您肯定会看到，这为矩阵 Y 的每个平面形成了一个由 2x2 行列式组成的向量。3x3 的情况也很简单，也可以写成六个 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，额外的开销是一个杀手。（我不会猜测盈亏平衡点落在哪里，但可以很容易地测试它。）

于 2012-06-02T11:43:10.917 回答
4

我会使用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.

结论：幼稚的方法是最快的。

于 2012-06-03T15:34:37.263 回答