如何移植np. linear .solve到R?

Question

我正试图将一些代码从Python移植到R，我在numpy.linalg.solve计算方面遇到了一点麻烦。

这是我的Python计算

import numpy as np
A = np.array([[1.21320066609265e-09, 7.41796679861689e-10, 3.88326978297399e-09, -8.0434635935612e-12],
[7.41796679861689e-10, 7.42222200630816e-09, -1.93295794328878e-09, -2.45873626472721e-10],
[3.88326978297399e-09, -1.93295794328878e-09, 1.52225096683626e-08, 9.93308495566108e-11],
[-8.0434635935612e-12, -2.45873626472721e-10, 9.93308495566108e-11, 1.27503537782411e-11]])
B = np.array([[4.87908846480423e-06, 5.95145387478624e-05, -2.52777639389024e-05, -8.97778789976655e-07],
[3.78363029085449e-06, 6.06099965807516e-05, -2.41823065607605e-05, -1.99323696392639e-06],
[3.78363029085449e-06, 6.06099965807516e-05, -2.41823065607605e-05, -1.99323696392639e-06],
[2.29254385430977e-05, 3.54386564822562e-05, 6.04902019176734e-05, -9.56686605491086e-07],
[2.29254385430977e-05, 3.54386564822562e-05, 6.04902019176734e-05, -9.56686605491086e-07]]).T
np.linalg.solve(A, B).T
Out[501]: 
array([[  3225.26,  15022.19,  -2128.05, 237883.89],
[  4483.82,   7277.74,  -1814.67,    979.25],
[  4483.82,   7277.74,  -1814.67,    979.25],
[ 20073.38,   1939.59,   -777.31, -18911.1 ],
[ 20073.38,   1939.59,   -777.31, -18911.1 ]])

在R中，我创建了等价矩阵:

A <- as.matrix(structure(list(V1 = c(1.21320066609265e-09, 7.41796679861689e-10,
3.88326978297399e-09, -8.0434635935612e-12),
V2 = c(7.41796679861689e-10, 7.42222200630816e-09,
-1.93295794328878e-09, -2.45873626472721e-10), 
V3 = c(3.88326978297399e-09, -1.93295794328878e-09,
1.52225096683626e-08, 9.93308495566108e-11),
V4 = c(-8.0434635935612e-12, -2.45873626472721e-10, 
9.93308495566108e-11, 1.27503537782411e-11)),
class = "data.frame", row.names = c(NA, -4L)))
B <- as.matrix(structure(list(V1 = c(4.87908846480423e-06, 5.95145387478624e-05, 
-2.52777639389024e-05, -8.97778789976655e-07),
V2 = c(3.78363029085449e-06, 6.06099965807516e-05,
-2.41823065607605e-05, -1.99323696392639e-06),
V3 = c(3.78363029085449e-06, 6.06099965807516e-05,
-2.41823065607605e-05, -1.99323696392639e-06),
V4 = c(2.29254385430977e-05, 3.54386564822562e-05,
6.04902019176734e-05, -9.56686605491086e-07),
V5 = c(2.29254385430977e-05, 3.54386564822562e-05,
6.04902019176734e-05, -9.56686605491086e-07)),
class = "data.frame", row.names = c(NA, -4L)))   

从阅读这篇文章numpy. linalgl .solve()和R solve()的区别我想我应该试试backsolve:

> t(backsolve(A, B, transpose = TRUE))
[,1]     [,2]       [,3]      [,4]
[1,]  4021.666 7616.490 -1719.3371 92393.246
[2,]  3118.718 7854.325 -1386.8310  7903.683
[3,]  3118.718 7854.325 -1386.8310  7903.683
[4,] 18896.658 2886.087  -480.3371 -3714.961
[5,] 18896.658 2886.087  -480.3371 -3714.961

…但这显然是非常不同的答案。

编辑:首先得到A的cholesky分解也不起作用:

> t(backsolve(chol(A), B, transpose = FALSE))
[,1]      [,2]      [,3]     [,4]
[1,]  842372.6 -185848.5 -230277.0 -1258487
[2,] 1870215.2 -412617.6 -511255.2 -2794078
[3,] 1870215.2 -412617.6 -511255.2 -2794078
[4,]  897618.3 -198038.6 -245379.0 -1341063
[5,]  897618.3 -198038.6 -245379.0 -1341063

有人能给我指一下正确的方向吗?

Answer 1

简短的回答:你不需要移植这个函数，已经有了R求解函数，以及许多其他方法，例如使用Cholesky分解(对称矩阵)，QR分解或矩阵包中的函数，它们可以处理稀疏矩阵。

---

正如user20650在下面的评论中指出的那样，你在R和Python中没有相同的数据。当A(3,4)上的错字被纠正时，R检测到矩阵A是奇异的，拒绝给出结果。

下面是R代码:

