Pregunta

I am using the choleski decomposition to compute the inverse of a matrix that is positive semidefinite. However, when my matrix becomes extremely large and has zeros in it I have that my matrix is no longer (numerically from the computers point of view) positive definite. So to get around this problem I use the pivot = TRUE option in the choleski command in R. However, (as you will see below) the two return the same output but with the rows and columns or the matrix rearranged. I am trying to figure out is there a way (or transformation) to make them the same. Here is my code:

X = matrix(rnorm(9),nrow=3)
A = X%*%t(X)

inv1 = function(A){
    Q = chol(A)
    L = t(Q)
    inverse = solve(Q)%*%solve(L)
    return(inverse)
}


inv2 = function(A){
    Q = chol(A,pivot=TRUE)
    L = t(Q)
    inverse = solve(Q)%*%solve(L)
    return(inverse)
}

Which when run results in:

 > inv1(A)
              [,1]      [,2]      [,3]
    [1,]  9.956119 -8.187262 -4.320911
    [2,] -8.187262  7.469862  3.756087
    [3,] -4.320911  3.756087  3.813175
    > 
    > inv2(A)
              [,1]      [,2]      [,3]
    [1,]  7.469862  3.756087 -8.187262
    [2,]  3.756087  3.813175 -4.320911
    [3,] -8.187262 -4.320911  9.956119

Is there a way to get the two answers to match? I want inv2() to return the answer from inv1().

¿Fue útil?

Solución

That is explained in ?chol: the column permutation is returned as an attribute.

inv2 <- function(A){
  Q <- chol(A,pivot=TRUE)
  Q <- Q[, order(attr(Q,"pivot"))]
  Qi <- solve(Q)
  Qi %*% t(Qi)
}
inv2(A)
solve(A)  # Identical

Otros consejos

Typically

M = matrix(rnorm(9),3)
M
           [,1]        [,2]       [,3]
[1,]  1.2109251 -0.58668426 -0.4311855
[2,] -0.8574944  0.07003322 -0.6112794
[3,]  0.4660271 -0.47364400 -1.6554356
library(Matrix)
pm1 <- as(as.integer(c(2,3,1)), "pMatrix")
M %*% pm1
           [,1]       [,2]        [,3]
[1,] -0.4311855  1.2109251 -0.58668426
[2,] -0.6112794 -0.8574944  0.07003322
[3,] -1.6554356  0.4660271 -0.47364400
Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top