파데 근사

1. 개요[편집]

파데 근사 Padé approximation
대상	유리함수 $r_{L/M} = p_L/q_M$ 로 함수 근사
정합 조건	테일러 계수를 $z^{L+M}$ 차까지 일치
계수 결정	$M \times M$ 퇴플리츠 선형계 한 번
대표 응용	스케일링-제곱 expm · A-안정 암시적 RK · 모멘트 정합 축소
대표 병리	프루아사르 이중항 (spurious pole–zero pair)
기원	Frobenius(1881) 정리화, Padé 박사학위논문(1892)

다항식은 무한대에서 반드시 발산한다. 근사하려는 함수가 그렇지 않다면, 애초에 도구를 잘못 고른 것이다.

파데 근사(Padé approximation)는 함수의 테일러 계수만 가지고 분자 차수 $L$, 분모 차수 $M$ 의 유리함수를 만들어, 그 유리함수의 테일러 전개가 원래 급수와 $z^{L+M}$ 차까지 일치하도록 하는 근사법이다. 기호로는 $[L/M]$ 로 쓴다. 정보량은 테일러 계수 $L+M+1$ 개로 똑같은데, 그 정보를 다항식이 아니라 분수에 담는다는 것 하나만 다르다.

그 한 끗 차이의 대가가 크다. 다항식은 $z \to \infty$ 에서 반드시 발산하고, 수렴반경도 가장 가까운 특이점까지로 못 박힌다. 반면 유리함수는 분모의 근으로 특이점을 흉내 낼 수 있다. 근사하려는 함수가 극점을 갖거나, 무한대에서 유한한 값으로 가거나, 절단면(branch cut)을 갖는다면 유리 근사가 같은 계수 예산으로 몇 자릿수 더 정확해진다. 테일러 급수가 “근처에서만 잘 맞는 국소 근사”라면 파데는 “그 급수가 알고 있는 정보를 최대한 짜내 멀리까지 밀고 가는 근사”다.¹

CAE 하는 사람이 파데를 만나는 자리는 대체로 셋이다. expm 안쪽의 스케일링-제곱, 암시적 시간적분법의 안정함수, 그리고 대형 선형계의 축소차수모델. 셋 다 아래에서 다룬다.

2. 정의와 계수 방정식[편집]

$f(z) = \sum_{k \ge 0} c_k z^k$ 가 주어졌다고 하자. $[L/M]$ 파데 근사는

$$r_{L/M}(z) = \frac{p(z)}{q(z)} = \frac{a_0 + a_1 z + \cdots + a_L z^L}{1 + b_1 z + \cdots + b_M z^M}, \qquad f(z) - r_{L/M}(z) = O\!\left(z^{L+M+1}\right)$$

로 정의한다. $q(0) = 1$ 은 정규화일 뿐이다(분자·분모를 동시에 상수배 해도 같은 함수이므로). 미지수는 $a_0,\dots,a_L$ 과 $b_1,\dots,b_M$ 해서 $L+M+1$ 개, 조건도 $L+M+1$ 개 — 개수가 딱 맞는다.

문제를 선형으로 만드는 요령이 있다. 위 식을 그대로 두면 비선형이지만, 양변에 $q$ 를 곱해 $f(z)q(z) - p(z) = O(z^{L+M+1})$ 로 바꾸면 계수에 대해 선형이다. $z^{L+1}, \dots, z^{L+M}$ 항을 모으면 분모 계수만의 퇴플리츠 선형계가 나온다($k < 0$ 이면 $c_k = 0$).

$$\sum_{j=1}^{M} b_j\, c_{L+i-j} = -\,c_{L+i}, \qquad i = 1, \dots, M$$

이걸 풀어 $b_j$ 를 얻고 나면 분자는 대입만 하면 된다.

$$a_i = \sum_{j=0}^{\min(i,M)} b_j\, c_{i-j}, \qquad i = 0, \dots, L$$

계산량은 $O(M^3)$ (퇴플리츠 구조를 이용하면 $O(M^2)$)이라 거의 공짜다. 다만 이 선형계가 특이할 수 있다는 점이 중요하다. 선형화된 조건은 항상 해를 갖지만, 그 해가 원래의 $O(z^{L+M+1})$ 조건을 만족한다는 보장은 없다. 이 결함(defect)이 바로 다음 절의 블록 구조를 낳는다.

