LU Decomposition Algorithm

LU Decomposition

정칙 정방 행렬 (Non-singular Square Matrix) 은 LU 분해가 가능하다. LU 분해는 Ax = b 형태의 Linear System 을 직접적으로 (directly) 푸는 효율적인 방법이다. 일단 LU 분해를 해놓으면 b 의 값에 상관없이 두번의 후 대입 (back substitution) 만으로 해를 구할 수 있다. LU 분해는 유일한 형태가 아니며, 그 형태에 따라 Doolittle, Crout, Cholesky 분해 알고리즘으로 나눠진다. 먼저 일반적인 형태를 알아보고 각각의 알고리즘을 자세히 살펴보자.

A 가 정칙 정방 행렬이고 행 교환이 필요없다고 하면 아래와 같은 일반적인 형태로 LU 분해가 가능하다.

    \[A = LU = \begin{bmatrix} l_{11} & 0 & \dots & 0 \\ l_{21} & l_{22} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \dots & l_{nn} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n} \\ 0 & u_{22} & \dots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & u_{nn} \end{bmatrix}\]

    \[a_{ij} = \sum_{k=1}^{\text{min}(i, j)} l_{ik} u_{kj}\]

이것을 상삼각 성분과 하삼각 성분으로 나누면

    \[\begin{aligned} a_{ij} &= \sum_{k=1}^i l_{ik} u_{kj}\quad \text{s.t. } i \le j \\ a_{ij} &= \sum_{k=1}^j l_{ik} u_{kj}\quad \text{s.t. } i \ge j \\ \end{aligned}\]

약간 고치면,

(1)   \[\begin{aligned}  a_{ij} &= \sum_{k=1}^{i-1} l_{ik} u_{kj} + l_{ii} u_{ij} \quad \text{s.t. } i \le j \\ u_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}}{l_{ii}} \\ \end{aligned}\]

(2)   \[\begin{aligned}  a_{ij} &= \sum_{k=1}^{j-1} l_{ik} u_{kj} + l_{ij} u_{jj} \quad \text{s.t. } i \ge j \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj}}{u_{jj}} \\ \end{aligned}\]

여기서 L 의 대각성분이 모두 1 이면 Doolittle 분해, U 의 대각성분이 모두 1 이면 Crout 분해라고 한다. 쉽게 생각하면 A = LDU 로 분해가 되는데, 대각행렬 D 가 왼쪽에 붙는지 오른쪽에 붙는지에 따라 알고리즘이 달라진다고 보면 된다. 그리고 만약 A 가 대칭행렬 (또는 hermitian 행렬) 이라면 A = LDL^T 로 분해가 가능하다. 이때 D 의 대각 성분이 양수라면 A = LL^T 형태로 분해가 가능하고 이것을 Cholesky 분해라고 부른다. 세가지 알고리즘 모두 in-place 로 구할 수 있으므로 메모리를 절약할 수 있다.

Doolittle Algorithm

Doolittle 알고리즘은 식 (1), (2) 에서 l_{ii} = 1 이므로,

    \[\begin{aligned} u_{ij} &= a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj} \quad &\text{s.t. } i \le j \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj}}{u_{jj}} \quad &\text{s.t. } i > j \\ \end{aligned}\]

u_{ij} 를 구하기 위해서는 l_{i1}, l_{i2}, \dots, l_{i(i-1)} 성분들이 필요하다. 또 l_{ij} 를 구하려면 u_{1j}, u_{2j}, \dots, u_{jj} 성분들이 필요하다. 그림으로 나타내면,

doolittle-lu-order

파란색으로 표시된 성분을 구하려면 빨간색 성분들이 필요하다. 잘 보면 좌상단부터 차례로 계산해 나갈 수 있다.
따라서 in-place 로 l_{ij}, u_{ij} 들을 계산해 나가려면 좌상단 성분부터 차례로 한 행씩 순서대로 구하면 된다. 의사 코드는 다음과 같다.

