선형 연립방정식 - 가우스 소거와 삼중대각

이 장에서 배우는 것

앞 장에서 실수 계산의 정밀도와 오차를 살펴보았다. 이제 계산 결과를 얻는 절차와 그 결과를 확인하는 절차를 함께 설계한다. 선형 연립방정식(linear system)은 여러 미지수가 서로 영향을 주는 계산에서 반복해서 나타난다. 금속판의 온도를 구할 때도 한 지점의 온도가 이웃 지점의 온도와 연결되므로, 여러 식을 한꺼번에 만족하는 값을 찾아야 한다.

이 장에서는 행렬의 구조에 따라 풀이 방법을 고른다. 일반적인 작은 행렬에는 부분 피벗을 사용하는 가우스 소거를 적용하고, 같은 행렬로 여러 문제를 풀 때는 분해 결과를 재사용한다. 이웃한 미지수끼리만 연결된 경우에는 삼중대각 구조를 이용한다. 마지막에는 구한 해를 원래 식에 대입하여 계산이 실제로 무엇을 만족하는지 확인한다.

• 부분 피벗 가우스 소거에서 행 교환과 소거의 역할을 설명한다.
• LU 분해 결과를 저장하고 서로 다른 우변에 재사용한다.
• 삼중대각 행렬을 세 배열로 표현하고 토머스 알고리즘으로 푼다.
• 원래 행렬과 우변으로 잔차를 계산하고 그 의미와 한계를 구분한다.

문제 상황

2차원 금속판 시뮬레이터를 개발하면서 먼저 작은 온도 계산 블록을 점검한다고 하자. 전체 격자를 다루기 전에, 판 위의 한 줄에서 선택한 세 지점의 온도를 구한다. 주변에서 주어진 온도와 열원은 이미 우변에 모았으며, 세 지점 사이의 연결은 다음 식으로 정리되어 있다. 이 식을 만드는 공간 이산화 과정은 여기서 다루지 않는다. 이번에는 주어진 식을 정확하게 구현하고 검증하는 데 집중한다.

 4 T1 -   T2        = 30
-  T1 + 4 T2 -   T3 = 10
        - T2 + 3 T3 = 50

계수는 온도 사이의 연결 강도를 나타낸다. 우변은 알려진 주변 조건과 공급 항을 합친 값이다. 마지막 대각 성분이 다른 것은 그 지점의 연결 조건이 다르기 때문이다. 이 예제의 수치는 풀이 절차를 손으로 확인하기 쉽게 정한 값이며, 특정 금속의 물성치를 나타내지는 않는다.

미지수 순서를 T1, T2, T3로 두면 행렬 A, 해 벡터 x, 우변 b를 사용하여 Ax = b로 쓸 수 있다. 정답은 10, 10, 20이다. 예를 들어 두 번째 식에는 -10 + 40 - 20 = 10이 들어간다. 이렇게 해를 알고 있는 작은 문제는 새 풀이 함수를 점검하는 기준이 된다.

실무에서는 같은 연결 구조에 여러 경계 조건을 적용하기도 한다. A는 그대로 두고 b만 바뀌는 경우다. 매번 처음부터 행렬을 소거하면 같은 계산을 반복한다. 반대로 행렬이 삼중대각인데도 모든 성분을 저장하면 대부분이 0인 공간과 연산에 비용을 쓴다. 계산기를 만들 때는 해를 얻는 것뿐 아니라, 무엇을 저장하고 무엇을 재사용할지도 결정해야 한다.

행렬의 구조와 반복 작업에 따라 풀이 방법을 선택한다

상황	선택	저장할 정보	주요 확인 사항
작은 일반 행렬	부분 피벗 가우스 소거	행렬과 우변	피벗 선택과 행 교환
같은 행렬, 여러 우변	부분 피벗 LU 분해	분해 행렬과 교환 기록	새 우변에도 교환 적용
삼중대각 행렬	토머스 알고리즘	세 대각선과 우변	소거 중 분모의 크기

부분 피벗 가우스 소거

가우스 소거(Gaussian elimination)는 아래쪽 계수를 차례로 없애서 상삼각 형태를 만드는 방법이다. k번째 단계에서는 k열의 대각 아래 성분을 0으로 만든다. i번째 행에서 k번째 행의 m배를 빼며, 배수는 m = a(i,k) / a(k,k)다. 같은 연산을 우변에도 적용해야 원래 연립방정식과 같은 해를 유지한다.