3. 파데 표와 블록 구조[편집]

$L$ 을 행, $M$ 을 열로 놓고 $[L/M]$ 을 전부 늘어놓은 무한 배열이 파데 표(Padé table)다. 0열($M=0$)은 그냥 테일러 부분합이고, 대각선 $L=M$ 이 실무에서 가장 많이 쓰인다.

$[L/M]$ 이 “정상”(normal)인지는 한켈 행렬식 $\det\!\left[c_{L-M+i+j}\right]_{i,j=0}^{M-1}$ 이 0이 아닌지로 판정된다. 이 값이 0인 자리에서는 표의 항목이 중복되는데, 중복 항목이 흩어지지 않고 정사각형 블록을 이룬다는 것이 파데의 블록 정리다. 알기 쉬운 예가 $f(z) = e^{z^2}$ 처럼 홀수차 계수가 전부 0인 함수 — 표 전체가 $2\times2$ 블록으로 타일링된다.

대각선이나 계단 모양으로 표를 훑어 내려가는 수열은 연분수의 근사분수(convergent)와 일대일로 대응한다. 파데 표·연분수·직교다항식의 3항 점화식이 사실상 같은 대상을 세 가지 언어로 부르는 것이며, 직교다항식의 분모가 곧 파데 분모라는 사실이 뒤에 나올 파데-란초스 연결의 뿌리다.

4. 유리함수가 이기는 이유, 그리고 지지 않는 방법[편집]

$f(z) = \log(1+z)$ 를 보자. 테일러 급수의 수렴반경은 $z = -1$ 의 특이점 때문에 정확히 1이다. $z = 3$ 을 넣으면 급수는 그냥 발산한다. 그런데 대각 파데 $[M/M]$ 은 극점과 영점을 절단면 $(-\infty, -1]$ 위에 촘촘히 늘어놓아 절단면을 흉내 내고, 그 결과 잘린 복소평면 전체에서 수렴한다. 계수는 똑같이 $2M+1$ 개만 썼는데 수렴 영역이 원반에서 평면 전체로 바뀐 것이다.

$e^{-x}$ 를 $[0,\infty)$ 에서 근사하는 문제는 더 극적이다. 차수 $n$ 다항식으로는 애초에 불가능하고(다항식은 발산한다), $[n/n]$ 유리함수의 최량 근사 오차는 $H^{-n}$ 꼴로 기하급수적으로 줄어든다. 여기서 $H \approx 9.28903$ 은 알펜 상수(Halphen constant)로, 한때 ”$1/9$ 추측”으로 불리다가 곤차르-라흐마노프(1989)가 정확한 값을 확정했다.² 대각 파데는 이 최량 근사만큼은 아니지만 같은 지수 속도를 낸다 — 레머즈 알고리즘으로 최량 근사를 따로 구할 필요가 거의 없는 이유다.

물론 공짜는 아니다. 파데 표 전체가 수렴한다는 정리는 없다. 너톨-폼머렌케 정리는 특이점 집합의 용량(capacity)이 0인 함수에 대해 대각 수열이 용량 수렴한다고만 말한다 — 예외 집합이 작아질 뿐 사라지지는 않는다. “적어도 부분수열은 국소균등 수렴한다”는 베이커-감멜-윌스 추측마저 2003년 루빈스키의 반례로 무너졌다. 파데는 강력하지만 수렴은 여전히 신에게 맡기는 구석이 있다.

4.1. 프루아사르 이중항[편집]

실무에서 진짜 문제가 되는 건 이론적 수렴이 아니라 잡음이다. 계수 $c_k$ 에 반올림 수준의 잡음 $\varepsilon$ 만 섞여도, 파데 근사는 그 잡음을 설명하려고 가짜 극점을 만들어낸다. 그런데 그 극점 바로 옆에는 거의 같은 자리에 영점이 붙어 나와 서로 상쇄한다 — 이 극점-영점 쌍이 프루아사르 이중항(Froissart doublet)이다. 잔차만 보면 근사는 멀쩡해 보이는데, 극점 위치를 물리적으로 해석하려는 순간(공진 주파수, 시스템 극점) 존재하지 않는 모드가 목록에 끼어든다.

