LU Decomposition Algorithm

Posted on 2016-01-202017-04-07 by Lee Ju Hyung

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}$ 들을 계산해 나가려면 좌상단 성분부터 차례로 한 행씩 순서대로 구하면 된다. 의사 코드는 다음과 같다.

decomp_doolittle(A)
  for i = 1, ..., n
    % Compute a line of L
    for j = 1, ..., i - 1
      sum = 0
      for k = 1, ..., j - 1
        sum += a(i, k) * a(k, j)
      end
      a(i, j) = (a(i, j) - sum) / a(j, j)
    end
    % Compute a line of U
    for j = i, ..., n
      sum = 0
      for k = 1, ..., i - 1
        sum += a(i, k) * a(k, j)
      end
      a(i, j) = a(i, j) - sum
    end
  end
end

decomp_doolittle(A)

for i = 1, ..., n

% Compute a line of L

for j = 1, ..., i - 1

sum = 0

for k = 1, ..., j - 1

sum += a(i, k) * a(k, j)

end

a(i, j) = (a(i, j) - sum) / a(j, j)

end

% Compute a line of U

for j = i, ..., n

sum = 0

for k = 1, ..., i - 1

sum += a(i, k) * a(k, j)

end

a(i, j) = a(i, j) - sum

end

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

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

solve_doolittle(A, b)
  % Solve Ly = b
  for i = 1, ..., n
    sum = 0
    for j = 1, ..., i - 1
      sum += a(i, j) * y(j)
    end
    y(i) = b(i) - sum
  end
  % Solve Ux = y
  for i = n, ..., 1
    sum = 0
    for j = i + 1, ..., n
      sum += a(i, j) * x(j)
    end
    x(i) = (y(i) - sum) / a(i, i)
  end
end

solve_doolittle(A, b)

% Solve Ly = b

for i = 1, ..., n

sum = 0

for j = 1, ..., i - 1

sum += a(i, j) * y(j)

end

y(i) = b(i) - sum

end

% Solve Ux = y

for i = n, ..., 1

sum = 0

for j = i + 1, ..., n

sum += a(i, j) * x(j)

end

x(i) = (y(i) - sum) / a(i, i)

end

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}$ 들을 계산해 나가려면 좌상단 성분부터 한 열씩 구하면 된다. 의사 코드는 다음과 같다.

decomp_crout(A)
  for j = 1, ..., n
    % Compute a column of U
    for i = 1, ..., j - 1
      sum = 0
      for k = 1, ..., i - 1
        sum += a(i, k) * a(k, j)
      end
      a(i, j) = (a(i, j) - sum) / a(i, i)
    end
    % Compute a column of L
    for i = j, ..., n
      sum = 0
      for k = 1, ..., j - 1
        sum += a(i, k) * a(k, j)
      end
      a(i, j) = a(i, j) - sum
    end 
  end
end

decomp_crout(A)

for j = 1, ..., n

% Compute a column of U

for i = 1, ..., j - 1

sum = 0

for k = 1, ..., i - 1

sum += a(i, k) * a(k, j)

end

a(i, j) = (a(i, j) - sum) / a(i, i)

end

% Compute a column of L

for i = j, ..., n

sum = 0

for k = 1, ..., j - 1

sum += a(i, k) * a(k, j)

end

a(i, j) = a(i, j) - sum

end

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

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

solve_crout(A, b)
  % Solve Ly = b
  for i = 1, ..., n
    sum = 0
    for j = 1, ..., i - 1
      sum += a(i, j) * y(j)
    end
    y(i) = (b(i) - sum) / a(i, i)
  end
  % Solve Ux = y
  for i = n, ..., 1
    sum = 0
    for j = i + 1, ..., n
      sum += a(i, j) * x(j)
    end
    x(i) = y(i) - sum
  end
end

solve_crout(A, b)

% Solve Ly = b

for i = 1, ..., n

sum = 0

for j = 1, ..., i - 1

sum += a(i, j) * y(j)

end

y(i) = (b(i) - sum) / a(i, i)

end

% Solve Ux = y

for i = n, ..., 1

sum = 0

for j = i + 1, ..., n

sum += a(i, j) * x(j)

end

x(i) = y(i) - sum

end

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 알고리즘은 다음과 같다.