한 행씩 계산하므로 row major 행렬에 적합하다.

다음은 LUx = b 를 푸는 의사 코드다. in-place 로 계산되었기 때문에 AL, U 가 모두 포함되어있다고 가정한다.

Crout Algorithm

Crout 알고리즘은 식 (1), (2) 에서 u_{jj} = 1 이므로,

    \[\begin{aligned} u_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}}{l_{ii}} \quad &\text{s.t. } i < j \\ l_{ij} &= a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \quad &\text{s.t. } i \ge j \\ \end{aligned}\]

u_{ij} 를 구하기 위해서는 l_{i1}, l_{i2}, \dots, l_{ii} 성분들이 필요하다. 또 l_{ij} 를 구하려면 u_{1j}, u_{2j}, \dots, u_{j(j-1)} 성분들이 필요하다. 그림으로 나타내면,

crout-lu-order

따라서 in-place 로 l_{ij}, u_{ij} 들을 계산해 나가려면 좌상단 성분부터 한 열씩 구하면 된다. 의사 코드는 다음과 같다.

한 열씩 계산하므로 column major 행렬에 적합하다.

다음은 LUx = b 를 푸는 의사 코드다. in-place 로 계산되었기 때문에 AL, U 가 모두 포함되어있다고 가정한다.

Cholesky Algorithm

Cholesky 알고리즘은 설명했듯이 A 가 대칭행렬인 것을 가정한다. l_{ij} = u_{ji} 이므로 식 (1), (2) 는 다음과 같이 바꿀 수 있다.

    \[\begin{aligned} u_{ij} &= \frac{a_{ij} - \sum_{k=1}^{i-1} u_{ki} u_{kj}}{u_{ii}} \quad &\text{s.t. } i \le j  \\ l_{ij} &= \frac{a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk}}{l_{jj}} \quad &\text{s.t. } i \ge j  \\ \end{aligned}\]

in-place 로 L 을 구할 것인지 U 를 구할 것인지에 따라 두가지 버젼으로 나눌 수 있다. L 을 구할 때는 열 우선으로, U 를 구할 때는 행 우선으로 구 할 수 있다.

먼저 L 을 구하는 Cholesky 알고리즘은 다음과 같다.

AL 성분만 in-place 로 구했기 때문에 대각성분을 제외한 상삼각성분은 그대로 남아있다.

다음은 LL^Tx = b 를 푸는 의사 코드다.

다음으로 U 를 구하는 Cholesky 알고리즘은 다음과 같다.

마찬가지로 AU 성분만 in-place 로 구했기 때문에 대각성분을 제외한 하삼각성분은 그대로 남아있다.

다음은 U^TUx = b 를 푸는 의사 코드다.

Schur / Real Schur Decomposition

Schur Decomposition

(Complex) Schur Decomposition 은 정사각행렬 A \in \mathbb{C}^{n\times n} 를 unitary 행렬 U \in \mathbb{C}^{n\times n} 와 상삼각행렬 T \in \mathbb{C}^{n\times n} 로 분해하는 방법이다. 즉, A = UTU^H 로 분해된다. 이것은 Hermitian 행렬의 고유값 분해 A = U\Lambda U^H 의 일반화된 형태다. 즉, A 가 실수 대칭행렬이거나 복소수 Hermitian 행렬이라면 T 는 대각행렬 \Lambda 가 된다.

행렬 A \in \mathbb{C}^{n\times n} 는 대각화 가능 여부와는 별도로 하나의 고유값 \lambda_{1(A)} 와 그 고유벡터 \mathbf{x}_{1(A)} 를 가진다. 그러면 \mathbf{x}_{1(A)} 를 이용해서 임의의 unitary 행렬 Q_1 을 아래와 같이 구성할 수 있다.

    \[Q_1 =  \begin{bmatrix} \\ \mathbf{x}_{1(A)} & \mathbf{w_2} & \mathbf{w_3} & \dots & \mathbf{w_n} \\ \\ \end{bmatrix}\]

