윌킨슨 다항식

1. 개요

윌킨슨 다항식 Wilkinson's Polynomial
정의	$w(x) = (x-1)(x-2)\cdots(x-20)$
발표	윌킨슨 (1959, 『Numerische Mathematik』 창간호)
유명한 섭동	$x^{19}$ 계수 $-210$ 을 $-210 - 2^{-23}$ 으로
결과	근 10개가 켤레 복소수 5쌍으로, 최대 허수부 $\approx 2.81$
가장 민감한 근	$x = 16$, $\partial r/\partial a_{19} \approx -2.4\times10^{9}$
교훈	멱급수 계수는 근을 나타내기에 나쁜 표현이다

근이 1, 2, …, 20 으로 사이좋게 떨어져 있는데, 계수 하나를 1억분의 1 만 건드리면 절반이 복소평면으로 날아간다.

윌킨슨 다항식(Wilkinson’s polynomial)은

$$w(x) = \prod_{i=1}^{20}(x - i) = x^{20} - 210x^{19} + 20615x^{18} - \cdots + 20!$$

이다. 제임스 윌킨슨이 1959년 『Numerische Mathematik』 창간호의 두 편짜리 논문에서 「근이 서로 잘 떨어져 있어도 다항식 근 찾기는 지독하게 악조건일 수 있다」는 것을 보이려고 쓴 예제다. 그 전까지의 통념은 「근이 몰려 있거나 중근일 때만 위험하다」였는데, 이 다항식은 근 사이 간격이 전부 1 이라 그 통념이 틀렸다는 것을 한 방에 보여 준다.

이 다항식이 수치해석에 남긴 교훈은 근 찾기 알고리즘에 대한 것이 아니다. 문제를 멱급수 계수로 표현하는 순간 이미 망가진다는 것이다. 이 문서는 유명한 섭동 실험을 60자리로 재현하고, 근의 조건수를 계산하고, 배정밀도에서는 계수를 저장하는 것만으로도 무슨 일이 생기는지 보이고, 그래서 실무가 어떻게 바뀌었는지(동반 행렬 고유값, 직교 기저)를 정리한다. 고유값 쪽 이야기는 고유값 조건수가 같은 윌킨슨의 이중대각 행렬로 이미 다뤘다.

2. 계수 하나, $2^{-23}$

윌킨슨이 실제로 해 보인 실험은 $x^{19}$ 의 계수 $-210$ 을 $-210 - 2^{-23}$ 으로 바꾸는 것이다. $2^{-23} \approx 1.19\times10^{-7}$ 이고 상대 변화로는 $5.7\times10^{-10}$ 이다. 섭동된 다항식의 근을 60 자리 정밀도로 구하면:

원래 근	섭동 후 근
1 ~ 4	1.000000, 2.000000, 3.000000, 4.000000
5 ~ 9	4.999999928, 6.000006944, 6.999697234, 8.007267603, 8.917250249
10, 11	$10.0952662 \pm 0.6435009\,i$
12, 13	$11.7936339 \pm 1.6523297\,i$
14, 15	$13.9923581 \pm 2.5188301\,i$
16, 17	$16.7307375 \pm 2.8126249\,i$
18, 19	$19.5024394 \pm 1.9403303\,i$
20	20.8469081

윌킨슨 원 논문의 표와 자릿수까지 같다. 실수 근 10 개가 켤레 복소수 5 쌍으로 바뀌었고, 허수부가 2.8 까지 커졌다. 근 사이 간격이 1 인데 근이 2.8 만큼 옆으로 움직였으니 「근 16 이 어디 있느냐」는 질문은 이 정도 섭동에서 의미를 잃는다.

섭동을 줄이면 어떻게 되나. $\delta = 2^{-e}$ 로 바꿔 가며 같은 계산을 하면

$\delta$	$2^{-30}$	$2^{-32}$	$2^{-34}$	$2^{-38}$	$2^{-44}$
복소 근 개수	6	4	0	0	0
최대 근 이동	0.96	0.58	0.14	$8.8\times10^{-3}$	$1.4\times10^{-4}$

이다. $2^{-34} \approx 5.8\times10^{-11}$, 상대로는 $2.8\times10^{-13}$ 섭동인데도 근이 0.14 나 움직인다. 배정밀도의 반올림 단위가 $1.1\times10^{-16}$ 이니, 이 다항식을 계수로 들고 다니는 한 근 몇 개는 유효숫자 서너 자리도 못 건진다.

3. 근의 조건수

단순근 $r$ 에서 계수 $a_k$ ($x^k$ 의 계수)를 $\delta a_k$ 만큼 바꾸면, $w(r) = 0$ 을 음함수 미분해