decomp_lcholesky(A)
  for j = 1, ..., n
    % Compute a diagonal element of L
    sum = 0
    for k = 1, ..., j - 1
      sum += a(j, k) * a(j, k)
    end
    a(j, j) = sqrt(a(j, j) - sum) 
    % Compute a column which is below the diagonal of L
    for i = j + 1, ..., n
      sum = 0
      for k = 1, ..., j - 1
        sum += a(i, k) * a(j, k)
      end
      a(i, j) = (a(i, j) - sum) / a(j, j)
    end
  end
end

decomp_lcholesky(A)

for j = 1, ..., n

% Compute a diagonal element of L

sum = 0

for k = 1, ..., j - 1

sum += a(j, k) * a(j, k)

end

a(j, j) = sqrt(a(j, j) - sum)

% Compute a column which is below the diagonal of L

for i = j + 1, ..., n

sum = 0

for k = 1, ..., j - 1

sum += a(i, k) * a(j, k)

end

a(i, j) = (a(i, j) - sum) / a(j, j)

end

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

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

solve_lcholesky(A, b)
  % Solve Ly = b
  for i = 1, ..., n
    sum = 0
    for j = 1, ..., i - 1
      sum += a(i, j) * y(j)
    end
    y(i) = (b(i) - sum) / a(i, i)
  end
  % Solve L^Tx = y
  for i = n, ..., 1
    sum = 0
    for j = i + 1, ..., n
      sum += a(j, i) * x(j)
    end
    x(i) = (y(i) - sum) / a(i, i)
  end
end

solve_lcholesky(A, b)

% Solve Ly = b

for i = 1, ..., n

sum = 0

for j = 1, ..., i - 1

sum += a(i, j) * y(j)

end

y(i) = (b(i) - sum) / a(i, i)

end

% Solve L^Tx = y

for i = n, ..., 1

sum = 0

for j = i + 1, ..., n

sum += a(j, i) * x(j)

end

x(i) = (y(i) - sum) / a(i, i)

end

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

decomp_ucholesky(A)
  for i = 1, ..., n
    % Compute a diagonal element of U
    sum = 0
    for k = 1, ..., i - 1
      sum += a(k, i) * a(k, i)
    end
    a(i, i) = sqrt(a(i, i) - sum)
    % Compute a row which is right the diagonal of U
    for j = i + 1, ..., n
      sum = 0
      for k = 1, ..., i - 1
        sum += a(k, i) * a(k, j)
      end
      a(i, j) = (a(i j) - sum) / a(i, i)
  end
end

decomp_ucholesky(A)

for i = 1, ..., n

% Compute a diagonal element of U

sum = 0

for k = 1, ..., i - 1

sum += a(k, i) * a(k, i)

end

a(i, i) = sqrt(a(i, i) - sum)

% Compute a row which is right the diagonal of U

for j = i + 1, ..., n

sum = 0

for k = 1, ..., i - 1

sum += a(k, i) * a(k, j)

end

a(i, j) = (a(i j) - sum) / a(i, i)

end

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

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

solve_ucholesky(A, b)
  % Solve U^Ty = b
  for i = 1, ..., n
    sum = 0
    for j = 1, ..., i - 1
      sum += a(j, i) * y(j)
    end
    y(i) = (b(i) - sum) / a(i, i)
  end
  % Solve Ux = y
  for i = n, ..., 1
    sum = 0
    for j = i + 1, ..., n
      sum += a(i, j) * x(j)
    end
    x(i) = (y(i) - sum) / a(i, i)
  end
end

solve_ucholesky(A, b)

% Solve U^Ty = b

for i = 1, ..., n

sum = 0

for j = 1, ..., i - 1

sum += a(j, i) * y(j)

end

y(i) = (b(i) - sum) / a(i, i)

end

% Solve Ux = y

for i = n, ..., 1

sum = 0

for j = i + 1, ..., n

sum += a(i, j) * x(j)

end

x(i) = (y(i) - sum) / a(i, i)

end

Schur / Real Schur Decomposition

Posted on 2016-01-142016-01-16 by Lee Ju Hyung

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_1$ 은 $A$ 의 첫번째 열을 소거한다.

$\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_2$ 는 $n-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

$A$ 가 $k$ 개의 켤례쌍 고유값을 가진다면, 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_1$ 은 $A$ 의 1번째 켤례쌍 고유값을 갖는 $\mathbb{R}^{2\times 2}$ 행렬이고, $A_2$ 는 $(n-2)\times(n-2)$ 크기의 정사각행렬이다. $A_2$ 는 $k-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

Posted on 2016-01-062016-01-17 by Lee Ju Hyung

QR Decomposition

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

