N7 · 순환·Toeplitz·띠 구조를 이용한 계산

1. 같은 규칙이 반복되는 측정값을 어떻게 계산할까

원형으로 배치한 네 측정점에서 각 값과 양옆 값을 결합한다고 합시다. 값은 같은 기준단위로 기록하고 입력 \(x_j\)는 기준상태에서의 변화량입니다. 출력은

\[ y_j=3x_j+x_{j-1}+x_{j+1},\qquad j=0,1,2,3 \]

입니다. 첨자는 4로 나눈 나머지로 읽어 \(x_{-1}=x_3\), \(x_4=x_0\)로 둡니다. 끝과 처음이 연결된다는 경계 가정이 행렬에 포함되어 있습니다.

\[\begin{split} C=\begin{pmatrix}3&1&0&1\\1&3&1&0\\0&1&3&1\\1&0&1&3\end{pmatrix},\qquad c=(3,1,0,1)^T. \end{split}\]

각 열은 첫 열 \(c\)를 한 칸씩 아래로 순환 이동한 것입니다. \(C_{ij}=c_{(i-j)\bmod4}\)처럼 첫 열만으로 전체를 정하는 행렬을 순환행렬이라고 합니다. \(x=(1,0,0,0)\)이면 출력이 \(c\)이고, 모든 입력이 1이면 출력은 모두 5입니다. 부호가 번갈아 드는 \((1,-1,1,-1)\)을 넣으면 자기 자신이 나와 고윳값이 1입니다. 완만한 변화와 번갈아 드는 변화에 서로 다른 배율을 주는 계산입니다.

그림 106 점은 네 개의 실제 좌표이고 연결선은 순서를 보여 준다. 좌우는 같은 세로축을 사용해 배율 차이를 비교했다. 처음과 끝을 연결하는 규칙은 본문의 순환 경계조건이다.

2. 푸리에 좌표에서는 곱셈이 성분별로 바뀐다

\(n\)개 좌표에 대해 단위 DFT 행렬을 \(F_{kj}=n^{-1/2}e^{-2\pi i kj/n}\)로 정합니다. 이 장에서는 순환행렬을 첫 열로 정의합니다. 첫 행 규약을 사용하면 주파수 부호가 바뀔 수 있어 두 규약을 섞지 않습니다.

\[ C=F^*\operatorname{diag}(\widehat c_k)F,\qquad \widehat c_k=\sum_{j=0}^{n-1}c_je^{-2\pi i kj/n}. \]

위 \(c=(3,1,0,1)\)에서는 \(\widehat c=(5,3,1,3)\)입니다. 예를 들어 \(k=1\)은 \(3+e^{-\pi i/2}+e^{-3\pi i/2}=3-i+i=3\)입니다. 허수 성분들이 상쇄되는 이유는 \(c_1=c_3\)인 실 대칭 구조입니다.

순환곱은 \(Cx=F^*(\widehat c\odot Fx)\)가 됩니다. NumPy의 기본 fft는 \(F\)보다 \(\sqrt n\)배 큰 정방향 변환이고 ifft에 \(1/n\) 정규화가 들어 있습니다. 따라서 코드에서는

\[ Cx=\operatorname{ifft}(\operatorname{fft}(c)\odot\operatorname{fft}(x)) \]

로 쓰며 여기에 \(\sqrt n\)을 추가로 곱하지 않습니다. NumPy FFT 규약의 부호와 정규화를 확인할 수 있습니다.

모든 \(\widehat c_k\ne0\)이면 역행렬도 순환입니다. 첫 열은 \(\operatorname{ifft}(1/\widehat c)\)입니다. 네 점 예의 역행렬 첫 열은 \((7/15,-1/5,2/15,-1/5)\)입니다. FFT가 나눗셈을 빠르게 해 주지만 작은 \(|\widehat c_k|\)의 민감도는 없애지 않습니다. 순환행렬은 정규이므로 \(\kappa_2(C)=\max|\widehat c_k|/\min|\widehat c_k|\)입니다.

3. 처음과 끝이 연결되지 않는 Toeplitz 행렬

시간에 따른 관측에서는 첫 시점과 마지막 시점을 이웃으로 놓는 가정이 보통 다릅니다. 같은 시차에 같은 계수를 주되 감싸지 않으면 Toeplitz 행렬입니다. 예를 들어

