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 。