Q_1 을 이용한 A 의 닮음 변환 Q_1^H A Q_1A 의 첫번째 열을 소거한다.

    \[\begin{aligned} Q_1^H A Q_1 &=  \begin{bmatrix} & \mathbf{x}_{1(A)}^H & \\  & \mathbf{w_2}^H & \\  & \mathbf{w_3}^H & \\  & \vdots & \\  & \mathbf{w_n}^H & \\ \end{bmatrix} \begin{bmatrix} \\ A\mathbf{x}_{1(A)} & A\mathbf{w_2} & A\mathbf{w_3} & \dots & A\mathbf{w}_n \\ \\ \end{bmatrix} \\  &= \begin{bmatrix} \lambda_A & \dots \\ 0 & \boxed{A_2}     \end{bmatrix} = T_1 \end{aligned}\]

A_2n-1 크기의 정사각행렬이다. A_2 역시 한개의 고유벡터 \lambda_{1(A_2)} 와 그 고유벡터 \mathbf{x}_{1(A_2)} 를 구해서 Q_2 를 만들 수 있을 것이다. 그럼 T_1 은 unitary 행렬인 \begin{bmatrix} I_{1\times 1} & 0 \\ 0 & Q_2     \end{bmatrix} 을 이용해서 두번째 열을 소거할 수 있다.

    \[\begin{bmatrix} I_{1\times 1} & 0 \\ 0 & Q_2     \end{bmatrix}^H \begin{bmatrix} \lambda_{1(A)} & \dots \\ 0 & \boxed{A_2}     \end{bmatrix} \begin{bmatrix} I_{1\times 1} & 0 \\ 0 & Q_2     \end{bmatrix} =  \begin{bmatrix} \lambda_{1(A)} & \dots \\ 0 & \begin{bmatrix}    \lambda_{1(A_2)} & \dots \\   0 & \boxed{A_3}   \end{bmatrix} \end{bmatrix} = T_2\]

한번 더 반복해보자.

    \[\begin{bmatrix} I_{2\times 2} & 0 \\ 0 & Q_3     \end{bmatrix}^H \begin{bmatrix} \lambda_{1(A)} & \dots \\ 0 & \begin{bmatrix}    \lambda_{1(A_2)} & \dots \\   0 & \boxed{A_3}   \end{bmatrix} \end{bmatrix} \begin{bmatrix} I_{2\times 2} & 0 \\ 0 & Q_3     \end{bmatrix} =  \begin{bmatrix} \lambda_{1(A)} & \dots \\ 0 & \begin{bmatrix}    \lambda_{1(A_2)} & \dots \\   0 & \begin{bmatrix}      \lambda_{1(A_3)} & \dots \\     0 & \boxed{A_4}     \end{bmatrix}   \end{bmatrix} \end{bmatrix} = T_3\]

이렇게 n-1 번 반복하면,

    \[\begin{aligned} U_0^H A U_0 &= T_1 \\ U_1^H T_1 U_1 &= T_2 \\ U_2^H T_2 U_2 &= T_3 \\ &\vdots \\ U_{n-2}^H T_{n-2} U_{n-2} &= T_{n-1} \end{aligned} \quad \begin{aligned} U_k = \begin{bmatrix} I_{k\times k} & 0 \\ 0 & Q_{k+1} \end{bmatrix} \end{aligned}\]

계속 열들이 소거되어 T_{n-1} 은 결국 상삼각행렬이 된다.

이제 A 를 표현하기 위해서 식을 이항하면,

    \[\begin{aligned} A &= U_0 T_1 U_0^H \\ T_1 &= U_1 T_2 U_1^H \\ T_2 &= U_2 T_3 U_2^H \\ &\vdots \\ T_{n-2} &= U_{n-2} T_{n-1} U_{n-2}^H \end{aligned}\]

