힐베르트 행렬

1. 개요

힐베르트 행렬 Hilbert Matrix
정의	$H_{ij} = 1/(i+j-1)$, $i, j = 1, \dots, n$
출처	힐베르트 (1894)
정체	단항식 $1, t, \dots, t^{n-1}$ 의 $[0,1]$ 위 그람 행렬
구조	대칭 양정치 · 행켈 · 코시 · 완전 양행렬
역행렬	정수 성분, 닫힌꼴, 성분 합 $= n^2$
노름	$\lVert H_n\rVert_2 \nearrow \pi$ (무한 행렬의 노름이 정확히 $\pi$)
쓰임	악조건 시험 행렬, 반올림 오차 교보재

성분이 전부 단위분수인 순한 얼굴의 행렬. 12×12 만 돼도 배정밀도를 통째로 집어삼킨다.

힐베르트 행렬(Hilbert matrix)은

$$H_n = \left[\frac{1}{i+j-1}\right]_{i,j=1}^{n} = \begin{bmatrix} 1 & \tfrac12 & \tfrac13 & \cdots \\ \tfrac12 & \tfrac13 & \tfrac14 & \cdots \\ \tfrac13 & \tfrac14 & \tfrac15 & \cdots \\ \vdots & & & \ddots \end{bmatrix}$$

로 정의되는 $n\times n$ 대칭 행렬이다. 다비트 힐베르트가 1894년 르장드르 다항식에 관한 짧은 논문에서 이 행렬식을 계산하면서 등장했다. 반대각선을 따라 성분이 같으니 행켈 행렬이고, $1/(x_i - y_j)$ 꼴이니 코시 행렬이기도 하다.

수치해석에서 이 행렬이 유명한 이유는 하나다. 조건수가 $n$ 에 대해 지수적으로 폭발한다. 그런데 역행렬·행렬식이 정수와 유리수의 닫힌꼴로 나와서 「정답」을 부동소수점 없이 정확히 알 수 있다. 악조건이 계산을 어떻게 망치는지를 반올림 오차 없는 기준값과 대조해 볼 수 있는 거의 유일한 교보재다. 조건수 실측표는 조건수 문서에, 코시 행렬로서의 성질과 변위 구조는 코시 행렬 문서에 있으니, 이 문서는 역행렬 닫힌꼴, 촐레스키 인수가 감추고 있는 르장드르 다항식, 무한 행렬의 노름 $\pi$, 그리고 「솔버가 아니라 입력이 범인」인 반올림 실험을 다룬다.

2. 어디서 나오나 — 단항식의 그람 행렬

$[0,1]$ 위에서 함수 $f$ 를 $n-1$ 차 다항식 $p(t) = \sum_{j=1}^{n} c_j t^{j-1}$ 로 최소제곱 근사한다고 하자. $\int_0^1 (f - p)^2\,dt$ 를 $c_j$ 로 미분해 0 으로 두면 정규방정식

$$\sum_{j=1}^{n} \underbrace{\int_0^1 t^{i-1} t^{j-1}\,dt}_{=\,1/(i+j-1)}\, c_j = \int_0^1 t^{i-1} f(t)\,dt$$

이 나온다. 계수 행렬이 정확히 $H_n$ 이다. 즉 힐베르트 행렬은 단항식 기저 $\{1, t, t^2, \dots\}$ 의 $L^2[0,1]$ 그람 행렬이다. $t^{10}$ 과 $t^{11}$ 은 $[0,1]$ 위에서 거의 같은 모양(0 근처에서 납작하다가 1 근처에서 솟는다)이라 두 기저 벡터가 거의 평행하고, 그 결과 그람 행렬이 거의 특이해진다. 악조건의 원인은 행렬이 아니라 기저다. 이 사실이 뒤의 촐레스키 절에서 정확한 식으로 드러난다.

양정치성은 공짜로 나온다. 임의의 $\mathbf c \ne 0$ 에 대해 $\mathbf c^T H_n \mathbf c = \int_0^1 p(t)^2\,dt > 0$ 이기 때문이다.

3. 역행렬 — 정수의 닫힌꼴

힐베르트 행렬의 역행렬은 모든 성분이 정수다. 닫힌꼴은