소거가 끝나면 마지막 식에는 마지막 미지수만 남는다. 그 값을 먼저 구하고 바로 위 식으로 올라간다. 이미 구한 값을 우변에서 빼고 대각 성분으로 나누는 이 과정을 후진 대입(back substitution)이라고 한다. 구현에서는 행 번호를 n부터 1까지 감소시키는 반복문이 자연스럽다.

문제는 나눗셈의 분모인 피벗(pivot)이다. 대각 성분이 0이면 그 자리에서 나눌 수 없다. 0이 아니더라도 아래쪽 성분보다 매우 작으면 소거 배수가 커져 반올림 오차가 증폭될 수 있다. 부분 피벗(partial pivoting)은 현재 열의 아직 처리하지 않은 행 중 절댓값이 가장 큰 성분을 찾아 대각 위치로 옮긴다.

완성 코드에서는 원래 온도 식의 첫째 행과 둘째 행을 바꾼 행렬을 일반 풀이기에 넣는다. 첫 열은 -1, 4, 0이 된다. 첫 피벗으로 4가 있는 둘째 행을 선택하면 아래 행의 소거 배수는 -1/4이 된다. 행 교환이 실제로 실행되는 사례를 넣어, 교환 기록과 우변 처리를 함께 확인한다.

행 교환은 식의 순서를 바꾸는 연산이다. 미지수의 순서는 바뀌지 않으므로 최종 해를 다시 재배열할 필요가 없다. 열을 바꾸는 방식이라면 미지수 순서까지 관리해야 하지만, 이 장의 알고리즘은 행만 교환한다.

부분 피벗은 널리 쓰이는 안정화 방법이지만 모든 행렬에서 오차를 작게 보장하지는 않는다. 서로 거의 종속인 식에서는 입력의 작은 변화가 해의 큰 변화로 이어질 수 있다. 또한 소거 과정에서 성분이 커지는 행렬도 있다. 피벗 선택은 계산 절차의 안정성을 돕는 장치이며, 문제 자체의 민감도를 없애는 장치는 아니다.

LU 분해로 소거 결과 재사용하기

LU 분해(LU factorization)는 소거 과정을 두 삼각행렬로 저장한다. 아래삼각행렬 L에는 소거 배수가 들어가고, 위삼각행렬 U에는 소거 후 계수가 들어간다. 부분 피벗의 행 교환을 P로 나타내면 관계는 PA = LU다. 여기서 P는 실제로 큰 행렬로 만들 필요가 없다. 각 단계에서 어느 행과 교환했는지 정수 배열에 기록하면 된다.

가우스 소거로 우변까지 동시에 처리하는 대신, 계수행렬만 먼저 분해한다고 생각하면 이해하기 쉽다. 분해 후에는 Ly = Pb를 풀고, 이어서 Ux = y를 푼다. 첫 과정은 위에서 아래로 값을 구하는 전진 대입(forward substitution)이며, 두 번째 과정은 앞에서 설명한 후진 대입이다.

완성 코드는 L과 U를 하나의 배열에 겹쳐 저장한다. 대각 위와 대각에는 U를 저장하고, 대각 아래에는 L의 소거 배수를 저장한다. L의 대각은 모두 1이므로 별도 공간을 쓰지 않는다. 따라서 대각 아래 성분은 소거 후의 0이 아니라, 나중에 우변을 처리할 때 사용할 배수다.

k번째 단계에서 행을 바꿀 때는 전체 행을 바꾼다. 오른쪽의 아직 처리하지 않은 계수뿐 아니라 왼쪽에 저장해 둔 L의 배수도 함께 움직여야 한다. 이 규칙을 빠뜨리면 첫 분해는 그럴듯해 보여도 뒤의 전진 대입이 잘못된다.

정수 배열 piv(k)는 k번째 단계에서 k행과 교환한 행 번호다. 이는 최종 행 순서를 직접 적은 배열이 아니다. 새 우변에 교환을 적용할 때도 k = 1부터 n까지 같은 교환을 차례로 수행한다. 교환들을 임의의 순서로 적용하면 다른 순열이 될 수 있다.