식 전체를 A 로 나타내면

    \[\begin{aligned} A &= U_0 U_1 U_2 \dots U_{n-2} T_{n-1} U_{n-2}^H \dots U_2^H U_1^H U_0^H \\ &= (U_0 U_1 U_2 \dots U_{n-2}) T_{n-1} (U_0 U_1 U_2 \dots U_{n-2})^H \\ &= UTU^H \end{aligned}\]

unitary 행렬끼리의 곱은 역시 unitary 행렬이므로 A = UTU^H 형태로 분해가 된다는 것을 알 수 있다.

Real Schur Decomposition

Real Schur Decomposition 은 실수 행렬만으로 이뤄진 Schur Decomposition 이라고 보면 된다. 다만 T 의 모양이 상삼각행렬이 아닐 수 있다.

    \[A = QTQ^T,\quad T = \begin{bmatrix} T_{11} & T_{12} & \dots & T_{1m} \\ & T_{22} & \dots & T_{2m} \\ & & \ddots & \vdots \\ & & & T_{mm} \\ \end{bmatrix}\]

T_{ii} 의 크기는 1\times 1 일 수도, 2\times 2 일 수도 있다. 1\times 1 블럭은 실수 고유값이고, 2\times 2 블럭은 켤례쌍 (conjugate pair) 복소수 고유값을 가진다.

Conjugate Pair Lambda

A \in \mathbb{R}^{n\times n} 의 고유값 중 하나가 실수라면 대응되는 고유벡터 역시 실수 벡터다. 반대로 고유값이 복소수라면 대응되는 고유벡터 역시 복소수 벡터가 된다. 또 복소수 고유값을 가진다면 반드시 켤례복소수 고유값도 가진다.

정말 켤례 복소수 고유값을 가지는지 증명해보자. 고유값이 \lambda = \alpha + i\beta, \beta \neq 0 라면 대응되는 고유벡터는 \mathbf{x} = \mathbf{u} + i\mathbf{v} 형태다.

    \[\begin{aligned} A\mathbf{x} &= A (\mathbf{u} + i\mathbf{v}) = A\mathbf{u} + iA\mathbf{v} \\ \lambda \mathbf{x} &= (\alpha + i\beta) (\mathbf{u} + i\mathbf{v}) = (\underbrace{\alpha\mathbf{u} - \beta\mathbf{v}}_{A\mathbf{u}}) + i(\underbrace{\beta\mathbf{u} + \alpha\mathbf{v}}_{A\mathbf{v}}) \end{aligned}\]

즉,

(1)   \[\begin{aligned}  A\mathbf{u} &= \alpha\mathbf{u} - \beta\mathbf{v} \\ A\mathbf{v} &= \beta\mathbf{u} + \alpha\mathbf{v} \end{aligned}\]

이번엔 A\bar{\mathbf{x}} 를 구해보자.

    \[\begin{aligned} A\bar{\mathbf{x}} &= A (\mathbf{u} - i\mathbf{v}) = A\mathbf{u} - iA\mathbf{v} \\ &= (\alpha\mathbf{u} - \beta\mathbf{v}) - i(\beta\mathbf{u} + \alpha\mathbf{v}) \\ &= (\alpha - i\beta)\mathbf{u} - i(\alpha - i\beta)\mathbf{v} \\ &= (\alpha - i\beta)(\mathbf{u} - i\mathbf{v}) = \bar{\lambda}\bar{\mathbf{x}} \end{aligned}\]

A\bar{\mathbf{x}} = \bar{\lambda}\bar{\mathbf{x}} 이므로 \bar{\lambda} = \alpha - i\beta, \beta\neq 0 역시 고유값, \bar{\mathbf{x}} 는 대응되는 고유벡터다.

Conjugate Pair Matrix

이번에는 A 의 켤례쌍 고유값 \lambda = \alpha + i\beta, \bar{\lambda} = \alpha - i\beta 와 같은 고유값을 갖는 닮음행렬 S \in \mathbb{R}^{2\times 2} 을 정의해보자.

식 (1) 을 다시 쓰면,

