Frage

I have a matrix in R that is supposed to be symmetric, however, due to machine precision the matrix is never symmetric (the values differ by around 10^-16). Since I know the matrix is symmetric I have been doing this so far to get around the problem:

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

Is there a better one line command for this?

War es hilfreich?

Lösung 3

Is the workaround really necessary if the values only differ by that much?

Someone pointed out that my previous answer was wrong. I like some of the other ones better, but since I can't delete this one (accepted by a user who left), here's yet another solution using the micEcon package:

symMatrix(s[upper.tri(s, TRUE)], nrow=nrow(s), byrow=TRUE)

Andere Tipps

s<-matrix(1:25,5)
s[lower.tri(s)] = t(s)[lower.tri(s)]

You can force the matrix to be symmetric using forceSymmetric function in Matrix package in R:

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
 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

I was curious to compare all the methods, so ran a quick microbenchmark. Clearly, the simplest 0.5 * (S + t(S)) is the fastest.

The specific function Matrix::forceSymmetric() is sometimes slightly faster, but it returns an object of a different class (dsyMatrix instead of matrix), and converting back to matrix takes a lot of time (although one might argue that it is a good idea to keep the output as dsyMatrix for further gains in computation).

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

Created on 2021-02-08 by the reprex package (v1.0.0)

as.dist() will overwrite the upper triangle of a matrix with the lower one and replace the diagonal with zeros. This method only works on numeric matrices.

mat <- matrix(1:25, 5)

unname(`diag<-`(as.matrix(as.dist(mat)), diag(mat)))

#      [,1] [,2] [,3] [,4] [,5]
# [1,]    1    2    3    4    5
# [2,]    2    7    8    9   10
# [3,]    3    8   13   14   15
# [4,]    4    9   14   19   20
# [5,]    5   10   15   20   25

Inspired by user3318600

    s<-matrix(1:25,5)
    s[lower.tri(s)]<-s[upper.tri(s)]
Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top