요약

가우스 소거를 할 때 “몇 배를 뺐는지”(곱수)를 버리지 않고 아래쪽 삼각 모양의 행렬에 적어 두면, 원래 행렬이 그 행렬과 소거 결과(위쪽 삼각 모양)의 곱으로 쪼개진다. 한 번 쪼개 두면 방정식의 우변이 바뀔 때마다 비싼 소거를 다시 하지 않고, 삼각행렬 두 개로 빠르게 푼다. 같은 회로에 입력만 바꿔 여러 번 시뮬레이션하는 경우처럼 반복 풀이가 많을 때 특히 이득이다. 다만 피벗 자리에 0이 오면 행을 바꿔야 해서, 일반적으로는 행 순서를 바꾸는 행렬까지 붙은 꼴이 된다.

예시로 보기

가우스 소거법의 예 \(A = \begin{pmatrix}1 & 2 & 1\\ 3 & 8 & 1\\ 0 & 4 & 1\end{pmatrix}\)을 소거할 때 곱수는 \(\ell_{21} = 3\)(2행 − 3 × 1행), \(\ell_{31} = 0\), \(\ell_{32} = 2\)(3행 − 2 × 2행)였다. 곱수를 제자리에 적으면

\[L = \begin{pmatrix}1 & 0 & 0\\ 3 & 1 & 0\\ 0 & 2 & 1\end{pmatrix}, \qquad U = \begin{pmatrix}1 & 2 & 1\\ 0 & 2 & -2\\ 0 & 0 & 5\end{pmatrix}, \qquad LU = A.\]

\(\mathbf{b} = (2, 12, 2)\)를 풀 때는 먼저 \(L\mathbf{c} = \mathbf{b}\)에서 위부터 \(c_1 = 2\), \(c_2 = 12 - 3 \cdot 2 = 6\), \(c_3 = 2 - 0 - 2 \cdot 6 = -10\)을 얻는다. 이것은 소거가 \(\mathbf{b}\)에 한 일과 같다. 이어서 \(U\mathbf{x} = \mathbf{c}\)를 아래부터 풀면 \(\mathbf{x} = (2, 1, -2)\)다. 곱수가 아래 정의의 \(L\)의 성분이다.

정의

입력: \(n \times n\) 행렬 \(A\). 출력: 순열 \(P\), 대각이 1인 아래삼각행렬 \(L\), 위삼각행렬 \(U\)로 \(PA = LU\).

LU(A)                               # 인덱스 1부터
  U ← A, L ← 0, P ← I
  for k = 1 to n
      p ← k..n 중 |U[i][k]|가 가장 큰 행        # 부분 피벗팅
      U, P, 그리고 L의 1..k−1열에서 k행과 p행을 바꾼다
      for i = k+1 to n
          L[i][k] ← U[i][k] / U[k][k]          # 곱수를 기록
          U[i] ← U[i] − L[i][k]·U[k]
  L의 대각을 1로
  return P, L, U

SOLVE(P, L, U, b)
  L c = P b 를 위에서부터 (전진 대입)
  U x = c   를 아래에서부터 (후진 대입)
곱수가 그대로 \(L\)이 된다

행 바꾸기가 필요 없으면(피벗이 모두 0이 아니면) \(A = LU\)이고, \(L\)의 \((i, k)\) 성분은 소거에서 쓴 곱수 \(\ell_{ik}\)이다1.

증명 스케치
  1. 소거 = 기본 행렬 곱: “\(i\)행 \(- \ell_{ik} \times k\)행”은 단위행렬의 \((i, k)\) 자리에 \(-\ell_{ik}\)를 넣은 \(E_{ik}\)를 왼쪽에 곱하는 것이다. 소거 전체는 \((\cdots E_{32}E_{31}E_{21})A = U\)다.
  2. 되돌리기: \(E_{ik}^{-1}\)은 같은 자리에 \(+\ell_{ik}\)를 넣은 행렬이다(빼기를 되돌리는 더하기). 그래서 \(A = (E_{21}^{-1}E_{31}^{-1}E_{32}^{-1}\cdots)U\).
  3. 곱수가 섞이지 않는다: 역행렬들을 이 순서로 곱하면 각 \(\ell_{ik}\)가 서로 곱해지지 않고 제자리에 그대로 놓인다. 뒤쪽 행렬의 1이 아닌 성분은 앞쪽 행렬의 1이 아닌 성분과 만나지 않기 때문이다. 그 곱이 \(L\)이다. ∎