(2)   \[A \begin{bmatrix} \mathbf{u} & \mathbf{v} \end{bmatrix} =  \begin{bmatrix} \mathbf{u} & \mathbf{v} \end{bmatrix} \begin{bmatrix} \alpha & \beta \\ -\beta & \alpha \end{bmatrix} \]

여기서 \mathbf{u}\mathbf{v} 는 선형독립이다. 만약 \mathbf{u}\mathbf{v} 가 선형종속이라면 그의 선형결합인 고유벡터 \mathbf{x}\bar{\mathbf{x}} 도 선형종속이 되고, 그러면 고유값 \lambda\bar{\lambda} 는 같아질 수 있으므로 \beta = 0 이 된다. 따라서 \mathbf{u}\mathbf{v} 는 선형독립이다.

\begin{bmatrix} \mathbf{u} & \mathbf{v} \end{bmatrix}QR 분해하면,

(3)   \[\begin{bmatrix}  \mathbf{u} & \mathbf{v}  \end{bmatrix} = \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix} R \]

이제 식 (3) 을 식 (2) 에 대입하면,

    \[\begin{aligned} A \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2      \end{bmatrix} R &=  \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix} R \begin{bmatrix} \alpha & \beta \\ -\beta & \alpha \end{bmatrix} \\ A \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2      \end{bmatrix} &=  \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix} R \begin{bmatrix} \alpha & \beta \\ -\beta & \alpha \end{bmatrix} R^{-1} \\ &= \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix} S \end{aligned}\]

S\begin{bmatrix} \alpha & \beta \\ -\beta & \alpha \end{bmatrix} 과 닮음 행렬이다.

    \[\begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix}^T  A \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix} = S\]

Decomposition

Ak 개의 켤례쌍 고유값을 가진다면, 1번째 켤례쌍 고유값과 대응되는 켤례쌍 고유벡터로부터 아래의 직교 행렬을 만들 수 있다.

    \[Q_1 = \begin{bmatrix} \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix} & \mathbf{w}_3 & \dots & \mathbf{w}_n     \end{bmatrix}\]

Q_1 을 이용한 A 의 닮음 변환 Q_1^T A Q 는 다음과 같이 좌상단 2\times 2 블럭 아래를 0 으로 만든다.

    \[\begin{aligned} Q_1^T A Q &= \begin{bmatrix}   \begin{bmatrix}   \mathbf{q}_1 & \mathbf{q}_2   \end{bmatrix}^T \\ \mathbf{w}_3^T \\ \vdots \\ \mathbf{w}_n^T      \end{bmatrix} \begin{bmatrix} A \begin{bmatrix} \mathbf{q}_1 & \mathbf{q}_2 \end{bmatrix} & A\mathbf{w}_3 & \dots & A\mathbf{w}_n  \end{bmatrix} \\ =& \begin{bmatrix} S_1 & \dots \\ 0 & \boxed{A_2} \end{bmatrix} = T_1 \end{aligned}\]

S_1A 의 1번째 켤례쌍 고유값을 갖는 \mathbb{R}^{2\times 2} 행렬이고, A_2(n-2)\times(n-2) 크기의 정사각행렬이다. A_2k-1 개의 켤례쌍 고유값을 가진다. 따라서 Q_2 를 만들 수 있을 것이다. 그럼 T_1 은 직교행렬 \begin{bmatrix} I_{2\times 2} & 0 \\ 0 & Q_2 \end{bmatrix} 를 이용해서 T_2 를 만들 수 있다.

    \[\begin{bmatrix}  I_{2\times 2} & 0 \\  0 & Q_2  \end{bmatrix}^T \begin{bmatrix} S_1 & \dots \\ 0 & \boxed{A_2} \end{bmatrix} \begin{bmatrix}  I_{2\times 2} & 0 \\  0 & Q_2  \end{bmatrix} =  \begin{bmatrix} S_1 & \dots \\ 0 &    \begin{bmatrix}   S_2 & \dots \\   0 & \boxed{A_3}   \end{bmatrix} \end{bmatrix} = T_2\]

