1 2## testing code is currently in 3## /u/maechler/R/MM/MISC/posdefify.R 4 5## TODO: probaby add 'rescale.kind = c("diag", "trace")' 6## ---- they would differ only when some EV's were negative 7 8### Higham's code by Ravi 9 10posdefify <- function(m, method = c("someEVadd", "allEVadd"), 11 symmetric = TRUE, 12 eigen.m = eigen(m, symmetric= symmetric), 13 eps.ev = 1e-7) 14{ 15 ## Purpose: From a matrix m, make a "close" positive definite one 16 ## ------------------------------------------------------------------------- 17 ## Arguments: m: numeric matrix (n x n), usually symmetric 18 ## ------------------------------------------------------------------------- 19 ## Author: Martin Maechler, Date: 19 Dec 1997; 7 Jul 2004 20 stopifnot(is.numeric(m) && is.matrix(m)) 21 method <- match.arg(method) 22 n <- length(lam <- eigen.m $values) 23 Eps <- eps.ev * abs(lam[1])# lam[1] is largest EV; "small" is *relative* 24 ## lam[n] is the SMALLEST eigenvalue 25 if(lam[n] < Eps) { # fix up small or negative values 26 switch(method, 27 "someEVadd" = lam[lam < Eps] <- Eps, 28 "allEVadd" = lam <- lam + Eps-lam[n] 29 ) 30 Q <- eigen.m $vectors 31 o.diag <- diag(m)# original one - for rescaling 32 m <- Q %*% (lam * t(Q)) ## == Q %*% diag(lam) %*% t(Q) 33 ## rescale to the original diagonal values 34 ## D <- sqrt(o.diag/diag(m)) 35 ## where they are >= Eps : 36 D <- sqrt(pmax(Eps,o.diag)/diag(m)) 37 m[] <- D * m * rep(D, each = n) ## == diag(D) %*% m %*% diag(D) 38 } 39 m 40} 41