비용. 분해는 소거와 같은 약 \(\frac23 n^3\)번의 연산이다. 한 번의 풀이(전진 + 후진 대입)는 약 \(2n^2\)번이다. 같은 \(A\)로 우변 \(k\)개를 풀면 \(\frac23 n^3 + 2kn^2\)번이라, 매번 소거하는 \(k \cdot \frac23 n^3\)보다 훨씬 싸다.

flowchart LR
    A["행렬 A"] --> F["분해 PA = LU, 한 번만, 약 2/3 n³"]
    B["새 우변 b"] --> PB["행 순서 바꾸기 Pb"]
    F -->|"P"| PB
    PB --> FW["전진 대입 Lc = Pb, 위에서부터"]
    F -->|"L"| FW
    FW --> BW["후진 대입 Ux = c, 아래에서부터"]
    F -->|"U"| BW
    BW --> X["해 x"]

윗줄의 분해는 행렬 A 하나에 한 번만 한다. 새 우변 b가 올 때마다 아랫줄의 두 대입만 다시 돌고, 둘을 합쳐 약 2n²번이다2.

예제

행 바꾸기가 필요한 경우. \(A = \begin{pmatrix}0 & 1\\ 1 & 1\end{pmatrix}\)은 첫 피벗 자리가 0이라 \(A = LU\)로 쓸 수 없다.

  1. 피벗 찾기: 첫 열에서 0이 아닌 수가 있는 2행을 올린다. \(P = \begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix}\).
  2. 분해: \(PA = \begin{pmatrix}1 & 1\\ 0 & 1\end{pmatrix}\)이 이미 위삼각이라 \(L = I\), \(U = PA\).
  3. 풀이: \(\mathbf{b} = (1, 3)\)이면 \(P\mathbf{b} = (3, 1)\). \(L\mathbf{c} = (3, 1)\)에서 \(\mathbf{c} = (3, 1)\), \(U\mathbf{x} = \mathbf{c}\)에서 \(y = 1\), \(x = 2\).
검증: 예시의 \(L\), \(U\), \(\mathbf{c} = (2, 6, -10)\), \(\mathbf{x} = (2, 1, -2)\), 행 바꾸기 예제, 무작위 행렬 200개에서 \(PA = LU\)와 곱수 크기 1 이하, 풀이의 잔차, 우변이 여러 개일 때 연산 수 비교, 카드의 값 — 08_lu-decomposition_impl.py, 08_lu-decomposition_verify.py

활용

  • 라이브러리. LAPACK의 getrf(부분 피벗팅 LU)와 getrs(풀이)가 표준이고, SciPy의 scipy.linalg.lu_factor, lu_solve가 이것을 부른다3. numpy.linalg.solve도 내부에서 LU를 쓴다.
  • 반복 풀이. 회로 시뮬레이션, 유한요소 해석, 뉴턴 방법의 반복처럼 같은(또는 거의 같은) 행렬로 여러 번 푸는 곳.
  • 행렬식. \(\det A = \pm u_{11}u_{22}\cdots u_{nn}\)(\(\pm\)는 행 바꾼 횟수의 홀짝)라, 행렬식도 분해의 부산물로 \(O(n^3)\)에 얻는다(행렬식).

연결

과목별 관점

수치해석 (2-2학기). 같은 분해를 소거 대신 공식으로 직접 계산한다. 대각이 1인 \(L\)을 쓰는 방식을 둘리틀 방법이라 하고, \(L\)과 \(U\)를 한 행렬에 겹쳐 저장한다(\(L\)의 대각 1은 적지 않는다). 열 \(j\)마다 위에서 아래로, 먼저 \(U\)를 채우고 그 아래 \(L\)을 채운다.

\[U_{ij} = M_{ij} - \sum_{k=1}^{i-1}L_{ik}U_{kj}\ (i \le j), \qquad L_{ij} = \frac{1}{U_{jj}}\left(M_{ij} - \sum_{k=1}^{j-1}L_{ik}U_{kj}\right)\ (i > j)\]

예를 들어 \(M = \begin{pmatrix}6 & 4 & 2\\ 3 & -2 & -1\\ 3 & 4 & 1\end{pmatrix}\), \(L = \begin{pmatrix}1 & 0 & 0\\ 1/2 & 1 & 0\\ 1/2 & -1/2 & 1\end{pmatrix}\), \(U = \begin{pmatrix}6 & 4 & 2\\ 0 & -4 & -2\\ 0 & 0 & -1\end{pmatrix}\)로 나뉜다. 삼각 연립방정식은 전진 대입 \(x_i = \frac{1}{L_{ii}}\big(r_i - \sum_{k<i}L_{ik}x_k\big)\)과 후진 대입 \(x_i = \frac{1}{U_{ii}}\big(r_i - \sum_{k>i}U_{ik}x_k\big)\)으로 푼다.

