R basic function eigen has problem with symmetric matrices

user_7614251 0
Viewed 200

eigen is a base R function that returns eigenvectors and eigenvalues for a given matrix. I just found that it can be not robust for the symmetric matrices when you specify this explicitly using symmetric=TRUE. I show the issue below. Are there libraries that fix the issue and find correct eigenvectors using information about symmetry?

Example

Let us firstly find an eigen vector for a symmetric matrix without specifying the symmetry.

set.seed(1)
I = 10; A = matrix( rnorm(I*I), ncol = I );
## make the matrix symmetric
A[lower.tri(A)] = A[upper.tri(A)]
ev = eigen(A, symmetric=F); ve = ev$vectors; va = ev$values;

Test the result taking the difference between A applied to its eigenvector and eigenvector scaled by the eigenvalue. That must be a zero vector. Hence, the sum of its elements must be zero:

sum(A %*% ve[,1] - va[1]*ve[,1] )

-3.663736e-15+0i

Now, let's do the same specifying symmetry:

ev = eigen(A, symmetric=TRUE); ve = ev$vectors; va = ev$values;
sum(A %*% ve[,1] - va[1]*ve[,1] )

-0.3534416

As it was mentioned in the answers, my example was bad since the matrix was not symmetric. The approach I used for the example

A[lower.tri(A)] = A[upper.tri(A)]

will actually give an asymmetric matrix. The more accurate way of getting a symmetric matrix would be

A[lower.tri(A)] = t(A)[lower.tri(A)]

or

A <- A+t(A)
eigenvalue eigenvector r
2 Answers

The example matrix was not made symmetric. Try the following:

set.seed(1)
I = 10; A = matrix(rnorm(I*I), ncol = I)

## original posters attempt
A1 <- A
A1[lower.tri(A1)] = A1[upper.tri(A1)]
isSymmetric(A1)
# [1] FALSE

## solution
ndx <- lower.tri(A)
A[ndx] <- t(A)[ndx]
isSymmetric(A)
# [1] TRUE

#ev = eigen(A, symmetric=FALSE); ve = ev$vectors; va = ev$values

ev = eigen(A, symmetric=TRUE); ve = ev$vectors; va = ev$values
sum(A %*% ve[,1] - va[1]*ve[,1])
# [1] -2.63678e-15
user_36775760

If we prend our matrix wasn't symmetric, R seems to use another algorithm resulting in a different output order but the eigen values are the same.

#?eigen


# random matrix
set.seed(1)
n = 3; 
A = matrix(rnorm(n*n), ncol = n, nrow = n)
A[lower.tri(A)] = t(A)[lower.tri(A)]

cat("symmetric matrix A\n")
A


cat("\nsymmetric = FALSE\n")
eig1 = eigen(A, symmetric = FALSE)
eig1$vectors
eig1$values
prod(eig1$values)

cat("\nsymmetric = TRUE\n")
eig2 = eigen(A, symmetric = TRUE)
eig2$vectors
eig2$values
prod(eig2$values)

cat("\nsymmetric = ...\n")
eig3 = eigen(A)
eig3$vectors
eig3$values
prod(eig3$values)
user_105658550
Related