증상도 특징적이다. 차수 $M$ 을 올릴수록 정확도가 좋아지다가 어느 지점부터 극점 개수만 늘고 정확도는 정체한다.³ 표준 처방은 곤네-귀텔-트레페텐(2013)의 강건 파데 — 퇴플리츠 계를 그냥 푸는 대신 특이값 분해로 수치적 계수(rank)를 판정해 잡음 차원을 잘라내고 공통 인수를 소거한다. 요즘은 파데 표를 아예 우회해 보간점을 탐욕적으로 고르는 AAA 알고리즘이 표준으로 자리 잡는 추세다. 교훈은 같다 — 부동소수점 연산에서 파데의 고차 항목은 데이터가 아니라 잡음을 적합한다.

5. $e^z$ 의 파데 근사와 A-안정성[편집]

$f(z) = e^z$ 의 파데 근사는 닫힌 형이 있다.

$$p_{L}(z) = \sum_{j=0}^{L} \frac{(L+M-j)!\,L!}{(L+M)!\,j!\,(L-j)!} z^{j}, \qquad q_{M}(z) = \sum_{j=0}^{M} \frac{(L+M-j)!\,M!}{(L+M)!\,j!\,(M-j)!} (-z)^{j}$$

이고 오차의 선두항은 $(-1)^{M}\dfrac{L!\,M!}{(L+M)!\,(L+M+1)!}z^{L+M+1}$ 이다. 가장 작은 대각 항목 $[1/1]$ 은 $\dfrac{1+z/2}{1-z/2}$ — 그렇다, 크랭크-니콜슨법(사다리꼴 법칙)의 증폭인자이자 케일리 변환 그 자체다.

여기서 시간적분 이론과 정면으로 만난다. 한 스텝 방법의 안정함수 $R(z)$ 가 좌반평면 전체에서 $|R(z)| \le 1$ 이면 A-안정이라 부르는데, $e^z$ 의 파데 근사 중 어느 것이 A-안정인지가 완전히 알려져 있다.

엘 정리. $e^z$ 의 $[L/M]$ 파데 근사가 A-안정일 필요충분조건은 $M-2 \le L \le M$ 이다.

즉 파데 표의 대각선과 그 아래 두 개의 부대각선만 A-안정이다. 엘(Ehle)이 1969년 박사논문에서 충분성을 보이고 “그 밖에는 없다”고 추측했으며, 1978년 바너-하이러-뇌르셋이 오더 스타(order star)라는 기하학적 도구를 발명해 필요성까지 증명했다.⁴ 이 정리가 암시적 룽게-쿠타법의 지형도를 그대로 결정한다.

방법	안정함수	차수	성질
$s$-단 가우스-르장드르	대각 $[s/s]$	$2s$	A-안정, 심플렉틱
$s$-단 라다우 IIA	부대각 $[s-1/s]$	$2s-1$	L-안정
$s$-단 로바토 IIIC	부대각 $[s-2/s]$	$2s-2$	L-안정

대각 파데는 $z \to -\infty$ 에서 $|R| \to 1$ 이라 L-안정이 아니다. 극도로 강성한 성분이 감쇠하지 않고 스텝마다 부호만 뒤집으며 살아남는다는 뜻이고, 강성 방정식에서 사다리꼴 법칙이 “안정한데 진동한다”는 악명을 얻는 이유가 정확히 이것이다. 그래서 진짜 강성 문제에는 라다우 IIA 같은 부대각 쪽을 쓴다.

행렬 지수함수 계산의 표준 알고리즘도 여기서 나온다. $e^{A} = \left(e^{A/2^{s}}\right)^{2^{s}}$ 로 노름을 줄인 뒤 대각 파데 $r_{m}(X) = q_m(X)^{-1}p_m(X)$ 를 적용하고 다시 제곱하는 것. 하이엄(2005)의 구현은 차수 13 대각 파데에 $\theta_{13} \approx 5.37$ 문턱값을 쓴다. 분모가 행렬이므로 “나눗셈”이 선형계 풀이 한 번이라는 점만 스칼라와 다르다 — 자세한 것은 해당 문서 참고.

6. 모멘트 정합과 모델 차수 축소[편집]