\[\begin{split} T=\begin{pmatrix}2&3&4\\1&2&3\\0&1&2\end{pmatrix} \end{split}\]

는 각 대각선의 값이 일정하지만 순환행렬은 아닙니다. 첫 열 \(a=(2,1,0)\)와 첫 행 \(b=(2,3,4)\)로 정해집니다.

이 행렬의 곱셈은 더 큰 순환행렬 안에 넣어 계산할 수 있습니다. 첫 열 \(c=(2,1,0,0,4,3)\)인 \(6\times6\) 순환행렬의 좌상단 \(3\times3\)이 정확히 \(T\)입니다. 입력을 \((x_1,x_2,x_3,0,0,0)\)으로 늘려 순환곱을 하고 앞 세 성분만 읽습니다. 아래쪽 좌표에 원래 입력을 넣지 않아 인위적인 끝 연결의 영향을 차단한 것입니다.

import numpy as np
import scipy.linalg as la
def toeplitz_product(first_col, first_row, x):
    a = np.asarray(first_col); b = np.asarray(first_row); x = np.asarray(x)
    n = len(a)
    if len(b) != n or len(x) != n or a[0] != b[0]:
        raise ValueError("크기와 공통 대각값을 확인하세요")
    c = np.r_[a, 0, b[:0:-1]]
    padded = np.r_[x, np.zeros(n)]
    return np.fft.ifft(np.fft.fft(c)*np.fft.fft(padded))[:n]

a = np.array([2., 1., 0.]); b = np.array([2., 3., 4.]); x = np.array([1., -2., 3.])
answer = toeplitz_product(a,b,x)
assert np.allclose(answer, [8., 6., 4.])
assert np.allclose(answer, la.toeplitz(a,b)@x)
print("패딩 FFT와 직접 Toeplitz 곱:", np.real_if_close(answer))

패딩 FFT와 직접 Toeplitz 곱: [8. 6. 4.]

이 방법은 곱셈을 빠르게 합니다. 큰 순환행렬의 역을 곱한 뒤 앞부분을 읽는다고 \(T^{-1}x\)가 되지는 않습니다. 블록 역행렬의 좌상단에는 Schur 보원이 들어가기 때문입니다. Toeplitz라는 구조만으로 모든 풀이가 즉시 \(O(n\log n)\)이 되었다고 말할 수 없습니다.

4. 양의 정부호 Toeplitz에서 예측 계수를 순서대로 구하기

실수 수열 \(\gamma_0,\gamma_1,\ldots\)로 \(\Gamma_n=(\gamma_{|i-j|})_{i,j=1}^n\succ0\)를 만든다고 합시다. 확률을 아직 사용하지 않고, 이 행렬을 같은 단위 변수들 사이의 주어진 양의 정부호 Gram 행렬로 생각합니다. 뒤에서 공분산으로 해석할 수 있지만 지금 필요한 계산은 선형계뿐입니다.

현재 값에서 앞 \(k\)개 값을 선형결합으로 빼는 오차의 제곱크기는

\[ v(a)=\gamma_0-2a^Tg_k+a^T\Gamma_ka,\qquad g_k=(\gamma_1,\ldots,\gamma_k)^T. \]

완전제곱을 만들면 최소계수는 \(\phi^{(k)}=\Gamma_k^{-1}g_k\), 최소값은 \(v_k=\gamma_0-g_k^T\phi^{(k)}\)입니다. 여기에서 계수 \(a_j\)는 \(j\)시점 전 변수에 곱하는 무차원 배율입니다.

매 \(k\)마다 새 행렬을 소거하지 않고 Levinson–Durbin 재귀를 사용합니다. \(v_0=\gamma_0\)에서 시작하고

\[ \alpha_k=\frac{\gamma_k-\sum_{j=1}^{k-1}\phi_j^{(k-1)}\gamma_{k-j}}{v_{k-1}},\qquad \phi_j^{(k)}=\phi_j^{(k-1)}-\alpha_k\phi_{k-j}^{(k-1)},\quad \phi_k^{(k)}=\alpha_k, \]

\[ v_k=v_{k-1}(1-\alpha_k^2). \]