$$(H_n^{-1})_{ij} = (-1)^{i+j}(i+j-1)\binom{n+i-1}{n-j}\binom{n+j-1}{n-i}\binom{i+j-2}{i-1}^2$$

이다. $n = 4$ 를 대입하면

$$H_4^{-1} = \begin{bmatrix} 16 & -120 & 240 & -140 \\ -120 & 1200 & -2700 & 1680 \\ 240 & -2700 & 6480 & -4200 \\ -140 & 1680 & -4200 & 2800 \end{bmatrix}$$

이고, 유리수 산술로 $H_4 H_4^{-1} = I$ 를 검산하면 정확히 맞는다($n = 10$ 까지 같은 방식으로 확인했다). 부호가 체스판처럼 번갈아 나오는 것은 $H_n$ 이 완전 양행렬(모든 소행렬식이 양수)이라서 생기는 일이다.

성분 합이 $n^2$ 이다. $H_4^{-1}$ 의 16 개 성분을 다 더하면 정확히 16, $H_{10}^{-1}$ 은 100 이다. 각 성분은 $10^{12}$ 급인데 부호가 번갈아 상쇄되어 100 만 남는다. 크기가 그렇게 큰 정수들의 합이 $n^2$ 으로 떨어진다는 것 자체가 「역행렬에 숨은 엄청난 상쇄」의 증거이고, 이 상쇄가 바로 부동소수점이 버티지 못하는 지점이다. 최대 성분 크기는 $n = 4$ 에서 6480, $n = 10$ 에서 $3.48\times10^{12}$, $n = 12$ 에서 $3.66\times10^{15}$ 다.

행렬식. $c_n = \prod_{k=1}^{n-1} k!$ 로 두면

$$\det H_n = \frac{c_n^4}{c_{2n}}$$

이다. $\det H_3 = 1/2160$, $\det H_5 = 1/266\,716\,800\,000$, $\det H_{10} \approx 2.16\times10^{-53}$. 행렬식이 작다는 것 자체는 악조건의 증거가 아니지만(대각 행렬 $0.1I$ 를 생각해 보라) 여기서는 고유값이 실제로 그만큼 작다. 가장 작은 고유값이 $n = 10$ 에서 $1.09\times10^{-13}$, $n = 20$ 에서 $7.8\times10^{-29}$ 다.

4. 촐레스키 인수 속의 르장드르 다항식

$H_n = LDL^T$ ($L$ 은 대각이 1 인 하삼각, $D$ 는 대각)로 분해하면 피벗이 깔끔한 닫힌꼴로 나온다. 유리수 산술로 $n = 10$ 까지 계산해 보면

$$D_k = \frac{1}{(2k-1)\binom{2k-2}{k-1}^2},\qquad D = \mathrm{diag}\!\left(1,\ \tfrac{1}{12},\ \tfrac{1}{180},\ \tfrac{1}{2800},\ \dots\right)$$

이 모든 $k$ 에서 정확히 성립한다. $n = 4$ 의 단위 하삼각 인수는

$$L = \begin{bmatrix} 1 & & & \\ \tfrac12 & 1 & & \\ \tfrac13 & 1 & 1 & \\ \tfrac14 & \tfrac{9}{10} & \tfrac32 & 1 \end{bmatrix}$$

이다. 이게 무슨 뜻인가. $LDL^T$ 분해는 그람 행렬 위에서 그람-슈미트 직교화를 하는 것과 같다. 단항식 $1, t, t^2, \dots$ 를 $[0,1]$ 위에서 차례로 직교화하면 이동 르장드르 다항식(shifted Legendre polynomial)의 최고차항 계수 1 짜리 판이 나오고, $D_k$ 는 그 $k-1$ 차 다항식의 노름 제곱 $\int_0^1 \tilde P_{k-1}(t)^2\,dt$ 다. 실제로 최고차항 계수가 1 인 $m$ 차 이동 르장드르 다항식의 노름 제곱은 $1/\bigl((2m+1)\binom{2m}{m}^2\bigr)$ 이고, $m = k-1$ 을 넣으면 위 식과 똑같다.

