此行为在的文档的限制部分中进行了解释,并在查找奇异矩阵的行列式小节中进行了举例说明:det
的行列式A
很大,尽管它A
是单数的。事实上, 的行列式A
应该正好为零!的不准确d
是由于 LU 分解的 MATLAB® 实现中舍入误差的聚合,该 LU 分解det
用于计算行列式。
也就是说,在这种情况下,您可以通过使用同一页面上给出的m 代码实现来产生您想要的结果,但对角元素U
按升序排序。考虑示例脚本:
clc();
clear();
A = rand(500,1500);
b = rand(500,1);
c = (A.')*A;
[L,U] = lu(c);
% Since det(L) is always (+/-)1, it doesn't impact anything
diagU = diag(U);
detU1 = prod(diagU);
detU2 = prod(sort(diagU,'descend'));
detU3 = prod(sort(diagU,'ascend'));
fprintf('Minimum: %+9.5e\n',min(abs(diagU)));
fprintf('Maximum: %+9.5e\n',max(abs(diagU)));
fprintf('Determinant:\n');
fprintf('\tNo Sort: %g\n' ,detU1);
fprintf('\tDescending Sort: %g\n' ,detU2);
fprintf('\tAscending Sort: %g\n\n',detU3);
这将产生输出:
Minimum: +1.53111e-13
Maximum: +1.72592e+02
Determinant:
No Sort: Inf
Descending Sort: Inf
Ascending Sort: 0
请注意,排序的方向很重要,并且由于对角线上不存在Inf
真值,因此无排序给出。0
降序排序首先看到最大值相乘,显然,它们超过realmax
并且永远不会乘以 true 0
,这将生成 a NaN
。升序排序将所有接近零的对角线值与极少数大的负值聚集在一起(实际上,更稳健的方法会根据幅度进行排序,但这里没有这样做),它们的乘法生成一个真值0
(意味着该值低于 IEEE-754 算术中可用的最小非规格化数)产生“正确”结果。
所有这些,正如其他人所暗示的那样,我将引用原始 Matlab 开发人员和 Mathworks 联合创始人 Cleve Moler 的话:
[行列式] 在理论考虑和手工计算中很有用,但不能为稳健的数值软件提供良好的基础。