Matrix Inverse

보통 역행렬을 계산할 때 solve() 함수를 사용하는데 solve 함수는 computational cost가 높다. solve 함수를 사용하지 않고 역행렬을 구하는 방법 중 하나는 Cholesky decomposition을 이용하는 것이다.

Cholesky decomposition

Cholesky decomposition은 어떤 행렬 \(A\)를 하삼각 행렬과 상삼각 행렬의 곱으로 분해하는 것이다.

\[A = L L^T\]

• \(L\)은 하삼각 행렬

chol()함수를 사용하면 상삼각 행렬을 만들어 주고, transpose를 하면 하삼각 행렬이 된다.

A <- matrix(c(
  4,2,2,4,
  2,5,7,0,
  2,7,19,11,
  4,0,11,25), 4, 4, byrow = T)
A

##      [,1] [,2] [,3] [,4]
## [1,]    4    2    2    4
## [2,]    2    5    7    0
## [3,]    2    7   19   11
## [4,]    4    0   11   25

L <- t(chol(A))
L

##      [,1] [,2] [,3] [,4]
## [1,]    2    0    0    0
## [2,]    1    2    0    0
## [3,]    1    3    3    0
## [4,]    2   -1    4    2

L%*%t(L)

##      [,1] [,2] [,3] [,4]
## [1,]    4    2    2    4
## [2,]    2    5    7    0
## [3,]    2    7   19   11
## [4,]    4    0   11   25

같은 기능을 하는 사용자 정의 함수는 다음과 같다.

my.chol <- function(A)
{
  n <- ncol(A)
  # initalization
  Lk <- sqrt(A[1,1])
  k <- 2
  for (k in 2:n) {
    ak <- A[1:(k-1),k]
    ell.k <- forwardsolve(Lk, ak)
    ell.kk <- sqrt(A[k,k] - crossprod(ell.k, ell.k))
    Lk <- rbind(cbind(Lk, rep(0, k-1)), c(ell.k, ell.kk))
  }
  L <- Lk
  return(L)  
}

• crossprod() 함수는 \(x^Ty\)를 계산할 때 t(x) %*% y 대신 사용할 수 있고, 속도가 조금 더 빠르다.

my.chol(A)

##      Lk       
## [1,]  2  0 0 0
## [2,]  1  2 0 0
## [3,]  1  3 3 0
## [4,]  2 -1 4 2

Matrix inverse using Cholesky decomposition

\(x^TA^{-1}x\)을 solve 함수를 사용해서 계산해보자.

x <- 1:4

# naive
xAx1 <- x %*% solve(A) %*% x
xAx1

##          [,1]
## [1,] 3.737847

같은 계산을 Cholesky decomposition을 이용해서 해보자.

1) \(A\)를 \(LL^T\)로 바꿈

2) \(Ly=x\)에서 \(y=L^{-1}x\)를 forwardsolve()함수를 이용해서 구함

3) \(y^Ty=x^{T}(L^{-1})^TL^{-1}x=x^T(LL^T)^{-1}x=x^TA^{-1}x\)

# use Cholesky
y <- forwardsolve(L, x)
xAx2 <- t(y) %*% y
xAx2

##          [,1]
## [1,] 3.737847

• forwardsolve() 함수는 \(L\)이 하삼각 행렬일 때 \(Lx=b\)형태의 문제를 풀기 위해 사용

• forwardsolve(L,b)를 쓰면 x를 구할 수 있음

\(B\)가 행렬일 때 \(B^TA^{-1}B\)도 같은 방식으로 풀 수 있다.

# compute t(B) %*% inv(A) %*% B
B <- matrix(12:1, 4, 3)

# naive
BXB1 <- t(B) %*% solve(A) %*% B
BXB1

##          [,1]     [,2]    [,3]
## [1,] 47.95139 30.63194 13.3125
## [2,] 30.63194 19.78472  8.9375
## [3,] 13.31250  8.93750  4.5625

# use Cholesky
C <- forwardsolve(L, B)
BXB2 <- t(C) %*% C
BXB2

##          [,1]     [,2]    [,3]
## [1,] 47.95139 30.63194 13.3125
## [2,] 30.63194 19.78472  8.9375
## [3,] 13.31250  8.93750  4.5625

Calculate OLS estimator using Cholesky decomposition

\(y=X\beta+e, \quad \quad e \sim N_p(0, \sigma^2 I_p)\)

\(\underset{\beta}{min}(y-X\beta)^T(y-X\beta)\)

\(X^TX\beta=X^Ty\)

\(\hat{\beta_{OLS}}=(X^TX)^{-1}X^Ty\)

(1) \(X^TX, X^Ty\) 계산

(2) \(X^TX=LL^T\)로 놓고 Cholesky decomposition을 통해 L을 구함

(3) \(LL^T\beta=X^Ty\)에서 \(L^T\beta=w\)라고 놓으면 \(Lw=X^Ty\)된다.

• forwardsolve(\(L\), \(X^Ty\))=\(w\)로 \(w\)를 구함

(4) \(L^T\beta=w\)에서 \(L^T=R\)(Right triangular = upper triangular)이므로 \(R\beta=w\)가 된다.

• backsolve(\(R\), \(w\))=\(\beta\)로 \(\beta\)를 구함

X <- cbind(rep(1, 6), 1:6, (1:6)^2)
y <- c(4.5, 5.5, 6.5, 8, 10, 12)

# use lm
obj <- lm(y ~ -1 + X) # 이미 X에 constant term이 포함되어 있어서 빼줌
obj$coefficients

##        X1        X2        X3 
## 4.0000000 0.3750000 0.1607143

L <- t(chol(t(X) %*% X))
L

##           [,1]    [,2]     [,3]
## [1,]  2.449490  0.0000 0.000000
## [2,]  8.573214  4.1833 0.000000
## [3,] 37.150594 29.2831 6.110101

w <- forwardsolve(L, t(X)%*%y)
beta <- backsolve(t(L), w)
beta

##           [,1]
## [1,] 4.0000000
## [2,] 0.3750000
## [3,] 0.1607143