R: Examples for 'Matrix::nearPD'

### ** Examples

 ## Higham(2002), p.334f - simple example
 A <- matrix(1, 3,3); A[1,3] <- A[3,1] <- 0
 n.A <- nearPD(A, corr=TRUE, do2eigen=FALSE)
 n.A[c("mat", "normF")]

$mat
3 x 3 Matrix of class "dpoMatrix"
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.7606899 0.1572981
[2,] 0.7606899 1.0000000 0.7606899
[3,] 0.1572981 0.7606899 1.0000000

$normF
[1] 0.5277903

 n.A.m <- nearPD(A, corr=TRUE, do2eigen=FALSE, base.matrix=TRUE)$mat
 stopifnot(exprs = {                           #=--------------
   all.equal(n.A$mat[1,2], 0.760689917)
   all.equal(n.A$normF, 0.52779033, tolerance=1e-9)
   all.equal(n.A.m, unname(as.matrix(n.A$mat)), tolerance = 1e-15)# seen rel.d.= 1.46e-16
 })
 set.seed(27)
 m <- matrix(round(rnorm(25),2), 5, 5)
 m <- m + t(m)
 diag(m) <- pmax(0, diag(m)) + 1
 (m <- round(cov2cor(m), 2))

      [,1]  [,2]  [,3]  [,4]  [,5]
[1,]  1.00  0.65 -0.46 -1.15 -0.76
[2,]  0.65  1.00  0.58  0.50 -0.90
[3,] -0.46  0.58  1.00 -0.45 -0.32
[4,] -1.15  0.50 -0.45  1.00  0.25
[5,] -0.76 -0.90 -0.32  0.25  1.00

 str(near.m <- nearPD(m, trace = TRUE))

iter   1 : #{p}=4, ||Y-X|| / ||Y||=    0.268591
iter   2 : #{p}=4, ||Y-X|| / ||Y||=           0
List of 7
 $ mat        :Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
  .. ..@ Dim     : int [1:2] 5 5
  .. ..@ Dimnames:List of 2
  .. .. ..$ : NULL
  .. .. ..$ : NULL
  .. ..@ x       : num [1:25] 1.313 0.406 -0.242 -0.852 -0.753 ...
  .. ..@ uplo    : chr "U"
  .. ..@ factors : list()
 $ eigenvalues: num [1:5] 2.80 1.83 1.23 7.70e-02 2.80e-08
 $ corr       : logi FALSE
 $ normF      : num 0.938
 $ iterations : int 2
 $ rel.tol    : num 0
 $ converged  : logi TRUE
 - attr(*, "class")= chr "nearPD"

 round(near.m$mat, 2)

5 x 5 Matrix of class "dpoMatrix"
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,]  1.31  0.41 -0.24 -0.85 -0.75
[2,]  0.41  1.19  0.41  0.27 -0.91
[3,] -0.24  0.41  1.15 -0.24 -0.32
[4,] -0.85  0.27 -0.24  1.28  0.26
[5,] -0.75 -0.91 -0.32  0.26  1.00

 norm(m - near.m$mat) # 1.102 / 1.08

[1] 1.079735

 if(require("sfsmisc")) {
    m2 <- posdefify(m) # a simpler approach
    norm(m - m2)  # 1.185, i.e., slightly "less near"
 }

Loading required package: sfsmisc

Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
logical.return = TRUE, : there is no package called 'sfsmisc'

 round(nearPD(m, only.values=TRUE), 9)

[1] 2.800681404 1.831722441 1.229003616 0.076994641 0.000000028

## A longer example, extended from Jens' original,
## showing the effects of some of the options:

pr <- Matrix(c(1,     0.477, 0.644, 0.478, 0.651, 0.826,
               0.477, 1,     0.516, 0.233, 0.682, 0.75,
               0.644, 0.516, 1,     0.599, 0.581, 0.742,
               0.478, 0.233, 0.599, 1,     0.741, 0.8,
               0.651, 0.682, 0.581, 0.741, 1,     0.798,
               0.826, 0.75,  0.742, 0.8,   0.798, 1),
             nrow = 6, ncol = 6)

nc.  <- nearPD(pr, conv.tol = 1e-7) # default
nc.$iterations  # 2

[1] 2

nc.1 <- nearPD(pr, conv.tol = 1e-7, corr = TRUE)
nc.1$iterations # 11 / 12 (!)

[1] 12

ncr   <- nearPD(pr, conv.tol = 1e-15)
str(ncr)# still 2 iterations

List of 7
 $ mat        :Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
  .. ..@ Dim     : int [1:2] 6 6
  .. ..@ Dimnames:List of 2
  .. .. ..$ : NULL
  .. .. ..$ : NULL
  .. ..@ x       : num [1:36] 1.006 0.485 0.643 0.487 0.646 ...
  .. ..@ uplo    : chr "U"
  .. ..@ factors : list()
 $ eigenvalues: num [1:6] 4.213 0.771 0.515 0.385 0.178 ...
 $ corr       : logi FALSE
 $ normF      : num 0.0626
 $ iterations : int 2
 $ rel.tol    : num 0
 $ converged  : logi TRUE
 - attr(*, "class")= chr "nearPD"

ncr.1 <- nearPD(pr, conv.tol = 1e-15, corr = TRUE)
ncr.1 $ iterations # 27 / 30 !

[1] 30

ncF <- nearPD(pr, conv.tol = 1e-15, conv.norm = "F")
stopifnot(all.equal(ncr, ncF))# norm type does not matter at all in this example

## But indeed, the 'corr = TRUE' constraint did ensure a better solution;
## cov2cor() does not just fix it up equivalently :
norm(pr - cov2cor(ncr$mat)) # = 0.09994

[1] 0.09994443

norm(pr -       ncr.1$mat)  # = 0.08746 / 0.08805

[1] 0.08804689

### 3) a real data example from a 'systemfit' model (3 eq.):
(load(system.file("external", "symW.rda", package="Matrix"))) # "symW"

[1] "symW"

dim(symW) #  24 x 24

[1] 24 24

class(symW)# "dsCMatrix": sparse symmetric

[1] "dsCMatrix"
attr(,"package")
[1] "Matrix"

if(dev.interactive())  image(symW)
EV <- eigen(symW, only=TRUE)$values
summary(EV) ## looking more closely {EV sorted decreasingly}:

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
        0        37       999  11031812      5536 177467251

tail(EV)# all 6 are negative

[1] -1.615418e-05 -4.803894e-05 -9.448841e-05 -8.402689e-04 -1.658025e-03
[6] -5.476273e-03

EV2 <- eigen(sWpos <- nearPD(symW)$mat, only=TRUE)$values
stopifnot(EV2 > 0)
if(require("sfsmisc")) {
       plot(pmax(1e-3,EV), EV2, type="o", log="xy", xaxt="n",yaxt="n")
       eaxis(1); eaxis(2)
} else plot(pmax(1e-3,EV), EV2, type="o", log="xy")

Loading required package: sfsmisc

Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
logical.return = TRUE, : there is no package called 'sfsmisc'

abline(0,1, col="red3",lty=2)

Examples for 'Matrix::nearPD'

Nearest Positive Definite Matrix