Solution for R basic function eigen has problem with symmetric matrices [closed]
is Given Below:
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)
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
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 An")
A
cat("nsymmetric = FALSEn")
eig1 = eigen(A, symmetric = FALSE)
eig1$vectors
eig1$values
prod(eig1$values)
cat("nsymmetric = TRUEn")
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)