일반 행렬의 분해에는 대략 n의 세제곱에 비례하는 연산이 필요하다. 이미 분해한 뒤 우변 하나를 푸는 데는 n의 제곱에 비례하는 연산이 필요하다. 우변이 여러 개일수록 분해를 재사용하는 효과가 커진다. 다만 계수나 연결 조건이 바뀌어 A가 달라지면 기존 분해도 다시 계산해야 한다.

삼중대각 풀이와 잔차 검증

세 대각선만 남기는 토머스 알고리즘

삼중대각 행렬(tridiagonal matrix)은 주대각선과 바로 아래·위 대각선에만 값이 있는 행렬이다. 앞의 온도 식이 이 구조다. 각 식이 자신의 온도와 바로 이웃한 온도만 포함하므로, 행 번호와 열 번호가 두 칸 이상 떨어진 성분은 0이다.

토머스 알고리즘(Thomas algorithm)은 이 구조에 맞춘 가우스 소거다. 아래 대각선 lower, 주대각선 diag, 위 대각선 upper를 길이 n의 배열로 저장한다. 이 장에서는 인덱스를 식의 행 번호와 맞추기 위해 lower(1)과 upper(n)을 사용하지 않는 자리로 두고 0을 넣는다.

lower = [ 0, -1, -1 ]
diag  = [ 4,  4,  3 ]
upper = [-1, -1,  0 ]

i번째 행을 처리할 때 배수는 lower(i) / d(i-1)이다. 여기서 d는 소거하면서 갱신하는 대각선이다. d(i)에서 배수와 upper(i-1)의 곱을 빼고, 우변 q(i)에서도 배수와 q(i-1)의 곱을 뺀다. 다른 성분을 훑을 필요가 없으므로 전체 소거와 대입이 n에 비례한다.

이 간단한 토머스 구현에는 행 교환이 없다. 원래 행렬이 가역이어도 중간 분모가 0이 되면 계산할 수 없다. 예를 들어 대각선이 모두 0이고 위·아래 성분이 1인 2행 행렬은 가역이지만 첫 분모가 0이다. 그런 문제에는 부분 피벗 일반 풀이기를 사용할 수 있다.

예제 행렬은 각 행의 대각 성분 절댓값이 나머지 성분 절댓값의 합보다 크다. 이처럼 엄격한 대각 우세를 갖는 삼중대각 행렬은 소거 중 0 피벗을 피하는 충분조건을 제공한다. 모든 삼중대각 행렬이 이런 조건을 만족하는 것은 아니므로, 구조만 보고 무조건 토머스 풀이기를 선택해서는 안 된다.

원래 식으로 돌아가 확인하기

잔차(residual)는 r = b - Ax다. 구한 x를 원래 식에 넣었을 때 얼마나 어긋나는지 나타낸다. 잔차의 각 성분은 대응하는 식의 불일치량이다. 소거로 바뀐 행렬에 해를 대입하면 원래 문제를 검사하지 못하므로, 원래 A와 b를 보존해야 한다.

온도나 우변의 단위와 크기가 달라지면 같은 잔차 값의 의미도 달라진다. 여기서는 무한대 노름(infinity norm)을 사용하여 다음과 같이 크기를 정규화한다. 벡터의 노름은 성분 절댓값의 최댓값이고, 행렬의 노름은 각 행의 절댓값 합 중 최댓값이다.

rho = ||b - A x||inf / (||A||inf ||x||inf + ||b||inf)

분모가 0이면 이 식으로 나누지 않는다. 완성 코드에서는 그 경우 잔차 노름 자체를 반환한다. 유한한 입력에서 분모가 0이면 b가 0이고 Ax도 0이므로 잔차 역시 0이다. 실제 응용에서는 입력과 중간 결과가 유한한지도 별도로 확인해야 한다.

정규화 잔차가 작다는 것은 원래 식을 수치적으로 잘 만족한다는 뜻이다. 해 자체의 오차가 작다는 뜻과는 다르다. 민감한 행렬에서는 큰 해 오차도 작은 잔차를 만들 수 있다. 이 장의 작은 예제는 알려진 정답과의 비교도 함께 수행하므로, 풀이 구현과 잔차 계산을 두 방향에서 점검한다.

출력에서는 잔차를 고정 소수점으로 보여 준다. 다만 화면에 0.000000으로 보인다고 해서 내부 값이 정확히 0이라는 뜻은 아니다. 출력보다 더 작은 오차는 반올림되어 사라진다. 그래서 프로그램은 내부 값으로 허용 기준을 검사하고, 화면에는 그 검사를 통과한 결과를 출력한다.