이전 계수들을 역순으로 사용하게 되는 것은 Toeplitz 행렬이 시간 순서를 뒤집어도 같은 형태를 갖기 때문입니다. 마지막 절에서 블록 방정식의 첫 부분과 마지막 부분에 각각 대입하여 재귀를 증명합니다. 양의 정부호인 선행 행렬들에서는 \(v_k>0\)이므로 \(|\alpha_k|<1\)입니다. 이 부등식은 유한 Gram 행렬의 조건이며 임의의 추정 시계열이 정상이라는 통계적 결론을 대신하지 않습니다.

\(\gamma_j=\rho^j\), \(|\rho|<1\)을 예로 들면 \(\alpha_1=\rho\), \(\phi^{(1)}=(\rho)\), \(v_1=1-\rho^2\)입니다. 다음 분자는 \(\rho^2-\rho\rho=0\)이므로 \(\alpha_2=0\)입니다. 이후도 첫 계수만 \(\rho\), 나머지는 0이고 \(v_k=1-\rho^2\)입니다. 앞의 가장 최근 값 하나가 이 Gram 구조의 선형 예측을 결정합니다.

5. 분해 없이도 로그 행렬식과 이차형식을 구하기

관측 수치 \(y=(y_1,\ldots,y_n)\)에 대해

\[ e_t=y_t-\sum_{j=1}^{t-1}\phi_j^{(t-1)}y_{t-j},\qquad e_1=y_1 \]

를 정의합니다. 앞에서 설명 가능한 부분을 순서대로 뺀 값입니다. \(e=Ly\)인 단위하삼각 \(L\)을 쓰면 \(L\Gamma_nL^T=\operatorname{diag}(v_0,\ldots,v_{n-1})\)이므로

\[ \log\det\Gamma_n=\sum_{t=1}^n\log v_{t-1},\qquad y^T\Gamma_n^{-1}y=\sum_{t=1}^n\frac{e_t^2}{v_{t-1}}. \]

각 계수 재귀와 잔여 계산에 \(O(t)\)가 들어가므로 전체는 \(O(n^2)\)입니다. 이는 정규방정식을 만들지 않는다는 일반 구호가 아니라 Toeplitz의 역순 구조를 사용한 구체적인 비용 절약입니다. 모든 일반 우변의 Toeplitz 풀이를 이 예측계수 재귀와 혼동하지 않습니다.

def innovations(gamma, y):
    gamma = np.asarray(gamma, dtype=float); y = np.asarray(y, dtype=float)
    n = len(y)
    if len(gamma) < n or gamma[0] <= 0: raise ValueError("양의 초기 Gram 값이 필요합니다")
    phi = np.empty(0); variance = gamma[0]
    residual = np.empty(n); variances = np.empty(n); alphas = []
    residual[0] = y[0]; variances[0] = variance
    for k in range(1,n):
        alpha = (gamma[k]-phi@gamma[1:k][::-1])/variance
        phi = np.r_[phi-alpha*phi[::-1], alpha]
        variance *= 1-alpha*alpha
        if variance <= 0: raise ValueError("양의 정부호성이 유지되지 않았습니다")
        residual[k] = y[k]-phi@y[:k][::-1]
        variances[k] = variance; alphas.append(alpha)
    return residual, variances, np.array(alphas)

rho = .5; y = np.array([1., 0., 2., -1.]); gamma = rho**np.arange(len(y))
e, v, alpha = innovations(gamma, y)
Gamma = la.toeplitz(gamma)
assert np.allclose(e, [1., -.5, 2., -2.])
assert np.allclose(v, [1., .75, .75, .75])
assert np.isclose(np.sum(e*e/v), y@la.solve(Gamma,y))
assert np.isclose(np.sum(np.log(v)), np.linalg.slogdet(Gamma)[1])
print("순차 잔여, 양의 피벗, 로그 행렬식과 이차형식 확인")

순차 잔여, 양의 피벗, 로그 행렬식과 이차형식 확인

확률을 도입한 뒤 \(\Gamma_n\)을 평균 0인 Gaussian 관측의 공분산으로 두면 \(\frac12(\log\det\Gamma_n+y^T\Gamma_n^{-1}y)\)는 상수 \(n\log(2\pi)/2\)를 제외한 음의 로그우도입니다. 이 해석에는 분포 가정이 필요합니다. 현재까지의 행렬식·이차형식 항등식은 그 가정을 사용하지 않았습니다.

빠른 재귀도 반올림에서 무조건 안정하지는 않습니다. \(v_k\)가 매우 작으면 나눗셈이 민감해집니다. SciPy의 Toeplitz 풀이 설명 역시 Levinson 방식의 속도와 수치 안정성을 별개로 다룹니다.

