mgcv (version 1.8-25)

slanczos: Compute truncated eigen decomposition of a symmetric matrix

Description

Uses Lanczos iteration to find the truncated eigen-decomposition of a symmetric matrix.

Usage

slanczos(A,k=10,kl=-1,tol=.Machine$double.eps^.5,nt=1)

Arguments

A

A symmetric matrix.

k

Must be non-negative. If kl is negative, then the k largest magnitude eigenvalues are found, together with the corresponding eigenvectors. If kl is non-negative then the k highest eigenvalues are found together with their eigenvectors and the kl lowest eigenvalues with eigenvectors are also returned.

kl

If kl is non-negative then the kl lowest eigenvalues are returned together with their corresponding eigenvectors (in addition to the k highest eignevalues + vectors). negative kl signals that the k largest magnitude eigenvalues should be returned, with eigenvectors.

tol

tolerance to use for convergence testing of eigenvalues. Error in eigenvalues will be less than the magnitude of the dominant eigenvalue multiplied by tol (or the machine precision!).

nt

number of threads to use for leading order iterative multiplication of A by vector. May show no speed improvement on two processor machine.

Value

A list with elements values (array of eigenvalues); vectors (matrix with eigenvectors in its columns); iter (number of iterations required).

Details

If kl is non-negative, returns the highest k and lowest kl eigenvalues, with their corresponding eigenvectors. If kl is negative, returns the largest magnitude k eigenvalues, with corresponding eigenvectors.

The routine implements Lanczos iteration with full re-orthogonalization as described in Demmel (1997). Lanczos iteraction iteratively constructs a tridiagonal matrix, the eigenvalues of which converge to the eigenvalues of A, as the iteration proceeds (most extreme first). Eigenvectors can also be computed. For small k and kl the approach is faster than computing the full symmetric eigendecompostion. The tridiagonal eigenproblems are handled using LAPACK.

The implementation is not optimal: in particular the inner triadiagonal problems could be handled more efficiently, and there would be some savings to be made by not always returning eigenvectors.

References

Demmel, J. (1997) Applied Numerical Linear Algebra. SIAM

See Also

cyclic.p.spline

Examples

# NOT RUN {
 require(mgcv)
 ## create some x's and knots...
 set.seed(1);
 n <- 700;A <- matrix(runif(n*n),n,n);A <- A+t(A)
 
 ## compare timings of slanczos and eigen
 system.time(er <- slanczos(A,10))
 system.time(um <- eigen(A,symmetric=TRUE))
 
 ## confirm values are the same...
 ind <- c(1:6,(n-3):n)
 range(er$values-um$values[ind]);range(abs(er$vectors)-abs(um$vectors[,ind]))
# }