[통계계산방법론] Gaussian Elimination Algorithm과 Cholesky Algorithm

Inverse Computing  

우리는 $Ax = b$ 라는 선형 모형에서의 solution을 얘기할 때, 행렬 A는 non-singularity의 성질을 가지고 있어야 하며 그때의 solution은 $A^{-1}b$으로 구한다. 하지만 computing에 있어서 A의 inverse를 구하는 것은 너무 복잡하고 많은 시간이 소요된다. 예를 들자면, $n \times n$ 의 행렬 A는 computing 시 $O(n^2)$ 의 flop이 필요하다. 이때 행렬 A를 $I + uv^T$의 형태로 만들어줌으로써 계산에 필요한 flop은 $O(n)$로 줄어든다. 따라서 행렬의 structure를 잘 이용하면 computation의 속도를 빨리할 수 있다. 

 

R 프로그램에서는 solve() 함수를 이용하여 inverse 함수를 구할 수 있지만, 행렬 A가 다음과 같은 형태의 행렬들이라면 solve() 함수를 사용하지 않으면서 더 빠르게 solution을 도출할 수 있다. 

 

• A : diagonal(a), solution x  = b/a
• A : upper-triangular, solution x = backsolve(A, b)
• A : lower-triangular, solution x = forwardsolve(A, b)

 

Gaussian Elimination

주어진 행렬 A가 위에서 언급된 diagonal, upper or lower triangular 이 아니더라도 이러한 형태로 만들어주어 computation 속도를 빠르게 할 수 있다. 그러한 방법들 중 하나가 Gauss Elimination이다. 이 알고리즘은 행렬을 Lower Triangular와 Upper triangular로 decompostion시킨다. 

$$ Ax = b \Leftrightarrow BAx = Bb$$ 

위의 식에서 full-rank인 행렬 B를 양변에 곱해줌으로써 BA는 triangular의 형태로 만들어주는 것의 아이디어이다. 

 

Gaussian Elimination은 각 열의 diagonal 원소 아래 모든 원소를 0으로 만드는 과정을 반복한다. 그 과정에서의 스텝은 다음과 같이 $3 \times 3$ 행렬로 보여준다.  

 

Step1. Let $A^{(0)} = A $

Step2. Then

  Step3. $M^{(2)}M^{(1)}A$의 행렬은 upper-triangular 행렬이 되고, R 프로그램에서 backsolve() 함수를 이용하여 solutino을 구할 수 있다.

 

하지만 위와 같은 과정에서 만약 diangonal의 원소값이 0인 경우에는 error가 발생할 수 있다. 이러한 이유 때문에 permutation matrix를 이용하여 partial pivoting을 해줘야 한다. 이때 permutation function을 $P$라고 하면, 다음과 같이 A를 LU decomposing 할 수 있다.

위의 과정에 의해 $Ax = PAx = LUx = Pb$ 가 성립되며 $Ux = y$ 라고 했을 때 $Ly = (Pb)$에서 y를 먼저 구한다 (forward.solve() 함수 이용). 그리고 $Ux = y$에서 x를 구하여 solution을 얻을 수 있다(backsolve() 함수 이용).  

 

Cholesky Decomposition (Square-root Algorithm)

Cholesky Decompostion도 하나의 LU분할로 볼 수 있는데, 이는 행렬 A를 $L$과 그의 transpose $L^T$로 분할하는 것이다. 이때 행렬 A는 positive definite의 property를 만족해야 한다. 

※ positive definite : $y^TAy >0$가 성립되는 행렬 A

 

Cholesky decomposition의 방법은 다음과 같다. 

1. $L^{(1)}$의 $L_{11}$을 $\sqrt{A_{11}}$로 initializing한다. 

2. k번 째의 과정에서 행렬 A가 다음과 같이 decomposed 되었을 때, 

$L^{(k-1)}l^{(k)} = a^{(k)}$를 이용하여  $l^{(k)} = (L_{k1}, ..., L_{kk-1})^T $를 업데이트 한다.

또한 $L^2_{kk} = A_{kk} - l^{(k)T}$를 이용하여 $L_{kk}$를 업데이트 한다. 

 

 

Reference 

• 고려대학교 신승준 교수님 ST509 강의자료