토머스 알고리즘

수치해석 시뮬레이션 마지막 수정: 2026-08-03 04:12:07

1. 개요[편집]

토머스 알고리즘 Thomas algorithm (TDMA)
다른 이름	TriDiagonal Matrix Algorithm, 삼중대각 알고리즘
분류	직접해법 — 피벗 없는 LU 분해의 삼중대각 특수화
비용	연산 약 $8N$, 메모리 $O(N)$
안전 조건	대각우세 또는 대칭 양정부호
주 고객	ADI 법, 크랭크-니콜슨법, 스플라인, 압력보정

밀집행렬을 $O(N^3)$ 에 푸는 게 국룰인 세상에서, 삼중대각만은 $N$ 에 선형으로 풀린다. 심지어 피벗도 안 한다.

토머스 알고리즘(Thomas algorithm, TDMA)은 삼중대각 행렬 연립방정식 $A\mathbf{x} = \mathbf{d}$ 를 피벗팅 없이 전진소거·후진대입 두 번의 훑기로 $O(N)$ 에 푸는 직접해법이다. 본질적으로는 가우스 소거법을 대역폭 1짜리 희소행렬에 특수화한 것으로, 새로 채워지는(fill-in) 원소가 하나도 없다는 성질 덕분에 밀집 소거의 $O(N^3)$ 이 통째로 $O(N)$ 으로 내려앉는다.¹

이름은 1949년경 이 절차를 정리한 르웰린 토머스(Llewellyn H. Thomas)에서 왔지만, 절차 자체는 삼중대각 행렬에 소거를 적용하면 누구나 도달하는 형태라 “재발견의 역사”가 길다. 중요한 건 이름이 아니라 왜 피벗 없이도 안전한가, 그리고 그 안전 조건이 깨지면 무슨 일이 벌어지는가다. 시뮬레이션 코드에서 이 알고리즘은 ADI 법·크랭크-니콜슨법·압력보정·스플라인 보간의 최하층에 깔려 있어서, 여기서 조용히 어긋나면 위층의 모든 결과가 함께 어긋난다.

2. 알고리즘[편집]

$i$ 번째 방정식이 이웃 셋만 묶는다고 하자.

a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i, \qquad i = 1,\dots,N,\quad a_1 = c_N = 0

전진소거는 아래 대각을 없애며 위 대각과 우변만 갱신한다. 수정된 대각을 $\beta_i$ 로 두면

\beta_1 = b_1,\qquad \beta_i = b_i - a_i c'_{i-1},\qquad c'_i = \frac{c_i}{\beta_i}

