R nearPD 最近正定矩阵

说明

计算最接近的正定矩阵，通常是相关矩阵或方差-协方差矩阵。

用法

nearPD(x, corr = FALSE, keepDiag = FALSE, base.matrix = FALSE,
       do2eigen = TRUE, doSym = FALSE,
       doDykstra = TRUE, only.values = FALSE,
       ensureSymmetry = !isSymmetric(x),
       eig.tol = 1e-06, conv.tol = 1e-07, posd.tol = 1e-08,
       maxit = 100, conv.norm.type = "I", trace = FALSE)

参数

x	数值n \times n 近似正定矩阵，通常是相关或协方差矩阵的近似值。如果x 不对称(且ensureSymmetry 不为假)，则使用symmpart(x)。
corr	逻辑指示矩阵是否应该是相关矩阵。
keepDiag	逻辑，概括 corr ：如果 TRUE ，结果矩阵应具有与输入矩阵相同的对角线 (diag(x) )。
base.matrix	逻辑指示生成的 mat 组件是否应该是 base matrix 或(默认情况下)类 dpoMatrix 的 Matrix 。
do2eigen	逻辑指示是否应将 posdefify() 特征步骤应用于 Higham 算法的结果。
doSym	逻辑指示X <- (X + t(X))/2是否应该在X <- tcrossprod(Qd, Q)之后完成；有人怀疑这是否有必要。
doDykstra	逻辑指示是否应使用 Dykstra 修正；默认为 true。如果为 false，则该算法本质上是直接固定点迭代 Y_k = P_U(P_S(Y_{k-1})) 。
only.values	逻辑性；如果 TRUE ，结果只是近似矩阵的特征值向量。
ensureSymmetry	逻辑性；默认情况下，只要 isSymmetric(x) 不为 true，就使用 symmpart(x)。用户可以显式地将其设置为 TRUE 或 FALSE ，保存对称性测试。但请注意，为非对称输入 x 设置 FALSE 通常是无意义的！
eig.tol	定义特征值与最大特征值 \lambda_1 相比的相对正值。当 \lambda_k / \lambda_1 \le eig.tol 时，特征值 \lambda_k 被视为零。
conv.tol	Higham 算法的收敛容差。
posd.tol	强制执行正定性的公差(在最终 posdefify 步骤中，当 do2eigen 为 TRUE 时)。
maxit	允许的最大迭代次数。
conv.norm.type	用于 Higham 算法的收敛范数类型 (norm(*, type))。出于速度(和向后兼容性)的原因，默认值为"I"(无穷大)；使用"F"更符合Higham的提议。
trace	逻辑或整数指定是否应跟踪收敛监控。

细节

这实现了 Higham (2002) 的算法，然后(如果 do2eigen 为真)使用 posdefify 中的代码强制正定性。 Knol 和 ten Berge (1989) 的算法(此处未实现)更通用，因为它允许约束 (1) 固定矩阵的某些行(和列)以及 (2) 强制最小特征值具有特定值。

请注意，设置 corr = TRUE 只是在算法内设置 diag(.) <- 1。

Higham (2002) 使用了 Dykstra 的修正，但 Jens Oehlschlägel 的版本没有使用它(无意中)，但仍然给出了合理的结果；这种简化现在仅在 doDykstra = FALSE 在 nearPD() 直至 Matrix 版本 0.999375-40 中处于活动状态时使用。

值

如果是 only.values = TRUE ，则为逼近矩阵特征值的数值向量；否则，默认情况下， class "nearPD" 的 S3 对象，本质上是一个包含组件的列表

mat	类 dpoMatrix 的矩阵，计算的正定矩阵。
eigenvalues	mat 的特征值的数值向量。
corr	逻辑上，只是参数 corr 。
normF	原始矩阵与结果矩阵之差的 Frobenius 范数 (norm(x-X, "F"))。
iterations	所需的迭代次数。
converged	逻辑指示迭代是否收敛。

例子

 ## 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")]
 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))

 str(near.m <- nearPD(m, trace = TRUE))
 round(near.m$mat, 2)
 norm(m - near.m$mat) # 1.102 / 1.08

 if(requireNamespace("sfsmisc")) {
    m2 <- sfsmisc::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(requireNamespace("sfsmisc")) {
    plot(pmax(1e-3,EV), EV2, type="o", log="xy", xaxt="n", yaxt="n")
    for(side in 1:2) sfsmisc::eaxis(side)
} else
    plot(pmax(1e-3,EV), EV2, type="o", log="xy")
abline(0, 1, col="red3", lty=2)

作者

Jens Oehlschlägel donated a first version. Subsequent changes by the Matrix package authors.

参考

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. Psychometrika 54, 53-61.

Higham, Nick (2002) Computing the nearest correlation matrix - a problem from finance; IMA Journal of Numerical Analysis 22, 329-343.

也可以看看

该版本的第一个版本(带有非可选的 corr=TRUE )已作为 nearcor() 提供；以及具有类似用途的更简单版本 posdefify() ，均来自包 sfsmisc 。