행켈 행렬

오류 제보
수치해석 마지막 수정: 2026-10-22 04:33:47

○ 공개 검수 대기 AI 보조 작성 표시 기준 보기

1. 개요

행켈 행렬
Hankel Matrix
정의Hij = hi+j — 반대각선마다 값이 일정
자유도n × n 에 2n − 1 개
토플리츠와의 관계행 순서를 뒤집으면(JH) 토플리츠
랭크 정리크로네커: 무한 행켈의 랭크 r ⇔ 생성함수가 r 차 유리함수
악명양정치 행켈은 조건수가 지수적으로 커진다 (힐베르트 행렬)
등장처모멘트 문제, 시스템 실현(ERA), 행렬 펜슬, 파데·섕크스, 모델 축소

토플리츠를 거울에 비춘 것. 그런데 거울 속에서는 신호가 몇 개의 지수함수로 이뤄졌는지가 랭크로 보인다.

행켈 행렬(Hankel matrix)은 (i,j)(i, j) 성분이 i+ji + j 에만 의존하는 행렬, 곧 각 반대각선(오른쪽 위에서 왼쪽 아래로 내려가는 선) 위의 값이 일정한 행렬이다. 수열 h0,h1,h2,…h_0, h_1, h_2, \dots 로

H=(h0h1h2⋯hn−1h1h2h3⋯hnh2h3h4⋯hn+1⋮⋱⋮hn−1hnhn+1⋯h2n−2)H = \begin{pmatrix} h_0 & h_1 & h_2 & \cdots & h_{n-1}\\ h_1 & h_2 & h_3 & \cdots & h_n\\ h_2 & h_3 & h_4 & \cdots & h_{n+1}\\ \vdots & & & \ddots & \vdots\\ h_{n-1} & h_n & h_{n+1} & \cdots & h_{2n-2} \end{pmatrix}

로 쓴다. 정사각이면 자동으로 대칭이다. 독일 수학자 헤르만 행켈(Hermann Hankel, 1839–1873)의 이름을 땄다1.

토플리츠 행렬이 i−ji - j 에 의존하는 「시간 차」 구조라면 행켈은 i+ji + j 에 의존하는 「시간 합」 구조다. 행 순서를 뒤집는 교환 행렬 JJ 를 곱하면 JHJH 가 토플리츠가 되므로 두 구조는 한 몸이다. 그래서 FFT 행렬-벡터 곱, 레빈슨형 빠른 풀이 같은 토플리츠 도구가 그대로 넘어온다. 그런데 응용에서 행켈이 나오는 자리는 토플리츠와 꽤 다르다. 토플리츠가 정상 과정의 공분산과 합성곱에서 나온다면, 행켈은 모멘트와 신호를 몇 개의 지수함수로 쪼개는 문제 — 곧 랭크 — 에서 나온다.

2. 모멘트와 힐베르트 행렬

측도 μ\mu 의 모멘트 mk=∫xk dμ(x)m_k = \int x^k\,d\mu(x) 로 만든 행켈 행렬 Hij=mi+jH_{ij} = m_{i+j} 는 단항식 1,x,x2,…1, x, x^2, \dots 의 그람 행렬이다.

cTH c=∫(∑icixi)2dμ(x)≥0\mathbf c^{\mathsf T} H\,\mathbf c = \int \Bigl(\sum_i c_i x^i\Bigr)^2 d\mu(x) \ge 0

이니까 양반정치이고, 지지집합이 무한하면 양정치다. 거꾸로 수열이 어떤 양의 측도의 모멘트가 되려면 모든 크기의 행켈 행렬이 양반정치여야 한다(함부르거 모멘트 문제). 행켈 행렬식은 직교다항식의 정규화 상수, 연분수 계수, 파데 근사의 존재 조건에도 그대로 등장한다. 섕크스 변환이 행켈 행렬식의 비로 정의되는 것도 같은 뿌리다.

[0,1][0, 1] 위 르베그 측도의 모멘트 mk=1/(k+1)m_k = 1/(k+1) 로 만든 행켈이 바로 힐베르트 행렬 Hij=1/(i+j+1)H_{ij} = 1/(i+j+1) 이다. 60자리로 2-노름 조건수를 재면:

nn4681012
κ2(Hn)\kappa_2(H_n)1.55×1041.55\times10^{4}1.50×1071.50\times10^{7}1.53×10101.53\times10^{10}1.60×10131.60\times10^{13}1.71×10161.71\times10^{16}
베커만 하한0.523.6272261.9×1031.9\times10^{3}

크기가 2 늘 때마다 약 1000배, 즉 한 칸에 약 33배씩 나빠진다. n=12n = 12 면 배정밀도로는 아무 자리도 믿을 수 없다.

2.1. 피할 수 없는 악조건

힐베르트 행렬이 유난히 나쁜 것일까? 베커만(2000)은 모든 실수 양정치 행켈 행렬 HnH_n (n≥3n \ge 3)에 대해

κ2(Hn)≥γ n−116n,γ=e4G/π≈3.210\kappa_2(H_n) \ge \frac{\gamma^{\,n-1}}{16n}, \qquad \gamma = e^{4G/\pi} \approx 3.210

(GG 는 카탈랑 상수)을 증명했다. 위 표의 둘째 줄이 이 하한이다. 힐베르트는 하한보다 훨씬 나쁘지만, 요점은 어떤 측도를 골라도 지수적 악화를 피할 수 없다는 것이다. 단항식 xkx^k 가 고차로 갈수록 서로 거의 평행해지는 탓이고, 실수 반데르몽드 행렬·크리로프 행렬이 악조건인 것과 같은 이유다.

실무적 교훈은 분명하다. 모멘트에서 직교다항식이나 가우스 구적 절점을 계산하지 마라. 행켈을 촐레스키 분해해 다항식 계수를 뽑는 교과서 방법은 배정밀도에서 십몇 차만 가도 무너진다. 대신 수정 모멘트(체비쇼프 다항식 같은 잘 조건화된 기저에 대한 모멘트), 스틸체스 절차, 또는 란초스 알고리즘으로 3항 점화식 계수를 직접 만든다. 이것이 가우치(Gautschi)가 평생 강조한 「모멘트 문제의 조건화」 이야기다2.

3. 크로네커 정리 — 랭크가 지수함수의 개수다

행켈 행렬의 핵심 정리는 크로네커(1881)의 것이다.

무한 행켈 행렬 (hi+j)i,j≥0(h_{i+j})_{i,j\ge0} 의 랭크가 유한한 rr 이다   ⟺  \iff 생성함수 ∑khkzk\sum_k h_k z^k 가 분모 차수 rr 인 유리함수다.

유리함수를 부분분수로 쪼개면, 랭크 rr 은 수열이

hk=∑iPi(k) zi k,∑i(deg⁡Pi+1)=rh_k = \sum_{i} P_i(k)\, z_i^{\,k}, \qquad \sum_i (\deg P_i + 1) = r

꼴, 곧 지수함수(감쇠 진동) rr 개의 합이라는 말과 같다. 이유는 간단하다 — 이런 수열은 rr 차 선형 점화식을 만족하고, 점화식은 행켈의 열들 사이의 선형 종속이다.

확인해 보자. hk=2(0.9)k+(−0.5)kh_k = 2(0.9)^k + (-0.5)^k 로 10×1010\times10 행켈을 만들면 특이값이 9.3659.365, 1.2151.215, 그다음 4×10−164\times10^{-16} — 정확히 랭크 2다. 반면 hk=1/(k+1)h_k = 1/(k+1) (힐베르트)은 유리 생성함수가 아니라 랭크가 꽉 차지만, 특이값이

1.75, 0.343, 3.6×10−2, 2.5×10−3, 1.3×10−4, 4.7×10−6, …1.75,\ 0.343,\ 3.6\times10^{-2},\ 2.5\times10^{-3},\ 1.3\times10^{-4},\ 4.7\times10^{-6},\ \dots

로 한 단계에 10~30배씩 떨어진다. 랭크는 꽉 찼는데 수치적으로는 저랭크인 셈이다. 이 「특이값 급감」이 다음 두 절, 시스템 식별과 모델 축소의 근거다.

4. 시스템 실현과 행켈 특이값

이산 선형 시스템 xk+1=Axk+Buk\mathbf x_{k+1} = A\mathbf x_k + B u_k, yk=Cxky_k = C\mathbf x_k 의 충격 응답(마르코프 파라미터) hk=CAkBh_k = CA^{k}B 로 행켈을 만들면

