슈어-파레 방법

1. 개요

슈어-파레 방법 Schur–Parlett algorithm
목적	일반 f 에 대한 행렬함수 f(A) 계산
1단계	슈어 분해 A = QTQH
2단계	상삼각 T 에 파레 점화식 (Parlett, 1974)
급소	분모 tii − tjj — 가까운 고유값
현대 처방	데이비스-하이엄(2003) 블록화 + 블록 내부 테일러
구현	MATLAB funm · 비용 ≈ 28n³ flops (블록이 작을 때)

$e^{A}$ 를 구하는 방법은 열아홉 가지쯤 알려져 있다. 그런데 $f$ 가 아무거나일 때는?

슈어-파레 방법(Schur–Parlett algorithm)은 정사각행렬 $A$ 와 임의의 스칼라 함수 $f$ 에 대해 행렬함수 $f(A)$ 를 계산하는 범용 알고리즘이다. 뼈대는 두 줄로 요약된다.

$$A = QTQ^{H} \;\Longrightarrow\; f(A) = Q\,f(T)\,Q^{H}$$

즉 슈어 분해로 문제를 상삼각행렬의 함수로 환원한 뒤, $f(T)$ 를 파레 점화식으로 채운다. 이 경로를 고르는 이유는 명확하다 — 유니터리 변환은 조건수가 1이라 오차를 증폭하지 않고, 조르당 표준형과 달리 결함 행렬에서도 안전하게 계산되며, 무엇보다 $f$ 가 무엇인지 묻지 않는다.

$e^{A}$ 에는 스케일링-제곱, $\log A$ 에는 역스케일링-제곱, $A^{1/2}$ 에는 뵈르크-함마를링 같은 전용 알고리즘이 있고 그쪽이 항상 더 낫다. 슈어-파레는 전용 알고리즘이 없는 나머지 전부를 담당하는 최후의 보루다. $\cos A$, $\operatorname{sinc}(A)$, 지수 적분기의 $\varphi_{k}$, 사용자가 방금 정의한 이상한 $f$ — 이런 것들을 계산할 수 있는 유일한 일반 절차다.

2. 파레 점화식

$F = f(T)$ 도 상삼각이고, 행렬함수의 기본 성질에 따라 $F$ 는 $T$ 와 교환한다.

$$FT = TF$$

이 한 줄이 알고리즘 전부다. 대각 원소는 $f_{ii} = f(t_{ii})$ 로 즉시 나오고, $(i,j)$ 성분 ($i<j$)에서 위 등식을 풀어 쓰면

$$f_{ij} = \frac{t_{ij}\,(f_{ii} - f_{jj}) + \displaystyle\sum_{k=i+1}^{j-1}\big(f_{ik}t_{kj} - t_{ik}f_{kj}\big)}{t_{ii} - t_{jj}}$$

가 된다(파레, 1974). 우변에 등장하는 $f_{ik}, f_{kj}$ 는 전부 $(i,j)$ 보다 더 짧은 대각선 위의 성분이므로, 주대각에서 시작해 첫 번째 초대각선, 두 번째 초대각선, … 순서로 한 겹씩 채워 나가면 된다. 비용은 $2n^{3}/3$ flops 정도로 슈어 분해($\approx 25n^{3}$)에 비하면 거저다.

이론적으로는 여기서 끝이다. 실제로는 여기서 시작이다.

3. 급소 — 분모가 죽는다

문제는 분모 $t_{ii} - t_{jj}$ 다. 두 고유값이 가까우면 이 값이 0에 접근하고 점화식이 폭발한다. 고유값이 정확히 중복이면 아예 0으로 나누기다.

여기서 짚어야 할 것은, 분모가 작다는 사실이 문제 자체가 나쁘다는 뜻이 아니라는 점이다. $A = \begin{bmatrix}1 & 1\\ 0 & 1+\epsilon\end{bmatrix}$ 의 $f(A)$ 는 $\epsilon \to 0$ 에서 멀쩡히 수렴한다 — 분자의 $f(t_{11})-f(t_{22})$ 도 함께 0으로 가서 극한이 $f'(1)$ 이 되기 때문이다. 즉 참값은 잘 정의되어 있는데 계산 경로가 0/0 을 통과하도록 짜여 있는 상황이고, 이것은 전형적인 자리수 상쇄다. 부동소수점에서는 분자의 유효숫자가 이미 대부분 날아간 뒤 작은 수로 나누므로 상대오차가 그대로 증폭된다.

게다가 이 오염은 국소적이지 않다. 한 성분이 오염되면 그 값이 뒤이은 대각선의 합에 들어가 오른쪽 위로 계속 번진다. 삼각행렬의 왼쪽 아래에서 시작된 작은 실수가 우상단 구석에서 자리수 전체를 잡아먹는 그림이 나온다.

4. 데이비스-하이엄의 블록화

현대 처방은 데이비스와 하이엄(2003)의 것이다. 발상은 단순하다 — 가까운 고유값끼리는 갈라놓지 말고 한 블록에 몰아넣고, 블록 안은 점화식 대신 다른 방법으로 처리한다.