그러니 처방도 여기서 바로 나온다. 단항식 대신 르장드르(또는 체비쇼프) 기저로 최소제곱을 풀면 그람 행렬이 대각 행렬이 되어 조건수 문제가 사라진다. 이것이 보간과 근사와 직교다항식 문서가 「고차 다항식 맞춤에는 직교 기저」를 국룰로 박아 두는 이유다. 피벗 $D_k$ 가 대략 $16^{-k}$ 로 줄어드는 것이 곧 단항식 기저가 얼마나 빨리 퇴화하는지의 정량이다.¹

5. 힐베르트 부등식과 노름 $\pi$

유한 $H_n$ 의 가장 큰 고유값(= 2-노름)은 $n$ 이 커져도 발산하지 않는다. 직접 재 보면

$n$	5	10	20	100	1000	4000
$\lVert H_n\rVert_2$	1.567	1.752	1.907	2.183	2.443	2.553

로 아주 느리게 증가하며 $\pi \approx 3.1416$ 아래에 머문다. 이것이 힐베르트 부등식이다. 제곱합이 유한한 수열 $a_1, a_2, \dots$ 에 대해

$$\sum_{i=1}^{\infty}\sum_{j=1}^{\infty} \frac{a_i a_j}{i+j-1} \le \pi \sum_{i=1}^{\infty} a_i^2$$

이 성립하고, 상수 $\pi$ 는 더 줄일 수 없다(슈어 1911). 연산자 언어로 말하면 무한 힐베르트 행렬은 $\ell^2$ 위의 유계 연산자이고 노름이 정확히 $\pi$ 다. 마그누스(1950)는 더 나아가 그 스펙트럼이 고유값 없이 연속적으로 $[0, \pi]$ 를 채운다는 것을 보였다.

표에서 보듯 수렴은 로그 속도라 끔찍하게 느리다. $n = 4000$ 에서도 $\pi$ 까지 0.59 가 남는다. $\pi$ 와의 차이가 $\pi^5/(2\ln n)^2$ 꼴로 줄어든다는 점근 결과가 알려져 있는데, 이 식을 위 표의 $n$ 에 대입하면 음수나 엉뚱한 값이 나온다 — 로그 점근식이 쓸 만해지려면 $\ln n$ 이 커야 하고, 그 $n$ 은 메모리에 올릴 수 있는 크기가 아니다.²

수치해석 쪽에서 이 사실이 주는 교훈은 명쾌하다. $\lVert H_n\rVert_2$ 는 3 을 넘지 않으므로 조건수 $\kappa_2 = \lVert H_n\rVert\lVert H_n^{-1}\rVert$ 의 폭발은 전부 역행렬 쪽, 즉 가장 작은 고유값 쪽에서 온다. 큰 쪽은 얌전하고 작은 쪽이 바닥으로 떨어진다.

6. 실험 — 범인은 솔버가 아니라 입력이다

$H_n$ 을 배정밀도 배열에 담는 순간 이미 $1/3$, $1/6$, $1/7$ 같은 성분이 반올림된다. 그러면 컴퓨터가 받은 행렬은 $H_n$ 이 아니라 상대오차 $10^{-16}$ 급 섭동이 들어간 $\tilde H_n$ 이다. 이 $\tilde H_n$ 의 역행렬을 60 자리 정밀도로 정확히 구한 것과, NumPy 의 inv 로 배정밀도에서 구한 것을 각각 참 역행렬 $H_n^{-1}$ 과 비교하면(1-노름 상대오차):

$n$	$\kappa_1(H_n)$	반올림된 입력의 정확한 역행렬	배정밀도 inv
6	$7.9\times10^{7}$	$7.8\times10^{-11}$	$1.3\times10^{-10}$
8	$9.9\times10^{10}$	$2.7\times10^{-9}$	$3.4\times10^{-9}$
10	$1.2\times10^{14}$	$9.0\times10^{-5}$	$3.4\times10^{-5}$
12	$1.5\times10^{17}$	$1.8\times10^{-2}$	$8.5\times10^{-3}$