만약 실수 고유값을 가진다면 Schur 분해에서 했던대로 좌상단 성분 아래를 0 으로 만든다. 결과적으로 A 는 다음과 같이 분해된다.

    \[A = QTQ^T,\quad T = \begin{bmatrix} T_{11} & T_{12} & \dots & T_{1m} \\ & T_{22} & \dots & T_{2m} \\ & & \ddots & \vdots \\ & & & T_{mm} \\ \end{bmatrix}\]

T_{ii} 는 실수 고유값이거나, 켤례 복소수 고유값을 갖는 2\times 2 행렬이다.

3 ways of QR Decomposition and Hessenberg matrix

QR Decomposition

QR 분해는 세가지 방법이 있다. 가장 기본적인 Gram-Schmidt 방법은 A 의 각 열들을 첫번째 열에 맞추어 차례로 직교화하는 방법을 사용한다. 결과적으로 R 은 Gram-Schmidt 과정에 따라 상삼각행렬이 된다. 나머지 두 방법 (Givens, Householder) 은 Gram-Schmidt 보다 효율적이다. 이 방법들은 가우스 소거 과정에서 나타났던 행 소거행렬 E 대신에 직교행렬 Q 를 이용해 열 소거를 해서 최종적으로 상삼각행렬 R 로 변환한다.

GivensHouseholder 방법은 Q 의 형태가 다르다. 여기서는 Gram-Schmidt 방법은 생략하고 나머지 두 QR 분해 방법들에 대해서 설명한다. Gram-Schmidt 방법에 대해서는 아래의 포스트로 대신한다.

선형대수 note 17: Orthogonal matrices and Gram-Schmidt

Givens Rotation

Givens rotation 행렬은 특정 ij 평면을 기준으로 반시계 방향으로 \theta 만큼 회전시키는 행렬이다. 일반적인 형태는 다음과 같다.

    \[G(i, j, \theta) = \begin{bmatrix} 1 & \dots & 0 & \dots & 0 & \dots & 0 \\ \vdots & \ddots & \vdots & & \vdots & & \vdots\\  0 & \dots & c & \dots & -s & \dots & 0 \\ \vdots & & \vdots & \ddots & \vdots & & \vdots\\  0 & \dots & s & \dots & c & \dots & 0 \\ \vdots & & \vdots & & \vdots & \ddots & \vdots\\  0 & \dots & 0 & \dots & 0 & \dots & 1 \\ \end{bmatrix} \quad \text{s.t. } 1 \leq i < j \leq n\]

    \[\begin{aligned} g_{ii} &= \cos\theta \\ g_{ij} &= -\sin\theta \\ g_{ji} &= \sin\theta \\ g_{jj} &= \cos\theta \\ g_{kl} &= \begin{cases}     0 & \text{if } k \neq l \\     1 & \text{if } k = l  \end{cases} \quad \text{s.t. } k, l \neq i, j \end{aligned}\]

\theta 를 직접 구할 필요없이 \cos\theta, \sin\theta 만 구하면 되므로,

    \[\cos\theta = \pm\frac{\alpha}{\sqrt{\alpha^2 + \beta^2}}\quad \sin\theta = \mp\frac{\beta}{\sqrt{\alpha^2 + \beta^2}}\]

여기서 \alpha 는 회전시켜 소거할 행렬 A 의 성분 a_{ii} 이고, \beta0 으로 만들 성분 a_{ji} 이다.

예를 들어 3\times 3 행렬의 (2, 1) 성분을 0 으로 만들 xy 평면 회전행렬 G_{21} 은 다음과 같다.

    \[G_{21} = G(1, 2, \theta_1) = \begin{bmatrix} \cos\theta_1 & -\sin\theta_1 & 0 \\ \sin\theta_1 & \cos\theta_1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}\]