6. Whittle의 순환 근사는 경계를 어떻게 바꾸는가

Toeplitz를 DFT로 대각화한 것처럼 처리하면 일반적으로 근사입니다. \(\gamma_j=\rho^{|j|}\)의 무한 푸리에 합은

\[ f(\omega)=\sum_{j\in\mathbb Z}\rho^{|j|}e^{-ij\omega} =\frac{1-\rho^2}{1-2\rho\cos\omega+\rho^2}>0. \]

이 장의 \(f\)에는 \(1/(2\pi)\)를 붙이지 않습니다. \(\omega_k=2\pi k/n\)에서 표본한 값을 고윳값으로 하는 순환행렬 \(C_n=F^*\operatorname{diag}(f(\omega_k))F\)를 만들면, 그 역 이차형식은

\[ y^TC_n^{-1}y=\frac1{1-\rho^2}\sum_{t=1}^n(y_t-\rho y_{t-1})^2,\qquad y_0=y_n. \]

정확한 Toeplitz 식에서는 첫 항이 \(y_1^2\)이고, 순환 근사에서는 \((y_1-\rho y_n)^2/(1-\rho^2)\)입니다. 따라서 차이는 정확히

\[ y^T\Gamma_n^{-1}y-y^TC_n^{-1}y =\frac{2\rho y_1y_n-\rho^2(y_1^2+y_n^2)}{1-\rho^2}. \]

경계의 두 값만으로 차이를 계산할 수 있습니다. \(n\)이 커진다는 사실만으로 이 차이가 0으로 가는 것은 아닙니다. 끝점 값들이 유계이면 차이는 \(O(1)\)이고 \(n\)으로 나눈 차이는 0으로 가지만, 끝점이 커지는 수열에서는 같은 결론을 자동으로 쓸 수 없습니다.

로그 행렬식도 정확히 비교할 수 있습니다.

\[ \log\det\Gamma_n=(n-1)\log(1-\rho^2),\qquad \log\det C_n=n\log(1-\rho^2)-2\log(1-\rho^n). \]

그 차이는 \(-\log(1-\rho^2)+2\log(1-\rho^n)\)여서 일반적으로 0이 아닌 상수로 갑니다. 유한표본의 정확한 우도와 순환 근사 우도를 같다고 써서는 안 됩니다. 더 일반적인 Szegő 극한과 Whittle 통계이론은 추가적인 수열·분포 가정을 갖는 외부 결과이며, 여기서는 위 AR(1)형 수열의 차이를 직접 증명합니다.

그림 107 \(\rho=0.8\)로 고정했다. 왼쪽은 전체 로그 행렬식 차이, 오른쪽은 \(n\)으로 나눈 차이이다. 양의 정부호 Toeplitz를 정확히 순환 대각화했다는 그림이 아니다.

7. 순환 임베딩으로 제곱근을 만들 때의 조건

곱셈용 임베딩은 대칭이어도 양의 준정부호일 필요가 없습니다. 목표 Toeplitz를 좌상단에 포함하는 실 대칭 순환행렬 \(C\)의 모든 FFT 고윳값이 비음수일 때에만

\[ B=F^*\operatorname{diag}(\sqrt{\widehat c_k})F,\qquad BB^T=C \]

로 실 제곱근을 만들 수 있습니다. 앞 \(n\)행을 \(B_n\)이라 하면 \(B_nB_n^T=\Gamma_n\)입니다. 여기까지는 결정론적인 Gram 인자분해입니다. 나중에 평균 0·공분산 \(I\)인 Gaussian \(z\)를 도입하면 \(B_nz\)가 목표 공분산을 갖는 표본이 됩니다. 실 입력 \(z\)에 실 대칭 \(B\)를 적용하면 복소 난수의 켤레대칭을 별도로 맞출 필요도 없습니다.

음의 고윳값을 그냥 0으로 자르면 목표 Gram 행렬도 바뀝니다. 더 큰 임베딩과 실제 수열의 추가 시차를 사용해 양의 준정부호성을 다시 검사해야 하며, 임의의 채움이나 패딩 확대가 항상 성공한다는 보장은 없습니다. FFT의 작은 허수부와 작은 음의 반올림값에 대한 허용오차도 스케일과 함께 기록합니다.

