矩阵计算 R 与 Stata 的精度答案

【问题标题】：Precision of matrix calculation R vs. Stata矩阵计算 R 与 Stata 的精度
【发布时间】：2018-02-01 07:41:18
【问题描述】：

我正在尝试确定矩阵是否为负半定义。出于这个原因，我检查是否所有特征值都小于或等于零。一个示例矩阵是：

              [,1]          [,2]          [,3]          [,4]
[1,] -1.181830e-05  0.0001576663 -2.602332e-07  1.472770e-05
[2,]  1.576663e-04 -0.0116220027  3.249607e-04 -2.348050e-04
[3,] -2.602332e-07  0.0003249607 -2.616447e-05  3.492998e-05
[4,]  1.472770e-05 -0.0002348050  3.492998e-05 -9.103073e-05

stata计算的特征值为1.045e-12，-0.00001559，-0.00009737，-0.01163805。但是，R计算的特征值是-1.207746e-20、-1.558760e-05、-9.737074e-05、-1.163806e-02。所以最后三个特征值非常相似，但第一个非常接近零的特征值不是。使用 stata 获得的特征值，矩阵不是半定义的，但使用 R 获得的特征值是半定义的。有没有办法找出哪个计算更精确？或者甚至可以重新调整矩阵以避免无限小的特征值？

提前非常感谢您。每一个提示都将受到高度赞赏。

【问题讨论】：

我使用了 options(scipen = 999) 并得到了第一个特征值 (0.000000000003696385) 的不同结果。 R 使用有限精度算术，因此答案在小数点后 15 或 16 位后不准确。
我在 R 中得到了不同的结果。您能否在矩阵上给出dput 的结果，并说明您使用的是什么版本的 R 以及计算特征值的函数？
感谢你们！ James，输入是structure(c(-1.18182990902922e-05, 0.000157666259977248, -2.60233161187324e-07, 1.47276984589202e-05, 0.000157666259977248, -0.0116220026732496, 0.000324960673501042, -0.000234805002079655, -2.60233161187324e-07, 0.000324960673501042, -2.61644706608917e-05, 3.49299778208411e-05, 1.47276984589202e-05, -0.000234805002079655, 3.49299778208411e-05, -9.10307332256831e-05), .Dim = c(4L, 4L))我使用的是RStudio 1.0.44版，计算特征值的函数是eigen()。

标签： r matrix stata zero eigenvalue

【解决方案1】：

您不能期望使用双精度浮点数的数值算法具有如此高的精度。

您可以期望不超过 17 个十进制数字，并且零附近的相对精度损失并不少见。也就是说，给定数值误差，1e-12 和 -1e-20 几乎都与 0 无法区分。

例如，对于最小的特征值（使用您在评论中给出的系数），我得到：

R 3.4.1：5.929231e-21，
MATLAB R2017a：3.412972022812169e-19
Stata 15：3.2998e-20 (matrix eigenvalues) 或 4.464e-19 (matrix symeigen)
带有 MKL 的英特尔 Fortran（DSYEV 函数）：2.2608e-19

您可以选择一个阈值，例如 1e-10，并在特征值与最大特征值的比率小于 1e-10 时将其强制为零。

不管怎样，你的 1e-12 看起来有点大。在 Stata 和 R 之间传输数据时，您可能会丢失一些精度：矩阵中的小相对误差可能会导致特征值为零附近的大相对误差。使用 Stata 和您问题中的数据（不在评论中），我得到例如 3.696e-12 的最小特征值。

但是，即使使用相同的矩阵，由于以下方面的差异，可能仍然存在差异（以上有）：

解析器，如果您以文本形式输入数字
用于特征值计算的算法
同一算法的实现细节（例如，浮点运算符不具有关联性）
用于编译计算例程或编译器选项的编译器
浮点硬件

这类问题的传统建议阅读：

What Every Computer Scientist Should Know About Floating-Point Arithmetic

【讨论】：

python (np.linalg.eigvals) 给了我5.96768236e-20 以增加示例数