완성 코드

다음 프로그램은 하나의 파일에 모듈과 실행 프로그램을 함께 둔다. 일반 풀이기에는 행 순서를 바꾼 온도 식을 넣고, 삼중대각 풀이기에는 원래 순서의 식을 넣는다. 두 풀이가 같은 온도를 내는지 검사한다. 또한 같은 LU 분해에 새 우변을 적용하여 20, 10, 10을 얻는지 확인한다.

풀이 함수는 크기가 맞는 유한한 배열을 받는다는 전제를 둔다. 오류 번호 info가 0이면 계산을 끝냈으며, 양수이면 해당 단계에서 0 피벗을 만났다는 뜻이다. 이것은 조건수 추정이나 근접 특이성 진단을 대신하지 않는다. 예제를 짧게 유지하면서 성공 여부와 잔차 검증을 분리하는 인터페이스다.

main.f90

module linear_solvers
  use iso_fortran_env, only : real64
  implicit none
  private
  public :: lu_factor, lu_solve, thomas_solve, scaled_residual

contains

  subroutine lu_factor(a, lu, piv, info)
    real(real64), intent(in) :: a(:, :)
    real(real64), intent(out) :: lu(:, :)
    integer, intent(out) :: piv(:), info
    real(real64) :: row(size(a, 2))
    integer :: n, k, p, i

    n = size(a, 1)
    lu = a
    info = 0
    piv = 0

    do k = 1, n
      p = k - 1 + maxloc(abs(lu(k:n, k)), dim=1)
      piv(k) = p
      if (abs(lu(p, k)) <= 0.0_real64) then
        info = k
        return
      end if

      if (p /= k) then
        row = lu(k, :)
        lu(k, :) = lu(p, :)
        lu(p, :) = row
      end if

      do i = k + 1, n
        lu(i, k) = lu(i, k) / lu(k, k)
        lu(i, k+1:n) = lu(i, k+1:n) &
                      - lu(i, k) * lu(k, k+1:n)
      end do
    end do
  end subroutine lu_factor

  subroutine lu_solve(lu, piv, b, x)
    real(real64), intent(in) :: lu(:, :), b(:)
    integer, intent(in) :: piv(:)
    real(real64), intent(out) :: x(:)
    real(real64) :: temp
    integer :: n, k, p, i

    n = size(b)
    x = b

    do k = 1, n
      p = piv(k)
      if (p /= k) then
        temp = x(k)
        x(k) = x(p)
        x(p) = temp
      end if
    end do

    do i = 2, n
      x(i) = x(i) - dot_product(lu(i, 1:i-1), x(1:i-1))
    end do

    do i = n, 1, -1
      x(i) = (x(i) - dot_product(lu(i, i+1:n), x(i+1:n))) &
             / lu(i, i)
    end do
  end subroutine lu_solve

  subroutine thomas_solve(lower, diag, upper, b, x, info)
    real(real64), intent(in) :: lower(:), diag(:), upper(:), b(:)
    real(real64), intent(out) :: x(:)
    integer, intent(out) :: info
    real(real64) :: d(size(diag)), q(size(b)), m
    integer :: n, i

    n = size(diag)
    d = diag
    q = b
    x = 0.0_real64
    info = 0

    if (abs(d(1)) <= 0.0_real64) then
      info = 1
      return
    end if

    do i = 2, n
      m = lower(i) / d(i-1)
      d(i) = d(i) - m * upper(i-1)
      q(i) = q(i) - m * q(i-1)
      if (abs(d(i)) <= 0.0_real64) then
        info = i
        return
      end if
    end do

    x(n) = q(n) / d(n)
    do i = n - 1, 1, -1
      x(i) = (q(i) - upper(i) * x(i+1)) / d(i)
    end do
  end subroutine thomas_solve

  function scaled_residual(a, x, b) result(rho)
    real(real64), intent(in) :: a(:, :), x(:), b(:)
    real(real64) :: rho, anorm, denom

    anorm = maxval(sum(abs(a), dim=2))
    denom = anorm * maxval(abs(x)) + maxval(abs(b))
    rho = maxval(abs(b - matmul(a, x)))
    if (denom > 0.0_real64) rho = rho / denom
  end function scaled_residual

end module linear_solvers