8. 시간 추세를 부드럽게 만드는 띠행렬

시간순 관측 \(y_t\)에서 부드러운 추세 \(\tau_t\)를 정하려면 두 번째 차분의 제곱을 벌점으로 줄 수 있습니다.

\[ \min_\tau\sum_{t=1}^n(y_t-\tau_t)^2 +\lambda\sum_{t=1}^{n-2}(\tau_t-2\tau_{t+1}+\tau_{t+2})^2. \]

두 합의 단위는 관측 단위의 제곱이며 시점 간격을 1로 고정했습니다. 시간 간격을 바꾸면 차분과 \(\lambda\)의 해석도 바꿔야 합니다. 차분행렬 \(D\)를 쓰면 해는 \(\tau=(I+\lambda D^TD)^{-1}y\)입니다. \(I\) 때문에 \(\lambda\ge0\)에서 항상 양의 정부호이고 반폭 2인 띠행렬이므로 N2의 띠 Cholesky로 \(O(n)\)에 풀 수 있습니다.

\(n=5\), \(\lambda=1\)이면

\[\begin{split} I+D^TD=\begin{pmatrix} 2&-2&1&0&0\\-2&6&-4&1&0\\1&-4&7&-4&1\\0&1&-4&6&-2\\0&0&1&-2&2 \end{pmatrix}. \end{split}\]

가운데 관측에만 1을 넣은 \(y=(0,0,1,0,0)\)의 해는 \((1/24,1/4,5/12,1/4,1/24)\)입니다. 첫 행은 \(2/24-2/4+5/12=0\), 가운데 행은 \(2/24-8/4+35/12=1\)이어서 원래 우변을 재현합니다. 완만한 선형 추세는 \(D\tau=0\)이므로 벌점을 받지 않습니다.

그림 108 왼쪽은 본문의 유한 경계 문제이고 점 사이 선은 시점 연결이다. 오른쪽은 주기 경계 또는 무한 내부에서 유도한 주파수식이다. 유한 행렬 전체가 그 식으로 정확히 DFT 대각화된다고 해석하지 않는다.

주기 경계에서는 차분의 주파수 절댓값 제곱이 \(|1-e^{i\omega}|^4=16\sin^4(\omega/2)\)여서 추세 이득이 \(h(\omega)=1/(1+16\lambda\sin^4(\omega/2))\)입니다. \(\lambda=1600\)에서 진폭 이득 \(1/2\)인 주기는 약 39.7분기입니다. 32분기에서는 약 0.297입니다. 따라서 이 값이 정확히 8년의 경계를 뜻하지는 않습니다. 전력은 이득의 제곱이므로 반전력 기준 \(|h|^2=1/2\)와도 구별해야 합니다. 마지막 두 시점의 처리가 다른 유한 표본에서는 끝점 효과가 추가됩니다.

9. Bartlett 가중치도 Gram 구조로 확인하기

\(L\ge0\)인 정수에 대해 \(W_{ij}=\max(1-|i-j|/(L+1),0)\)를 쓰면 대칭 Toeplitz 띠행렬입니다. 각 시점에 길이 \(L+1\)인 상자벡터를 놓고 \(1/\sqrt{L+1}\)로 정규화하면 두 상자의 겹치는 길이가 \(L+1-|i-j|\)입니다. 따라서 이 벡터들을 행으로 모은 \(R\)에 대해 \(W=RR^T\succeq0\)입니다.

스칼라 점수열 \(g_t\)의 가중 이차값은

\[ g^TWg=\sum_tg_t^2+2\sum_{h=1}^L\left(1-\frac h{L+1}\right) \sum_{t=1}^{n-h}g_tg_{t+h}\ge0 \]

이고 \(O(nL)\)에 계산합니다. 여러 변수 점수를 행으로 모은 \(G\)에서는 \(G^TWG\)가 양의 준정부호입니다. \(W\)가 Toeplitz라고 최종 변수별 행렬 \(G^TWG\)도 Toeplitz인 것은 아닙니다. HAC 공분산으로의 해석에는 평균 제거·정규화·시계열 가정이 추가되지만, 비음수 이차값의 근거는 이 Gram 인자분해입니다.

10. 연습과 전체 풀이

1. 네 점 순환행렬에서 \(C(1,0,-1,0)^T\)를 구하세요.