d'_1 = \frac{d_1}{\beta_1},\qquad d'_i = \frac{d_i - a_i d'_{i-1}}{\beta_i}

후진대입은 끝에서부터 되짚는다.

x_N = d'_N,\qquad x_i = d'_i - c'_i x_{i+1}\quad (i = N-1,\dots,1)

곱셈·나눗셈·뺄셈을 전부 세면 대략 $8N$ 회, 임시 배열은 $c'$ 와 $d'$ 두 개면 충분하다.² 구조를 보면 $\beta_i$ 가 곧 $U$ 의 대각, $a_i/\beta_{i-1}$ 이 $L$ 의 부대각이므로 토머스 = 삼중대각 LU + 전방·후방 치환이다. 그래서 우변이 여러 개면 소거를 한 번만 하고 치환만 반복하면 된다 — ADI 법처럼 계수는 그대로인데 우변만 매 스텝 바뀌는 상황에서 이 재사용이 결정적이다.

3. 왜 피벗 없이 안전한가[편집]

일반 소거에서 피벗팅이 필요한 이유는 대각 원소가 0에 가까워지면 승수가 폭발하기 때문이다. 삼중대각에서는 대각우세가 이걸 구조적으로 막는다. 모든 행에서 $\lvert b_i \rvert > \lvert a_i \rvert + \lvert c_i \rvert$ 라 하고, $\lvert c'_i \rvert < 1$ 을 귀납으로 보이자. 우선 $\lvert c'_1 \rvert = \lvert c_1 \rvert / \lvert b_1 \rvert < 1$ 이다. $\lvert c'_{i-1} \rvert < 1$ 을 가정하면

\lvert \beta_i \rvert = \lvert b_i - a_i c'_{i-1} \rvert \ \ge\ \lvert b_i \rvert - \lvert a_i \rvert \lvert c'_{i-1} \rvert \ >\ \lvert b_i \rvert - \lvert a_i \rvert \ >\ \lvert c_i \rvert

따라서 $\beta_i \neq 0$ 이고 $\lvert c'_i \rvert = \lvert c_i \rvert / \lvert \beta_i \rvert < 1$ 이다. 즉 대각우세는 소거를 따라 전파된다. 이 한 줄에 세 가지가 동시에 들어 있다.

0으로 나눌 일이 없다. $\lvert \beta_i \rvert$ 는 항상 대각우세 여유분 $\lvert b_i \rvert - \lvert a_i \rvert$ 보다 크다. 피벗 탐색이 불필요한 이유가 이것이다.
수정 대각이 원래 대각의 부호와 크기를 유지한다. 빼는 양 $\lvert a_i c'_{i-1} \rvert$ 이 $\lvert a_i \rvert$ 보다 작으므로 $\beta_i$ 는 $b_i$ 의 부호를 그대로 물려받고, 크기도 여유분 아래로 내려가지 않는다. 확산 이산화에서 흔한 $a_i = c_i = -\mu$ , $b_i = 1+2\mu$ 를 넣어 보면 $\beta_i \ge 1 + \mu$ 로 고정된다. 자리수 손실이 일어날 여지 자체가 없다.
후진대입이 오차를 증폭하지 않는다. $\lvert c'_i \rvert < 1$ 이므로 $x_{i+1}$ 의 오차가 $x_i$ 로 넘어올 때 감쇠한다.

정량적으로도 결과가 좋다. 대각우세 삼중대각계에서 피벗 없는 소거의 성장인자는 2를 넘지 않으며, 따라서 토머스 알고리즘은 후진안정이다 — 계산된 해는 계수를 반올림 단위 정도만 흔든 문제의 정확한 해다. 자세한 의미는 후진 오차 해석 참고. 대칭 양정부호인 경우에도 같은 결론이 서며, 이때 토머스는 삼중대각 촐레스키 분해와 동치다.

4. 대각우세가 깨지는 순간[편집]

문제는 이 조건이 실제 시뮬레이션에서 생각보다 쉽게 깨진다는 것이다. 정상 대류-확산 방정식 $a u_x = \alpha u_{xx}$ 를 균일 격자에서 중심차분하고 셀 페클레 수 $Pe_h = ah/\alpha$ 로 정리하면

\left(1 + \tfrac{Pe_h}{2}\right)u_{i-1} - 2u_i + \left(1 - \tfrac{Pe_h}{2}\right)u_{i+1} = 0

이다. $Pe_h < 2$ 이면 $\lvert a_i \rvert + \lvert c_i \rvert = 2 = \lvert b_i \rvert$ 로 겨우 등호(약한 대각우세)를 유지하지만, $Pe_h > 2$ 를 넘는 순간 $c_i$ 의 부호가 뒤집히면서 $\lvert a_i \rvert + \lvert c_i \rvert = Pe_h > 2$ 가 되어 대각우세가 완전히 사라진다. 중심차분이 $Pe_h > 2$ 에서 격자 진동을 일으킨다는 그 유명한 임계값과, 토머스의 안전 보장이 무너지는 임계값이 정확히 같은 지점이라는 게 포인트다. 페클레 수 문서의 진동 현상은 이 선형대수적 사건의 물리적 얼굴인 셈이다.

이때의 대응은 셋 중 하나다. (1) 풍상 차분으로 대각우세를 인위적으로 회복시킨다 — 대신 수치 소산과 분산을 감수한다. (2) 문제되는 항을 우변으로 넘겨 명시적으로 처리한다(지연 보정, IMEX). (3) 그냥 피벗을 하는 삼중대각 루틴을 쓴다. LAPACK의 dgtsv가 정확히 이쪽으로, 부분 피벗팅을 하면 두 번째 상부 대각에 채움이 생기므로 배열 하나를 더 잡아야 하지만 비용은 여전히 $O(N)$ 이다. “삼중대각이니까 토머스” 라는 반사신경이 통하지 않는 유일한 구간이 여기다.

5. 순환 삼중대각과 블록 삼중대각[편집]

주기경계조건을 걸면 $(1,N)$ 과 $(N,1)$ 모서리에 원소가 하나씩 붙어 순환(cyclic) 삼중대각이 된다. 순수 삼중대각이 아니므로 토머스를 그대로 돌리면 경계 근처 한 줄만 조용히 틀린다.³ 표준 처방은 셔먼-모리슨 공식이다. 모서리 항을 계수 1짜리 랭크-1 수정으로 몰아내 $A = T + \mathbf{u}\mathbf{v}^{\mathsf{T}}$ 로 쓰고

\mathbf{x} = \mathbf{y} - \frac{\mathbf{v}^{\mathsf{T}}\mathbf{y}}{1 + \mathbf{v}^{\mathsf{T}}\mathbf{z}}\,\mathbf{z}, \qquad T\mathbf{y} = \mathbf{d},\ \ T\mathbf{z} = \mathbf{u}

로 푼다. 즉 토머스를 두 번 돌리면 끝이고 비용은 여전히 $O(N)$ 이다. 실무 팁 하나: $T$ 를 만들 때 $b_1$ 에서 빼는 보정 파라미터 $\gamma$ 는 $\gamma = -b_1$ 처럼 잡아 $b_1 - \gamma$ 에서 자리수 손실이 나지 않게 한다. 이건 부동소수점 연산의 기본기지만, 순환 TDMA 구현에서 가장 자주 빠뜨리는 부분이기도 하다.

미지수가 격자점마다 여러 개인 연성 문제에서는 $a_i, b_i, c_i$ 가 $m \times m$ 블록이 되어 블록 삼중대각이 된다. 점화식은 형태 그대로이고 나눗셈이 블록 역행렬로 바뀐다.

\mathbf{C}'_i = \left(\mathbf{B}_i - \mathbf{A}_i \mathbf{C}'_{i-1}\right)^{-1}\mathbf{C}_i

