그람-슈미트

1. 개요

그람-슈미트 직교화 Gram–Schmidt orthogonalization
산출물	축약 QR A = QR (Q 는 m×n, 명시적으로 나옴)
두 변형	고전(CGS) · 수정(MGS) — 대수적으로 동일
비용	2mn2 flops (둘 다 같다)
직교성 손실	CGS O(uκ2) · MGS O(uκ) · 하우스홀더 O(u)
주 서식지	아놀디/GMRES 내부 직교화
병렬 후예	블록 GS · TSQR · CholeskyQR2

대수적으로 완전히 같은 두 줄의 코드가 유효숫자 8자리 차이를 낸다. 수치해석이 별개의 학문인 이유.

그람-슈미트 직교화(Gram–Schmidt orthogonalization)는 선형독립인 벡터열 $a_1,\dots,a_n$ 에서, 매 단계 앞서 만든 정규직교 벡터 방향 성분을 빼고 정규화하는 것만으로 같은 부분공간을 생성하는 정규직교열 $q_1,\dots,q_n$ 을 순서대로 만들어내는 절차다.

$$\tilde q_k = a_k - \sum_{i<k}(q_i^{\top}a_k)\,q_i, \qquad q_k = \frac{\tilde q_k}{\|\tilde q_k\|_2}$$

빼낸 계수 $r_{ik} = q_i^{\top}a_k$ 와 $r_{kk} = \|\tilde q_k\|_2$ 를 상삼각행렬에 모으면 그대로 축약 QR 분해 $A = QR$ 이 된다. 즉 이 절차는 “직교기저를 만드는 법”이자 “QR을 계산하는 법”이며, $\mathrm{span}\{a_1,\dots,a_k\} = \mathrm{span}\{q_1,\dots,q_k\}$ 가 모든 $k$ 에서 성립한다는 중첩 부분공간 성질이 뒤에 나올 모든 이야기의 씨앗이다.

이름은 그람(1883)과 슈미트(1907)에게서 왔지만 절차 자체는 라플라스와 코시에게 이미 있었다.¹ 문제는 이름이 아니라, 같은 식을 코드로 옮기는 방법이 둘이고 그 둘의 수치적 운명이 완전히 다르다는 것이다.

2. CGS와 MGS — 같은 식, 다른 코드

고전 그람-슈미트(CGS) 는 위 식을 글자 그대로 옮긴다. 계수 $r_{ik}$ 를 전부 원본 $a_k$ 에 대해 계산해 두고, 마지막에 한꺼번에 뺀다.

수정 그람-슈미트(MGS) 는 순서를 바꾼다. 작업 벡터 $v \leftarrow a_k$ 를 두고, $i = 1,\dots,k-1$ 에 대해 $r_{ik} = q_i^{\top}v$ 를 그 시점의 $v$ 에서 재고 곧바로 $v \leftarrow v - r_{ik}q_i$ 로 갱신한다.

정확 산술에서는 $q_i \perp q_j$ 이므로 $q_i^{\top}v = q_i^{\top}a_k$ 이고 두 결과가 완전히 같다. 부동소수점에서는 $q_i$ 들이 정확히 직교하지 않으므로 두 값이 갈라진다. 결정적 차이는 MGS가 이미 오염된 상태를 다시 측정한다는 점이다. 앞 단계에서 생긴 오차가 다음 단계의 계수에 반영되어 일부 상쇄되는 반면, CGS는 모든 계수를 오차를 모르는 원본에서 뽑기 때문에 오차가 그대로 누적된다.

연산량은 둘 다 $2mn^2$ flops로 동일하다. 그럼에도 CGS가 사라지지 않는 이유는 메모리 접근 패턴에 있다. CGS의 한 열은 $Q_{k-1}^{\top}a_k$ 라는 행렬-벡터 곱 한 번(BLAS-2)으로 끝나지만, MGS는 내적 $k-1$ 개를 순차로 해야 한다(BLAS-1). 분산 메모리에서는 이것이 통신 횟수 $O(n)$ 대 $O(n^2)$ 의 차이가 된다.

3. 직교성 손실 — $\kappa$, $\kappa^2$, 그리고 $O(u)$

계산된 $\hat Q$ 가 얼마나 직교인지를 $\|I - \hat Q^{\top}\hat Q\|_2$ 로 재면, 세 방법의 성적표가 이렇게 갈린다($u \approx 1.1\times10^{-16}$ 은 배정도 단위 반올림, $\kappa = \kappa_2(A)$).

방법	직교성 손실	성립 조건	비용
고전 GS (CGS)	$O(u\,\kappa^2)$	$u\kappa^2 \lesssim 1$	$2mn^2$
수정 GS (MGS)	$O(u\,\kappa)$	$u\kappa \lesssim 1$	$2mn^2$
CGS2 / MGS2 (2회)	$O(u)$	$u\kappa \lesssim 1$	$4mn^2$
하우스홀더 QR	$O(u)$	무조건	$2n^2(m-n/3)$

