Matrix (version 0.3-24)

eigen: Spectral Decomposition of a Matrix

Description

eigen is a generic function. The default method eigen.default computes eigenvalues and eigenvectors by providing an interface to the EISPACK routines RS, RG, CH and CG. eigen.Matrix provides an interface to the Lapack functions DSYEV and DGEEV.

Usage

## S3 method for class 'Matrix':
eigen(x, vectors = TRUE, balance = FALSE, 
   rcond = FALSE, \dots)

Arguments

x

a matrix whose spectral decomposition is to be computed.

vectors

if TRUE, both the eigenvalues and eigenvectors are returned, otherwise only the eigenvalues are returned.

balance

logical. Not used at present.

rcond

logical. Not used at present.

...

further arguments passed to or from other methods.

Value

• The spectral decomposition of x is returned as components of a list.
• values
  a vector containing the $p$ eigenvalues of x, sorted in decreasing order, according to Mod(values) if they are complex.
• vectors
  a $p\times p$ matrix whose columns contain the eigenvectors of x, or NULL if only.values is TRUE.

References

Smith, B. T, Boyle, J. M., Dongarra, J. J., Garbow, B. S., Ikebe,Y., Klema, V., and Moler, C. B. (1976). Matrix Eigensystems Routines -- EISPACK Guide. Springer-Verlag Lecture Notes in Computer Science.

Anderson, E., et al. (1999). LAPACK User's Guide, 3rd edition, SIAM, Philadelphia.

See Also

eigen, svd, a generalization of eigen; qr, and chol for related decompositions.

To compute the determinant of a matrix, the qr decomposition is much more efficient: det.

Examples

eigen(cbind(c(1,-1),c(-1,1)))   # uses Eispack
eigen(cbind(c(1,-1),c(-1,1)), symmetric = FALSE)# Eispack (different algorithm).
eigen(as.Matrix(cbind(c(1,-1),c(-1,1))))  # uses Lapack

eigen(cbind(1,c(1,-1)), only.values = TRUE)
eigen(as.Matrix(cbind(1,c(1,-1))), vectors = FALSE)
eigen(cbind(-1,2:1)) # complex values
eigen(as.Matrix(cbind(-1,2:1))) # complex values from Lapack
eigen(print(cbind(c(0,1i), c(-1i,0))))# Hermite ==> real Eigen values
## 3 x 3:
eigen(cbind( 1,3:1,1:3))
eigen(as.Matrix(cbind( 1,3:1,1:3)))
eigen(cbind(-1,c(1:2,0),0:2)) # complex values
eigen(as.Matrix(cbind(-1,c(1:2,0),0:2))) # complex values

Meps <- .Machine$double.eps
m <- matrix(round(rnorm(25),3), 5,5)
sm <- m + t(m) #- symmetric matrix
eigen(as.Matrix(sm), vectors = FALSE) # ordered INcreasingly
em <- eigen(sm); V <- em$vect
print(lam <- em$values) # ordered DEcreasingly


stopifnot(
 abs(sm %*% V - V %*% diag(lam))          < 60*Meps,
 abs(sm       - V %*% diag(lam) %*% t(V)) < 60*Meps)

##------- Symmetric = FALSE:  -- different to above : ---

em <- eigen(sm, symmetric = FALSE); V2 <- em$vect
print(lam2 <- em$values) # ordered decreasingly in ABSolute value !
                         # and V2 is not normalized (where V is):
print(i <- rev(order(lam2)))
stopifnot(abs(1 - lam2[i] / lam) < 60 * Meps)

zapsmall(Diag <- t(V2) %*% V2) # orthogonal, but not normalized
print(norm2V <- apply(V2 * V2, 2, sum))
stopifnot( abs(1- norm2V / diag(Diag)) < 60*Meps)

V2n <- sweep(V2,2, STATS= sqrt(norm2V), FUN="/")## V2n are now Normalized EV
apply(V2n * V2n, 2, sum)
##[1] 1 1 1 1 1

## Both are now TRUE:
stopifnot(abs(sm %*% V2n - V2n %*% diag(lam2))            < 60*Meps,
          abs(sm         - V2n %*% diag(lam2) %*% t(V2n)) < 60*Meps)

## Re-ordered as with symmetric:
sV <- V2n[,i]
slam <- lam2[i]
all(abs(sm %*% sV -  sV %*% diag(slam))             < 60*Meps)
all(abs(sm        -  sV %*% diag(slam) %*% t(sV)) < 60*Meps)
## sV  *is* now equal to V  -- up to sign (+-) and rounding errors
all(abs(c(1 - abs(sV / V)))       <     1000*Meps) # TRUE (P ~ 0.95)