가우스 소거 행렬 E_{21}(2, 1) 성분을 0 으로 소거했던것 처럼, G_{21}(2, 1) 성분을 0 으로 소거한다.
이어서 G_{21} A(2, 1) 성분을 0 으로 만들 xz 평면 회전행렬 G_{31}G_{31} G_{21} A(3, 2) 성분을 0 으로 만들 yz 평면 회전행렬 G_{32} 는 각각 다음과 같다.

    \[G_{31} = G(1, 3, \theta_2) = \begin{bmatrix} \cos\theta_2 & 0 & -\sin\theta_2 \\ 0 & 1 & 0 \\ \sin\theta_2 & 0 & \cos\theta_2 \\ \end{bmatrix}, \quad G_{32} = G(2, 3, \theta_3) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\theta_3 & -\sin\theta_3 \\ 0 & \sin\theta_3 & \cos\theta_3 \\ \end{bmatrix}\]

이들을 이용해서 3\times 3 행렬인 A 의 열들을 차례로 모두 소거하면

    \[(G_{32} G_{31} G_{21}) A = R\]

G_{ij} 는 직교행렬이다. 따라서 식을 정리하면,

    \[A = (G_{32} G_{31} G_{21})^T R  = QR\]

아래는 직접 짜본 MATLAB 코드다.

Householder Reflection

Householder reflection 행렬 H \in \mathbb{R}^{m\times n} 는 어떤 벡터에 그 벡터를 \mathbf{u} 로 투영한 벡터의 두배 만큼 뺀다. 즉, \mathbf{u} 가 직선의 단위 직교 벡터 (normal) 라면 H 는 임의의 벡터를 직선에 반사시키는 행렬이다. QR 분해에서는 반사시켜서 축에 정렬시킨다.

householder-reflection-graph

행렬 A 의 첫번째 열 \mathbf{a}_1 은 아래와 같이 소거되어 \mathbf{r}_1 을 만든다.

    \[H_1 = I - 2\mathbf{u}_1\mathbf{u}_1^T \qquad H_1 \mathbf{a}_1 =  \begin{bmatrix} ||\mathbf{a}_1|| \\ 0 \\ \vdots \\ 0 \\ \end{bmatrix} \ \text{or}\  \begin{bmatrix} -||\mathbf{a}_1|| \\ 0 \\ \vdots \\ 0 \\ \end{bmatrix} = \mathbf{r}_1\]

\mathbf{u}_1\mathbf{a}_1 - \mathbf{r}_1 방향이다. H_1 A 의 오른쪽 아래 서브 행렬에 대해서 계속 서브 반사 행렬 Q_k = \begin{bmatrix}I_{k-1} & 0 \\ 0 & H_k\end{bmatrix} 를 곱하면 결과적으로 상삼각화된다.

    \[(Q_{n} \dots Q_3 Q_2 Q_1) A = R\]

H_k, Q_k 는 직교행렬이면서 대칭행렬이다. 따라서 식을 정리하면,

    \[\begin{aligned} A &= (Q_1 Q_2 Q_3 \dots Q_{n}) R = QR \\ A_{m\times n} &= Q_{m\times m}R_{m\times n} \end{aligned}\]

전체 과정은 \mathbb{C}^{m\times n} 에서도 성립한다.

아래는 직접 짜본 MATLAB 코드다.

Hessenberg Decomposition

Hessenberg 행렬이란 sub diagonal 또는 super diagonal 성분까지 포함한 삼각행렬을 말한다. 예를 들어 5\times 5 upper Hessenberg 행렬은 다음과 같은 형태다.

    \[H = \begin{bmatrix} \times & \times & \times & \times & \times \\ \times & \times & \times & \times & \times \\ 0 & \times & \times & \times & \times \\ 0 & 0 & \times & \times & \times \\ 0 & 0 & 0 & \times & \times\\ \end{bmatrix}\]