비용은 $O(N m^3)$ . 블록 안쪽 소거에는 피벗팅을 마음껏 써도 되고(블록 크기가 작아 비싸지 않다), 안전 조건은 $\lVert \mathbf{B}_i^{-1}\rVert(\lVert \mathbf{A}_i\rVert + \lVert \mathbf{C}_i\rVert) < 1$ 꼴의 블록 대각우세로 자연스럽게 일반화된다. 3차원 압축성 전산유체역학 코드의 라인 암시(line-implicit) 풀이가 대표적으로 $m = 5$ 짜리 블록 TDMA를 돌린다.

6. 병렬화가 어려운 이유[편집]

토머스의 약점은 딱 하나, 완전한 순차 재귀라는 것이다. $\beta_i$ 는 $\beta_{i-1}$ 을 알아야 계산되고 $x_i$ 는 $x_{i+1}$ 을 알아야 나온다. 임계 경로 길이가 $2N$ 이라 프로세서를 $N$ 개 붙여도 얻는 게 없다. 대안은 두 갈래다.

순환 축소(cyclic reduction, CR). 홀수 번호 미지수를 이웃 두 식으로 소거하면 짝수 번호만 남은 절반 크기의 삼중대각계가 나온다. 이걸 $\log_2 N$ 번 재귀하면 병렬 깊이가 $O(\log N)$ 으로 떨어진다. 총 연산량은 토머스의 두 배 남짓이지만 각 레벨 안은 완전 병렬이다.⁴ 다행히 각 축소 단계가 대각우세를 보존하므로 안정성 논거도 그대로 물려받는다.
PCR(parallel cyclic reduction). 축소를 홀수 행만이 아니라 모든 행에 동시에 적용한다. $\log_2 N$ 단계 뒤에는 모든 방정식이 완전히 분리되어 후진대입 자체가 사라진다. 연산량은 $O(N \log N)$ 으로 늘지만 단계마다 통신 패턴이 균일해 GPU에 잘 맞는다. 실제 GPU 라이브러리의 삼중대각 루틴은 대개 CR과 PCR을 섞은 하이브리드다.