$$\delta r \approx -\frac{r^k}{w'(r)}\,\delta a_k$$

를 얻는다. $w'(j) = \prod_{i \ne j}(j - i)$ 는 정확히 계산된다. $k = 19$ 에 대해 근마다 민감도를 재면

근 $j$	1	5	10	14	15	16	20
$\partial r_j/\partial a_{19}$	$8.2\times10^{-18}$	0.61	$-7.6\times10^{6}$	$-1.3\times10^{9}$	$2.1\times10^{9}$	$-2.4\times10^{9}$	$-4.3\times10^{7}$

로, 가장 민감한 근은 끝(20)이 아니라 가운데보다 조금 위인 15, 16 이다. 분자 $j^{19}$ 는 $j$ 와 함께 폭증하지만, 분모 $\lvert w'(j)\rvert = (j-1)!\,(20-j)!$ 는 가운데에서 가장 작기 때문이다. 두 효과가 겹쳐 15~16 에서 최대가 된다.

선형 이론이 맞는지 확인해 보자. $\delta = 2^{-44} \approx 5.7\times10^{-14}$ 에 $2.4\times10^{9}$ 를 곱하면 $1.37\times10^{-4}$ 이고, 위 표의 실측 최대 이동 $1.4\times10^{-4}$ 와 맞는다. 반면 $2^{-23}$ 에 대입하면 근 16 이 287 만큼 움직여야 하는데 실제로는 그렇지 않다 — 섭동이 이 정도로 크면 선형 근사의 수렴 반경을 한참 넘어, 이웃 근끼리 충돌해 복소 쌍으로 갈라지는 비선형 영역이다. 근 두 개가 실수축에서 만나 위아래로 갈라지는 순간이 정확히 이중근이고, 고유값 조건수 문서의 「조건수가 크다 = 결함 행렬이 바로 옆에 있다」와 같은 이야기다.

상대 조건수 $\lvert a_{19}\,r^{18}/w'(r)\rvert$ 로 바꾸면 근 16 에서 $3.2\times10^{10}$ 이다. 즉 $x^{19}$ 계수 하나의 상대오차가 근 16 의 상대오차로 300억 배 증폭된다. 다른 계수까지 동시에 흔들리면 더 나빠진다.

4. 배정밀도에서는 저장만 해도 깨진다

$w(x)$ 의 계수는 $x^2$ 의 계수에서 최대 $1.38\times10^{19}$ 까지 커져 $2^{53} \approx 9.0\times10^{15}$ 를 한참 넘는다. 그렇다고 큰 정수가 전부 깨지는 건 아니다 — $20!$ 처럼 2 의 거듭제곱을 많이 품은 정수는 정확히 표현된다. 실제로 세어 보면 21 개 정수 계수 중 5 개($x^3$ ~ $x^7$ 의 계수)가 배정밀도로 정확히 표현되지 않아 각각 16~512 만큼 반올림된다. 이 반올림된 계수로 만든 다항식의 근을 정확히(60 자리로) 구하면, 근이 최대 $6.2\times10^{-4}$ 움직인다. 아직 실수 근이긴 하지만 14.000548, 14.999627, 16.000192 처럼 소수 넷째 자리부터 틀린다. 알고리즘은 아직 한 줄도 돌지 않았다.

여기에 실제 근 찾기 알고리즘을 얹으면 더 나빠진다. NumPy 의 np.roots 는 동반 행렬을 만들어 QR 알고리즘으로 고유값을 구하는데, 같은 반올림 계수로 돌리면 근 15 를 15.0755 로, 13 을 13.0743 으로 내놓는다. 오차가 정확 근 대비 100 배 더 크다. 에델만-무라카미(1995)가 이유를 정리했다 — 동반 행렬 QR 은 행렬에 대해서는 후진 안정이지만, 그 후진 오차를 다항식 계수 쪽으로 옮기면 계수별 상대오차가 아니라 「가장 큰 계수 크기 × 기계 엡실론」 급 절대오차가 된다. 계수 크기가 1 에서 $10^{18}$ 까지 퍼진 다항식에서는 작은 계수가 상대적으로 크게 흔들린다는 뜻이다.

평가도 마찬가지다. $x = 15.5$ 에서 호너 방법으로 $w$ 를 배정밀도로 평가하면 $-5.5415\times10^{12}$, 정확한 값은 $-5.5788\times10^{12}$ 로 상대오차 0.7% 다. 값의 크기가 $10^{12}$ 인데 중간 합은 $10^{24}$ 급까지 커졌다가 상쇄되기 때문이다. 근 근처에서는 값 자체가 0 에 가까워져 상쇄가 더 심해지고, 뉴턴-랩슨법이 반복마다 계산하는 $w(x)/w'(x)$ 가 사실상 잡음이 된다.¹

5. 무엇이 바뀌었나 — 실무의 처방

윌킨슨 본인이 정리한 결론은 「근 찾기 알고리즘을 더 잘 만들어라」가 아니라 **「애초에 문제를 계수로 만들지 마라」**였다.

• 고유값 문제를 특성다항식으로 풀지 않는다. 1950년대까지는 행렬의 특성다항식 계수를 구한 다음 근을 찾는 방법이 쓰였다. 윌킨슨 다항식은 $\mathrm{diag}(1, \dots, 20)$ 의 특성다항식이다. 행렬로 들고 있으면 고유값은 대각 성분 그대로 정확한데, 계수로 바꾸는 순간 위와 같은 지옥이 된다. 현대 고유값 문제 해법(QR 알고리즘)은 특성다항식을 한 번도 만들지 않는다.
• 다항식이 주어지면 동반 행렬로, 가능하면 다른 기저로. 문제 자체가 다항식이면 동반 행렬 고유값이 표준이고(MATLAB roots, NumPy roots), 균형화(balancing)를 켜면 계수 스케일 문제가 많이 완화된다. 더 근본적으로는 다항식을 단항식 대신 체비쇼프 다항식 기저로 표현하고 「동료 행렬」(colleague matrix)의 고유값을 구한다. 구간 위에서 정의된 함수의 근을 찾을 때는 이쪽이 훨씬 잘 조건화되어 있다.
• 근을 알고 있다면 곱 형태 그대로. $\prod(x - r_i)$ 를 전개하지 말고 곱으로 평가하면 각 인수의 상대오차가 그대로 유지된다. 다항식 계수를 「정보를 압축한 표현」이라 생각하기 쉽지만, 근에 관해서는 오히려 정보를 잃는 표현이다.
• 의사근(pseudozero)으로 생각한다. 토-트레페슨(1994)처럼 「계수 상대오차 $\varepsilon$ 이내의 모든 다항식이 가질 수 있는 근의 집합」을 그리면, 근 15·16 근처에 커다란 덩어리가 보인다. 이게 의사스펙트럼의 다항식 판이다. 계산된 근 하나보다 그 덩어리가 문제의 진짜 답이다.

6. 왜 하필 1 ~ 20 인가

근을 $1, \dots, n$ 으로 잡으면 $n$ 이 커질수록 민감도가 폭증하고, $n = 20$ 은 당시 기계에서 「사람이 보기엔 순하고 컴퓨터에겐 치명적인」 경계였다. 같은 개수의 근을 $[-1, 1]$ 에 체비쇼프 점으로 놓으면 조건수가 훨씬 작고, 근을 단위원 위에 균등하게 놓은 $x^{20} - 1$ 은 아주 잘 조건화되어 있다.² 즉 문제는 「근이 많다」가 아니라 근이 실수축 위 한쪽에 몰려 있고 원점에서 멀다는 배치다. 단항식 $x^k$ 이 $[1, 20]$ 위에서 서로 거의 평행해져 표현이 퇴화하는 것이고, 힐베르트 행렬이 $[0,1]$ 위 단항식의 퇴화를 보여 주는 것과 같은 현상의 다른 얼굴이다.

윌킨슨은 훗날 이 발견을 회고하며 『The perfidious polynomial』(배신의 다항식, 1984)이라는 글을 남겼고 이 글로 미국수학협회(MAA)의 쇼브네 상을 받았다. 영국 국립물리연구소(NPL)의 초기 컴퓨터에서 다항식 근 찾기 프로그램을 시험하다 이 현상을 만났고, 그 글에서 이것을 「수치해석가로서 겪은 가장 충격적인 경험」이라고 회고했다.³

7. 관련 문서

• 고유값 조건수 · 바우어-파이크 정리 · 의사스펙트럼 · 조건수
• 후진 오차 해석 · 부동소수점 연산 · 힐베르트 행렬
• 뉴턴-랩슨법 · 뮬러 방법 · QR 알고리즘 · 체비쇼프 다항식
• 동반 행렬 · 호너 방법 · 젠킨스-트라우브 알고리즘

8. Footnotes

1. 근 근처에서 계산된 $w(x)$ 의 부호가 반올림 잡음으로 마구 바뀌는 구간이 생긴다. 이분법을 돌리면 그 구간 어딘가에서 「수렴」하는데, 그 근은 정확히 「잡음 띠의 폭」만큼만 믿을 수 있다. 수렴했다고 정답이라는 보장은 없다 — 이 바닥의 영원한 교훈. ↩

2. 단위원 위 균등 배치는 고속 푸리에 변환의 노드이기도 하다. 잘 조건화된 다항식 문제는 대개 원 위에 살고, 실수축 위에 줄 세운 근은 대개 지옥에 산다. ↩

3. 근이 1, 2, …, 20 이라는 걸 뻔히 아는 다항식에서 엉뚱한 근이 나오면 누구라도 코드부터 의심한다. 그런데 범인은 코드가 아니라 수학이었다. 오늘날 힐베르트 행렬로 솔버를 시험하다 같은 경험을 하는 대학원생과 70년을 사이에 두고 동지다. ↩