H=(CCACA2⋮)(BABA2B⋯)=O CH = \begin{pmatrix} C \\ CA \\ CA^2 \\ \vdots \end{pmatrix} \begin{pmatrix} B & AB & A^2B & \cdots \end{pmatrix} = \mathcal O\,\mathcal C

로 관측 가능성 행렬 × 제어 가능성 행렬로 쪼개진다. 따라서 rank⁡H\operatorname{rank} H 는 시스템의 최소 차수이고, HH 를 특이값 분해해 O\mathcal O, C\mathcal C 를 복원하면 (A,B,C)(A, B, C) 가 나온다. 이것이 호-칼만(1966)의 실현 이론이다.

구조 동역학에서 이것을 잡음 섞인 실측에 쓰도록 만든 것이 ERA(eigensystem realization algorithm, 준-파파 1985)다. 가속도계로 잰 충격 응답으로 행켈을 만들고, 특이값의 급감 지점에서 잘라 차수를 정한 뒤, 한 칸 밀린 행켈 H(1)H^{(1)} (원소 hi+j+1h_{i+j+1})과의 관계로 AA 를 뽑는다. AA 의 고유값이 고유진동수와 감쇠비를 준다 — 실험 모드 해석의 표준 도구 중 하나이고, NASA 랭글리에서 우주 구조물 시험용으로 나왔다. 넓게 보면 시스템 식별의 부분공간 방법 전체가 이 행켈 분해의 변주다.

시스템의 행켈 특이값 σ1≥σ2≥⋯\sigma_1 \ge \sigma_2 \ge \cdots (행켈 연산자의 특이값, 연속시간이면 두 그람 행렬 곱 PQPQ 고유값의 제곱근)은 각 상태가 입력→출력 에너지 전달에 얼마나 기여하는지를 잰다. 작은 행켈 특이값에 대응하는 상태를 버리는 것이 균형 절단이고, 오차가 2∑i>rσi2\sum_{i>r}\sigma_i 로 묶인다. 글로버(1984)는 행켈 노름에서 최적인 차수 rr 근사를 모두 구성했고 그 오차가 정확히 σr+1\sigma_{r+1} 임을 보였다(AAK 정리의 행렬판). 이것이 축소차수모델의 시스템 이론 쪽 기둥이다.

5. 신호에서 지수함수 뽑기 — 행렬 펜슬

크로네커 정리를 데이터에 거꾸로 쓰면 「신호를 감쇠 진동 몇 개로 분해」하는 알고리즘이 된다. 고전판이 1795년의 프로니 방법이고, 잡음에 강한 현대판이 행렬 펜슬 방법(후아-사르카르, 1990)이다.

  1. 표본 x0,…,xN−1x_0, \dots, x_{N-1} 로 (N−L)×(L+1)(N-L)\times(L+1) 행켈을 만든다.
  2. SVD 로 신호 부분공간(큰 특이값 rr 개)의 오른쪽 특이벡터 VV 를 얻는다.
  3. VV 의 마지막 행을 뺀 것 V1V_1, 첫 행을 뺀 것 V2V_2 사이의 관계 V2≈V1ZV_2 \approx V_1 Z 를 최소제곱으로 풀고, ZZ 의 고유값을 극 ziz_i 로 읽는다.

3단계가 동작하는 이유는 행켈 구조의 이동 불변성이다 — 한 칸 밀린 창은 각 지수 성분에 ziz_i 를 곱한 것과 같다.

실측: xk=2(0.98)kcos⁡(0.6k)+2(0.95)kcos⁡(1.1k)+x_k = 2(0.98)^k\cos(0.6k) + 2(0.95)^k\cos(1.1k) + 표준편차 0.01 잡음, N=64N = 64, L=24L = 24. 행켈 특이값은

18.38, 17.59, 10.06, 9.08 ∣ 0.098, 0.090, 0.082, …18.38,\ 17.59,\ 10.06,\ 9.08\ \big\vert\ 0.098,\ 0.090,\ 0.082,\ \dots

로 4번째와 5번째 사이에 100배 낭떠러지가 있어 r=4r = 4 (켤레쌍 2개)가 자명하다. 펜슬로 뽑은 극은 크기 0.98000.9800, 0.95000.9500, 각도 ±0.5999\pm0.5999, ±1.0999\pm1.0999 로 참값 0.980.98, 0.950.95, 0.60.6, 1.11.1 을 넷째 자리까지 맞힌다. 64개 표본으로 주파수를 10−410^{-4} rad 정밀도로 잡는 것인데, 같은 길이의 FFT 빈 간격은 2π/64≈0.0982\pi/64 \approx 0.098 rad 다.