Givens 과 Householder 방법은 $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}$ 이고, $\beta$ 는 $0$ 으로 만들 성분 $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 코드다.

function [c, s] = givens (A, i, j)
  a = A(i, i);
  b = A(j, i);
  % check if zeroing element is already zero
  if (b == 0)
    c = 1.0;
    s = 0.0;
  else
    % divide by bigger number for the stability
    if (abs(b) > abs(a))
      r = a / b;
      s = 1.0 / sqrt(1 + r^2);
      c = -s * r;
    else
      r = b / a;
      c = 1.0 / sqrt(1 + r^2);
      s = -c * r;
    end
  end
endfunction

function [Q, R] = qr_givens (A)
  [m, n] = size(A);
  R = A;
  Q = eye(m);
  for i = 1:n
    for j = i+1:m
      [c, s] = givens(R, i, j);
      G = eye(m);
      G(i, i) = c;
      G(i, j) = s;
      G(j, i) = -s;
      G(j, j) = c;
      
      Q = G * Q;
      R = G * R;
    end
  end
  Q = conj(transpose(Q));
endfunction

function [c, s] = givens (A, i, j)

a = A(i, i);

b = A(j, i);

% check if zeroing element is already zero

if (b == 0)

c = 1.0;

s = 0.0;

else

% divide by bigger number for the stability

if (abs(b) > abs(a))

r = a / b;

s = 1.0 / sqrt(1 + r^2);

c = -s * r;

else

r = b / a;

c = 1.0 / sqrt(1 + r^2);

s = -c * r;

end

endfunction

function [Q, R] = qr_givens (A)

[m, n] = size(A);

R = A;

Q = eye(m);

for i = 1:n

for j = i+1:m

[c, s] = givens(R, i, j);

G = eye(m);

G(i, i) = c;

G(i, j) = s;

G(j, i) = -s;

G(j, j) = c;

Q = G * Q;

R = G * R;

end

Q = conj(transpose(Q));

endfunction

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 코드다.

function [H] = householder (a)
  m = rows(a);
  H = eye(m);
  an = conj(transpose(a)) * a; % squared norm of a
  if (an > a(1) * a(1)) % check if column of a is already eliminated
    % numerically stable with large |u|
    u = a;
    if (sign(a(1)) > 0)
      u(1) += sqrt(an);
    else
      u(1) -= sqrt(an);
    end
    H -= 2 * (u * u') / (u' * u);
  end
endfunction

function [Q, R] = qr_householder (A)
  [m, n] = size(A);
  Q = eye(m);
  R = A;
  for i = 1:min(m, n)
    H = householder(R(i:m, i));
    R(i:m, i:n) = H * R(i:m, i:n);
    Q(i:m, :) = H * Q(i:m, :);
  end
  Q = conj(transpose(Q));
endfunction

function [H] = householder (a)

m = rows(a);

H = eye(m);

an = conj(transpose(a)) * a; % squared norm of a

if (an > a(1) * a(1)) % check if column of a is already eliminated

% numerically stable with large |u|

u = a;

if (sign(a(1)) > 0)

u(1) += sqrt(an);

else

u(1) -= sqrt(an);

end

H -= 2 * (u * u') / (u' * u);

end

endfunction

function [Q, R] = qr_householder (A)

[m, n] = size(A);

Q = eye(m);

R = A;

for i = 1:min(m, n)

H = householder(R(i:m, i));

R(i:m, i:n) = H * R(i:m, i:n);

Q(i:m, :) = H * Q(i:m, :);

end

Q = conj(transpose(Q));

endfunction

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 코드다.

function [T, Q] = hessenberg_similarity_transform (A)
  m = rows(A);
  T = A;
  Q = eye(m);
  for i = 2:m - 1
    H = householder(T(i:m, i - 1));
    T(i:m, i - 1:m) = H * T(i:m, i - 1:m);
    T(:, i:m) = T(:, i:m) * H;
    Q(:, i:m) = Q(:, i:m) * H;
  end
endfunction

function [T, Q] = hessenberg_similarity_transform (A)

m = rows(A);

T = A;

Q = eye(m);

for i = 2:m - 1

H = householder(T(i:m, i - 1));

T(i:m, i - 1:m) = H * T(i:m, i - 1:m);

T(:, i:m) = T(:, i:m) * H;

Q(:, i:m) = Q(:, i:m) * H;

end

endfunction

JUDIS.me

notes-to-self

Linear Algebra

LU Decomposition Algorithm