program main
  use iso_fortran_env, only : real64
  use linear_solvers, only : lu_factor, lu_solve, thomas_solve, &
                            scaled_residual
  implicit none
  integer, parameter :: n = 3
  real(real64), parameter :: tol = 1.0e-12_real64
  real(real64) :: a(n, n), lu(n, n), tri(n, n)
  real(real64) :: b(n), b2(n), bt(n), x(n), x2(n), xt(n)
  real(real64) :: lower(n), diag(n), upper(n)
  real(real64) :: rho, rho2, rhot
  integer :: piv(n), info

  ! The first two equations are exchanged.
  a = reshape([ -1.0_real64,  4.0_real64,  0.0_real64, &
                 4.0_real64, -1.0_real64, -1.0_real64, &
                -1.0_real64,  0.0_real64,  3.0_real64 ], [n, n])
  b = [10.0_real64, 30.0_real64, 50.0_real64]

  call lu_factor(a, lu, piv, info)
  if (info /= 0) error stop 'LU factorization failed'
  call lu_solve(lu, piv, b, x)

  ! Reuse the same factors for a different right-hand side.
  b2 = [10.0_real64, 70.0_real64, 20.0_real64]
  call lu_solve(lu, piv, b2, x2)

  lower = [ 0.0_real64, -1.0_real64, -1.0_real64]
  diag  = [ 4.0_real64,  4.0_real64,  3.0_real64]
  upper = [-1.0_real64, -1.0_real64,  0.0_real64]
  bt = [30.0_real64, 10.0_real64, 50.0_real64]
  call thomas_solve(lower, diag, upper, bt, xt, info)
  if (info /= 0) error stop 'Thomas elimination failed'

  tri = reshape([ 4.0_real64, -1.0_real64,  0.0_real64, &
                 -1.0_real64,  4.0_real64, -1.0_real64, &
                  0.0_real64, -1.0_real64,  3.0_real64 ], [n, n])
  rho = scaled_residual(a, x, b)
  rho2 = scaled_residual(a, x2, b2)
  rhot = scaled_residual(tri, xt, bt)

  if (maxval(abs(x - [10.0_real64, 10.0_real64, &
                     20.0_real64])) > tol) error stop 'Wrong LU solution'
  if (maxval(abs(x2 - [20.0_real64, 10.0_real64, &
                      10.0_real64])) > tol) error stop 'Wrong reused solution'
  if (maxval(abs(x - xt)) > tol) error stop 'Solvers disagree'
  if (max(rho, rho2, rhot) > tol) error stop 'Residual too large'

  write(*, '(A,I2)') 'First pivot row:', piv(1)
  write(*, '(A,3F8.2)') 'LU temperature:    ', x
  write(*, '(A,3F8.2)') 'Thomas temperature:', xt
  write(*, '(A,3F8.2)') 'Reused LU solution:', x2
  write(*, '(A,F10.6)') 'LU residual:       ', rho
  write(*, '(A,F10.6)') 'Thomas residual:   ', rhot
  write(*, '(A,F10.6)') 'Reused residual:   ', rho2
  write(*, '(A)') 'Validation: PASS'
end program main

줄별 해설

use iso_fortran_env, only : real64는 계산에 사용할 실수 종류를 가져온다. 모든 실수 상수에 같은 종류 접미사를 붙여 상수가 먼저 더 낮은 정밀도로 해석되는 일을 피한다. 모듈은 계산 절차를 제공하고, 실행 프로그램은 입력 준비와 결과 검사를 담당한다.

lu_factor의 a는 입력 전용이며, lu = a로 작업용 복사본을 만든다. 출력 인수 lu와 piv는 호출자가 올바른 크기로 준비한다. 이 예제의 계약은 n이 1 이상이고, 행렬이 n행 n열이며, 교환 배열 길이가 n이라는 것이다. 크기가 맞지 않는 호출을 허용하는 범용 인터페이스로 확장한다면 진입 시 크기 검사도 추가해야 한다.

maxloc(abs(lu(k:n, k)), dim=1)은 현재 열의 남은 구간에서 최대 절댓값의 위치를 반환한다. 반환 위치는 부분 배열 안에서 1부터 시작한다. 따라서 원래 행 번호로 옮기기 위해 k - 1을 더한다. 같은 최댓값이 여러 개면 이 호출은 처음 나타난 위치를 선택한다.

