젠킨스-트라우브 알고리즘

1. 개요

젠킨스-트라우브 알고리즘(Jenkins–Traub algorithm)¹은 다항식 $P(z)$ 의 근을 하나씩 찾아 수축(deflation)해 나가는 3단계 가변 시프트 반복법이다. 마이클 젠킨스의 스탠퍼드 박사논문(1969)과 조지프 트라우브와 함께 낸 1970년 두 논문에서 나왔고, 복소 계수판 CPOLY 와 실계수판 RPOLY 가 각각 ACM 알고리즘 419·493 으로 공개되면서 1970~90년대 범용 다항식 근 찾기의 사실상 표준이 됐다. IMSL 같은 상용 라이브러리도 실계수 다항식 근을 이 방법으로 푼다.

핵심 착상은 $P$ 대신 그 보조 다항식 $H(z)$ 를 반복적으로 갱신하는 것이다. $H$ 는 차수 $n-1$ 이고, 반복이 진행되면 $P(z)/(z-\alpha_1)$ — 즉 근 $\alpha_1$ 하나만 뺀 나머지 — 의 상수배로 수렴한다. 그러면 $\alpha_1 \approx s - P(s)/\bar H(s)$ 로 근이 바로 읽힌다. 이 구조가 동반 행렬에 대한 역반복·레일리 몫 반복과 정확히 같다는 것이 원 논문 제목의 “generalized Rayleigh iteration” 이다.

오늘날 NumPy·MATLAB 의 roots 는 동반 행렬 고유값으로 가지만, 젠킨스-트라우브는 여전히 $O(n^2)$ 이고, 근을 작은 것부터 순서대로 내주고, 실계수 다항식을 실수 연산만으로 푸는 실용적 장점 때문에 신호처리·제어 코드 곳곳에 남아 있다.

2. H 다항식 — 근 하나만 빼기

$P$ 를 모닉 차수 $n$, 단순근 $\alpha_1,\dots,\alpha_n$ 이라 하자. 출발은 $H^{(0)} = P'$ 이다. 갱신식은 시프트 $s_\lambda$ 에 대해

$$H^{(\lambda+1)}(z) = \frac{1}{z - s_\lambda}\left[ H^{(\lambda)}(z) - \frac{H^{(\lambda)}(s_\lambda)}{P(s_\lambda)}\,P(z) \right]$$

이다. 괄호 안은 $z = s_\lambda$ 에서 0이므로 나눗셈이 정확히 떨어지고(호너식 조립제법 한 번), 결과는 다시 차수 $n-1$ 이다.

왜 이게 통하는지는 $H$ 를 부분분수 기저로 펼치면 보인다. $P'(z) = \sum_i P(z)/(z-\alpha_i)$ 이므로 $H^{(0)} = \sum_i c_i^{(0)}\,\frac{P(z)}{z-\alpha_i}$ ($c_i^{(0)} = 1$) 이고, 갱신 한 번은 각 계수를

$$c_i^{(\lambda+1)} = \frac{c_i^{(\lambda)}}{\alpha_i - s_\lambda}$$

로 바꾼다. 즉 시프트에 가장 가까운 근의 계수가 가장 빨리 커진다. 이것은 행렬 $C$(동반 행렬)에 $(C - s I)^{-1}$ 을 반복 적용하는 역반복법과 같은 식이다. 계수 하나가 나머지를 압도하면 $H$ 를 최고차 계수로 정규화한 $\bar H$ 가 $P(z)/(z-\alpha_1)$ 에 가까워지고,

$$t = s - \frac{P(s)}{\bar H(s)}$$

가 $\alpha_1$ 의 추정이 된다. $\bar H = P'$ 이면 이것은 그냥 뉴턴-랩슨법이니, 젠킨스-트라우브는 “뉴턴법의 분모를 점점 똑똑한 것으로 바꿔 가는 방법” 이라고도 읽을 수 있다.

3. 세 단계

3.1. 1단계 — 무시프트

$s = 0$ 으로 몇 번(CPOLY 는 5번) 갱신한다. 계수가 $c_i \propto \alpha_i^{-\lambda}$ 로 바뀌므로 절댓값이 작은 근 쪽 성분이 키워진다. 목적은 수렴이 아니라 $H$ 에서 큰 근 성분을 걸러내 2단계를 준비하는 것이다.

3.2. 2단계 — 고정 시프트

시프트를 $s = \beta e^{i\theta}$ 로 고정해 반복한다. $\beta$ 는 코시 다항식 $|a_n|z^n + \cdots + |a_1|z - |a_0|$ 의 양의 근으로, 모든 근의 절댓값의 하한이다. 그러니 $s$ 는 원점 근처 근들 바로 안쪽 원 위에 놓이고, 2단계는 그 원 위에서 $s$ 에 가장 가까운 근 하나를 골라낸다. 추정값 $t$ 가 연속으로 안정되면 3단계로 넘어가고, 정해진 횟수 안에 안정되지 않으면 $\theta$ 를 돌려(CPOLY 는 94°씩) 다시 시도한다.