① 고유값 군집화. $T$ 의 대각 원소를 분리 파라미터 $\delta$(기본값 $0.1$)로 묶는다. $|t_{ii} - t_{jj}| \le \delta$ 이면 같은 군집으로 보내고, 이 관계를 연쇄적으로 확장한다(사슬처럼 이어지면 한 군집). 결과적으로 서로 다른 군집의 고유값은 $\delta$ 이상 떨어져 있고, 같은 군집 안은 가까울 수도 있다.

② 순서화 슈어 분해. 같은 군집의 고유값이 대각에서 연속으로 붙어 있도록 $T$ 를 재정렬한다. 인접 블록을 맞바꾸는 유니터리 상사변환을 버블 정렬하듯 반복하는 표준 연산이고(LAPACK dtrsen), 이 단계 자체는 후진 안정하다.

③ 대각 블록은 테일러로. 블록 $T_{ii}$ 의 고유값은 전부 반경 $\delta$ 안에 모여 있으므로, 중심 $\sigma = \operatorname{tr}(T_{ii})/m_{i}$ 를 빼면 $T_{ii} - \sigma I$ 의 노름이 작다. 그래서

$$f(T_{ii}) = \sum_{k\ge 0} \frac{f^{(k)}(\sigma)}{k!}\,(T_{ii} - \sigma I)^{k}$$

테일러 급수를 직접 평가한다. 점화식을 아예 쓰지 않으므로 분모 문제가 사라진다.

④ 블록 사이는 실베스터 방정식으로. 블록 단위로 다시 쓴 파레 점화식은

$$T_{ii}F_{ij} - F_{ij}T_{jj} = F_{ii}T_{ij} - T_{ij}F_{jj} + \sum_{k=i+1}^{j-1}\big(F_{ik}T_{kj} - T_{ik}F_{kj}\big)$$

이라는 실베스터 방정식 이 된다. 이 방정식은 $T_{ii}$ 와 $T_{jj}$ 가 공통 고유값을 갖지 않을 때 유일해를 가지며, 해의 민감도는 두 스펙트럼 사이의 분리도 $\operatorname{sep}(T_{ii}, T_{jj})$ 가 결정한다. ①에서 군집을 $\delta$ 만큼 떼어 놓은 것이 정확히 이 $\operatorname{sep}$ 을 확보하기 위한 조치였다. 양변이 이미 삼각이므로 크로네커 곱을 만들 필요 없이 블록 후진대입으로 푼다.

$\delta$ 는 두 오차 사이의 저울이다. 키우면 군집이 커져 블록 크기가 늘고 $\|T_{ii}-\sigma I\|$ 가 커져 테일러가 나빠진다. 줄이면 블록은 작아지지만 $\operatorname{sep}$ 이 작아져 실베스터 풀이가 나빠진다. $\delta = 0.1$ 은 실험적으로 잡은 타협점이지 이론적 최적값이 아니다.¹

5. 비용과 한계

블록이 전부 $1\times1$ 로 갈라지는 좋은 경우 총 비용은 $\approx 28n^{3}$ flops — 슈어 분해 $25n^{3}$ + 점화식 $2n^{3}/3$ + 잡비. 대각화 가능하고 고유값이 잘 흩어진 행렬에서는 $f$ 가 무엇이든 이 값이다.

문제는 군집이 커질 때다. 극단적으로 모든 고유값이 한 점에 몰린 행렬(예: 단일 조르당 블록)이면 군집이 하나, 블록이 $n\times n$ 이 되어 $n\times n$ 행렬의 테일러 급수를 통째로 돌려야 한다. 비용이 $n^{4}$ 급으로 뛰고, $\|T-\sigma I\|$ 가 크면 급수 수렴도 느려지고 정확도도 떨어진다. 비정규성이 강한 행렬이 슈어-파레의 진짜 적이라는 뜻이다.

실무적 제약이 하나 더 있다. $f$ 의 고차 도함수를 사용자가 제공해야 한다. MATLAB funm(A, fun) 은 fun(x,k) 형태로 $k$ 차 도함수를 요구하며, 이걸 손으로 못 넣으면 알고리즘이 성립하지 않는다.² 이 불편을 없애려는 시도가 근래의 다중정밀도 무도함수(derivative-free) 변형으로, 대각 블록을 고정밀 산술에서 처리해 도함수 요구를 우회한다.

정확도 목표도 겸손하게 잡혀 있다. 일반 $f$ 에 대해 후진 안정한 알고리즘은 알려져 있지 않고 존재 여부도 불분명하다. 데이비스-하이엄이 내건 목표는 “전진 오차가 문제의 조건수에 비례하는 수준” — 즉 $\|\hat F - f(A)\| \lesssim \operatorname{cond}(f,A)\,u\,\|f(A)\|$ 이고, 이것이 이 분야에서 현실적으로 요구할 수 있는 최선이다.

6. 조건수와 프레셰 도함수

“이 답을 몇 자리까지 믿어도 되는가”에 답하려면 행렬함수의 프레셰 도함수 $L_{f}(A,E)$ 가 필요하다. 상대 조건수는

