[Rd] eigen(symmetric=TRUE) for complex matrices

[Berend Hasselman]

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