풀이. 성분을 계산하면 \((3,0,-3,0)^T\)이므로 배율은 3입니다. 같은 고윳값을 갖는 \((0,1,0,-1)\) 방향도 있습니다. 두 방향을 복소수로 묶으면 주파수 1과 3의 DFT 방향이 됩니다.

2. \(c=(1,1,1,1)\)의 순환행렬은 역행렬을 갖나요?

풀이. FFT 값은 \((4,0,0,0)\)이고 행렬의 모든 성분은 1입니다. 성분합이 0인 모든 벡터를 0으로 보내므로 계수 1이며 역행렬이 없습니다. FFT에서 0으로 나눌 수 없는 것이 정확한 영공간과 대응합니다.

3. Toeplitz 행렬 \(\begin{pmatrix}2&1\\1&2\end{pmatrix}\)를 첫 열 \((2,1,0,1)\)의 순환행렬에 넣으면 큰 행렬도 양의 정부호인가요?

풀이. 작은 행렬의 고윳값은 3과 1로 양수입니다. 큰 순환행렬의 고윳값은 \((4,2,0,2)\)이므로 특이한 양의 준정부호입니다. 작은 주부분행렬의 양의 정부호가 큰 임베딩의 엄격한 양의 정부호를 보장하지 않습니다.

4. \(\gamma=(1,1/2,1/5)\)의 두 단계 Levinson 계수와 오차값을 구하세요.

풀이. \(\alpha_1=1/2\), \(v_1=3/4\)입니다. \(\alpha_2=(1/5-1/4)/(3/4)=-1/15\)입니다. 첫 계수는 \(1/2-(-1/15)(1/2)=8/15\), 둘째는 \(-1/15\)입니다. \(v_2=(3/4)(1-1/225)=56/75\)입니다. 직접 \(\gamma_0-(\gamma_1,\gamma_2)(8/15,-1/15)^T=1-(4/15-1/75)=56/75\)로 확인합니다.

5. \(\rho=1/2\), \(y=(1,0,2,-1)\)의 정확한 Toeplitz 이차형식을 구하세요.

풀이. 잔여는 \((1,-1/2,2,-2)\)이고 분모는 \((1,3/4,3/4,3/4)\)입니다. 따라서 \(1+(1/4+4+4)/(3/4)=12\)입니다. 행렬식은 \((3/4)^3=27/64\)입니다.

6. 유한 HP 평활이 \(y_t=a+bt\)를 그대로 보존함을 보이세요.

풀이. \(y_t-2y_{t+1}+y_{t+2}=a+bt-2(a+b(t+1))+a+b(t+2)=0\)입니다. 따라서 \(Dy=0\), \((I+\lambda D^TD)y=y\)이고 유일한 해는 \(\tau=y\)입니다.

7. \(L=1\), \(g=(1,-1,1)\)의 Bartlett 이차값을 계산하세요.

풀이. 대각합은 3이고 시차 1의 곱합은 \(-2\)입니다. 두 배와 가중치 \(1/2\)를 곱해 더하면 \(3-2=1\)입니다. 개별 시차합은 음수일 수 있어도 전체 Gram 이차값은 비음수입니다.

11. 지금까지의 내용을 수학의 언어로 정리해 봅시다

이번 계산의 핵심은 같은 대각값·순환 경계·띠폭을 서로 다른 구조로 구분하는 것입니다. 확률적 해석에 앞서 성립하는 행렬 항등식을 먼저 완성합니다.

정리 1. 순환 DFT 대각화와 곱셈

\(C_{ij}=c_{(i-j)\bmod n}\)이면 2절의 DFT 분해가 성립합니다. 순환행렬들은 가환하며 곱과, 존재하는 경우 역행렬도 순환입니다.

증명. \(F^*\)의 \(k\)열을 \(v^{(k)}_j=n^{-1/2}e^{2\pi i jk/n}\)라 쓰겠습니다. \(l=i-j\pmod n\)로 합의 첨자를 바꾸면

\[ (Cv^{(k)})_i=\frac1{\sqrt n}\sum_lc_l e^{2\pi i(i-l)k/n} =v_i^{(k)}\sum_lc_le^{-2\pi i lk/n}=\widehat c_kv_i^{(k)}. \]

