Find a Close Positive Definite Matrix

From a matrix m, construct a "close" positive definite one.

posdefify(m, method = c("someEVadd", "allEVadd"),
          symmetric = TRUE, eigen.m = eigen(m, symmetric= symmetric),
          eps.ev = 1e-07)

Arguments

m: a numeric (square) matrix.
method: a string specifying the method to apply; can be abbreviated.
symmetric: logical, simply passed to eigen (unless eigen.m is specified); currently, we do not see any reason for not using TRUE.
eigen.m: the eigen value decomposition of m, can be specified in case it is already available.
eps.ev: number specifying the tolerance to use, see Details below.

Details

We form the eigen decomposition $$m = V \Lambda V'$$ where $\Lambda$ is the diagonal matrix of eigenvalues, $\Lambda_{j,j} = \lambda_j$, with decreasing eigenvalues $\lambda_1 \ge \lambda_2 \ge \ldots \ge \lambda_n$.

When the smallest eigenvalue $\lambda_n$ are less than Eps <- eps.ev * abs(lambda[1]), i.e., negative or “almost zero”, some or all eigenvalues are replaced by positive (>= Eps) values, $\tilde\Lambda_{j,j} = \tilde\lambda_j$. Then, $\tilde m = V \tilde\Lambda V'$ is computed and rescaled in order to keep the original diagonal (where that is >= Eps).

Value

a matrix of the same dimensions and the “same” diagonal (i.e. diag) as m but with the property to be positive definite.

Author

Martin Maechler, July 2004

Note

As we found out, there are more sophisticated algorithms to solve this and related problems. See the references and the nearPD() function in the Matrix package. We consider nearPD() to also be the successor of this package's nearcor().

References

Section 4.4.2 of Gill, P.~E., Murray, W. and Wright, M.~H. (1981) Practical Optimization, Academic Press.

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.

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

Lucas (2001) Computing nearest covariance and correlation matrices. A thesis submitted to the University of Manchester for the degree of Master of Science in the Faculty of Science and Engeneering.

Examples

 set.seed(12)
 m <- matrix(round(rnorm(25),2), 5, 5); m <- 1+ m + t(m); diag(m) <- diag(m) + 4
 m
#>       [,1]  [,2]  [,3]  [,4]  [,5]
#> [1,]  2.04  2.31 -0.74 -0.62 -0.78
#> [2,]  2.31  4.36 -0.92  2.08  3.44
#> [3,] -0.74 -0.92  3.44  1.35  1.86
#> [4,] -0.62  2.08  1.35  6.02  0.41
#> [5,] -0.78  3.44  1.86  0.41  2.94
 posdefify(m)
#>        [,1]   [,2]   [,3]   [,4]   [,5]
#> [1,]  2.040  1.560 -0.887 -0.408 -0.246
#> [2,]  1.560  4.360 -0.533  1.706  2.345
#> [3,] -0.887 -0.533  3.440  1.195  1.353
#> [4,] -0.408  1.706  1.195  6.020  0.578
#> [5,] -0.246  2.345  1.353  0.578  2.940
 1000 * zapsmall(m - posdefify(m))
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0  750  147 -212 -534
#> [2,]  750    0 -387  374 1095
#> [3,]  147 -387    0  155  507
#> [4,] -212  374  155    0 -168
#> [5,] -534 1095  507 -168    0