이 단계의 수렴은 1차다. 대신 전역적이다 — 고정 시프트 역반복은 가장 가까운 고유값이 유일하기만 하면 출발점과 무관하게 수렴한다. 뉴턴법이 출발점을 잘못 고르면 엉뚱한 데서 헤매는 문제를 이 단계가 막아 준다.

3.3. 3단계 — 가변 시프트

이제 매 반복마다 시프트를 최신 추정으로 바꾼다.

$$s_{\lambda+1} = s_\lambda - \frac{P(s_\lambda)}{\bar H^{(\lambda+1)}(s_\lambda)}$$

시프트를 바꾸는 역반복이 레일리 몫 반복이 되듯, 수렴이 초선형으로 바뀐다. 다중 정밀도(300자리)로 직접 구현해 오차 $e_k$ 의 비 $\ln e_{k+1}/\ln e_k$ 를 재 보면 이렇다.

반복	근 다섯 개짜리 복소 다항식	윌킨슨 $W_{10}$
1	3.26	2.91
2	2.78	2.66
3	2.67	2.624
4	2.637	2.619
5	2.625	— (300자리 바닥)

두 다항식 모두 $1 + \varphi = \varphi^2 = 2.618\ldots$ 로 다가간다($\varphi$ 는 황금비). 오차가 $8\times10^{-4} \to 2.6\times10^{-9} \to 1.1\times10^{-23} \to 3.1\times10^{-61}$ 으로 떨어지니, 배정밀도라면 3단계 진입 후 3~4번이면 끝난다. 이 차수는 뉴턴(2차)보다 높지만 반복당 $P$ 평가·$\bar H$ 평가·$H$ 갱신으로 $O(n)$ 작업이 세 번 들므로, 작업량 대비 효율($2.618^{1/3} \approx 1.38$)은 뉴턴($2^{1/2} \approx 1.41$, $P$·$P'$ 두 번)과 비슷한 급이다. 젠킨스-트라우브의 진짜 이득은 속도보다 출발점을 고를 필요가 없다는 데 있다.

실험에서 하나 더 보인 것: 첫 번째 다항식에서 2단계가 고른 근은 절댓값이 가장 작은 $-0.8+0.2i$ ($|\cdot| = 0.82$)가 아니라 $1+0.5i$ ($|\cdot| = 1.12$)였다. 시프트 $0.38+0.32i$ 에서 거리가 각각 1.19, 0.65 였기 때문이다. “작은 근부터” 는 경향이지 보장이 아니고, 정확한 규칙은 “시프트에서 가장 가까운 근부터” 다.

4. 수축과 안정성

근 하나를 찾으면 $P(z) = (z - \alpha)Q(z) + r$ 로 나누어 $Q$ 에 대해 처음부터 다시 한다. 이 수축이 근 찾기의 고질병이다. 찾은 근에 오차가 있으면 $Q$ 의 계수가 오염되고, 그 오차가 뒤에 찾는 근으로 누적된다.

윌킨슨이 보인 규칙은 절댓값이 작은 근부터 찾아 앞쪽(최고차부터의) 조립제법으로 나누면 수축이 안정하다는 것이다. 젠킨스-트라우브의 1·2단계가 작은 근 쪽으로 기울어 있는 것이 바로 이 규칙과 맞물린다 — 알고리즘 설계가 “빠른 수렴” 과 “안정한 수축” 을 같은 장치로 얻는다. 수축된 다항식에서 찾은 근을 원래 $P$ 에 대해 뉴턴 한두 번으로 다듬는(polishing) 것도 표준 마무리다.

그래도 근 자체가 나쁜 조건이면 어떤 알고리즘도 구할 수 없다. 윌킨슨 다항식에서 계수 하나를 $2^{-23}$ 만큼 건드리면 근 10개가 복소평면으로 튀어나가는데, 이건 알고리즘의 결함이 아니라 문제의 조건수다. 중근 $m$ 개짜리는 정확도가 $\varepsilon^{1/m}$ 으로 제한되고 수렴도 느려진다.

5. RPOLY — 실계수는 실수로

실계수 다항식의 복소근은 켤레쌍으로 나온다. 복소 시프트로 하나씩 찾으면 복소 연산이 들고, 수축할 때 켤레까지 정확히 맞추기도 번거롭다. RPOLY 는 선형 인수 $(z - s)$ 대신 이차 인수 $\sigma(z) = z^2 + uz + v$ 로 시프트한다. $H$ 갱신을 이차식으로 나누는 조립제법으로 바꾸면 모든 연산이 실수로 돌고, 한 번에 켤레쌍 하나를 꺼낸다. 3단계에서 이차 인수가 수렴하지 않으면 실근 하나를 노리는 선형 반복으로 바꿔 시도한다. 실계수 문제(제어 시스템의 특성다항식, 필터 설계)에서는 복소 연산을 피하는 만큼 싸고, 결과가 정확히 켤레쌍으로 나온다는 실용적 장점이 있다.

6. 동반 행렬 고유값과 비교

	젠킨스-트라우브	동반 행렬 + QR
비용	근당 $O(n)$ × 반복, 전체 $O(n^2)$	일반 QR 은 $O(n^3)$, 구조 이용 시 $O(n^2)$
메모리	$O(n)$	일반 $O(n^2)$
안정성	경험적으로 우수, 수축 순서에 의존	행렬 의미의 후진 안정(에델만·무라카미 1995)
근 일부만	가능(작은 것부터)	전부 계산
구현	휴리스틱 많음(수렴 판정, 시프트 회전)	LAPACK 호출 한 줄

에델만·무라카미(1995)는 동반 행렬의 QR 이 계수 노름 대비 작은 섭동을 받은 다항식의 정확한 근을 준다는 것을 보였다. 그래서 “검증된 LAPACK 한 줄” 이 휴리스틱 덩어리를 이기고 범용 roots 의 기본값이 됐다. 다만 계수의 크기가 수십 자릿수 차이 나면 노름 기준의 후진 안정은 작은 계수를 지켜 주지 못하니, 이럴 때는 계수 스케일링이나 다른 방법을 고려한다. 2015년 오렌츠·마흐·반데브릴·왓킨스는 동반 행렬의 구조를 살린 $O(n^2)$ 후진 안정 QR 을 내놓아 비용 차이도 지웠다.

다른 경쟁자도 있다. 근 전체를 동시에 갱신하는 애버스 방법(에얼리히-애버스)은 3차 수렴에 병렬화가 쉽고, 임의 정밀도 패키지 MPSolve 의 엔진이다. 실근이 많은 경우 라게르 방법은 거의 어디서 출발해도 수렴하는 성질로 『Numerical Recipes』 의 선택이 됐다. 낮은 차수(≤4)는 닫힌꼴 공식보다 이런 반복법이 수치적으로 더 안전한 경우가 많다.²

7. 쓸 때 알아 둘 것

1. 계수 정규화. 최고차 계수로 나누고, 필요하면 $z \to \gamma z$ 로 스케일해 계수 크기를 고르게 한다. 오버플로·언더플로를 막는 가장 싼 방법이다.
2. 0 근 먼저. 상수항이 0이면 $z$ 를 먼저 인수분해해 빼낸다. 1단계가 $P(0)$ 으로 나누기 때문이다.
3. 중근·근접근은 기대치를 낮춰라. 수렴 판정이 오래 걸리거나 시프트 회전이 반복되면 대개 근이 뭉쳐 있다는 신호다.
4. 고유값으로 오는 문제라면 다항식으로 가지 마라. 행렬 $A$ 의 고유값을 특성다항식의 근으로 구하는 것은 조건수를 일부러 망치는 일이다. 다항식은 원래 다항식으로 주어진 문제(필터 극점, 특성방정식, 프로니 방법의 예측 다항식)에만 쓴다.³

8. 관련 문서

• 동반 행렬 · 역반복법 · 레일리 몫 · QR 알고리즘
• 뉴턴-랩슨법 · 뮬러 방법 · 라게르 방법 · 애버스 방법
• 윌킨슨 다항식 · 호너 방법 · 프로니 방법
• 고유값 문제 · 고유값 조건수

9. Footnotes

1. 트라우브는 이후 정보 기반 복잡도(information-based complexity) 분야를 세운 사람이기도 하다. 반복법의 수렴 차수와 함수 평가 비용을 저울질하는 감각이 이 알고리즘에서 이미 보인다. ↩

2. 3차·4차 근의 공식은 이론적으로는 우아하지만 소거 오차를 사랑한다. 카르다노 공식으로 실근 세 개를 구하려다 복소수 세제곱근을 거쳐야 하는 “기약 경우” 를 만나면 대부분 조용히 반복법으로 갈아탄다. ↩

3. 고유값 → 특성다항식 → 근 찾기는 학부 선형대수 시험의 국룰이지만, 수치해석에서는 금기 1순위다. 반대 방향(다항식 → 동반 행렬 → 고유값)은 오히려 표준이라는 게 웃픈 포인트. ↩