I asked a question earlier on Inverse of $A+B$ and $A+BCD$?
I got a very good hint on using the Sherman–Morrison–Woodbury formula. I used R to code the idea and using some arbitrary matrices, I tested the formula.
With low dimensions, the formula works fine but once the dimension increases (about $20\times 20$ or more), the part $(C^{-1} + V A^{-1} U)$ becomes singular and therefore $(C^{-1} + V A^{-1} U)^{-1}$ does not exist!
Does anyone know how to fix this? Wikipedia says that Sherman–Morrison–Woodbury formula is an exact formula which means its exactly equivalent to $(A + UCV)^{-1}$. It it's equivalent, whenever the inverse of $A + UCV$ exists, the Sherman–Morrison–Woodbury formula should also exist, right? then why I get error?
For your reference, below is the code I wrote in R:
lambda2 <- 0.5
eta2 <- 1
rho2 <- 1
sigma2 <- 0.25
alpha <- 2
x <- seq(from = -2, to = 2, by = 0.2) # try by = 0.1, 0.05, 0.01
lengh.x <- length(x)
distanceMat.alpha <- ( abs(matrix(x, nrow = length(x), ncol = length(x)) - matrix(x,
nrow = length(x), ncol = length(x), byrow = TRUE)) )^alpha
OneMat <- rep(1,lengh.x)%*%t(rep(1,lengh.x))
H <- lambda2*OneMat
G <- sigma2*diag(length(x))
A <- G + H
B <- eta2*exp(-rho2*(distanceMat.alpha))
CovMat.test <- A + B
# Getting Inverse using solve():
CovMat.test.Inv.SLOW <- solve(CovMat.test)
# Getting inverse using Sherman-Morrison-Woodbury formula:
A.inv <- (1/sigma2)*(diag(lengh.x) - (lambda2/(sigma2 + lengh.x*lambda2)) * OneMat)
B.eigen <- eigen(B, symmetric = TRUE)
Delta <- diag(B.eigen$values)
Delta.inv <- diag(1/B.eigen$values)
P <- B.eigen$vectors
# B can be decomposed as: P%*%Delta%*%t(P)
CovMat.test.Inv.FAST <- A.inv - A.inv%*%P%*%solve(Delta.inv +
t(P)%*%A.inv%*%P)%*%t(P)%*%A.inv
# Test:
TwoMatAreEq <- sum(round(CovMat.test.Inv.FAST, 10) == round(CovMat.test.Inv.SLOW, 10))
== (dim(CovMat.test)[1])^2
TwoMatAreEq