다만 실무에서 병렬 TDMA가 필요한 상황은 생각보다 드물다. ADI 법에서 반스텝마다 푸는 라인은 서로 독립인 계가 $N$ 개이므로, 계 하나를 쪼갤 게 아니라 계들을 나눠 맡기고 각자 평범한 토머스를 돌리면 된다. CR/PCR이 실제로 필요한 건 계가 하나뿐이거나 한 라인이 여러 랭크에 걸쳐 분산된 경우다. 병렬 컴퓨팅에서 늘 그렇듯, 알고리즘을 병렬화하기 전에 문제에 이미 있는 병렬성을 먼저 쓰는 게 이긴다.

7. 어디에 쓰이나[편집]

ADI 법과 라인 완화. 방향별 반스텝마다 격자선 하나가 삼중대각계 하나다. 다중격자법의 평활자로 쓰이는 라인 완화도 정확히 같은 구조이며, 종횡비가 큰 경계층 격자에서 점 완화보다 압도적으로 강하다.
1차원 암시 시간적분. 크랭크-니콜슨법으로 푼 1D 확산·슈뢰딩거 방정식, 금융의 옵션 가격 차분법이 전부 매 스텝 TDMA 한 번이다.
삼차 스플라인 보간. 자연 경계 조건의 스플라인은 대각 4, 부대각 1인 삼중대각계로 귀결된다. $4 > 1 + 1$ 이라 강한 대각우세가 공짜로 따라온다 — 보간과 근사에서 스플라인이 “안전하다”고 불리는 근거 중 하나.
컴팩트 차분법(Padé 도식). 도함수를 암시적으로 정의하는 고차 차분은 매 방향마다 삼중대각계를 푼다. 스펙트럴 방법에 준하는 분해능을 구조격자에서 얻는 대신, 값을 얻는 통로가 TDMA다.
압력보정과 전처리기. SIMPLE 알고리즘의 보정식을 라인 단위로 풀거나, 삼중대각 근사를 통째로 전처리기로 얹는 용법.
대칭 고유값 문제. 란초스 알고리즘이 만들어내는 삼중대각 $T_m$ 의 고유값을 이분법으로 찾을 때 쓰는 스텀 수열 계산이 사실상 토머스의 전진소거와 같은 점화식이다.

8. 관련 문서[편집]

대역행렬 일반으로 보면 대역폭 $p$ 인 계를 푸는 비용이 $O(Np^2)$ 인데, 삼중대각은 $p = 1$ 이라 $O(N)$ 이 된다. 즉 토머스가 특별히 영리한 게 아니라 삼중대각이 특별히 착한 것이다. 이 착함을 지키려고 사람들이 격자를 굳이 구조화된 라인으로 자르는 것이기도 하다. ↩
” $8N$ 연산”이라는 수치를 곧이곧대로 성능으로 번역하면 안 된다. 연산 대비 메모리 접근 비율이 극도로 낮아서 실제로는 완전한 메모리 대역폭 바운드다. 라인을 캐시에 맞게 묶어 도는 것(라인 블로킹)이 flop 수 줄이기보다 훨씬 이득인 이유. ↩
학부 과제에서 순환 삼중대각을 만나면 열에 아홉은 그냥 토머스를 돌리고 “경계 한 점만 살짝 이상한데 괜찮겠지” 하고 넘어간다. 괜찮지 않다. 주기 방향으로 전역 보존이 깨지므로 장시간 적분에서 질량이 야금야금 새는 형태로 되돌아온다. ↩
순환 축소는 1965년 호크니가 포아송 방정식 고속 풀이용으로 도입했다. 벡터 컴퓨터 시대에 한 번 각광받았다가 캐시 지역성이 나빠 CPU에서는 밀려났고, GPU 시대에 다시 부활했다. 하드웨어가 바뀌면 죽었던 알고리즘이 살아 돌아오는 이 바닥의 국룰. ↩