我们如何才能知道 4x4 矩阵是否为奇异矩阵?
我们可以在不使用单位矩阵扩充给定矩阵然后进行行操作的情况下知道这一点吗?
那么如何确定一个矩阵是否真的是奇异的呢?(在 MATLAB 中,不使用纸笔或符号计算,或手写行操作?)教科书有时会告诉学生使用行列式,所以我将从那里开始。
理论上,可以简单地测试矩阵的行列式是否为零。因此
M = magic(4)
M =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
det(M)
ans =
-1.4495e-12
事实证明,这个矩阵确实是奇异的,所以有一种方法可以将 M 的一行写成其他行的线性组合(对于列也是如此。)但是我们从 det 得到的值不完全为零. 真的是零吗,MATLAB 刚刚搞糊涂了?这是为什么?符号行列式确实为零。
det(sym(M))
ans =
0
事实证明,行列式的计算对于较大的数组来说是一件非常低效的事情。因此,一个不错的选择是使用我们的方阵的特定矩阵分解的对角元素的乘积。事实上,这就是 MATLAB 在 det 内部对非符号输入所做的事情。
[L,U,P] = lu(M)
L =
1 0 0 0
0.25 1 0 0
0.3125 0.76852 1 0
0.5625 0.43519 1 1
U =
16 2 3 13
0 13.5 14.25 -2.25
0 0 -1.8889 5.6667
0 0 0 3.5527e-15
P =
1 0 0 0
0 0 0 1
0 1 0 0
0 0 1 0
看到 L 的对角元素是 1,但 U 有非零对角元素。矩阵乘积的行列式以及上三角矩阵或下三角矩阵的行列式也有很好的性质。
prod(diag(U))
ans =
-1.4495e-12
看到我们得到了相同的答案。由于 LU 分解即使对于大型数组也相当快,因此它是计算行列式的好方法。问题是,它使用浮点运算。所以 U 的那些对角元素是实数,浮点数。当我们取产品时,我们得到的结果并不完全为零。最后一个元素只是略微非零。
det还有其他问题。例如,我们知道矩阵 eye(100) 的条件非常好。毕竟它是一个单位矩阵。
det(eye(100))
ans =
1
现在,如果我们将矩阵乘以一个常数,这不会改变矩阵作为奇异矩阵的状态。但是行列式会发生什么?
det(.1*eye(100))
ans =
1e-100
那么这个矩阵是奇异的吗?毕竟,如果 det(magic(4)) 给我们 1e-12,那么 1e-100 必须对应一个奇异矩阵!但事实并非如此。更糟,
det(.0001*eye(100))
ans =
0
事实上,行列式按 0.0001^100 缩放,在 matlab 中为 1e-400。至少如果 matlab 可以使用双精度表示一个那么小的数字,那将是正确的。它不能这样做。这个数字将下溢。或者同样容易,我们可以让它溢出。
det(10000*eye(100))
ans =
Inf
显然,所有这些缩放的单位矩阵都是同样非奇异的,但是可以使 det 给我们任何我们想看到的答案!因此,我们必须得出结论,计算行列式对矩阵来说是一件可怕的事情。我不在乎你的教科书很久以前告诉你什么,或者你的老板或老师告诉你什么。如果有人告诉您为此目的使用计算机计算行列式,那是一个糟糕的建议。时期。行列式只是有太多的问题。
我们可以做其他事情来测试奇点。最好的工具是使用排名。因此,如果 NxM 矩阵的秩小于 min(N,M),则该矩阵是奇异的。这里有几个测试:
rank(M)
ans =
3
rank(.0001*eye(100))
ans =
100
所以 rank 能够告诉我们 4x4 幻方是奇异的,但我们的缩放单位矩阵不是奇异的。(一件好事是秩可以测试非方阵的奇异性。)
我们还可以使用 cond 来测试数值奇异性。最小可能的条件数是 1.0,它对应于一个表现非常好的矩阵。大的条件数是不好的。在双精度中,这意味着当条件数大于大约 1e15 时,您的矩阵非常有问题。
cond(M)
ans =
8.148e+16
cond(.0001*eye(100))
ans =
1
事实上,cond 认为缩放的单位矩阵毕竟条件良好。大的条件数是不好的。对于双精度矩阵,接近 1e15 左右的条件数表示可能是数值奇异的矩阵。所以我们看到 M 显然是奇异的。同样, cond 能够处理非方阵。
或者,我们可以使用 rcond,一个估计条件数倒数的工具。对于非常大的数组,这是一个很好的工具。当 rcond 给出的数字接近 eps 时,请注意!
rcond(M)
ans =
1.3061e-17
rcond(.0001*eye(100))
ans =
1
最后,对于数学齿轮头(像我一样),我们可能会退出 svd。毕竟,svd 是 cond 和 rank 所依据的工具。当矩阵的一个或多个奇异值与最大奇异值相比很小时,我们再次具有奇异性。
svd(M)
ans =
34
17.889
4.4721
4.1728e-16
在这里,我们看一下与矩阵的最大奇异值相比,奇异值何时较小。一件好事是 svd 可以告诉我们矩阵与奇异点的接近程度,如果有多个小的奇异值,则可以告诉我们有关矩阵秩的信息。
好消息是我展示的所有工具都不需要用户进行基本的行操作或任何花哨的操作。
但请不要使用 DET!是的,它出现在教科书中。是的,也许你的导师或你的老板告诉你使用它。他们是错的,因为像这样的工具在使用浮点运算的计算机上应用时会失败。而且您根本不想计算符号行列式,这将非常低效。
很抱歉,如果这是一个长篇阅读。我现在要离开我的肥皂盒。
计算排名并与维度进行比较。如果秩低于维度,则矩阵是奇异的。
最可靠的方法是对矩阵进行奇异值分解。最大奇异值与最小奇异值之比应在一些合理的公差范围内。这个比率是矩阵的条件数。对于双精度值,当条件数超过一百万或更多时,双精度值会变得非常危险,这是一个相当高的限制。请注意,一旦你有了 SVD,它对很多其他事情都很有用,而不仅仅是计算条件数。
奇异值分解是数值分析的瑞士军用电锯;如果您知道矩阵不是奇异/病态的,那么它可能会有点笨拙。但如果你不知道,这是一个很好的了解工具。Matlab 尤其如此,因为它是一个内置工具。
我会用cond
. 这为您提供了一个矩阵与奇异矩阵的接近程度的数值估计(其中 Inf 是奇异矩阵)。
例如:
m = randn(4);
cond(m) %Well conditioned, usually in the 10's
m = diag([1e-6 1 2 1e6]);
cond(m) %Less well conditioned, 1e12
m = diag([0 1 2 3]);
cond(m) %Singular: Inf