36

我在 R 中有一个应该是对称的矩阵，但是，由于机器精度，该矩阵永远不是对称的（值相差大约 10^-16）。因为我知道矩阵是对称的，所以到目前为止我一直在这样做以解决这个问题：

s.diag = diag(s)
s[lower.tri(s,diag=T)] = 0
s = s + t(s) + diag(s.diag,S)

有没有更好的单行命令呢？

4

6 回答 6

69 
s<-matrix(1:25,5)
s[lower.tri(s)] = t(s)[lower.tri(s)]
于 2014-02-17T10:03:07.030 回答
15

您可以使用R 中的包中的forceSymmetric函数强制矩阵对称：Matrix

library(Matrix)
x<-Matrix(rnorm(9), 3)
> x
3 x 3 Matrix of class "dgeMatrix"
           [,1]       [,2]       [,3]
[1,] -1.3484514 -0.4460452 -0.2828216
[2,]  0.7076883 -1.0411563  0.4324291
[3,] -0.4108909 -0.3292247 -0.3076071

A <- forceSymmetric(x)
> A
3 x 3 Matrix of class "dsyMatrix"
           [,1]       [,2]       [,3]
[1,] -1.3484514 -0.4460452 -0.2828216
[2,] -0.4460452 -1.0411563  0.4324291
[3,] -0.2828216  0.4324291 -0.3076071
于 2013-08-10T19:20:59.517 回答
10

如果值仅相差这么多，解决方法真的有必要吗？

有人指出我之前的回答是错误的。我更喜欢其他一些，但由于我无法删除这个（已离开的用户接受），这是使用该micEcon包的另一种解决方案：

symMatrix(s[upper.tri(s, TRUE)], nrow=nrow(s), byrow=TRUE)
于 2013-08-10T19:07:30.723 回答
8 
 s<-matrix(1:25,5)
 pmean <- function(x,y) (x+y)/2
 s[] <- pmean(s, matrix(s, nrow(s), byrow=TRUE))
 s
#-------
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    4    7   10   13
[2,]    4    7   10   13   16
[3,]    7   10   13   16   19
[4,]   10   13   16   19   22
[5,]   13   16   19   22   25
于 2013-08-10T19:08:08.437 回答
2

我很想比较所有的方法，所以跑了一个快速的microbenchmark. 显然，最简单0.5 * (S + t(S))的就是最快的。

特定的函数Matrix::forceSymmetric()有时会稍微快一些，但它返回一个不同类的对象（dsyMatrix而不是matrix），并且转换回 tomatrix需要很多时间（尽管有人可能会争辩说保留输出以dsyMatrix获得进一步的收益是一个好主意在计算中）。

S <-matrix(1:50^2,50)
pick_lower <- function(M) M[lower.tri(M)] = t(M)[lower.tri(M)]

microbenchmark::microbenchmark(micEcon=miscTools::symMatrix(S[upper.tri(S, TRUE)], nrow=nrow(S), byrow=TRUE),
                               Matri_raw =Matrix::forceSymmetric(S),
                               Matri_conv =as.matrix(Matrix::forceSymmetric(S)),
                               pick_lower = pick_lower(S),
                               base =0.5 * (S + t(S)),
                               times=100) 
#> Unit: microseconds
#>        expr    min      lq       mean   median       uq        max neval cld
#>     micEcon 62.133 74.7515  136.49538 104.2430 115.6950   3581.001   100   a
#>   Matri_raw 14.766 17.9130   24.15157  24.5060  26.6050     63.939   100   a
#>  Matri_conv 46.767 59.8165 5621.96140  66.3785  73.5380 555393.346   100   a
#>  pick_lower 27.907 30.7930  235.65058  48.9760  53.0425  12484.779   100   a
#>        base 10.771 12.4535   16.97627  17.1190  18.3175     47.623   100   a

由reprex 包于 2021-02-08 创建(v1.0.0)

于 2021-02-08T22:53:49.613 回答
-2

灵感来自 user3318600

    s<-matrix(1:25,5)
    s[lower.tri(s)]<-s[upper.tri(s)]
于 2020-11-11T22:05:16.830 回答