마지막 무대는 축소차수모델이다. 대형 선형계 $E\dot x = Ax + bu$, $y = c^{\top}x$ 의 전달함수 $H(s) = c^{\top}(sE-A)^{-1}b$ 를 $s_0$ 근방에서 전개한 계수를 모멘트라 부른다. 차원 $q$ 짜리 축소 모델이 모멘트를 $2q$ 개까지 맞추면, 그 축소 전달함수는 정확히 $H$ 의 $[q-1/q]$ 파데 근사다.

문제는 모멘트를 직접 계산하면 안 된다는 것. $m_k$ 는 본질적으로 $(A^{-1}E)^k$ 를 벡터에 반복 작용시키는 거듭제곱 반복이라, $k$ 가 조금만 커져도 모든 모멘트가 지배 고유벡터 방향으로 붕괴한다. 1990년대 초 회로 시뮬레이션에서 유행한 AWE(asymptotic waveform evaluation)가 이 벽에 부딪혀 차수 10 언저리에서 무너졌다.

해법이 파데-란초스(Padé via Lanczos, 펠트만-프로인트 1995)다. 모멘트를 명시적으로 만들지 않고 란초스 알고리즘의 양측(two-sided) 버전으로 삼중대각 $T_q$ 를 뽑으면, $c^{\top}(sE-A)^{-1}b$ 의 $[q-1/q]$ 파데 근사가 $T_q$ 의 분해행렬에서 곧바로 나온다. 같은 답을 크리로프 부분공간법의 직교화로 얻으니 수치적으로 안정하고, 차수 수십~수백까지 밀 수 있다. SPICE 계열 상용 회로 시뮬레이터의 상호연결선(interconnect) 축소, 반도체 소자 시뮬레이션의 소신호 해석, MEMS·구조 진동의 주파수 응답 축소가 전부 이 계열이다. 전개점을 여러 개 잡아 각 점에서 몇 개씩 모멘트를 맞추는 다점 파데(유리 크릴로프)는 넓은 주파수 대역을 한 모델로 덮는 표준 기법이 됐다.

7. 관련 문서[편집]

• 테일러 급수 · 보간과 근사 · 직교다항식 · 연분수
• 행렬 지수함수 · 행렬함수 · 지수 적분기
• 강성 방정식 · 룽게-쿠타법 · 크랭크-니콜슨법 · 암시적 룽게-쿠타법
• 축소차수모델 · 란초스 알고리즘 · 크리로프 부분공간법
• 특이값 분해 · 조건수 · 부동소수점 연산

8. Footnotes[편집]

1. 앙리 파데(Henri Padé)의 1892년 박사논문 지도교수가 에르미트였다. 다만 표의 구조 자체는 1881년 프로베니우스가 먼저 정리했고, 근사분수 아이디어는 오일러·야코비까지 거슬러 올라간다. 수학사의 국룰대로 이름은 정리한 사람이 아니라 표를 예쁘게 그린 사람에게 갔다. ↩

2. ”$1/9$ 추측”은 수치실험으로 얻은 오차비가 $0.1076\ldots$ 여서 “대충 $1/9 = 0.111\ldots$ 아니냐”고 부른 데서 나왔다. 실제 값은 $1/9.28903\ldots = 0.10765\ldots$ 로, 9와 9.28903 사이의 그 미세한 틈을 메우는 데 20년 가까이 걸렸다. 수치해석자가 상수 하나에 얼마나 집착할 수 있는지 보여주는 사례. ↩

3. 파데 근사를 처음 써보면 열에 아홉은 “차수를 올렸더니 그래프에 뾰족한 못이 박혔다”는 경험을 한다. 그게 프루아사르 이중항이 상쇄에 실패해 극점만 노출된 상태다. 차수를 낮추거나 강건 파데를 쓰면 사라지는데, 초심자는 보통 그 전에 파데를 포기한다. ↩

4. 오더 스타는 $\left|R(z)e^{-z}\right| > 1$ 인 영역을 복소평면에 칠해서 나오는 별 모양 그림이다. 근사의 차수는 원점에 모인 꽃잎 개수로, 안정성은 꽃잎이 좌반평면을 침범하는지로 읽힌다. 해석적으로 40년 안 풀리던 추측이 “그림을 그렸더니 보였다”로 끝난, 이 바닥에서 손꼽히는 사건이다. ↩