📂행렬대수

가역행렬의 LU 분해

빌드업

행렬 $A \in \mathbb{R}^{m \times m}$ 의 왼쪽에 곱해졌을 때 $(i, j)$ 성분을 $0$ 이 되도록 하는 행렬 $E_{ij}$ 를 $A$ 에 대한 $ij$-소거연산자라고 정의해보자.

구체적으로 정방행렬 $(a_{ij}) \in \mathbb{R}^{m \times m}$ 에 대한 $E_{ij}$ 는 대각성분이 $1$ 이고 $(i,j)$-성분이 $\displaystyle -m_{ij} = -{{a_{ij}} \over {a_{jj}}}$, 나머지 성분이 $0$으로 구해진다. 이는 연립 방정식의 풀이에서 같은 변수끼리 계수를 맞춰서 소거하는 연산을 행렬로 나타낸 것이다. 이러한 소거는 연립 방정식의 해가 단 하나만 존재할 때 의미가 있으므로, 우리가 다룰 $A$ 는 가역행렬로 가정한다.

예시

간단한 예로써 $$ B := \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix} $$ 라고 하면 $E_{21} B = \begin{bmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix}$가 된다. 그리고 $B = \begin{bmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{bmatrix}$이므로 $B$ 에 대한 $(2,1)$-소거연산자는 $$ E_{21} = \begin{bmatrix} 1 & 0 & 0 \\ -m_{ij} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ -{{1} \over {2}} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ 으로 구할 수 있다. 한편 이러한 연산을 여러번 취한 $$ U_{B} : = E_{32} E_{31} E_{21} B = \begin{bmatrix} 2 & 1 & 0 \\ 0 & 2 & 1 \\ 0 & 0 & 2 \end{bmatrix} $$ 는 상삼각행렬이고, $E_{21} , E_{31} , E_{32}$ 은 모두 하삼각행렬이다. 하삼각행렬끼리의 곱은 하삼각행렬이므로 $L_{B} : = E_{21} E_{31} E_{32}$ 도 하삼각행렬이고 $L^{-1} B = U$ 으로 나타낼 수 있다. 다시 $B$ 에 대해 정리하면 $B = L_{B} U_{B}$ 즉 하^lower삼각행렬과 상^upper삼각행렬의 곱으로 분해한 꼴이 된다.

이러한 분해법은 너무나 단순한 덕분에 $A \in \mathbb{R}^{n \times n}$ 인 경우에 대해서도 별다른 문제 없이 일반화가 가능하다. 이를 이용한 연립방정식 $A \mathbf{x} = \mathbf{b}$ 의 풀이법을 생각해보면 $A$ 가 $A = LU$ 로 분해되므로 $$ \begin{align*} LU \mathbf{x} =& \mathbf{b} \\ \mathbf{y} :=& U \mathbf{x} \end{align*} $$ 라 둘 수 있다. 여기서 $U$ 가 상삼각행렬이므로 $\mathbf{y}$ 는 후진대입법을 통해 쉽게 구할 수 있고, $$ \begin{align*} L (U \mathbf{x} ) =& L \\ \mathbf{y} =& \mathbf{b} \end{align*} $$ 에서 $L$ 이 하삼각행렬이므로 $\mathbf{b}$ 는 전진대입법을 통해 쉽게 구할 수 있다. 따라서 LU 분해를 실제로 구할 수 있다는 것은 상당히 많은 문제를 상당히 쉽고 빠르게 풀 수 있다는 뜻이 된다.

유도

이제 직접적으로 $L$ 과 $U$ 를 찾아내기 위해, $A := L U$ 를 만족하는 하삼각행렬 $L$ 과 상삼각행렬 $U$ 를 다음과 같이 정의하자. $$ \begin{align*} L : =& \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ l_{21} & 1 & 0 & \cdots & 0 \\ l_{31} & l_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \cdots & 1 \end{bmatrix} \\ U : =& \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ 0 & u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & u_{nn} \end{bmatrix}\end{align*} $$ 여기서 $L$ 은 소거연산자들의 곱이므로 대각성분들이 모두 $1$ 임을 쉽게 확인할 수 있다. 위 정의에 따라 직접 $A = (a_{ij})$ 의 성분 $a_{ij}$ 를 계산해보면 $$ a_{ij} = \sum_{s = 1}^{n} l_{is} u_{sj} = \sum_{s = 1}^{min(i , j)} l_{is} u_{sj} $$ 이 된다. 이제 $a_{ij}$ 를 대각선 상에 있는 경우, 대각선 윗쪽에 있는 경우, 대각성 아랫쪽에 있는 경우로 나누어서 생각해보자.

$i = j = k$이면

$$ \begin{align*} a_{kk} =& \sum_{s = 1}^{\min(i, j)} l_{ks} u_{sk} \\ =& \sum_{s = 1}^{k} l_{ks} u_{sk} \\ =& \sum_{s = 1}^{k-1} l_{ks} u_{sk} + l_{kk} u_{kk} \\ =& \sum_{s = 1}^{k-1} l_{ks} u_{sk} + u_{kk} \end{align*} $$

이를 $u_{kk}$ 에 대해서 정리하면

$$ u_{kk} = a_{kk} - \sum_{s = 1}^{k-1} l_{ks} u_{sk} $$

를 얻을 수 있다. 비슷하게 $k < j$ 일 땐 $$ a_{kj} = \sum_{s = 1}^{k-1} l_{ks} u_{sj} + u_{kj} $$ 를, $k < i$ 일 땐 $$ a_{ik} = \sum_{s = 1}^{k-1} l_{is} u_{sk} + l_{ik} u_{kk} $$ 를 얻을 수 있다. 이를 각각 $u_{kj}$ 와 $l_{ik}$ 에 대해서 정리하면 다음을 얻는다. $$ \begin{align*} u_{kj} =& a_{kj} - \sum_{s = 1}^{k-1} l_{ks} u_{sj} \\ l_{ik} =& {{ 1 } \over {u_{kk} }} \left\{ a_{ik} - \sum_{s = 1}^{k-1} l_{is} u_{sk} \right\} \end{align*} $$

■

알고리즘

$(a_{ij}) \in \mathbb{R}^{n \times n}$ 가 가역행렬이라고 하자.

Step 1. $k = 1$

$u_{1j} = a_{1j}$ 을 대입하고 $\displaystyle l_{i1} = {{1} \over {u_{11}}} a_{i1}$ 을 계산한다.

---

Step 2. $k = 2, 3, \cdots , n-1$

• Step 2-1. 다음을 계산한다. $$ u_{kk} = a_{kk} - \sum_{s = 1}^{k-1} l_{ks} u_{sk} $$
• Step 2-2. $j = k+1 , k+2 , \cdots , n-1$ 에 대해 다음을 계산한다. $$ u_{kj} = a_{kj} - \sum_{s = 1}^{k-1} l_{ks} u_{sj} $$
• Step 2-3. $i = k+1 , k+2 , \cdots , n-1$ 에 대해 다음을 계산한다. $$ l_{ik} = {{ 1 } \over {u_{kk} }} \left\{ a_{ik} - \sum_{s = 1}^{k-1} l_{is} u_{sk} \right\} $$

---

Step 3. $k = n$ 에 대해 다음을 계산한다. $$ u_{nn} = a_{nn} - \sum_{s = 1}^{n-1} l_{ns} u_{sn} $$

---

• 알고리즘을 잘 살펴보면 $k$ 가 증가함에 따라 $L$ 은 왼쪽에서 오른쪽으로, $U$ 는 위에서 아래로 구해짐을 알 수 있다. 한 단계를 거칠 때 필요한 값들은 적어도 그 전 단계에서 모두 구해졌으므로 단계마다 계산하는 단위는 한 행이나 한 열이 된다.

코드

이제 앞에서 예시로 쓰였던 행렬 $B$ 를 실제로 분해해보자.

아래의 코드는 위의 알고리즘을 실제로 R 에서 구현한 것이다.

B <- matrix(c(2,1,0,1,2,1,0,1,2), ncol= 3); B
 
LU<-function(A)
{
  n=length(A[,1])
  L<-diag(1,n)
  U<-diag(0,n)
 
  k=1
  j=1
  i=1
  U[1,]<-A[1,]
  L[,1]<-A[,1]/A[1,1]
  
  for(k in 2:(n-1))
  {
    U[k,k] = A[k,k] - sum(L[k,1:(k-1)]*U[1:(k-1),k])
    for(j in (k+1):n)
    {
      U[k,j] = A[k,j] - sum(L[k,1:(k-1)]*U[1:(k-1),j])
    }
    for(i in (k+1):n)
    {
      L[i,k] = ( A[i,k] - sum(L[i,1:(k-1)]*U[1:(k-1),k]) ) / U[k,k]
    }
  }
  
  U[n,n] = A[n,n] - sum(L[n,1:(n-1)]*U[1:(n-1),n])
  
  return(list(L=L,U=U))
}
 
LU(B)
require(Matrix)
expand(lu(B))

실행 결과는 다음과 같다.

R 자체에서 구현된 LU 분해와 결과가 정확히 일치함을 확인 할 수 있다.

같이보기

• PLU 분해