MGS의 상계는 비외르크(Å. Björck, 1967)가 증명했고, CGS의 $\kappa^2$ 는 오랫동안 “경험적으로 그렇다”였다가 지로-랑구-로즐로즈니크(2005)가 정리로 못 박았다. 숫자로 보면 실감이 난다. $\kappa = 10^4$ 이면 CGS의 손실이 $10^{-8}$, MGS는 $10^{-12}$ — 둘 다 쓸 만하다. $\kappa = 10^8$ 이 되는 순간 MGS는 $10^{-8}$ 로 버티지만 CGS는 $O(1)$, 즉 직교성을 완전히 잃는다. 학부 과제에서 CGS로 만든 벡터 두 개의 내적이 0.3쯤 나와 학생을 좌절시키는 그 현상이다.

하우스홀더가 $O(u)$ 인 것은 곱해지는 것이 전부 정확한 노름 보존 연산이기 때문이고, 여기에는 $\kappa$ 가 아예 등장하지 않는다. 대신 하우스홀더는 $Q$ 를 반사의 곱으로만 들고 있어서, 열을 하나 만들 때마다 그 열이 필요한 상황에는 쓰기 어렵다. 두 방법은 대체재가 아니라 용도가 다르다.

4. 재직교화 — “두 번이면 충분하다”

망가진 직교성을 되살리는 처방은 어이없을 만큼 단순하다. 한 번 더 직교화한다. 첫 통과가 남긴 오염 성분은 원본 대비 $O(u\kappa)$ 규모인데, 두 번째 통과는 그 작아진 벡터를 기준으로 다시 같은 비율을 깎으므로 결과가 $O(u)$ 로 떨어진다. “Twice is enough”는 케이헌의 구전 격언이고 파를레의 교과서를 통해 정설이 됐다.² $u\kappa \lesssim 1$ 만 만족하면 CGS2도, MGS2도 하우스홀더와 같은 등급의 직교성을 준다.

비용을 두 배로 낼 필요는 없다. DGKS 판정(Daniel–Gragg–Kaufman–Stewart, 1976)은 직교화 후 노름이 이전의 $1/\sqrt{2}$ 배 미만으로 줄었을 때만 — 즉 상쇄가 심하게 일어난 열에서만 — 두 번째 통과를 건다. 대부분의 열은 한 번으로 끝나므로 실측 오버헤드는 20% 수준이다. 병렬 환경에서는 오히려 CGS2를 일부러 고른다. 통신 2회짜리 CGS2가 통신 $n$ 회짜리 MGS보다 빠르고, 안정성은 더 좋기 때문이다.

5. 크릴로프 안에서는 왜 MGS인가

밀집 QR만 놓고 보면 하우스홀더가 이기는데, 크리로프 부분공간법의 심장부에서는 여전히 MGS가 기본값이다. 이유는 벡터가 만들어지는 순서 때문이다. 아놀디 알고리즘은 $Av_j$ 를 계산해야 다음 열이 생기고, 그러려면 $v_j$ 가 이미 완성돼 있어야 한다. 직교화할 행렬 전체를 미리 손에 쥐어야 하는 하우스홀더와는 궁합이 나쁘다.

여기에 더 강력한 변론이 있다. 페이지-로즐로즈니크-스트라코시(2006)는 MGS 기반 GMRES가 후진 안정임을 증명했다. $\hat Q$ 의 직교성이 무너지는데도 그렇다. 핵심은 직교성 손실이 잔차 노름에 반비례해서 커진다는 것이다. 즉 직교성이 죽는 시점은 잔차가 이미 반올림 수준까지 내려간 뒤이고, 그때는 더 이상 직교성이 필요하지 않다. 알고리즘이 망가지는 속도와 문제가 풀리는 속도가 정확히 맞물려 있다.³

비외르크와 페이지(1992)의 다른 결과도 같은 방향이다. $A$ 에 MGS를 돌린 것은 $\binom{O}{A}$ 라는 $(m+n)\times n$ 확대행렬에 하우스홀더 QR을 돌린 것과 정확히 같다. MGS는 열등한 방법이 아니라 변장한 하우스홀더이고, 그래서 MGS로 푼 최소자승법 해는 하우스홀더와 같은 등급으로 후진 안정하다. 망가지는 것은 $\hat Q$ 이지 $x$ 가 아니다.

대칭 행렬에서 이 절차가 3항 점화식으로 붕괴하는 것이 란초스 알고리즘이고, 거기서 벌어지는 직교성 붕괴는 무작위가 아니라 “수렴한 리츠벡터 방향”으로 일어난다는 페이지의 정리로 관리된다.

6. 블록과 TSQR — 병렬 기계 위에서

노드 수천 개짜리 기계에서 비용은 flops가 아니라 동기화 횟수다. 그래서 현대적 후예들은 전부 통신을 줄이는 방향으로 진화했다.