둘째 열과 셋째 열이 같은 자릿수다. $n = 12$ 에서 역행렬이 1~2% 틀리는 것은 LU 분해 탓이 아니라, 성분을 배정밀도에 저장한 것만으로 이미 생긴 일이다. 아무리 완벽한 솔버를 가져와도 이 이상 좋아질 수 없다. 이것이 후진 오차 해석이 말하는 「후진 안정한 알고리즘의 오차는 조건수 × 입력 오차를 넘지 않고, 그 이상 줄일 방법도 없다」는 명제의 가장 깔끔한 실연이다. (정확한 비교를 원하면 MATLAB 의 invhilb, SciPy 의 scipy.linalg.invhilbert 처럼 닫힌꼴로 정수 역행렬을 바로 주는 함수를 쓴다. 다만 $n$ 이 커지면 그 정수가 $2^{53}$ 을 넘어 배정밀도에서 정확히 표현되지 않는다.)

그래서 힐베르트 행렬로 자작 솔버를 시험할 때는 두 가지를 구분해야 한다. 「참 $H_n$ 에 대한 오차」는 입력 반올림이 지배하니 솔버 품질을 말해 주지 않는다. 솔버를 평가하려면 「컴퓨터가 실제로 받은 $\tilde H_n$ 에 대한 잔차·후진 오차」를 봐야 한다. 많은 교재 연습문제가 이 둘을 섞어서 「가우스 소거가 망가졌다」는 잘못된 결론을 유도한다.³

7. 그 밖의 얼굴들

• 완전 양행렬. 모든 소행렬식이 양수다. 코시 행렬식 공식에서 $x_i = i$, $y_j = 1 - j$ 로 두면 차의 곱이 전부 양수가 되기 때문이다. 이 덕분에 고유값이 모두 서로 다르고, $k$ 번째로 큰 고유값의 고유벡터가 정확히 $k-1$ 번 부호를 바꾼다(간트마허-크레인의 진동 행렬 성질).
• 높은 상대 정확도. 성분이 차의 곱으로 정확히 표현되는 구조를 이용하면 가장 작은 특이값까지 상대 정확도로 계산할 수 있다(데멜의 구조 SVD, 코시 행렬 참고). 일반 SVD 는 $\lambda_{\min} \approx 10^{-28}$ 을 $10^{-16}\lVert H\rVert$ 급 절대오차로만 주니 자릿수가 하나도 안 맞는다.
• 시험 행렬로서의 한계. 대칭 양정치이고 고유벡터 구조가 규칙적이라, 비대칭 문제나 피벗팅 전략을 시험하기엔 부족하다. 가우스 소거법의 성장 인자를 키우는 행렬, 비정규 행렬(의사스펙트럼) 같은 다른 시험 행렬과 섞어 써야 한다.
• 연속판. 무한 행렬의 연속 유사물은 $[0,\infty)$ 위의 적분 연산자 $(Kf)(x) = \int_0^\infty f(y)/(x+y)\,dy$ 이고 이것도 노름이 $\pi$ 다. 힐베르트 부등식의 적분형이 바로 이 연산자의 유계성이다.

8. 관련 문서

• 조건수 · 코시 행렬 · 행켈 행렬 · 반데르몽드 행렬
• 르장드르 다항식 · 직교다항식 · 그람-슈미트 · 보간과 근사
• 최소자승법 · 촐레스키 분해 · 후진 오차 해석 · 부동소수점 연산
• 힐베르트 공간 · 그람 행렬 · 완전 양행렬

9. Footnotes

1. 같은 그람-슈미트를 부동소수점에서 단항식 벡터에 직접 돌리면 직교성이 순식간에 무너진다. 그러니 「르장드르 기저를 쓰라」는 말은 「단항식을 직교화하라」가 아니라 「처음부터 3항 점화식으로 르장드르를 만들어라」는 뜻이다. 직교화를 수치적으로 하는 순간 힐베르트 행렬의 저주를 고스란히 물려받는다. ↩

2. 점근식은 맞는데 숫자를 넣으면 틀리는 사태. 로그가 들어간 점근식의 국룰이다. 「$n \to \infty$」의 $\infty$ 가 생각보다 훨씬 멀다. ↩

3. 「n=12 에서 해가 엉망 → 솔버 버그다 → 밤새 디버깅」 루트를 탄 대학원생이 한둘이 아니다. 솔버는 무죄, 범인은 1/3 을 저장한 순간이었다. ↩