abs(lu(p, k)) <= 0.0_real64는 유한한 값에 대해 0 피벗을 검사한다. 작은 양수까지 일괄 실패로 처리하지는 않는다. 예를 들어 행렬 전체를 아주 작은 수로 배율 조정했다고 해서 해의 민감도가 반드시 나빠지는 것은 아니다. 실제 진단에는 행렬 규모와 조건 정보를 함께 봐야 한다.

row는 행 교환 중 값을 보관하는 임시 배열이다. 이후 lu(i, k)에는 나눗셈으로 얻은 배수를 저장한다. 바로 다음 문장은 k열 오른쪽만 갱신하므로 저장한 배수가 남는다. 마지막 단계에는 아래 행이 없어 내부 반복문이 실행되지 않는다.

lu_solve는 먼저 x = b로 우변을 복사한다. 해 배열을 작업 공간으로 사용하므로 우변 원본은 보존된다. 교환 기록을 순서대로 적용한 뒤, L의 대각이 1이라는 사실을 이용하여 전진 대입에서는 나누지 않는다. 이 함수는 성공한 분해 결과만 받으므로 피벗 오류 검사를 반복하지 않는다.

후진 대입의 i+1:n은 i = n일 때 길이가 0인 배열 구간이 된다. 길이가 같은 빈 배열들의 dot_product 결과는 0이다. 따라서 마지막 행만 따로 처리하는 분기 없이 같은 식으로 계산할 수 있다. 배열 구간의 시작과 끝이 뒤집혔다고 해서 항상 오류가 되는 것은 아니다.

thomas_solve는 대각선과 우변을 각각 d, q에 복사한다. 원래 세 대각선과 우변은 바뀌지 않는다. 첫 분모를 먼저 검사하고, 각 단계에서 새 대각 성분을 만든 직후 다시 검사한다. 오류가 있으면 호출자는 해 배열을 사용하지 않고 종료한다.

sum(abs(a), dim=2)는 열 방향으로 더해 각 행의 절댓값 합을 만든다. 그 최댓값이 행렬의 무한대 노름이다. maxval(abs(a))는 가장 큰 성분 하나만 구하므로 여기서 필요한 행렬 노름과 다르다.

reshape에 넣은 값은 열 단위로 채워진다. 일반 행렬의 첫 세 값 -1, 4, 0은 첫째 행이 아니라 첫째 열이다. 배열 생성자의 겉모양을 행렬처럼 읽지 말고, 작성 후 식의 계수와 행·열 위치를 대조해야 한다.

마지막 검사는 알려진 해, 두 풀이의 일치, 정규화 잔차를 각각 확인한다. tol은 이 작은 검증 문제의 기준이다. 다른 크기와 조건의 문제에 그대로 적용할 보편적인 정확도 목표는 아니다. 검사를 모두 통과한 뒤에만 Validation: PASS가 출력된다.

실행 결과

파일을 저장한 디렉터리에서 다음과 같이 컴파일하고 실행한다. 입력 파일이나 표준 입력은 필요하지 않다. 모듈이 같은 소스 파일에서 실행 프로그램보다 먼저 나오므로 별도 컴파일 순서를 정할 필요도 없다.

gfortran -std=f2018 -Wall main.f90 -o linear_demo
./linear_demo

예상 출력은 다음과 같다. 첫 피벗 행 번호가 2이므로 행 교환 경로가 사용되었음을 확인할 수 있다. 두 온도 풀이의 값은 같고, 새 우변을 적용한 결과는 별도의 알려진 해와 일치한다.

First pivot row: 2
LU temperature:       10.00   10.00   20.00
Thomas temperature:   10.00   10.00   20.00
Reused LU solution:   20.00   10.00   10.00
LU residual:         0.000000
Thomas residual:     0.000000
Reused residual:     0.000000
Validation: PASS

실수 값은 모두 폭과 소수 자릿수가 지정된 서식으로 출력한다. 이 예제의 잔차는 출력 해상도보다 작아서 0으로 표시된다. 출력된 여섯 자리만 보고 오차의 유무를 판정하지 않고, 내부 허용 기준 검사와 함께 해석한다.

실무에서 자주 틀리는 것

행 교환을 새 우변에 적용하지 않는다

LU 배열만 재사용하고 새 우변을 그대로 전진 대입에 넣으면 LUx = b를 푸는 셈이다. 필요한 식은 LUx = Pb다. 아래 코드는 lu_solve 내부의 우변 준비 부분을 비교한다. 그 뒤의 전진 대입과 후진 대입은 동일하다.