같은 아이디어가 이름을 바꿔 여러 곳에 있다.

  • 특이 스펙트럼 분석(SSA): 시계열을 시간 지연 행켈에 넣고 SVD 로 추세·주기·잡음을 가른다. 기후·경제 시계열에서 쓴다.
  • 행켈 DMD: 동적 모드 분해를 지연 좌표로 쌓은 행켈 위에서 하면 관측량 하나로도 쿠프만 연산자의 고유값을 근사할 수 있다.
  • 캐드조우 잡음 제거: 「저랭크로 자르기」와 「반대각선 평균으로 행켈 구조 되돌리기」를 번갈아 하는 교대 투영. 구조화된 저랭크 근사의 원조격이다.

6. 계산 요령

  • 곱셈은 FFT 로. HvH\mathbf v 는 수열 hh 와 뒤집힌 v\mathbf v 의 합성곱 일부이므로 고속 푸리에 변환으로 O(nlog⁡n)O(n\log n). 행렬을 저장할 필요도 없이 2n−12n-1 개 수만 있으면 된다. 큰 행켈 SVD 를 랜덤화 SVD나 란초스로 할 때 이 곱이 핵심 부품이다.
  • 풀이는 토플리츠로 바꿔서. Hx=bH\mathbf x = \mathbf b 는 (JH)x=Jb(JH)\mathbf x = J\mathbf b 로 토플리츠 계가 되어 레빈슨형 O(n2)O(n^2) 알고리즘을 쓸 수 있다. 다만 JHJH 는 일반적으로 대칭 토플리츠가 아니므로 대칭판 레빈슨-더빈 알고리즘이 아니라 일반 레빈슨(비대칭판)이 필요하다.
  • 행렬식·랭크는 믿지 마라. 위 악조건 때문에 행켈 행렬식을 직접 계산하는 알고리즘(교과서의 섕크스 정의, 파데 계수의 크래머 해)은 고차에서 망가진다. 실무에서는 ε-알고리즘·qd-알고리즘처럼 행렬식을 피하는 점화식이나, SVD 로 수치 랭크를 재는 방식을 쓴다3.

7. 관련 문서

8. Footnotes

  1. 같은 사람이 베셀 함수의 행켈 함수, 행켈 변환의 그 행켈이다. 서른넷에 요절했는데 이름이 붙은 물건이 셋 이상이라, 수치해석 시험 범위에서 마주칠 확률이 은근히 높다. ↩

  2. 가우치의 결론을 한 줄로 줄이면 “모멘트는 정보가 아니라 함정이다”. 같은 측도라도 어떤 기저의 모멘트를 들고 오느냐에 따라 조건수가 「천문학적」과 「한두 자리」 사이를 오간다. 격자가 아니라 기저가 다 했다. ↩

  3. 수치 랭크를 「0이 아닌 특이값의 개수」로 세면 힐베르트 행렬은 영원히 꽉 찬 랭크다. 「허용오차 위 특이값의 개수」로 세야 하고, 그 허용오차를 고르는 일이 사실상 모델 차수를 고르는 일이다. 낭떠러지가 선명하면 국룰, 완만하면 현생. ↩

검증 자료

이 문서의 주요 출처

  1. B. Beckermann (2000), The condition number of real Vandermonde, Krylov and positive definite Hankel matrices, Numer. Math. 85, 553-577
  2. J.-N. Juang & R. S. Pappa (1985), An eigensystem realization algorithm for modal parameter identification and model reduction, J. Guidance, Control, and Dynamics 8(5), 620-627
  3. K. Glover (1984), All optimal Hankel-norm approximations of linear multivariable systems and their L-infinity error bounds, Int. J. Control 39(6), 1115-1193
  4. Y. Hua & T. K. Sarkar (1990), Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise, IEEE Trans. ASSP 38(5), 814-824
  5. V. V. Peller, Hankel Operators and Their Applications, Springer 2003
이 문서가 문제를 푸는 데 도움이 됐나요? 한 번의 클릭이 다음 검수 순서를 정합니다.

의견을 반영할게요. 감사합니다.

다음으로 읽기 동적 모드 분해 약 5분 더 탐색 →