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?
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*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*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)]
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) #  FALSE ## solution ndx <- lower.tri(A) A[ndx] <- t(A)[ndx] isSymmetric(A) #  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*ve[,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)