틀린 코드다.

x = b
! Forward substitution starts here without applying piv.

고친 코드다.

x = b
do k = 1, n
  p = piv(k)
  if (p /= k) then
    temp = x(k)
    x(k) = x(p)
    x(p) = temp
  end if
end do

교환 기록은 분해와 한 묶음으로 보관한다. 다른 행렬에서 얻은 교환 기록을 사용하거나, 기록을 최종 행 순서로 오해하여 한 번에 재배열하는 오류도 피해야 한다.

행의 오른쪽 부분만 교환한다

두 번째 이후 피벗 단계에서는 왼쪽 열에 이미 L의 배수가 저장되어 있다. 오른쪽 계수만 바꾸면 L과 U가 서로 다른 행 순서를 나타낸다. 아래의 row는 완성 코드처럼 길이 n인 임시 배열이다.

틀린 코드다.

row(k:n) = lu(k, k:n)
lu(k, k:n) = lu(p, k:n)
lu(p, k:n) = row(k:n)

고친 코드다.

row = lu(k, :)
lu(k, :) = lu(p, :)
lu(p, :) = row

이 오류는 첫 단계에서만 교환하는 예제로는 드러나지 않을 수 있다. 풀이기를 확장할 때는 뒤 단계에서 교환이 필요한 알려진 해 문제도 검증 자료에 포함한다.

토머스 소거에서 원래 대각선으로 나눈다

토머스 알고리즘의 다음 단계는 이전 단계에서 갱신한 대각선에 의존한다. 원래 diag를 분모로 사용하면 첫 소거의 영향이 다음 식에 전달되지 않는다. 예제에서는 두 번째 작업 대각 성분이 4에서 3.75로 바뀐다.

틀린 코드다.

m = lower(i) / diag(i-1)
d(i) = d(i) - m * upper(i-1)

고친 코드다.

m = lower(i) / d(i-1)
d(i) = d(i) - m * upper(i-1)
q(i) = q(i) - m * q(i-1)

계수와 우변을 같은 소거 배수로 갱신해야 한다. 계수만 갱신하면 원래 식과 같은 해를 유지하지 못한다.

분해 배열로 잔차를 계산한다

겹쳐 저장한 LU 배열은 원래 A가 아니다. 특히 대각 아래에 들어 있는 값은 원래 계수가 아니라 소거 배수다. 이 배열 전체에 matmul을 적용해도 L과 U의 곱을 계산하는 것이 아니다.

틀린 코드다.

rho = scaled_residual(lu, x, b)

고친 코드다.

rho = scaled_residual(a, x, b)

큰 문제에서 원래 행렬을 따로 저장하기 어렵다면, 원래 계수로 Ax를 계산하는 절차를 별도로 보존할 수 있다. 삼중대각 문제도 세 대각선만으로 각 행의 곱을 계산할 수 있다. 핵심은 검사에 사용하는 연산이 원래 문제를 나타내야 한다는 점이다.

한눈에 보기

풀이 비용과 검증 책임은 서로 다른 항목이다

항목	핵심 계산	비용 규모	기억할 점
부분 피벗	현재 열의 최대 절댓값 선택	일반 분해에 포함	전체 행과 교환 기록을 관리한다
LU 분해	PA = LU	연산 O(n³), 저장 O(n²)	A가 같을 때 재사용한다
삼각 대입	Ly = Pb, Ux = y	우변마다 O(n²)	성공한 분해를 입력으로 받는다
토머스 알고리즘	대각선과 우변의 순차 갱신	연산·저장 O(n)	이 구현은 행 교환을 하지 않는다
정규화 잔차	원래 b - Ax의 크기 비교	밀집 O(n²), 삼중대각 O(n)	작은 잔차만으로 해 오차를 단정하지 않는다

온도 계산 블록의 계수가 고정되어 있으면 일반 행렬은 한 번 분해하고 우변만 교체할 수 있다. 삼중대각 계산도 대각선 소거 결과와 배수를 별도로 저장하는 방식으로 재사용할 수 있다. 완성 코드의 토머스 함수는 이해하기 쉬운 단일 호출 형태이며, 호출할 때마다 계수 소거를 다시 수행한다.