역행렬. \(M^{-1}\)의 \(j\)번째 열은 \(LU\mathbf d_j = \mathbf e_j\)(\(\mathbf e_j\)는 \(j\)번째 칸만 1)의 해다. 한 번 분해해 두고 우변 \(n\)개를 차례로 푼다. 예의 \(M\)이면 \(M^{-1} = \begin{pmatrix}1/12 & 1/6 & 0\\ -1/4 & 0 & 1/2\\ 3/4 & -1/2 & -1\end{pmatrix}\)이다4.

오차 줄이기(반복 개선). 반올림 오차 때문에 컴퓨터로 구한 해 \(\mathbf x_0\)은 참값 \(\mathbf x\)와 조금 다르다. \(\mathbf x_0 = \mathbf x + \Delta\mathbf x\)로 두면 \(M\Delta\mathbf x = M\mathbf x_0 - \mathbf r\)이다. 우변은 계산할 수 있으니 이 식을 같은 \(LU\)로 풀어 \(\mathbf x = \mathbf x_0 - \Delta\mathbf x\)로 고친다. 잔차 \(M\mathbf x_0 - \mathbf r\)을 더 정밀하게 계산할 때 효과가 크다. 검증 코드에서 조건이 나쁜 \(6 \times 6\) 힐베르트 행렬의 해가 한 번의 개선으로 참값에 더 가까워졌다4.

확인 문제

C1 \(A = \begin{pmatrix}2 & 1\\ 6 & 8\end{pmatrix}\)를 \(LU\)로 분해하라(행 바꾸기 없이).

답: 곱수 \(\ell_{21} = 6/2 = 3\). 2행 \(-\) 3 × 1행 \(= (0, 5)\). \(L = \begin{pmatrix}1 & 0\\ 3 & 1\end{pmatrix}\), \(U = \begin{pmatrix}2 & 1\\ 0 & 5\end{pmatrix}\).

C2 아래 코드가 하는 일을 쉬운 말로 설명하라. L은 대각이 1인 아래삼각행렬이다. 
for i in range(n):
    c[i] = b[i] - sum(L[i][j] * c[j] for j in range(i))

답: \(L\mathbf{c} = \mathbf{b}\)를 위에서부터 푸는 전진 대입이다. 이미 구한 \(c_0, \dots, c_{i-1}\)에 곱수를 곱해 빼면 \(c_i\)가 나온다. 소거가 우변 \(\mathbf{b}\)에 했던 일을 다시 하는 것과 같다.

C3 같은 \(A\)로 우변 100개를 풀 때 LU 분해를 먼저 해 두는 것이 유리한 이유를 연산 수로 설명하라.

답: 분해는 한 번 \(\frac23 n^3\), 풀이는 한 번에 \(2n^2\)이라 합이 \(\frac23 n^3 + 200n^2\)이다. 매번 소거하면 \(100 \times \frac23 n^3\)이다. \(n = 1000\)이면 약 77배, \(n\)이 더 커지면 100배에 가까워진다.

C4 위의 \(L\), \(U\)로 \(M\mathbf x = (12, 0, 8)\)을 풀라.

답: 전진 대입 \(L\mathbf y = (12, 0, 8)\): \(y_1 = 12\), \(y_2 = 0 - \frac12\cdot12 = -6\), \(y_3 = 8 - \frac12\cdot12 + \frac12\cdot(-6) = -1\). 후진 대입 \(U\mathbf x = \mathbf y\): \(x_3 = 1\), \(-4x_2 - 2 = -6\)이라 \(x_2 = 1\), \(6x_1 + 4 + 2 = 12\)라 \(x_1 = 1\). \(\mathbf x = (1, 1, 1)\).

  1. Strang, Introduction to Linear Algebra 5판, 2.6절 “Elimination = Factorization: A = LU”(곱수가 \(L\)에 그대로 놓이는 이유, 연산 수), 2.7절 “Transposes and Permutations”(\(PA = LU\)). ↩

  2. 이 문서 정의의 LU·SOLVE 의사코드와 비용 문단을 옮겼다(Strang 5판 2.6~2.7절). ↩

  3. SciPy 문서는 lu_factor가 LAPACK getrf를 쓴다고 밝힌다. 행렬식과 \(U\)의 대각 원소 곱의 관계는 08_lu-decomposition_verify.py에서 확인했다. ↩

  4. 해 \((1, 1, 1)\), \(M^{-1}\), 힐베르트 행렬 실험은 08_lu-decomposition_verify.py로 계산했다. ↩ ↩2