A <- matrix(c(1.21320066609265e-09, 7.41796679861689e-10, 3.88326978297399e-09, -8.0434635935612e-12,
7.41796679861689e-10, 7.42222200630816e-09, -1.93295794328878e-09, -2.45873626472721e-10,
3.88326978297399e-09, -1.93295794328878e-09, 1.52225096683626e-08, 9.93308495566108e-11,
-8.0434635935612e-12, -2.45873626472721e-10, 9.93308495566108e-11, 1.27503537782411e-11),
4, 4, byrow = T)
B <- matrix(c(4.87908846480423e-06, 5.95145387478624e-05, -2.52777639389024e-05, -8.97778789976655e-07,
3.78363029085449e-06, 6.06099965807516e-05, -2.41823065607605e-05, -1.99323696392639e-06,
3.78363029085449e-06, 6.06099965807516e-05, -2.41823065607605e-05, -1.99323696392639e-06,
2.29254385430977e-05, 3.54386564822562e-05, 6.04902019176734e-05, -9.56686605491086e-07,
2.29254385430977e-05, 3.54386564822562e-05, 6.04902019176734e-05, -9.56686605491086e-07),
4, 5, byrow = F)
> solve(A, B)
Error in solve.default(A, B) : 
system is computationally singular: reciprocal condition number = 1.46364e-17
Calls: solve -> solve.default

但是，请注意，Python为条件号: 生成此值。

>>> np.linalg.cond(A)
5.955813150735809e+17

即系统是非常病态的，即矩阵是近奇异的，并且结果对数值误差高度敏感。即使Python和R都使用LAPACK，库编译方式的微小差异也可能解释为什么它们不能产生相同的结果，甚至一个失败而另一个没有。特别是因为numpy使用的LAPACK显然是由f2py转换的。我用intel-numpy也没有得到相同的结果。您可能想知道为什么使用MKL的Intel Fortran不能产生与使用Intel -numpy的Intel Python相同的结果。的确，这有点令人惊讶。在Python/numpy上的第三个版本(安装了MSYS2)，我得到了另一个结果。

使用Cholesky分解，R产生一个答案(这个答案不再有任何价值，但无论如何，它是一个答案):

s <- chol(A)
backsolve(s, forwardsolve(t(s), B))
# or: backsolve(s, backsolve(s, B, transpose = T))
[,1]      [,2]      [,3]        [,4]        [,5]
[1,]   8432.584  7593.700  7593.700  25745.1517  25745.1517
[2,]  13873.313  6591.618  6591.618    688.2458    688.2458
[3,]  -3551.561 -2664.813 -2664.813  -2327.7813  -2327.7813
[4,] 230104.178 -3666.895 -3666.895 -27384.6857 -27384.6857

---

一个直接调用LAPACK的快速而肮脏的Fortran程序(Intel Fortran &MKL)。对DGESV的调用以错误3结束，即第三个主元完全为零，即矩阵A完全是奇异的，并且尝试解决这个系统真的不是一个好主意，无论你试图"移植"什么。

program test
implicit none
real(8) :: a(4, 4), b(4, 5)
integer :: ipiv(4), info, i

a(1,:) = [1.21320066609266d-09, 7.41796679861689d-10, 3.88326978297399d-09, -8.0434635935612d-12]
a(2,:) = [7.41796679861689d-10, 7.42222200630816d-09, -1.93295794328878d-09, -2.45873626472721d-10]
a(2,:) = [3.88326978297399d-09, -1.93295794328878d-09, 1.52225096683626d-08, 9.93308495566108d-11]
a(2,:) = [-8.0434635935612d-12, -2.45873626472721d-10, 9.93308495566108d-11, 1.27503537782411d-11]
b = reshape([4.87908846480423d-06, 5.95145387478624d-05, -2.52777639389024d-05, -8.97778789976655d-07, &
3.78363029085449d-06, 6.06099965807516d-05, -2.41823065607605d-05, -1.99323696392639d-06, &
3.78363029085449d-06, 6.06099965807516d-05, -2.41823065607605d-05, -1.99323696392639d-06, &
2.29254385430977d-05, 3.54386564822562d-05, 6.04902019176734d-05, -9.56686605491086d-07, &
2.29254385430977d-05, 3.54386564822562d-05, 6.04902019176734d-05, -9.56686605491086d-07], &
[4, 5])
call dgesv(4, 5, a, 4, ipiv, b, 4, info)
print *, info
do i = 1, 4
print *, b(i,:)
end do
end program

---

让我们看一个更好的例子:

import numpy as np
import numpy.linalg
a = np.matrix([[2, 9, 4], [7, 5, 3], [6, 1, 8]], dtype=float)
b = np.eye(3)
np.linalg.solve(a, b)
array([[-0.10277778,  0.18888889, -0.01944444],
[ 0.10555556,  0.02222222, -0.06111111],
[ 0.06388889, -0.14444444,  0.14722222]])

While in R:

a <- matrix(c(2, 9, 4, 7, 5, 3, 6, 1, 8), 3, 3, byrow = T)
b <- diag(3)
solve(a, b)
[,1]        [,2]        [,3]
[1,] -0.10277778  0.18888889 -0.01944444
[2,]  0.10555556  0.02222222 -0.06111111
[3,]  0.06388889 -0.14444444  0.14722222