앞 장의 정밀도 선택은 이 장에서도 계속 중요하다. 그러나 실수 종류를 넓히는 것만으로 잘못된 행 교환이나 우변 갱신을 고칠 수는 없다. 알고리즘의 계약을 지키고, 원래 식을 보존하고, 알려진 작은 문제로 점검하는 습관이 먼저 필요하다.

연습 문제

1. 원래 삼중대각 행렬에서 해가 5, 10, 15가 되도록 새 우변을 손으로 계산한다. 일반 풀이기에 사용할 우변의 순서도 구하고, 기존 LU 분해를 재사용하여 확인한다.
2. 토머스 알고리즘에서 소거 후 대각선 d의 세 성분을 계산한다. 원래 대각선으로 계속 나누는 잘못된 코드가 어느 단계부터 다른 결과를 만드는지 설명한다.
3. 2행 행렬 A = [0, 1; 1, 0]과 우변 b = [2, 3]을 생각한다. 토머스 풀이기의 오류 번호와 부분 피벗 LU 풀이기의 해를 구한다. 행렬이 가역인지도 설명한다.
4. 삼중대각 행렬을 밀집 배열로 만들지 않고 잔차 벡터와 행렬 무한대 노름을 계산하는 코드를 작성한다. 배열 길이는 n ≥ 1이며 lower(1), upper(n)은 0이라는 규칙을 사용한다.

정답과 해설

1. 원래 순서의 우변은 10, 20, 35다. 첫 식은 4×5 - 10 = 10이고, 둘째 식은 -5 + 4×10 - 15 = 20이며, 셋째 식은 -10 + 3×15 = 35다. 일반 풀이기의 A는 첫 두 식을 교환한 상태이므로 우변은 20, 10, 35로 준비한다.

   b2 = [20.0_real64, 10.0_real64, 35.0_real64]
   call lu_solve(lu, piv, b2, x2)

   이때 A는 바뀌지 않았으므로 다시 분해하지 않는다. 완성 코드의 기존 두 번째 정답 검사도 새 해 5, 10, 15에 맞게 바꿔야 한다. 입력만 바꾸고 검증 기대값을 그대로 두면 검사는 의도대로 실패한다.

2. 첫 작업 대각선은 4다. 두 번째 배수는 -1/4이므로 d(2) = 4 - 1/4 = 15/4다. 세 번째 배수는 -4/15이고, d(3) = 3 - 4/15 = 41/15다. 소수로는 4, 3.75, 약 2.733333이다.

   잘못된 코드도 두 번째 행에서는 같은 분모 4를 사용한다. 세 번째 행부터는 갱신된 15/4 대신 원래 값 4로 나누므로 d(3)을 11/4로 계산한다. 오류가 시작되는 위치와 최종 결과가 눈에 띄는 위치는 다를 수 있다.

3. 토머스 풀이기는 첫 대각 성분이 0이므로 info = 1을 반환한다. 부분 피벗은 첫째 행과 둘째 행을 교환하여 단위행렬을 얻는다. 우변도 3, 2로 바뀌므로 해는 x1 = 3, x2 = 2다.

   이 행렬의 행렬식은 -1이므로 가역이다. 따라서 토머스 구현의 실패를 곧바로 원래 문제의 해가 없다는 뜻으로 해석해서는 안 된다. 실패 이유는 이 구현이 요구하는 소거 경로에 0 분모가 있다는 것이다.

4. 아래 코드는 r(n)과 rowsum(n)을 실수 배열로, anorm을 실수로, i를 정수로 선언한 상황의 계산 부분이다. 주대각 기여를 먼저 넣고, 존재하는 이웃 항만 추가한다.

   r = b - diag * x
   rowsum = abs(diag)
   
   do i = 1, n
     if (i > 1) then
       r(i) = r(i) - lower(i) * x(i-1)
       rowsum(i) = rowsum(i) + abs(lower(i))
     end if
     if (i < n) then
       r(i) = r(i) - upper(i) * x(i+1)
       rowsum(i) = rowsum(i) + abs(upper(i))
     end if
   end do
   
   anorm = maxval(rowsum)

   잔차 노름은 maxval(abs(r))다. 이후 정규화 분모와 0 분모 처리는 완성 코드와 같다. 이 방식은 밀집 행렬을 만들지 않으므로 연산과 저장 공간이 모두 n에 비례한다. n = 1일 때도 두 이웃 분기가 실행되지 않아 같은 코드가 동작한다.