레빈슨-더빈 알고리즘

1. 개요

레빈슨-더빈 알고리즘 Levinson–Durbin Recursion
푸는 것	대칭 양정치 토플리츠 선형계
비용	더빈 약 2n2, 레빈슨 약 4n2 flops
메모리	O(n)
부산물	반사 계수(부분 자기상관), 차수별 예측 오차
안정성	양정치에서 약한 안정 (사이벤코 1980)
제안	레빈슨 1947, 더빈 1960

크기 $m$ 짜리 답을 알면 크기 $m+1$ 짜리 답은 내적 하나로 나온다. 그걸 $n$ 번 하면 끝.

레빈슨-더빈 알고리즘(Levinson–Durbin recursion)은 대칭 양정치 토플리츠 행렬을 계수로 하는 선형계를, 크기 1부터 한 차수씩 키워 가는 점화식으로 $O(n^2)$ 연산에 푸는 방법이다. 일반 촐레스키 분해의 $n^3/3$ 에 비해 한 차수가 떨어지고, 저장도 $O(n)$ 이면 된다.

역사는 둘로 갈린다. 노먼 레빈슨이 1947년 위너 필터의 이산판 — 일반 우변을 가진 토플리츠 계 — 을 풀려고 이 점화식을 냈고¹, 제임스 더빈이 1960년 시계열 모형 적합에서 우변이 행렬의 첫 열 자체인 특수한 경우(율-워커 방정식)를 더 짧게 정리했다. 오늘날 「레빈슨-더빈」은 보통 더빈판을, 「레빈슨 알고리즘」은 일반 우변판을 가리킨다.

이 알고리즘이 단순한 선형계 풀이 요령에 그치지 않는 이유는 중간 산물이 전부 의미를 갖기 때문이다. 각 단계의 계수는 차수 $m$ 의 최적 선형 예측기이고, 새로 생기는 한 개의 수(반사 계수)는 통계학의 부분 자기상관이자 음향관 모형의 반사율이며, 그 크기가 1 미만인지가 행렬의 양정치성과 필터의 안정성을 동시에 판정한다.

2. 율-워커 방정식과 더빈 점화식

정상 과정 $x_t$ 를 과거 $p$ 개 값으로 예측하는 $\hat x_t = \sum_{i=1}^{p} a_i x_{t-i}$ 의 최적 계수는 자기공분산 $r_k = \mathrm{E}[x_t x_{t-k}]$ 로 쓴 율-워커 방정식

$$\begin{pmatrix} r_0 & r_1 & \cdots & r_{p-1} \\ r_1 & r_0 & \cdots & r_{p-2} \\ \vdots & & \ddots & \vdots \\ r_{p-1} & r_{p-2} & \cdots & r_0 \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_p \end{pmatrix} = \begin{pmatrix} r_1 \\ r_2 \\ \vdots \\ r_p \end{pmatrix}$$

을 만족한다. 행렬은 대칭 토플리츠이고 우변은 그 첫 열을 한 칸 민 것이다.

더빈 점화식은 $E_0 = r_0$ 에서 시작해 $m = 1, 2, \dots, p$ 마다

$$k_m = \frac{r_m - \sum_{i=1}^{m-1} a_i^{(m-1)}\, r_{m-i}}{E_{m-1}},$$ $$a_i^{(m)} = a_i^{(m-1)} - k_m\, a_{m-i}^{(m-1)} \quad (i < m), \qquad a_m^{(m)} = k_m, \qquad E_m = (1 - k_m^2)\, E_{m-1}$$

를 반복한다. 단계 $m$ 의 비용은 내적 하나와 벡터 갱신 하나, 약 $4m$ flops 이므로 전체가 약 $2p^2$ 이다.

2.1. 왜 되는가 — 뒤집어도 같은 행렬

열쇠는 페르시메트리(persymmetry)다. 순서를 뒤집는 교환 행렬 $J$ 에 대해 대칭 토플리츠는 $JTJ = T$ 를 만족한다. 그래서 $T_m \mathbf{a} = \mathbf{r}$ 의 해를 뒤집은 $J\mathbf{a}$ 는 「우변을 뒤집은 계」의 해가 된다 — 시간을 거꾸로 돌린 후방 예측기다.

크기 $m+1$ 계를 $\begin{pmatrix} T_m & J\mathbf{r} \\ (J\mathbf{r})^\top & r_0\end{pmatrix}$ 꼴로 나눠 보면, 새 해는 이전 전방 해에 후방 해를 $k_m$ 배 섞고 마지막 칸에 $k_m$ 을 붙인 것이 된다. 위의 $a_i - k_m a_{m-i}$ 가 정확히 그 「뒤집어서 섞기」다. 새로 정할 미지수가 하나뿐이니 $k_m$ 은 스칼라 방정식 하나로 나오고, 그 분모 $E_{m-1}$ 은 차수 $m-1$ 예측기의 평균제곱 오차다.

