R、C和Python中伪逆的不同结果

Question

我在把一个算法从R转换成C，我需要得到矩阵的伪逆，但是C中得到的结果与R中得到的结果有些不同，这些差异改变了算法的行为。在

我用来在C中得到伪逆的代码是this。在

我读了一些书，得到伪逆的方法有很多种，C语言中使用的方法是Moore Penrose。R中使用的函数来自库corpcor。两者都使用“奇异值分解”。在

这是我想从中得到伪逆的矩阵

1                  0.920980394593472  0.996160973582776   0.996772980609752   0.997372221594439 0.999972797627027
0.920980394593472  1                  0.885601439824631   0.88878682654952    0.892173764646865 0.923738536637407
0.996160973582776  0.885601439824631  1                   0.999973383442349   0.999885329646229 0.99549326808266
0.996772980609752  0.88878682654952   0.999973383442349   1                   0.999969202115456 0.996158288591094
0.997372221594439  0.892173764646865  0.999885329646229   0.999969202115456   1                 0.996814694067663
0.999972797627027  0.923738536637407  0.99549326808266    0.996158288591094   0.996814694067663 1

我从R中的函数pseudo inverse（）得到的结果是：

^{pr2}$

我在C中得到的结果是：

1398795243.74255  79184.33844201    -9594022229.12525 28322858223.2099  -19819644215.1338 -305583186.690388
79166.91917247    3402.48426033     52556628.829717   -126546466.939768 76466567.769084   -2564764.38775363
-9594334089.78616 52556515.9039231  1004461808180.58  -2455360323666.24 1502806633291.96  -42481639977.8112
28323609294.95    -126546129.049526 -2455359143404.21 6020330778543.35  -3694093433789.59  101211765648.895
-19820098170.0141 76466329.4304944  1502805309171.23  -3694091962863.6   2271455511686.72  -60603547743.7687
-305568392.855205 -2564768.40243798 -42481807759.1065 101212225714.588   -60603854784.616  2185698311.36118

两者之间的区别是

-118562.671649933 415.6198192764    8247876.90765953  -20311027.5418015 12507618.5904007 -327309.978267968
424.5539327194    -1.3609798784     -27269.4164030999 67275.2745009959  -41490.291238904  1064.5449163299
8413314.25839043  -27227.552922301  -544970420.589966 1344206932.90039  -828869777.349854 21313146.4894028
-20709242.8262024 67140.0062440038  1343568061.89014  -3313879228.39941 2043353828.96973  -52563162.2870026
12748195.2462006  -41390.7198514938 -828153259.419922 2042558174.66016  -1259419586.19043 32408074.3294983
-335134.889266014 1067.29834637     21400799.0577011  -52804363.5690002 32569427.5587997  -834391.050109863

为了检查我在C语言中使用的算法是否有问题，我在python中用纽比.利纳格.皮涅夫（）使用“奇异值分解”。结果与C和R不同

1398224882.37767  81521.32618159    -9548319116.82994 28210636794.0452  -19750702778.4149 -307443670.558374
81576.67749763    3392.80756354     52367028.3401356  -126080750.377468 76180379.3995419  -2557069.77374461
-9547349936.09641 52367486.8455529  1000758728845.37  -2446264734953.02 1497217439225.67  -42331313003.6236
28208301799.8629  -126082060.163116 -2446268326785.52 5998001838415.43  -3680372478514.1  100842703532.378
-19749291055.22   76181277.4796568  1497221470187.79  -3680376958173.79 2263027785174.03  -60376849475.2803
-307489737.200422 -2557061.32729561 -42330783514.2789 100841257137.344  -60375886615.3659 2179570267.21681

• 如果使用的方法和数据是相同的，导致结果的原因是什么 与众不同？。在
• 哪个结果最准确？在

编辑我犯了一个错误，我没有把所有的数字都放在矩阵中重新创建结果，我已经用正确的矩阵更新了问题。在

Answer 1

Ageneralized inverseA^g对于A应该满足