서로 다른 두 열의 내적은 \(n^{-1}\sum_{j=0}^{n-1}e^{2\pi i j(k-l)/n}=0\)입니다. 비율이 1이 아닌 \(n\)제곱근의 등비합이기 때문입니다. 같은 열은 길이 1입니다. 따라서 \(F\)가 유니터리이고 대각화식이 성립합니다. 같은 \(F\)로 대각화되므로 곱은 대각값의 곱, 역은 역수이며 다시 어떤 첫 열의 순환행렬입니다. 대각행렬들의 곱이 가환하므로 원래 행렬들도 가환합니다. ∎

FFT의 \(O(n\log n)\) 비용은 \(n=2^p\)에서 짝수·홀수 첨자 합으로 나누어 설명할 수 있습니다. 길이 \(n\)의 변환은 길이 \(n/2\)의 변환 두 개와 \(n\)개 이하의 위상 곱·덧셈으로 조립됩니다. 따라서 \(T(n)=2T(n/2)+O(n)\)이고 각 깊이의 비용 \(O(n)\)을 \(p\)번 합해 \(O(n\log n)\)입니다. 일반 길이는 선택한 FFT 구현의 인수분해·패딩 전략을 사용합니다.

정리 2. Toeplitz 순환 임베딩

첫 열 \(a\), 첫 행 \(b\) (\(a_0=b_0\))인 \(n\times n\) Toeplitz를 첫 열 \((a_0,\ldots,a_{n-1},0,b_{n-1},\ldots,b_1)\)의 \(2n\) 순환에 넣으면 좌상단 블록이 원래 행렬입니다. 따라서 3절의 패딩 곱이 정확합니다.

증명. \(0\le i,j<n\)에서 \(i\ge j\)이면 순환 첨자 \(i-j\in[0,n-1]\)이므로 성분이 \(a_{i-j}\)입니다. \(i<j\)이면 첨자는 \(2n-(j-i)\)이고 뒤쪽에 넣은 \(b_{j-i}\)를 읽습니다. 이것이 Toeplitz의 두 경우입니다. 나머지 입력 성분이 모두 0이므로 앞 \(n\)개의 출력에는 이 좌상 블록의 곱만 남습니다. ∎

정리 3. Levinson–Durbin과 양의 오차 피벗

\(\Gamma_{k+1}\succ0\)인 실 대칭 Toeplitz에서 4절의 재귀는 \(\Gamma_k\phi^{(k)}=g_k\)를 정확히 풀고 \(v_k=\gamma_0-g_k^T\phi^{(k)}>0\)를 줍니다.

증명. \(k=1\)은 \(\gamma_0\phi_1=\gamma_1\)입니다. \(k>1\)에서 \(J\)를 길이 \(k-1\)의 역순 행렬이라 두겠습니다. Toeplitz 대칭성에서 \(J\Gamma_{k-1}J=\Gamma_{k-1}\)이고 \(u=(\gamma_{k-1},\ldots,\gamma_1)^T=Jg_{k-1}\)입니다. 이전 해 \(\phi\)에 대해 \(\Gamma_{k-1}J\phi=J\Gamma_{k-1}\phi=u\)입니다.

새 후보를 \((\phi-\alpha J\phi,\alpha)\)로 쓰고 블록 행렬 \(\Gamma_k=\begin{pmatrix}\Gamma_{k-1}&u\\u^T&\gamma_0\end{pmatrix}\)에 곱합니다. 첫 블록은 \(g_{k-1}-\alpha u+\alpha u=g_{k-1}\)입니다. 마지막 성분은 \(u^T\phi+\alpha(\gamma_0-u^TJ\phi)\)입니다. \(u^TJ\phi=g_{k-1}^T\phi\)이므로 괄호는 \(v_{k-1}\)입니다. 이것을 \(\gamma_k\)와 같게 놓으면 바로 \(\alpha_k\)의 식입니다.

새 오차값은

\[ v_k=\gamma_0-g_{k-1}^T(\phi-\alpha_kJ\phi)-\gamma_k\alpha_k =v_{k-1}-\alpha_k(\gamma_k-u^T\phi) =v_{k-1}(1-\alpha_k^2). \]

각 \(v_k\)는 \(\Gamma_{k+1}\)에서 이전 좌표를 소거한 Schur 보원입니다. H7의 양의 정부호 Schur 정리로 양수이므로 나눗셈이 가능하고 \(|\alpha_k|<1\)입니다. 모든 단계에서 양의 피벗이면 뒤에서 증명할 합동 대각화를 거꾸로 써서 해당 선행 Toeplitz도 양의 정부호임을 확인할 수 있습니다. ∎