$$\kappa_{f}(A) = \frac{\|L_{f}(A)\|\,\|A\|}{\|f(A)\|}, \qquad \|L_{f}(A)\| = \max_{E\ne 0}\frac{\|L_{f}(A,E)\|}{\|E\|}$$

이고, 이 값이 $10^{10}$ 이면 배정도로는 여섯 자리쯤만 살아남는 것이 정상이다. 알고리즘을 의심하기 전에 이 숫자부터 봐야 한다.

계산에는 크기 $2n$ 짜리 블록 삼각 항등식이 유용하다.

$$f\!\left(\begin{bmatrix} A & E \\ 0 & A\end{bmatrix}\right) = \begin{bmatrix} f(A) & L_{f}(A,E) \\ 0 & f(A)\end{bmatrix}$$

즉 $f$ 를 계산할 수 있으면 프레셰 도함수도 같은 루틴 한 번 더로 얻는다. 여기에 유도 노름 추정을 위한 하이엄-티슬 1-노름 거듭제곱법을 얹어 $\kappa_{f}$ 를 추정하는 것이 funm_condest1 이 하는 일이다. 다만 이 $2n\times2n$ 행렬은 고유값이 정확히 중복($A$ 의 스펙트럼이 두 번)이므로, 슈어-파레로 처리하면 반드시 큰 블록이 생긴다 — 조건수 추정이 원래 계산보다 비싸고 까다로운 이유다.

7. 언제 쓰고 언제 피하는가

• 쓴다. $f$ 가 전용 알고리즘이 없는 함수일 때. 삼각함수, $\varphi$ 함수, 사용자 정의 해석함수, 분수 거듭제곱의 일반형 등. 행렬이 정규에 가깝고 고유값이 흩어져 있으면 특히 잘 맞는다.
• 피한다. $f$ 가 $\exp$, $\log$, $\sqrt{\cdot}$, $\operatorname{sign}$ 중 하나라면 전용 알고리즘이 언제나 낫다. 행렬 지수함수의 스케일링-제곱은 후진 오차 해석까지 갖춰져 있고, 행렬 로그의 슈어-파데 역스케일링-제곱, 행렬 제곱근의 뵈르크-함마를링, 행렬 부호 함수의 뉴턴 반복이 각각 그렇다.
• 아예 다른 문제로 바꾼다. $A$ 가 대형 희소행렬이면 슈어 분해 자체가 불가능하다($O(n^{3})$ 시간, $O(n^{2})$ 저장, 게다가 $f(A)$ 는 조밀). 이때 필요한 것은 대개 $f(A)b$ 하나이므로 크리로프 부분공간법이나 등고선 적분으로 갈아탄다. 슈어-파레는 어디까지나 조밀·중소형 알고리즘이다.

이름이 붙은 두 사람 중 파레(Beresford Parlett)는 대칭 고유값 문제의 교과서로 더 유명하고, 이 점화식은 1974년의 짧은 논문에서 나왔다. 슈어는 1909년에 분해를 증명했을 뿐 행렬함수와는 무관했다 — 65년 시차의 두 결과가 한 알고리즘 이름으로 묶인 셈이다.³

8. 관련 문서

• 행렬함수 · 행렬 지수함수 · 행렬 제곱근 · 행렬 부호 함수 · 행렬 로그
• 슈어 분해 · 조르당 표준형 · 비정규 행렬 · QR 알고리즘
• 실베스터 방정식 · 바텔스-스튜어트 알고리즘 · 리아푸노프 방정식
• 조건수 · 후진 오차 해석 · 부동소수점 연산 · 테일러 급수 · 파데 근사
• LAPACK · MATLAB · 지수 적분기 · 크리로프 부분공간법

9. Footnotes

1. 이런 상수는 논문의 표에 딱 한 줄 나오고 마는데, 실제로는 알고리즘의 성패를 좌우한다. $\delta$ 를 바꿔 가며 같은 행렬을 돌려 보면 오차가 U자 곡선을 그리는 것이 보인다. “기본값을 건드리지 마세요”라는 조언이 대개 옳지만, 왜 그 값인지는 알고 있는 편이 좋다. ↩

2. 도함수를 넘겨야 한다는 요구는 생각보다 자주 사람을 막는다. $f$ 가 스플라인 보간이거나 다른 코드가 뱉어 주는 블랙박스면 $k$ 차 도함수라는 물건이 애초에 없다. 그럴 때 “수치미분으로 대충 넣자”는 유혹이 오는데, 테일러 급수의 고차 계수를 유한차분으로 만들면 상쇄오차가 그대로 답에 들어가므로 대개 안 하느니만 못하다. ↩

3. 파레의 점화식은 원래 “행렬함수를 삼각화해서 구하자”는 발상 자체가 새롭던 시절의 것이라, 논문에서 이미 “가까운 고유값이 있으면 곤란하다”고 스스로 경고하고 블록 버전까지 제시해 뒀다. 문제를 아는 것과 실용적으로 해결하는 것 사이에 30년이 걸렸다는 점이, 수치해석에서 “아이디어는 절반도 아니다”라는 말의 좋은 예시다. ↩