• 블록 그람-슈미트(BGS). 열을 블록 단위로 묶어 이전 블록 전체에 대해 한 번에 직교화한다. 블록 간 직교화가 행렬-행렬 곱(BLAS-3)이 되어 캐시가 살아나고, 통신은 블록 수만큼으로 줄어든다. 블록 내부는 안정한 국소 QR로 처리하고, 전체를 두 번 돌리는 BCGS2가 $O(u)$ 를 회복한다.
• TSQR(Tall-Skinny QR, Demmel–Grigori–Hoemmen–Langou 2012). $m \gg n$ 인 세로로 긴 행렬을 행 블록으로 쪼개 프로세서마다 국소 하우스홀더 변환 QR을 돌리고, 나온 $R$ 들을 둘씩 쌓아 다시 QR하는 이진 트리로 합친다. 통신이 $\log_2 p$ 회로 끝나 통신 최적이고, 전 과정이 하우스홀더라 안정성은 무조건 $O(u)$ 다. s-step 크릴로프와 통신 회피 알고리즘 계열의 기본 부품.

7. 촐레스키와의 뒷문 — CholeskyQR

$A = QR$ 의 양변에 $A^{\top}$ 를 곱하면 $A^{\top}A = R^{\top}Q^{\top}QR = R^{\top}R$ 이다. 즉 $R$ 은 그람 행렬 $A^{\top}A$ 의 촐레스키 분해 인자와 같다. 여기서 알고리즘 하나가 바로 나온다: $G = A^{\top}A$ 를 만들고($2mn^2$, 통신 1회), $R = \mathrm{chol}(G)$ 를 구하고($n^3/3$), $Q = AR^{-1}$ 로 끝낸다. 전부 BLAS-3에 동기화 한 번 — GPU에서 압도적으로 빠르다.

대가는 정직하다. $A^{\top}A$ 를 만드는 순간 조건수가 제곱되므로 직교성 손실이 $O(u\kappa^2)$, 정확히 CGS와 같은 등급이다. 게다가 $\kappa(A) \gtrsim u^{-1/2} \approx 10^8$ 이면 $G$ 가 수치적으로 양의 정부호가 아니게 되어 촐레스키 자체가 실패한다.⁴ 처방도 같다 — CholeskyQR2, 즉 두 번 돌리면 $\kappa \lesssim u^{-1/2}$ 범위에서 $O(u)$ 를 회복한다. 그 너머를 노린다면 $G$ 의 대각에 작은 시프트를 더해 촐레스키를 살려낸 뒤 세 번 돌리는 shifted CholeskyQR3가 있다.

결국 CGS의 $\kappa^2$, 정규방정식의 $\kappa^2$, CholeskyQR의 $\kappa^2$ 는 전부 같은 병이다. 그람 행렬을 명시적으로 만들었거나, 만든 것과 다름없는 계산을 했거나.

8. 관련 문서

• QR 분해 · 하우스홀더 변환 · 기븐스 회전 · 촐레스키 분해
• 최소자승법 · 조건수 · 후진 오차 해석 · 부동소수점 연산
• 크리로프 부분공간법 · 아놀디 알고리즘 · 란초스 알고리즘 · 삼중대각화
• 특이값 분해 · 쌍대각화 · LAPACK · 병렬 컴퓨팅

9. Footnotes

1. 요르겐 그람은 1883년 최소자승 논문에서, 에르하르트 슈미트는 1907년 적분방정식 논문에서 이 절차를 썼는데, 슈미트 본인이 각주에 “그람의 방법과 본질적으로 같다”고 적어 두었다. 그런데 라플라스가 1812년과 1816년에 이미 같은 계산을 하고 있었고, 심지어 라플라스 쪽 형태가 오늘날 “수정” GS라 부르는 그것이다. 원본이 수정판이고 고전판이 나중이라는 뜻. 스티글러의 명명 법칙 표본실에 전시해도 될 사례다. ↩

2. Parlett, The Symmetric Eigenvalue Problem(1980)에 “twice is enough”가 케이헌의 말로 실려 있다. 세 번째는 왜 안 하냐고? 필요 없어서다. 두 번 돌리고도 노름이 붕괴하면 그건 반올림 문제가 아니라 그 벡터가 애초에 이전 부분공간 안에 들어 있다는 뜻이므로, 난수 벡터로 갈아끼우고 랭크 결핍을 보고하는 것이 옳다. ↩

3. 이걸 처음 보면 사기처럼 들린다. “직교성이 깨지는데 답은 맞다”니. 하지만 GMRES가 최종적으로 요구하는 것은 $\hat Q$ 의 직교성 자체가 아니라 최소자승 잔차의 정확도이고, 그 둘이 같은 양의 역수 관계로 묶여 있다는 것이 정리의 내용이다. 반대로 말하면 $\hat Q$ 를 GMRES 밖으로 꺼내 다른 계산에 재활용하는 순간 이 보증은 사라진다. ↩

4. $u^{-1/2} \approx 10^8$ 이 마의 선인 이유는 산수 한 줄이다. $\kappa^2 > u^{-1}$ 이 되는 순간 $A^{\top}A$ 는 부동소수점 위에서 특이행렬과 구별되지 않는다. “정규방정식 쓰지 마라”는 잔소리와 정확히 같은 계산이며, 조건수 $10^8$ 짜리 행렬을 들고 오는 사람이 이 두 잔소리를 동시에 듣게 되는 이유다. ↩