3. 중간 산물의 의미

3.1. 반사 계수 = 부분 자기상관

$k_m$ 은 $x_t$ 와 $x_{t-m}$ 사이에서 중간의 $x_{t-1}, \dots, x_{t-m+1}$ 효과를 걷어낸 상관, 즉 부분 자기상관(PACF)이다. 자기회귀 모형 AR($p$) 에서는 $m > p$ 부터 $k_m = 0$ 이 되므로 PACF 그림이 끊기는 지점이 차수다.

$x_t = 1.5\,x_{t-1} - 0.75\,x_{t-2} + \varepsilon_t$ ($\mathrm{Var}\,\varepsilon = 1$) 의 정확한 자기공분산 $r_0, r_1, r_2, \dots = 112/13,\ 96/13,\ 60/13,\ 18/13, \dots$ 을 넣으면

$m$	1	2	3	4	5
$k_m$	$6/7 = 0.857143$	$-0.75$	$0$	$0$	$0$
$E_m$	$16/7 = 2.2857$	$1$	$1$	$1$	$1$

로, 2단계에서 계수가 정확히 $(1.5, -0.75)$ 가 되고 예측 오차가 잡음 분산 1에 닿은 뒤 더 줄지 않는다. 같은 과정을 19000개 시뮬레이션해 표본 자기공분산으로 돌리면 $k = (0.8571, -0.7524, -0.0038, 0.0171, 0.0053)$ — 3차부터 표본 잡음 수준($\pm 2/\sqrt{N} \approx \pm0.015$ 근처)으로 떨어진다.

3.2. $\lvert k_m \rvert < 1$ 의 삼중 의미

$E_m = (1-k_m^2)E_{m-1}$ 이고 $E_m$ 은 예측 오차 분산이니 음수가 될 수 없다. 실제로 다음이 동치다.

• 모든 $m$ 에서 $\lvert k_m \rvert < 1$
• $T_p$ 가 양정치 (모든 선행 주소행렬식이 양수, $\det T_m = \prod_{j<m} E_j$)
• 예측 오차 필터 $A(z) = 1 - \sum a_i z^{-i}$ 의 영점이 모두 단위원 안 — 최소 위상, 역필터 $1/A(z)$ 가 안정

그래서 반사 계수는 안정성이 보장된 필터 매개변수화다. 계수 $a_i$ 를 양자화하면 극점이 단위원 밖으로 튀어 나갈 수 있지만, $k_m$ 을 $(-1, 1)$ 안에서 양자화하면 어떻게 깎아도 필터는 안정하다.

4. 레빈슨 — 일반 우변

$T\mathbf{x} = \mathbf{b}$ 처럼 우변이 임의이면, 더빈 점화식으로 후방 예측기 $J\mathbf{a}^{(m)}$ 를 같이 키우면서 해를

$$\mathbf{x}^{(m+1)} = \begin{pmatrix} \mathbf{x}^{(m)} \\ 0 \end{pmatrix} + \mu_m \begin{pmatrix} -J\mathbf{a}^{(m)} \\ 1 \end{pmatrix}, \qquad \mu_m = \frac{b_{m+1} - (J\mathbf{r}^{(m)})^\top \mathbf{x}^{(m)}}{E_m}$$

로 확장한다. 내적이 하나 더 들어 비용은 약 $4n^2$ flops 다(골럽·반 론). 위너 필터의 정규방정식이 정확히 이 꼴이다.

4.1. 실측 — 비용

자기공분산 $r_k = (1+k)^{-1.5}$ 의 토플리츠 계를 풀었다.

$n$	더빈 (순수 파이썬 루프)	SciPy solve_toeplitz	조밀 LU
500	4.7 ms	1.0 ms	33 ms
1000	10 ms	3.2 ms	62 ms
2000	24 ms	12 ms	242 ms

인터프리터 루프로 짠 더빈조차 최적화된 LAPACK 조밀 풀이를 10배 앞선다. $n$ 이 두 배일 때 조밀 쪽은 $n^3$ 대로 늘고 레빈슨 쪽은 $n^2$ 대로 는다(작은 $n$ 에서는 오버헤드 때문에 비율이 깔끔하지 않다).

5. 수치 안정성

레빈슨-더빈은 피벗 없는 점화식이라 안정성이 오래 논쟁거리였다. 사이벤코(1980)는 양정치 토플리츠에서 더빈 알고리즘의 오차가 촐레스키 분해와 같은 크기 수준임을 보였다 — 이른바 약한 안정성이다. 다만 후진 오차 해석 의미의 후진 안정성까지 보장되는 것은 아니어서, 잔차가 촐레스키보다 조건수에 더 민감하게 커질 수 있다.

가우스형 자기공분산 $r_k = e^{-(k/L)^2}$ ($n = 60$) 로 조건수를 키워 가며 상대잔차를 촐레스키와 비교했다.