직교행렬 Q 를 이용해서 정방행렬 A 를 상삼각 형태로 변환했던 것 처럼, 직교행렬 Q 를 이용해서 A 를 upper Hessenberg 형태로 만들 수 있다. Householder reflection 행렬 Q_{k+1} = \begin{bmatrix}I_{k} & 0 \\ 0 & H_{k+1}\end{bmatrix} 을 사용할 수 있다. 원리는 Householder reflection 과 같지만 한단계 아래 성분들에 대해서 소거하는 것이다.

    \[\begin{aligned} Q_2 A &= \begin{bmatrix} 1 & 0^T \\  \mathbf{0} & H_2 \\ \end{bmatrix} \begin{bmatrix} a_{11} & \dots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \dots & a_{nn}   \end{bmatrix} \\ &=  \begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ \multicolumn{4}{c}{H_2 \begin{bmatrix}   a_{21} & \dots & a_{2n} \\   \vdots & \ddots & \vdots \\   a_{n1} & \dots & a_{nn}     \end{bmatrix}} \end{bmatrix} \\ &=  \begin{bmatrix} \times & \times & \dots \\  \times & \times & \dots \\ 0 & \times & \dots \\ \vdots & \vdots & \ddots \\ 0 & \times & \dots \end{bmatrix} \end{aligned}\]

Q_{k+1}k 번째 열을 upper Hessenberg form 으로 만든다. n-2 번 적용하면,

    \[(Q_{n-1} \dots Q_3 Q_2) A = H\]

여기서 H 는 upper Hessenberg 행렬이다.

Similar Hessenberg matrix

Q_{k+1}^T = Q_{k+1} 을 오른쪽에 곱하면 왼쪽 서브 행렬 (n\times k) 을 유지한다. 예를 들어 Q_2 A Q_2 는 결과 행렬의 형태 (왼쪽 서브 행렬 (n\times 1) 에 대해서만 upper Hessenberg form) 를 유지하는 닮음 변환을 나타낸다.

    \[\begin{aligned} Q_2 A Q_2 &= \begin{bmatrix} 1 & 0^T \\  \mathbf{0} & H_2 \\ \end{bmatrix} \begin{bmatrix} a_{11} & \dots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \dots & a_{nn}   \end{bmatrix}  \begin{bmatrix} 1 & 0^T \\  \mathbf{0} & H_2 \\ \end{bmatrix} \\ &=  \begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ H_2 \begin{bmatrix}   a_{21} \\ \vdots \\ a_{n1}   \end{bmatrix} & \multicolumn{3}{c}{H_2 \begin{bmatrix}   a_{22} & \dots & a_{2n} \\   \vdots & \ddots & \vdots \\   a_{n2} & \dots & a_{nn}     \end{bmatrix}} \end{bmatrix}  \begin{bmatrix} 1 & 0^T \\  \mathbf{0} & H_2 \\ \end{bmatrix} \\ &= \begin{bmatrix} a_{11} & \begin{bmatrix}a_{12} & \dots & a_{1n}\end{bmatrix} H_2 \\ H_2 \begin{bmatrix}   a_{21} \\ \vdots \\ a_{n1}   \end{bmatrix} & \multicolumn{3}{c}{H_2 \begin{bmatrix}   a_{22} & \dots & a_{2n} \\   \vdots & \ddots & \vdots \\   a_{n2} & \dots & a_{nn}   \end{bmatrix} H_2 } \end{bmatrix}  =  \begin{bmatrix} \times & \times & \dots \\  \times & \times & \dots \\ 0 & \times & \dots \\ \vdots & \vdots & \ddots \\ 0 & \times & \dots \end{bmatrix} \end{aligned}\]

따라서

    \[(Q_{n-1} \dots Q_3 Q_2) A (Q_2 Q_3 \dots Q_{n-1}) = H\]

upper Hessenberg 행렬 H 는 행렬 A 와 닮음이다. 만약 A 가 대칭행렬이라면 결과는 대칭 삼중 대각 행렬 (symmetric tridiagonal matrix) 이 된다는 것을 알 수 있다.

아래는 직접 짜본 MATLAB 코드다.