numeric n * n approximately positive
definite matrix, typically an approximation to a correlation or
covariance matrix. If x is not symmetric (and
ensureSymmetry is not false), symmpart(x) is used.
corr
logical indicating if the matrix should be a
correlation matrix.
keepDiag
logical, generalizing corr: if TRUE, the
resulting matrix should have the same diagonal
(diag(x)) as the input matrix.
do2eigen
logical indicating if a
posdefify() eigen step should be applied to
the result of the Higham algorithm.
doSym
logical indicating if X <- (X + t(X))/2 should be
done, after X <- tcrossprod(Qd, Q); some doubt if this is necessary.
doDykstra
logical indicating if Dykstra's correction should be
used; true by default. If false, the algorithm is basically the
direct fixpoint iteration
Y(k) = P_U(P_S(Y(k-1))).
only.values
logical; if TRUE, the result is just the
vector of eigen values of the approximating matrix.
ensureSymmetry
logical; by default, symmpart(x)
is used whenever isSymmetric(x) is not true. The user
can explicitly set this to TRUE or FALSE, saving the
symmetry test. Beware however that setting it FALSE
for an asymmetric input x, is typically nonsense!
eig.tol
defines relative positiveness of eigenvalues compared
to largest one, λ_1. Eigen values λ_k are
treated as if zero when λ_k / λ_1 ≤ eig.tol.
conv.tol
convergence tolerance for Higham algorithm.
posd.tol
tolerance for enforcing positive definiteness (in the
final posdefify step when do2eigen is TRUE).
maxit
maximum number of iterations allowed.
conv.norm.type
convergence norm type (norm(*,
type)) used for Higham algorithm. The default is "I"
(infinity), for reasons of speed (and back compatibility); using
"F" is more in line with Higham's proposal.
trace
logical or integer specifying if convergence monitoring
should be traced.
Details
This implements the algorithm of Higham (2002), and then (if
do2eigen is true) forces positive definiteness using code from
posdefify. The algorithm of Knol DL and ten
Berge (1989) (not implemented here) is more general in (1) that it
allows constraints to fix some rows (and columns) of the matrix and (2)
to force the smallest eigenvalue to have a certain value.
Note that setting corr = TRUE just sets diag(.) <- 1
within the algorithm.
Higham (2002) uses Dykstra's correction, but the version by Jens
Oehlschlaegel did not use it (accidentally), and has still lead to
good results; this simplification, now only via doDykstra = FALSE,
was active in nearPD() upto Matrix version 0.999375-40.
Value
If only.values = TRUE, a numeric vector of eigen values of the
approximating matrix;
Otherwise, as by default, an S3 object of class"nearPD", basically a list with components
mat
a matrix of class dpoMatrix, the
computed positive-definite matrix.
eigenvalues
numeric vector of eigen values of mat.
corr
logical, just the argument corr.
normF
the Frobenius norm (norm(x-X, "F")) of the
difference between the original and the resulting matrix.
iterations
number of iterations needed.
converged
logical indicating if iterations converged.
Author(s)
Jens Oehlschlaegel donated a first version. Subsequent changes
by the Matrix package authors.
References
Cheng, Sheung Hun and Higham, Nick (1998)
A Modified Cholesky Algorithm Based on a Symmetric Indefinite Factorization;
SIAM J. Matrix Anal. Appl., 19, 1097–1110.
Knol DL, ten Berge JMF (1989)
Least-squares approximation of an improper correlation matrix by a
proper one.
Psychometrika54, 53–61.
Higham, Nick (2002)
Computing the nearest correlation matrix - a problem from finance;
IMA Journal of Numerical Analysis22, 329–343.
See Also
A first version of this (with non-optional corr=TRUE)
has been available as nearcor(); and
more simple versions with a similar purpose
posdefify(), both from package sfsmisc.
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")]
stopifnot(all.equal(n.A$mat[1,2], 0.760689917),
all.equal(n.A$normF, 0.52779033, tolerance=1e-9) )
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))
str(near.m <- nearPD(m, trace = TRUE))
round(near.m$mat, 2)
norm(m - near.m$mat) # 1.102 / 1.08
if(require("sfsmisc")) {
m2 <- posdefify(m) # a simpler approach
norm(m - m2) # 1.185, i.e., slightly "less near"
}
round(nearPD(m, only.values=TRUE), 9)
## 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
nc.1 <- nearPD(pr, conv.tol = 1e-7, corr = TRUE)
nc.1$iterations # 11 / 12 (!)
ncr <- nearPD(pr, conv.tol = 1e-15)
str(ncr)# still 2 iterations
ncr.1 <- nearPD(pr, conv.tol = 1e-15, corr = TRUE)
ncr.1 $ iterations # 27 / 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
norm(pr - ncr.1$mat) # = 0.08746 / 0.08805
### 3) a real data example from a 'systemfit' model (3 eq.):
(load(system.file("external", "symW.rda", package="Matrix"))) # "symW"
dim(symW) # 24 x 24
class(symW)# "dsCMatrix": sparse symmetric
if(dev.interactive()) image(symW)
EV <- eigen(symW, only=TRUE)$values
summary(EV) ## looking more closely {EV sorted decreasingly}:
tail(EV)# all 6 are negative
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")
abline(0,1, col="red3",lty=2)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Matrix)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Matrix/nearPD.Rd_%03d_medium.png", width=480, height=480)
> ### Name: nearPD
> ### Title: Nearest Positive Definite Matrix
> ### Aliases: nearPD
> ### Keywords: algebra array
>
> ### ** 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
> stopifnot(all.equal(n.A$mat[1,2], 0.760689917),
+ all.equal(n.A$normF, 0.52779033, tolerance=1e-9) )
>
> 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
.. ..@ x : num [1:25] 1.313 0.406 -0.242 -0.852 -0.753 ...
.. ..@ Dim : int [1:2] 5 5
.. ..@ Dimnames:List of 2
.. .. ..$ : NULL
.. .. ..$ : NULL
.. ..@ 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 : num 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
[1] 1.184855
>
> 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
.. ..@ x : num [1:36] 1.006 0.485 0.643 0.487 0.646 ...
.. ..@ Dim : int [1:2] 6 6
.. ..@ Dimnames:List of 2
.. .. ..$ : NULL
.. .. ..$ : NULL
.. ..@ 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 : num 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 11030000 5536 177500000
> 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")
> abline(0,1, col="red3",lty=2)
>
>
>
>
>
>
> dev.off()
null device
1
>