$L$	조건수	더빈 잔차	촐레스키 잔차	$\max\lvert k_m \rvert$
1.0	$5.9$	$3.6\times10^{-16}$	$1.2\times10^{-16}$	0.368
1.5	$1.3\times10^{2}$	$2.0\times10^{-15}$	$2.5\times10^{-16}$	0.641
2.0	$9.1\times10^{3}$	$4.9\times10^{-15}$	$1.6\times10^{-15}$	0.779
2.5	$2.1\times10^{6}$	$1.5\times10^{-14}$	$6.0\times10^{-15}$	0.852

더빈의 잔차가 촐레스키보다 몇 배 크지만 조건수가 $10^6$ 이 되어도 $10^{-14}$ 수준에 머문다. 사이벤코의 결론 그대로다. 문제가 커지는 것은 $\lvert k_m \rvert$ 가 1에 바짝 붙을 때 — $E_m$ 이 상쇄로 작아져 다음 단계의 분모가 부정확해진다. 스펙트럼이 거의 선 스펙트럼인 신호(순수 사인파 몇 개)가 그렇다².

5.1. 양정치가 아니면

$r = (0, 1, 0)$ 이면 $T = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$ 는 정칙인데 첫 단계 분모가 $E_0 = r_0 = 0$ 이라 즉시 0으로 나눈다. 선행 주소행렬식 하나가 0이면 행렬 전체가 멀쩡해도 깨지는 것이 피벗 없는 소거의 숙명이다. 부정치 토플리츠에는 문제가 되는 단계를 건너뛰는 룩어헤드(look-ahead) 변형이나, 구조를 버리고 번치-카우프만 피벗 같은 일반 대칭 부정치 분해를 쓴다.

6. 친척과 쓰임

• 격자 필터(lattice filter). 전방 오차 $f_m$ 과 후방 오차 $b_m$ 을 $f_m(t) = f_{m-1}(t) - k_m b_{m-1}(t-1)$, $b_m(t) = b_{m-1}(t-1) - k_m f_{m-1}(t)$ 로 한 단씩 쌓는 구조다. 단을 추가해도 앞 단이 바뀌지 않아 차수를 실시간으로 바꿀 수 있다. 재귀 최소제곱의 격자형 구현도 같은 뼈대다.
• 버그 방법. 자기공분산을 먼저 추정하지 않고 데이터에서 전·후방 오차 에너지 합을 최소화하도록 $k_m$ 을 바로 구한다. 짧은 데이터에서 해상도가 좋고 $\lvert k_m \rvert \le 1$ 이 자동 보장된다.
• 선형 예측 음성 부호화. 성도(vocal tract)를 단면적이 다른 관의 연결로 모형화하면 관 접합부의 음향 반사율이 바로 $k_m$ 이다 — 「반사 계수」라는 이름이 여기서 왔다. 미국 연방 표준 LPC-10 음성 코덱은 10차 예측기의 반사 계수를 양자화해 보냈다³. 선형 예측 부호화 계열 코덱의 핵심 루프다.
• 슈어 알고리즘. 내적 없이 반사 계수를 구하는 대안이다. 내적은 병렬화에 걸림돌이라 하드웨어 구현에서는 슈어 쪽이 선호되기도 한다.
• 스펙트럼 추정. 율-워커 계수로 만든 $\sigma^2/\lvert A(e^{i\omega})\rvert^2$ 가 최대 엔트로피(자기회귀) 파워 스펙트럼 추정이다.

7. 관련 문서

• 토플리츠 행렬 · 순환 행렬 · 촐레스키 분해 · 번치-카우프만 피벗 · 후진 오차 해석
• 위너 필터 · LMS 필터 · 재귀 최소제곱 · 칼만 필터
• 자기회귀 모형 · 선형 예측 부호화 · 파워 스펙트럼 · 상관함수

8. Footnotes

1. 레빈슨은 이 논문을 위너의 「황화」 보고서를 공학자들이 읽을 수 있게 풀어 쓰려는 해설로 썼는데, 정작 수치해석사에는 해설 대상보다 해설에 딸린 점화식이 더 오래 남았다. 부록이 본편을 이긴 사례다. ↩

2. 순수 사인파의 자기공분산 행렬은 크기가 사인파 개수의 두 배를 넘는 순간 특이해진다. 레빈슨-더빈은 이걸 $\lvert k_m \rvert \to 1$, $E_m \to 0$ 으로 정직하게 보고하는데, 이걸 「알고리즘이 터졌다」로 읽으면 억울한 누명이다. 신호가 그 차수에서 완벽히 예측 가능하다는 뜻이다. ↩

3. LPC-10 은 2400 bit/s 로 음성을 보냈다. 사람 목소리가 로봇처럼 들리긴 했지만, 반사 계수 10개와 피치·음량만으로 「누가 무슨 말을 하는지」는 전달됐다. 냉전기 보안 전화의 소리가 대체로 그런 느낌이었다. ↩