A^gAA^g=A^g

AA^gA=A

（A^gA）^T=Ag

对于给定的矩阵，corpcor::pseudoinverse的结果不满足这些属性，而{}的结果满足：

check_pinv <- function(mat, fun, ...) {
    pinv <- fun(mat, ...)
    isTRUE(all.equal(mat %*% pinv %*% mat, mat)) &&
        isTRUE(all.equal(pinv %*% mat %*% pinv, pinv)) &&
        isTRUE(all.equal(pinv %*% mat, t(mat %*% pinv))) &&
        isTRUE(all.equal(mat %*% pinv, t(pinv %*% mat)))
}

mat <- matrix(c(                                                       
   1,                  0.920980394593472,  0.996160973582776,   0.996772980609752,   0.997372221594439, 0.999972797627027,
   0.920980394593472,  1,                  0.885601439824631,   0.88878682654952,    0.892173764646865, 0.923738536637407,
   0.996160973582776,  0.885601439824631,  1,                   0.999973383442349,   0.999885329646229, 0.99549326808266,
   0.996772980609752,  0.88878682654952,   0.999973383442349,   1,                   0.999969202115456, 0.996158288591094,
   0.997372221594439,  0.892173764646865,  0.999885329646229,   0.999969202115456,   1,                 0.996814694067663,
   0.999972797627027,  0.923738536637407,  0.99549326808266,    0.996158288591094,   0.996814694067663, 1), nrow = 6, ncol = 6)

check_pinv(mat, corpcor::pseudoinverse)
#> [1] FALSE
check_pinv(mat, MASS::ginv)
#> [1] TRUE

这两个函数之间的一个重要区别是用于确定是否应将奇异值视为零的默认公差级别。如果使用MASS::ginv（即sqrt(.Machine$double.eps)）中使用的值也用于corpcor::pseudoinverse，则满足伪逆属性：

^{pr2}$

注意，必须使用max(svd(mat)$d) * sqrt(.Machine$double.eps)，因为corpcor::pseudoinverse从绝对意义上解释公差，而{}将公差视为相对于最大单数值。在这个公差水平下，产生的伪逆矩阵是相同的。在

all.equal(corpcor::pseudoinverse(mat, max(svd(mat)$d) * sqrt(.Machine$double.eps)), 
          MASS::ginv(mat))
#> [1] TRUE

在python中，numpy.linalg.pinv和{}都不满足以下属性：

import numpy
mat = numpy.array([[1,            0.9209803946, 0.9961609736, 0.9967729806, 0.9973722216, 0.9999727976],
                   [0.9209803946, 1,            0.8856014398, 0.8887868265, 0.8921737646, 0.9237385366],
                   [0.9961609736, 0.8856014398, 1,            0.9999733834, 0.9998853296, 0.9954932681],
                   [0.9967729806, 0.8887868265, 0.9999733834, 1,            0.9999692021, 0.9961582886],
                   [0.9973722216, 0.8921737646, 0.9998853296, 0.9999692021, 1,            0.9968146941],
                   [0.9999727976, 0.9237385366, 0.9954932681, 0.9961582886, 0.9968146941, 1]])

pinv1 = numpy.linalg.pinv(mat)
print numpy.allclose(pinv1.dot(mat).dot(pinv1), pinv1)
# False
print numpy.allclose(mat.dot(pinv1).dot(mat), mat)
# True

from scipy import linalg
pinv2 = linalg.pinv(mat)
print numpy.allclose(pinv2.dot(mat).dot(pinv2), pinv2)
# False
print numpy.allclose(mat.dot(pinv2).dot(mat), mat)
# False

print numpy.allclose(pinv1, pinv2)
# True

注意：矩阵使用原始值。结果不受影响，因为只有最小的奇异值显示出显著的变化。在

同样，如果使用1e-8而不是默认的1e-15作为公差，则满足这些伪逆特性。C版本也是如此，可以从R和RcppGSL一起使用。在