정리 4. 순차 잔여의 대각화

5절의 단위하삼각 \(L\)은 \(L\Gamma_nL^T=\operatorname{diag}(v_0,\ldots,v_{n-1})\)를 만족합니다. 따라서 로그 행렬식과 역 이차형식의 두 식이 성립합니다.

증명. \(\Gamma_n\)을 좌표벡터들의 내적으로 정의하면 \(t\)번째 잔여는 앞 \(t-1\)개 좌표의 최적 선형결합을 뺀 것입니다. 정규방정식에 의해 그 잔여는 앞 좌표들 각각과 직교합니다. 이전 잔여들도 앞 좌표의 선형결합이므로 새 잔여와 직교합니다. 잔여의 제곱길이는 정리 3의 \(v_{t-1}\)입니다. 이 내적들을 행렬로 적으면 합동 대각화식입니다.

\(L\)은 대각 1인 삼각이므로 \(\det L=1\)이고 \(\det\Gamma_n=\prod v_{t-1}\)입니다. 또한 \(\Gamma_n^{-1}=L^T\operatorname{diag}(v^{-1})L\)이므로 \(y\)를 좌우로 곱하면 \(\sum e_t^2/v_{t-1}\)입니다. ∎

정리 5. AR(1)형 Toeplitz와 순환 근사의 정확한 차이

\(|\rho|<1\), \(n\ge3\)에서 6절의 기호·두 로그 행렬식·경계 이차형식 차이가 성립합니다.

증명. 양의 시차와 음의 시차를 나누어 등비급수를 합하면 \(1+\rho e^{-i\omega}/(1-\rho e^{-i\omega})+\rho e^{i\omega}/(1-\rho e^{i\omega})=(1-\rho^2)/|1-\rho e^{i\omega}|^2\)입니다. 앞 재귀의 \(v_0=1\), \(v_k=1-\rho^2\)에서 Toeplitz 행렬식과 비주기 이차형식을 얻습니다.

순환행렬은 \(f(\omega_k)\)가 고윳값이므로 행렬식은 그 곱입니다. \(n\)제곱근에 대해 \(\prod_k(1-\rho e^{2\pi ik/n})=1-\rho^n\)이므로 \(\det C_n=(1-\rho^2)^n/(1-\rho^n)^2\)입니다. 또한 순환 역행렬의 기호는 \(|1-\rho e^{i\omega}|^2/(1-\rho^2)\)입니다. DFT 대각화와 Parseval 등식에서 역 이차형식은 순환 차분의 제곱합입니다. Toeplitz 식과 공통인 \(t=2,\ldots,n\) 항을 빼면 첫 경계항만 남고 전개하면 6절의 차이입니다. ∎

정리 6. 유한 평활과 Bartlett 양의 준정부호성

8절 평활문제는 \(\lambda\ge0\)에서 유일한 해를 갖고 제시한 띠행렬식으로 계산됩니다. 9절의 \(W\)는 양의 준정부호입니다.

증명. 임의의 \(z\ne0\)에서 \(z^T(I+\lambda D^TD)z=\|z\|^2+\lambda\|Dz\|^2>0\)입니다. 목적함수의 도함수는 \(2((I+\lambda D^TD)\tau-y)\)이므로 유일한 최소점이 그 선형계의 해입니다. 각 \(D\)행은 연속 세 좌표에만 비영이므로 \(D^TD\)는 반폭 2이고, N2의 띠 보존 논법으로 \(O(n)\) 풀이를 얻습니다.

Bartlett에는 \(R\in\mathbb R^{n\times(n+L)}\)를 \(i\le t\le i+L\)이면 \(R_{it}=1/\sqrt{L+1}\), 그 밖에는 0으로 정의합니다. 두 행의 지지집합은 \(\max(L+1-|i-j|,0)\)개 좌표에서 겹치므로 \((RR^T)_{ij}=W_{ij}\)입니다. 따라서 모든 \(g\)에서 \(g^TWg=\|R^Tg\|^2\ge0\)입니다. ∎

구조를 이용하면 전체 행렬을 그대로 저장하거나 분해하지 않아도 됩니다. 다음 장에서는 행렬과 벡터의 곱만으로 작은 탐색공간을 만들고 그 안에서 해를 개선합니다. N8로 이어 읽기.