[R] Difference between eigs() and eigen()

[Pierrick Bruneau]

Dear R users and contributors,

I recently observed a difference between the outputs of the classic
eigen() function, and the Arnoldi variant eigs() that extracts only
the few first eigenpairs. Here is some sample code illustrating the
problem:

library(rARPACK)
library(speccalt)
set.seed(1)

# compute kernel matrix from rows of synth5
# then its Laplacian
kern <- local.rbfdot(synth5)
diag(kern) <- 0
deg <- sapply(1:(dim(synth5)[1]), function(i) {
return(sum(kern[i,]))
})
L <- diag(1/sqrt(deg)) %*% kern %*% diag(1/sqrt(deg))

eig1 <- eigs(L, 6)
eig2 <- eigen(L, symmetric=TRUE)

eig1$values then reads:
1.0000000 1.0000000 0.9993805 0.9992561 0.9985084 0.9975311

whereas eig2$values reads:
1.0000000 1.0000000 1.0000000 1.0000000 0.9993805 0.9992561
which is the correct result (eigenvalue 1 has multiplicity 4 in that example).

I guess there is an issue between Arnoldi methods and eigenvalues with
multiplicities greater than 1 (indeed at the end of the series the
unique eigenvals look identical), but as the problem is not documented
in the package PDF, I'm quite unclear if this is
implementation-specific or Arnoldi-general... The issue is quite
important in my case, as the associated eigenvectors then differ quite
significantly, and this impacts negatively my further operations.

I guess the next step is to dig into the mathematical literature, but
before this I wondered if someone already encountered this issue?

Any help would be appreciated,
Pierrick