行列式显示无穷大而不是零!为什么?
Determinant is showing infinity instead of zero! Why?
这是我为作业中遇到的问题编写的 matlab 代码。根据所有同学的代码(不同的代码),将 A 与其转置相乘后,所得方阵的行列式应为零。为什么我的代码没有给出 c 和 d 的行列式为无穷大
A = rand(500,1500);
b = rand(500,1);
c = (A.')*A;
detc = det(c);
cinv = inv((A.')*A);
d = A*(A.');
detd = det(d);
dinv = inv(A*(A.'));
x1 = (inv((A.')*A))*((A.')*b);
x2 = A.'*((inv(A*(A.')))*b);
我认为这几乎只是一个舍入问题,Inf 并不意味着您得到无穷大作为答案,只是您的行列式确实很大并且超过了 realmax
。正如阿迪尔所说,A*A。生成一个对称矩阵,并且它的行列式应该有一个数值。例如,设置:
A=rand(5,15)
你应该会发现 A*A 的数据。'只是一个数值。
所以你的朋友是怎么得到零的,好吧,对于 det
的大矩阵很容易得到 0 或 inf(你为什么首先要这样做我不知道)。所以我认为他们只是遇到了 same/similar 舍入问题。
好的。所以 det(A'*A) 不为零的事实并不能很好地表明 A'*A 的(非)奇异性。
行列式取决于缩放比例,显然非奇异矩阵可以有非常小的行列式。例如,矩阵
1/2 * I_n
其中 I_n 是 nxn 恒等式,其行列式为 (1/2)^n,随着 n 趋于无穷大,它会(快速)收敛到 0。但是 1/2 * I_n 根本不是单数。
因此,检查矩阵奇异性的最佳方法是条件数。
在你的情况下,在做了一些测试后
>> A = rand(500, 1500) ;
>> det(A'*A)
ans =
Inf
您可以看到(计算出的)行列式显然不为零。但这实际上并不奇怪,也不应该真正打扰您。行列式很难计算,所以是的,它只是舍入误差。如果你想要一个更好的近似值,你可以这样做
>> s = eig(A'*A) ;
>> prod(s)
ans =
0
在那里,你看到它更接近于零。
另一方面,条件数是矩阵(非)奇异性的更好估计。在这里,是
>> cond(A'*A)
ans =
1.4853e+20
而且,由于它比 1e+16 大得多,所以矩阵显然是奇异的。 1e+16的原因有点繁琐,但主要是计算机在做浮点运算时的精度问题。
此行为在 Limitations section of the det
's documentation and exemplified in the Find Determinant of Singular Matrix subsection 中有说明:
The determinant of A
is quite large despite the fact that A
is singular. In fact, the determinant of A
should be exactly zero! The inaccuracy of d
is due to an aggregation of round-off errors in the MATLAB® implementation of the LU decomposition, which det
uses to calculate the determinant.
就是说,在这种情况下,您 可以 通过使用同一页面上给出的 m-code implementation 但对 [=18 的对角线元素进行排序来产生您想要的结果=] 在一个上升的问题上。考虑示例脚本:
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
并且永远不会乘以真正的 0
,这将生成 NaN
。升序排序将所有接近零的对角线值与非常少的大负值聚集在一起(事实上,更稳健的方法会根据大小排序,但这里没有这样做),它们的乘法生成一个真正的 0
(意味着该值低于 IEEE-754 算术中可用的最小非规范化数字)产生 "correct" 结果。
写了这么多,正如其他人暗示的那样,我将 quote 原始 Matlab 开发人员和 Mathworks 联合创始人 Cleve Moler:
[The determinant] is useful in theoretical considerations and hand calculations, but does not provide a sound basis for robust numerical software.
这是我为作业中遇到的问题编写的 matlab 代码。根据所有同学的代码(不同的代码),将 A 与其转置相乘后,所得方阵的行列式应为零。为什么我的代码没有给出 c 和 d 的行列式为无穷大
A = rand(500,1500);
b = rand(500,1);
c = (A.')*A;
detc = det(c);
cinv = inv((A.')*A);
d = A*(A.');
detd = det(d);
dinv = inv(A*(A.'));
x1 = (inv((A.')*A))*((A.')*b);
x2 = A.'*((inv(A*(A.')))*b);
我认为这几乎只是一个舍入问题,Inf 并不意味着您得到无穷大作为答案,只是您的行列式确实很大并且超过了 realmax
。正如阿迪尔所说,A*A。生成一个对称矩阵,并且它的行列式应该有一个数值。例如,设置:
A=rand(5,15)
你应该会发现 A*A 的数据。'只是一个数值。
所以你的朋友是怎么得到零的,好吧,对于 det
的大矩阵很容易得到 0 或 inf(你为什么首先要这样做我不知道)。所以我认为他们只是遇到了 same/similar 舍入问题。
好的。所以 det(A'*A) 不为零的事实并不能很好地表明 A'*A 的(非)奇异性。 行列式取决于缩放比例,显然非奇异矩阵可以有非常小的行列式。例如,矩阵 1/2 * I_n 其中 I_n 是 nxn 恒等式,其行列式为 (1/2)^n,随着 n 趋于无穷大,它会(快速)收敛到 0。但是 1/2 * I_n 根本不是单数。
因此,检查矩阵奇异性的最佳方法是条件数。
在你的情况下,在做了一些测试后
>> A = rand(500, 1500) ;
>> det(A'*A)
ans =
Inf
您可以看到(计算出的)行列式显然不为零。但这实际上并不奇怪,也不应该真正打扰您。行列式很难计算,所以是的,它只是舍入误差。如果你想要一个更好的近似值,你可以这样做
>> s = eig(A'*A) ;
>> prod(s)
ans =
0
在那里,你看到它更接近于零。
另一方面,条件数是矩阵(非)奇异性的更好估计。在这里,是
>> cond(A'*A)
ans =
1.4853e+20
而且,由于它比 1e+16 大得多,所以矩阵显然是奇异的。 1e+16的原因有点繁琐,但主要是计算机在做浮点运算时的精度问题。
此行为在 Limitations section of the det
's documentation and exemplified in the Find Determinant of Singular Matrix subsection 中有说明:
The determinant of
A
is quite large despite the fact thatA
is singular. In fact, the determinant ofA
should be exactly zero! The inaccuracy ofd
is due to an aggregation of round-off errors in the MATLAB® implementation of the LU decomposition, whichdet
uses to calculate the determinant.
就是说,在这种情况下,您 可以 通过使用同一页面上给出的 m-code implementation 但对 [=18 的对角线元素进行排序来产生您想要的结果=] 在一个上升的问题上。考虑示例脚本:
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
并且永远不会乘以真正的 0
,这将生成 NaN
。升序排序将所有接近零的对角线值与非常少的大负值聚集在一起(事实上,更稳健的方法会根据大小排序,但这里没有这样做),它们的乘法生成一个真正的 0
(意味着该值低于 IEEE-754 算术中可用的最小非规范化数字)产生 "correct" 结果。
写了这么多,正如其他人暗示的那样,我将 quote 原始 Matlab 开发人员和 Mathworks 联合创始人 Cleve Moler:
[The determinant] is useful in theoretical considerations and hand calculations, but does not provide a sound basis for robust numerical software.