[Rd] eigen(symmetric=TRUE) for complex matrices
Berend Hasselman
bhh at xs4all.nl
Tue Jun 18 21:23:16 CEST 2013
On 18-06-2013, at 13:23, peter dalgaard <pdalgd at gmail.com> wrote:
> On some further digging, on the Lion machine, the problem seems absent from builds against the veclib framework. I strongly suspect that the issue is the same as PR#14964, which also affects Hermitian matrices of dimension 33x33 and above.
>
With the CRAN distribution of R-3.0.1 on OS X 10.8.4 I encountered the problem too.
I've done a small experiment to see if I could find a possible cause.
This version of R uses the libRblas and libRlapack.
For complex matrices I'm assuming that eigen uses the Lapack routine zheev and that R first calls zheev to determine the optimal workspace.
I've made an alternative eigen that also uses Lapack's zheev but calls it with minimal workspace to force zheev to use the unblocked algorithm. To determine if the blocked algorithm that Lapack uses could be the cause of the problem.
This will work when eigen is called first so that the Lapack library is dyn.loaded.
Code:
# quick and dirty
zeigen <- function(A,symmetric,only.values=FALSE) {
# SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
# * INFO )
# .. Scalar Arguments ..
# * CHARACTER JOBZ, UPLO
# * INTEGER INFO, LDA, LWORK, N
# * ..
# * .. Array Arguments ..
# * DOUBLE PRECISION RWORK( * ), W( * )
# * COMPLEX*16 A( LDA, * ), WORK( * )
n <- dim(A)[1]
# use minimal workspace
lwork <- as.integer(2*n-1)
# warning: "N" and "L" work on OSX but you may have to write an interface routine since some systems
# give Fortran runtime errors when passing character arguments to Fortran routines.
z <- .Fortran("zheev", "N", "L", n, A=A, n, w=numeric(n), work=complex(lwork), lwork,
numeric(3*n-2), info=integer(1L))
list(values=z$w, vectors=z$A)
}
N <- 100
jj <- matrix(0,N,N)
A <- exp(-0.1*(row(jj)-col(jj))^2)
min(eigen(A,only.values=TRUE,symmetric=TRUE)$values)
## [1] 2.521153e-10
# Coercing A to a complex matrix should make no difference, but makes
# eigen() return the wrong answer:
min(eigen(A+0i,only.values=TRUE,symmetric=TRUE)$values)
## [1] -0.359347
min(zeigen(A+0i,only.values=TRUE,symmetric=TRUE)$values)
##[1] 2.521154e-10
So it seems that the blocked algorithm is the cause of the error and that using the (possibly slow) unblocked algorithm gives the correct result.
Berend
More information about the R-devel
mailing list