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