N8 · 행렬곱으로 만드는 Krylov 반복법

1. 모든 계수를 저장하지 않고 평형을 구할 수 있을까

양끝을 고정한 막대의 세 내부 지점에서 평형 변위를 구한다고 합시다. 길이·강성·외력을 기준단위로 나누어 무차원화하면 교육용 차분 모형은

\[\begin{split} T_3x=b,\qquad T_3=\begin{pmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{pmatrix},\qquad b=e_1 \end{split}\]

입니다. 경계변위는 \(x_0=x_4=0\)이고 \(i\)번째 식은 \(2x_i-x_{i-1}-x_{i+1}=b_i\)입니다. 계수는 최근접 지점만 연결하며 변위와 외력의 부호를 허용합니다. 비선형 재료나 큰 변형은 이 선형 모형의 범위 밖입니다.

뒤 식부터 풀면 \(x_2=2x_3\), 가운데 식에서 \(x_1=3x_3\), 첫 식에서 \(4x_3=1\)입니다. 따라서 \(x_*=(3/4,1/2,1/4)^T\)입니다. 세 지점을 수만 개로 늘려도 행마다 필요한 값은 자기 자신과 두 이웃뿐입니다. \(Tx\)를 계산하는 함수를 제공하면 전체 \(n\times n\) 배열 없이 같은 연산을 할 수 있습니다.

이 장에서 반복의 품질은 원래 평형식의 잔차 \(r_k=b-Ax_k\)로 확인합니다. 변위오차 \(e_k=x_*-x_k\)와는 \(r_k=Ae_k\)로 연결되지만 같지 않습니다. 알고리즘 안에서 갱신하는 잔차와 원래 행렬곱으로 다시 계산한 잔차도 반올림 아래에서는 조금 다를 수 있습니다.

2. 이웃 값을 그대로 대입하는 정상반복

Jacobi는 \(x_i^{(k+1)}=(b_i+x_{i-1}^{(k)}+x_{i+1}^{(k)})/2\)처럼 모든 이웃의 이전 값을 씁니다. \(x_0=0\)에서 \(T_3\)의 첫 두 반복은 \((1/2,0,0)\), \((1/2,1/4,0)\)입니다. Gauss–Seidel은 같은 순회에서 이미 갱신한 왼쪽 값을 사용하므로 첫 반복이 \((1/2,1/4,1/8)\)입니다.

\(A=D+L+U\)를 대각·엄격 하삼각·엄격 상삼각으로 나누면 두 반복행렬은

\[ G_J=-D^{-1}(L+U),\qquad G_{GS}=-(D+L)^{-1}U. \]

일반 분할 \(A=M-N\)의 반복은 \(x_{k+1}=M^{-1}Nx_k+M^{-1}b\)입니다. \(A,M\)이 가역일 때 모든 초기값에서 같은 정확해로 수렴할 필요충분조건은 \(\rho(M^{-1}N)<1\)입니다. \(G=I\), \(c=0\)처럼 각 초기값에 머무는 수열도 수렴하는 수열이므로, 목표를 특정하지 않고 “수열이 수렴하면 \(\rho(G)<1\)”이라고 쓰면 틀립니다.

\(T_3\)에서는 \(\rho(G_J)=1/\sqrt2\), \(\rho(G_{GS})=1/2\)입니다. 더 큰 격자에서는 느린 공간 방향이 점점 잘 줄지 않아 반복 횟수가 커집니다. 한 번의 대입을 더 효율적으로 누적할 방법을 찾는 것이 다음 질문입니다.

3. 이미 계산한 벡터를 모두 이용하기

초기 잔차 \(r_0=b-Ax_0\)에서 \(r_0,Ar_0,A^2r_0,\ldots\)를 계산합니다. 앞 \(k\)개가 만드는 공간을

\[ \mathcal K_k(A,r_0)=\operatorname{span}\{r_0,Ar_0,\ldots,A^{k-1}r_0\} \]

라고 합니다. Krylov 부분공간입니다. 원소는 \(p_{k-1}(A)r_0\) 꼴의 다항식 작용으로 쓸 수 있습니다. 새 행렬곱 하나로 공간을 한 차원씩 늘리고 그 안에서 가장 좋은 보정을 선택합니다.

거듭제곱 벡터들은 거의 평행해질 수 있어 그대로 기저로 저장하지 않습니다. Arnoldi는 \(q_1=r_0/\|r_0\|\)에서 출발해 \(Aq_j\)를 기존 \(q_i\)들과 직교화합니다.

\[ Aq_j=\sum_{i=1}^jh_{ij}q_i+h_{j+1,j}q_{j+1},\qquad AQ_k=Q_{k+1}\overline H_k. \]

\(\overline H_k\)는 \((k+1)\times k\)의 상 Hessenberg 행렬입니다. \(h_{j+1,j}=0\)이면 새 방향이 없고 현재 공간이 \(A\)-불변입니다. 가역계의 GMRES에서는 이때 정확한 해를 공간 안에서 찾을 수 있는 경우를 행복한 중단이라고 부릅니다. 반올림에서 아주 작은 값은 수치적 종속인지 검사해야 합니다.

실 대칭 \(A\)에서는 \(H_k=Q_k^TAQ_k\)도 대칭이므로 Hessenberg와 합쳐 삼중대각이 됩니다. 이로써 Lanczos의 세 항 재귀 \(Aq_j=\beta_jq_{j-1}+\alpha_jq_j+\beta_{j+1}q_{j+1}\)를 얻습니다. 정확한 산술에서 이전 모든 방향과 직교하다는 사실에 의존합니다. 부동소수점에서는 오래된 방향 성분이 다시 생겨 이미 찾은 고윳값의 복사처럼 보이는 값이 나올 수 있습니다. 필요한 경우 N3의 재직교화를 사용하고 \(\|Q^TQ-I\|\)와 원래 고유쌍 잔차를 확인합니다.

4. 작은 행렬의 고윳값을 원래 공간에서 검증하기

\(H_ky=\theta y\), \(\|y\|=1\)이면 \(v=Q_ky\)를 원래 공간의 후보로 사용합니다. \(\theta\)를 Ritz 값이라고 합니다. Arnoldi 관계식에서

\[ Av-\theta v=h_{k+1,k}q_{k+1}e_k^Ty,\qquad \|Av-\theta v\|_2=|h_{k+1,k}|\,|e_k^Ty| \]

입니다. 따라서 작은 문제를 푼 뒤 원래 차원의 잔차를 싸게 인증할 수 있습니다. 실제 계산에서는 \(Q\)의 직교성 손실과 Arnoldi 관계의 잔차까지 확인합니다.

대칭행렬의 극단 고윳값에는 Rayleigh 몫의 최소·최대 성질이 있어 부분공간을 늘릴수록 극단 Ritz 값이 바깥쪽 정확한 극단값으로 향하는 단조성이 있습니다. 그러나 모든 행렬에서 “10번이면 양끝이 먼저 정확하다” 같은 고정 횟수는 보장되지 않습니다. 시작벡터의 성분, 간극과 다항식 근사 가능성에 따라 달라집니다.

그림 109 가로축은 부분공간 차수이며 세로축은 절대 크기의 로그 눈금이다. 참 최대 고윳값 80과의 오차와 재계산한 잔차를 함께 보인다. 잔차가 인증하는 것은 가까운 어떤 고윳값이며, 초기부터 최대 고윳값의 오차를 상계하는 것은 아니다. 두 번 직교화한 Arnoldi 구현을 사용했으며, \(10^{-15}\)보다 작은 결과는 그 표시 하한에 놓았다.

5. 켤레기울기법은 어떤 오차를 가장 작게 만드는가

\(A=A^T\succ0\)인 평형계는 에너지 \(\varphi(x)=\frac12x^TAx-b^Tx\)의 최소점입니다. 정확해 \(x_*=A^{-1}b\)에 대해

\[ \varphi(x)-\varphi(x_*)=\frac12(x-x_*)^TA(x-x_*) =\frac12\|x-x_*\|_A^2. \]

CG는 \(x_0+\mathcal K_k(A,r_0)\) 안에서 이 에너지 오차를 최소화합니다. 잔차의 유클리드 노름을 매번 최소화하는 방법은 아니므로 \(\|r_k\|_2\)가 항상 단조 감소한다고 주장하지 않습니다.

긴 공간 최적화를 짧은 재귀로 실행할 수 있습니다. \(p_0=r_0\)에서

\[ \alpha_k=\frac{r_k^Tr_k}{p_k^TAp_k},\quad x_{k+1}=x_k+\alpha_kp_k,\quad r_{k+1}=r_k-\alpha_kAp_k, \]

\[ \beta_k=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k},\qquad p_{k+1}=r_{k+1}+\beta_kp_k. \]

서로 다른 탐색방향은 \(p_i^TAp_j=0\)이고 잔차들은 보통 내적에서 직교합니다. 두 종류의 직교성을 혼동하지 않습니다. 마지막 절에서 두 성질과 재귀가 왜 공간 최소점을 만드는지 함께 증명합니다.

처음 \(T_3\)에서 \(x_0=0\)으로 직접 계산하면 다음과 같습니다.

단계	\(x_k\)	\(r_k\)	다음 탐색방향 \(p_k\)
0	\((0,0,0)\)	\((1,0,0)\)	\((1,0,0)\)
1	\((1/2,0,0)\)	\((0,1/2,0)\)	\((1/4,1/2,0)\)
2	\((2/3,1/3,0)\)	\((0,0,1/3)\)	\((1/9,2/9,1/3)\)
3	\((3/4,1/2,1/4)\)	\((0,0,0)\)	종료

첫 보폭은 \(1/2\)입니다. 둘째는 \((1/4)/(3/8)=2/3\), 셋째는 \((1/9)/(4/27)=3/4\)입니다. 세 번의 행렬곱으로 정확한 해에 도달했습니다. 정확한 산술의 유한종료와 실제 계산에서의 오차 허용 종료를 구별합니다.

6. 조건수와 고윳값 군집이 반복 횟수에 미치는 영향

\(\kappa=\lambda_{\max}(A)/\lambda_{\min}(A)\)이면 정확한 CG는

\[ \|e_k\|_A\le2\left(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\right)^k\|e_0\|_A \]

를 만족합니다. 실제 스펙트럼이 같은 구간 안에서도 몇 곳에 모여 있으면 이 구간 전체의 최악 상한보다 빨라질 수 있습니다. 그러나 “군집이 \(m\)개이면 \(m\)번에 끝난다”는 정확한 명제는 아닙니다. 서로 다른 정확한 고윳값이 \(m\)개이면 정확산술에서 \(m\)번 이내 종료하고, 폭을 가진 군집에서는 폭·간격·시작오차에 따라 추가 반복이 필요합니다.

전처리에서는 쉽게 풀 수 있는 \(M\succ0\)를 사용해 \(M^{-1/2}AM^{-1/2}\)의 조건수를 개선합니다. \(M^{-1}A\)는 보통 유클리드 내적에서 대칭이 아니므로 원래 CG의 대칭 가정을 그 행렬에 직접 붙이면 안 됩니다. 대칭 변환으로 분석하고 구현에서는 \(Mz=r\)만 풉니다.

Poisson 행렬 \(T_n\)의 대각은 모두 2여서 Jacobi 전처리 \(M=2I\)는 스칼라배일 뿐입니다. 이 경우 원래 CG와 반복 횟수를 개선하지 않습니다. 반면 \(A=DT_nD\), \(D\)가 매우 다른 양의 대각을 갖는 경우에는 \(M=2D^2\)로 단위 불균형을 제거할 수 있습니다. 다음 그림은 이 두 문제 중 스케일이 다른 후자를 사용합니다.

그림 110 행렬은 \(DT_nD\)이고 \(D\)의 대각 범위는 1에서 20이다. 가로축은 10회까지 선형이고 그 이후 로그 간격이며, 세로축은 로그 눈금이다. 모든 곡선은 원래 식 \(b-Ax\)로 다시 계산한 같은 잔차를 사용한다. 전처리된 잔차와 원래 잔차를 서로 다른 기준으로 비교하지 않는다.

7. 직접 실행하는 무행렬 CG

import numpy as np
def poisson_product(x):
    y = 2*np.asarray(x, dtype=float).copy()
    y[:-1] -= x[1:]; y[1:] -= x[:-1]
    return y

def cg_operator(matvec, b, solve_m=None, rtol=1e-11, atol=0., limit=None):
    b = np.asarray(b, dtype=float)
    x = np.zeros_like(b); r = b-matvec(x)
    threshold = atol+rtol*np.linalg.norm(b)
    history = [np.linalg.norm(r)]
    if history[-1] <= threshold: return x, np.array(history)
    if solve_m is None: solve_m = lambda r: r.copy()
    z = solve_m(r); p = z.copy(); rho = r@z
    if rho <= 0: raise ValueError("양의 정부호 전처리자가 필요합니다")
    if limit is None: limit = 10*len(b)
    for _ in range(limit):
        Ap = matvec(p); curvature = p@Ap
        if curvature <= 0: raise ValueError("양의 곡률이 아닙니다")
        alpha = rho/curvature
        x += alpha*p; r -= alpha*Ap
        # 교육용 검증: 매 단계 원래 식의 잔차를 별도로 계산합니다.
        true_residual = b-matvec(x); history.append(np.linalg.norm(true_residual))
        if history[-1] <= threshold: return x, np.array(history)
        z = solve_m(r); new_rho = r@z
        if new_rho <= 0: raise RuntimeError("재귀 잔차를 다시 점검해야 합니다")
        p = z+(new_rho/rho)*p; rho = new_rho
    raise RuntimeError("CG 반복 한도 내에 목표 잔차에 도달하지 못했습니다")

x, history = cg_operator(poisson_product, np.array([1.,0.,0.]))
assert np.allclose(x,[.75,.5,.25],atol=1e-14)
assert len(history)-1 == 3
print("세 점 평형의 CG 해와 반복 횟수:", x, len(history)-1)

세 점 평형의 CG 해와 반복 횟수: [0.75 0.5  0.25] 3

위 검증용 구현은 매번 잔차를 다시 계산하므로 단계당 행렬곱이 하나 더 듭니다. 일반 구현에서는 일정 간격으로 검사할 수 있습니다. 재귀 잔차가 아주 작아도 실제 잔차가 허용치보다 크면 종료를 선언하지 않습니다. 오차를 더 줄일 수 없는 경우에는 반복 횟수만 늘리지 말고 정밀도·전처리·문제의 조건수를 다시 봅니다.

8. 비대칭 또는 부정치에서는 어떤 최소화를 사용할까

GMRES는 \(x_k\in x_0+\mathcal K_k(A,r_0)\) 중 \(\|b-Ax_k\|_2\)를 최소화합니다. Arnoldi를 대입하면 큰 최적화가

\[ \min_y\|\beta e_1-\overline H_ky\|_2,\qquad \beta=\|r_0\|,quad x_k=x_0+Q_ky \]

로 줄어듭니다. 모든 이전 기저를 보관하므로 비용과 저장이 늘며, 재시작 GMRES는 일정 차수마다 공간을 버리는 대신 수렴이 달라질 수 있습니다.

고윳값만으로 이 수렴을 결정할 수 없다는 작은 반례를 보겠습니다.

\[\begin{split} C=\begin{pmatrix}0&0&0&-1\\1&0&0&4\\0&1&0&-6\\0&0&1&4\end{pmatrix},\qquad b=e_1. \end{split}\]

특성다항식은 \((z-1)^4\)라서 모든 고윳값이 1입니다. 그러나 \(Ce_1=e_2\), \(C^2e_1=e_3\), \(C^3e_1=e_4\)입니다. \(k<4\)에서 \(C\mathcal K_k(C,e_1)=\operatorname{span}(e_2,\ldots,e_{k+1})\)는 \(b\)와 직교합니다. 따라서 어떤 보정을 골라도 잔차 노름의 최소는 1이고 \(x_k=0\)입니다. 4차원 전체를 사용해야 정확해를 얻습니다. 같은 고윳값 목록을 갖는 \(I_4\)는 한 번이면 끝납니다.

그림 111 점은 작은 Krylov 최소제곱을 실제로 풀어 얻은 단계별 잔차이다. 정확산술의 마지막 잔차는 0이며, 계산값은 반올림 범위에서 0에 가깝다. 선형 세로축에 계산값을 그대로 표시했다. 동반행렬의 정체는 고윳값을 수치적으로 정렬한 결과가 아니라 본문의 기저 계산으로 증명했다.

대칭 부정치에서는 CG 분모 \(p^TAp\)가 양수라는 보장이 없어 MINRES처럼 잔차를 최소화하는 방법을 사용합니다. 직사각 최소제곱에는 \(Av\)와 \(A^Tu\)만 사용하는 Golub–Kahan 과정에서 작은 이중대각 최소제곱을 푸는 LSQR이 있습니다.

\(AV_k=U_{k+1}B_k\)라면 \(x=V_ky\)의 목적함수는 \(\|\beta e_1-B_ky\|\)입니다. 정규행렬을 명시적으로 만들지 않지만 원래 최소제곱의 조건수나 잡음 민감도가 없어지는 것은 아닙니다. LSMR은 관련 과정에서 정규방정식 잔차를 줄이는 기준을 사용합니다. 상세한 짧은 재귀와 각 방법의 가정은 Netlib의 반복법 교재를 외부 조망으로 두고, 이 장의 핵심 CG·Arnoldi·GMRES 최소화는 아래에서 증명합니다.

9. 희소 저장과 전처리는 각각 어떤 비용을 줄이는가

\(T_3\)의 CSR 저장은 값 목록 \((2,-1,-1,2,-1,-1,2)\), 열 첨자 \((0,1,0,1,2,1,2)\), 행 시작 위치 \((0,2,5,7)\)입니다. \(i\)행은 시작 위치 사이의 원소만 읽습니다. CSC는 같은 정보를 열 기준으로 저장합니다. 이러한 형식은 행렬·벡터곱을 \(O(\mathrm{nnz})\)로 실행하고, 무행렬 함수는 더 규칙적인 스텐실에서 첨자 배열도 생략합니다.

전처리자는 좋은 스펙트럼 또는 좋은 근사 작용과 싼 적용비용을 함께 가져야 합니다. \(M=A\)는 한 단계에 풀게 하지만 \(M\)을 푸는 일이 원래 문제입니다. \(M=I\)는 적용이 쉽지만 개선이 없습니다. IC는 대칭 양의 정부호에 맞춘 불완전 Cholesky이고, ILU는 일반적으로 비대칭이므로 이를 아무 확인 없이 PCG에 넣으면 가정을 잃습니다. IC가 양의 피벗으로 완료되는지도 확인해야 합니다. Nyström 등 저계수 전처리는 저계수 근사 품질과 적용비용을 함께 평가하는 후속 선택지입니다.

희소한 원래 행렬도 분해 중 fill-in으로 조밀해질 수 있습니다. RCM·AMD 같은 순열은 행렬의 고윳값을 바꾸기보다 저장과 분해 비용을 바꾸는 데 목적이 있습니다. 어떤 순열이든 항상 비영 성분을 일정 배수만큼 줄인다는 보장은 없습니다. N2의 별 모양 그래프 예와 연결해서 읽을 수 있습니다.

10. 대각을 전부 구하지 않고 자취와 로그 행렬식 계산하기

대칭 \(B\)에 대해 모든 부호벡터 \(z\in\{-1,1\}^n\)의 유한 평균을 취하면

\[ 2^{-n}\sum_z z^TBz=\operatorname{tr}B. \]

대각 성분은 \(z_i^2=1\)로 그대로 남고, 비대각은 \(z_iz_j\)의 부호가 절반씩 상쇄됩니다. 일부 부호만 고른 평균이 Hutchinson 자취 추정의 출발점입니다. 이 단원에서는 전체 유한 평균과 제곱편차를 증명하며, 무작위 표본의 신뢰구간은 확률 단원에서 추가합니다.

대칭 \(W\)와 \(|\rho|\|W\|_2=a<1\)이면 \(I-\rho W\succ0\)이고

\[ \log\det(I-\rho W)=-\sum_{j=1}^\infty\frac{\rho^j}{j}\operatorname{tr}(W^j). \]

각 프로브에 \(W\)를 반복해서 곱하면 \(z^TW^jz\)를 얻어 큰 거듭제곱 행렬을 만들지 않습니다. \(K\)차 절단의 절대오차는 \(na^{K+1}/((K+1)(1-a))\) 이하입니다. 일부 프로브만 쓰는 오차와 다항식 절단오차는 별개입니다. Chebyshev 근사나 Lanczos 구적은 필요한 스펙트럼 구간과 함수 근사오차를 바탕으로 이 계산을 개선하는 후속 방식입니다.

네 점 순환 그래프의 이웃 평균 \(W\)는 고윳값이 \(1,0,0,-1\)입니다. 따라서 정확한 로그 행렬식은 \(\log(1-\rho^2)\)입니다. 이 값을 기준으로 절단과 부호평균을 각각 검증할 수 있습니다.

그림 112 \(\rho=0.8\)이며 왼쪽은 16개 부호를 모두 열거했다. 오른쪽은 정확한 자취를 사용한 급수 절단오차이므로 프로브 표본오차가 포함되지 않는다. 오른쪽 세로축은 로그 눈금이다.

11. 연습과 전체 풀이

1. \(T_3\)의 첫 CG 보정이 왜 \(e_1/2\)인지 에너지에서 구하세요.

풀이. 첫 공간은 \(x=te_1\)입니다. \(\varphi(te_1)=t^2-t=(t-1/2)^2-1/4\)이므로 최소 \(t=1/2\)입니다. 이는 \(\alpha_0=(r_0^Tr_0)/(p_0^TAp_0)=1/2\)와 같습니다.

2. 같은 문제에서 \(p_0^TT_3p_1=0\)을 확인하세요.

풀이. \(p_1=(1/4,1/2,0)\)에 곱하면 \(T_3p_1=(0,3/4,-1/2)\)입니다. \(p_0=e_1\)와의 내적은 0입니다. 반면 \(p_0^Tp_1=1/4\)이므로 보통 내적에서는 직교하지 않습니다.

3. \(A=\operatorname{diag}(2,2,5,5)\)인 CG의 정확산술 종료 횟수를 제한하세요.

풀이. \(q(t)=(1-t/2)(1-t/5)\)는 \(q(0)=1\)이고 모든 고윳값에서 0입니다. 따라서 \(q(A)=0\)이고 두 단계 Krylov 보정으로 오차를 0으로 만들 수 있습니다. CG의 최소성에서 최대 두 번입니다. 시작오차가 한 고유공간에만 있으면 한 번입니다.

4. \(M=2I\)가 \(T_n\)의 CG 조건수를 바꾸지 않는 이유를 쓰세요.

풀이. 대칭 변환은 \(M^{-1/2}T_nM^{-1/2}=T_n/2\)입니다. 최대·최소 고윳값을 같은 수로 나누므로 비는 같습니다. 스칼라 전처리는 이 예의 스펙트럼 분포를 개선하지 않습니다.

5. 동반행렬 \(C\)의 GMRES가 \(k=2\)에서 잔차를 줄이지 못함을 직접 쓰세요.

풀이. 보정 \(x=ae_1+be_2\)이면 \(Cx=ae_2+be_3\)입니다. 잔차 제곱은 \(\|e_1-ae_2-be_3\|^2=1+a^2+b^2\)이며 최소는 \(a=b=0\), 값 1입니다.

6. \(B=\begin{pmatrix}2&1\\1&3\end{pmatrix}\)의 네 부호 프로브 값과 평균을 구하세요.

풀이. \(z^TBz=5+2z_1z_2\)이므로 같은 부호 두 벡터에서는 7, 반대 부호 두 벡터에서는 3입니다. 평균은 5로 자취와 같고 평균 제곱편차는 4입니다. 공식 \(2\sum_{i\ne j}b_{ij}^2=4\)와 일치합니다.

7. \(\rho=1/2\)인 네 점 그래프에서 로그 행렬식의 4차 절단을 계산하세요.

풀이. \(\operatorname{tr}W=\operatorname{tr}W^3=0\), \(\operatorname{tr}W^2=\operatorname{tr}W^4=2\)입니다. 따라서 \(-\rho^2-\rho^4/2=-1/4-1/32=-9/32\)입니다. 정확한 값은 \(\log(3/4)\)이고 남은 음의 짝수항들 때문에 절단값이 더 큽니다.

8. \(G=I\), \(c=0\)가 정상반복 수렴 명제에 주는 주의를 설명하세요.

풀이. 모든 수열은 \(x_k=x_0\)이므로 수렴하지만 초기값마다 극한이 다르고 \(\rho(G)=1\)입니다. “모든 초기값에서 정해진 유일한 해로 수렴한다”는 조건이어야 \(G^k\to0\)과 동치가 됩니다.

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

다항식이 오차에 작용하는 방식, 부분공간에서의 최적성, 짧은 재귀가 같은 알고리즘을 설명합니다. 증명에서는 정확한 산술을 사용하고 부동소수점 검증은 원래 잔차와 별도로 연결합니다.

정리 1. 정상반복과 대각우세·양의 정부호

\(A=M-N\)이고 \(A,M\)이 가역이면 모든 초기값에서 \(A^{-1}b\)로 수렴할 필요충분조건은 \(\rho(M^{-1}N)<1\)입니다. 엄격한 행 대각우세에서는 Jacobi와 Gauss–Seidel이 수렴합니다. 대칭 양의 정부호에서는 Gauss–Seidel과 \(0<\omega<2\)인 SOR가 수렴합니다.

증명. 정확해를 반복식에서 빼면 \(e_{k+1}=M^{-1}Ne_k\)입니다. 모든 \(e_0\)에서 0으로 가는 것은 H8의 거듭제곱 수렴에서 스펙트럼 반지름 1 미만과 동치입니다.

엄격한 행 대각우세에서는 \(\|G_J\|_\infty=\max_i\sum_{j\ne i}|a_{ij}|/|a_{ii}|<1\)입니다. Gauss–Seidel에 절댓값 1 이상인 고윳값 \(\mu\)와 고유벡터 \(z\)가 있다고 가정하자. \(|z_i|=\|z\|_\infty>0\)인 행에서 \(\mu(D+L)z=-Uz\)를 쓰고 \(|\mu|\)로 나누면

\[ |a_{ii}||z_i|\le\sum_{j<i}|a_{ij}||z_j|+|\mu|^{-1}\sum_{j>i}|a_{ij}||z_j| \le\sum_{j\ne i}|a_{ij}|\|z\|_\infty<|a_{ii}||z_i|, \]

모순입니다.

SPD에서는 오차 에너지 \(e^TAe\)를 사용합니다. 좌표 \(i\)를 \(e_i\leftarrow e_i-\omega(Ae)_i/a_{ii}\)로 바꾸면 제곱을 전개한 에너지 변화는 \(-\omega(2-\omega)(Ae)_i^2/a_{ii}\le0\)입니다. 한 순회에서 감소가 전혀 없다면 모든 갱신이 0이고 \(Ae=0\), 따라서 \(e=0\)입니다. 에너지 단위구면의 콤팩트성과 선형 반복의 연속성에서 한 순회의 최대 에너지 비가 1 미만입니다. 따라서 그 에너지 노름에서 축약이며 수렴합니다. \(\omega=1\)이 Gauss–Seidel입니다. ∎

정리 2. Arnoldi 관계, Ritz 잔차와 GMRES

직교화가 중단되지 않은 \(k\)단계 Arnoldi는 \(\mathcal K_k\)의 정규직교기저와 3절의 관계를 만듭니다. Ritz 잔차와 GMRES의 작은 최소제곱 식은 각각 4·8절의 식입니다.

증명. 첫 벡터는 \(r_0\)를 정규화한 것입니다. \(q_j\in\mathcal K_j\)이면 \(Aq_j\in\mathcal K_{j+1}\)이고 앞 기저를 뺀 잔여도 그 공간에 있습니다. 잔여가 비영이면 앞 공간과 직교하는 새 방향을 얻습니다. 생성 벡터 수와 차원을 귀납하면 같은 Krylov 공간의 기저입니다. 각 열의 사영식을 나란히 붙이면 Arnoldi 관계입니다. 대칭일 때 \(h_{ij}=q_i^TAq_j=h_{ji}\)이므로 위 Hessenberg는 삼중대각입니다.

\(H_ky=\theta y\)를 관계에 넣으면 첫 \(k\)행의 차는 0이고 마지막 항만 \(h_{k+1,k}q_{k+1}e_k^Ty\)로 남습니다. \(q_{k+1}\)의 길이가 1이므로 잔차 노름식이 나옵니다. GMRES에는 \(x=x_0+Q_ky\), \(r_0=\beta q_1\)을 대입하면 \(b-Ax=Q_{k+1}(\beta e_1-\overline H_ky)\)이고 \(Q_{k+1}\)가 길이를 보존하므로 두 최소화가 같습니다. ∎

정리 3. CG 재귀와 에너지 최적성

\(A=A^T\succ0\)에서 5절의 CG는 정확해를 만나기 전에는 양의 분모를 갖고 \(x_0+\mathcal K_k(A,r_0)\) 위의 유일한 에너지 최소점을 만듭니다. 잔차들은 서로 직교하고 탐색방향들은 \(A\)-직교합니다.

증명. 어느 부분공간 \(S\)에서 \(x\)가 에너지 최소일 필요충분조건은 모든 \(p\in S\)에 \(p^T(Ax-b)=0\)인 것입니다. 도함수를 0으로 놓아 필요성을 얻고, 임의 보정 \(p\)의 비용 차이가 \(-p^Tr+\frac12p^TAp\)임을 쓰면 충분성과 유일성도 얻습니다.

귀납적으로 \(p_0,\ldots,p_{k-1}\)가 \(\mathcal K_k\)의 \(A\)-직교기저이고 \(r_k\perp\mathcal K_k\)라 합시다. \(r_k\in\mathcal K_{k+1}\)이며 \(r_k\ne0\)이면 기존 공간 밖이므로 새 방향이 됩니다. \(r_k\)를 이전 \(p_i\)들에 대해 \(A\)-직교화합니다. \(i\le k-2\)에서는 \(Ap_i=(r_i-r_{i+1})/\alpha_i\)이고 두 잔차가 \(\mathcal K_k\)에 있으므로 \(r_k^TAp_i=0\)입니다. 마지막 방향에만 보정이 필요하며

\[ -\frac{r_k^TAp_{k-1}}{p_{k-1}^TAp_{k-1}} =\frac{\|r_k\|^2}{\alpha_{k-1}p_{k-1}^TAp_{k-1}} =\frac{\|r_k\|^2}{\|r_{k-1}\|^2}. \]

여기서 \(r_k\perp r_{k-1}\)와 이전 보폭식을 썼습니다. 따라서 새 방향은 정확히 재귀의 \(p_k\)입니다. \(p_k\ne0\)와 SPD에서 \(p_k^TAp_k>0\)입니다. 이전 방향 성분과 \(r_k\)가 직교하므로 \(p_k^Tr_k=\|r_k\|^2\)이고, 그 직선의 최소 보폭은 \(\alpha_k=\|r_k\|^2/(p_k^TAp_k)\)입니다.

새 잔차 \(r_{k+1}=r_k-\alpha_kAp_k\)는 \(p_k\)와 직교하고, 이전 \(p_i\)와도 \(r_k\perp p_i\), \(p_i^TAp_k=0\) 때문에 직교합니다. 따라서 확장 공간 \(\mathcal K_{k+1}\)에서의 최소점입니다. 이전 잔차들은 그 공간에 속하므로 새 잔차와도 직교합니다. \(k=0\)은 \(p_0=r_0\)로 시작해 같은 보폭 계산을 적용하면 됩니다. ∎

정리 4. 다항식 오차, 유한종료와 Chebyshev 상한

CG 오차는 \(\|e_k\|_A=\min_{q\in\mathcal P_k,\ q(0)=1}\|q(A)e_0\|_A\)입니다. 서로 다른 고윳값이 \(d\)개이면 최대 \(d\)단계에 정확히 종료하고 6절의 조건수 상한을 만족합니다.

증명. 허용 보정은 \(p_{k-1}(A)r_0=p_{k-1}(A)Ae_0\)입니다. 따라서 오차는 \((I-Ap_{k-1}(A))e_0=q(A)e_0\)이고 \(q(0)=1\)입니다. 역으로 이런 \(q\)는 \(1-q(t)\)가 \(t\)로 나누어져 같은 보정 형태입니다. 정리 3의 최소성을 적용합니다. \(q(t)=\prod_{j=1}^d(1-t/\lambda_j)\)이면 \(q(A)=0\)이므로 유한종료입니다.

\(0<a=\lambda_{\min}<b=\lambda_{\max}\)라 합시다. 고유기저에서 \(\|q(A)e_0\|_A\le\max_{t\in[a,b]}|q(t)|\|e_0\|_A\)입니다. \(T_k(\cos\theta)=\cos(k\theta)\)인 Chebyshev 다항식을 쓰고

\[ q_k(t)=\frac{T_k((a+b-2t)/(b-a))}{T_k((a+b)/(b-a))} \]

로 놓으면 \(q_k(0)=1\)이며 분자 절댓값은 구간에서 1 이하입니다. 코사인 덧셈식은 \(T_{k+1}(z)=2zT_k(z)-T_{k-1}(z)\)를 줍니다. 이 재귀와 초깃값을 확인하면 \(z>1\)에서 \(T_k(z)=\frac12((z+\sqrt{z^2-1})^k+(z-\sqrt{z^2-1})^k)\)입니다. \(z=(a+b)/(b-a)\)에서는 \(z+\sqrt{z^2-1}=(\sqrt b+\sqrt a)/(\sqrt b-\sqrt a)\)입니다. 둘째 항도 양수이므로 분모는 첫 항의 절반 이상이고, 역수를 취하면 \(2((\sqrt\kappa-1)/(\sqrt\kappa+1))^k\)를 얻습니다. \(a=b\)에서는 \(A=aI\)이고 한 번에 종료합니다. ∎

정리 5. 대칭 전처리의 의미

\(M\succ0\)인 PCG는 \(\widetilde A=M^{-1/2}AM^{-1/2}\)의 보통 CG와 동치이며 수렴 상한에는 \(\kappa_2(\widetilde A)\)가 들어갑니다.

증명. \(y=M^{1/2}x\), \(\widetilde b=M^{-1/2}b\)로 바꾸면 \(Ax=b\)는 \(\widetilde Ay=\widetilde b\)가 됩니다. 이 행렬은 대칭이고 \(v^T\widetilde Av=(M^{-1/2}v)^TA(M^{-1/2}v)>0\)입니다. 변환된 잔차는 \(M^{-1/2}r\)여서 그 제곱노름이 \(r^TM^{-1}r=r^Tz\)입니다. 변환된 방향을 다시 \(x\)좌표로 돌리면 보폭 분모는 \(p^TAp\)이고 재귀는 7절의 PCG 코드와 같습니다. 따라서 앞 CG 정리를 변환 문제에 적용합니다. ∎

정리 6. 유한 부호평균과 로그급수 오차

대칭 \(B\)에서 부호벡터 전체의 평균은 자취이고 평균 제곱편차는 \(2\sum_{i\ne j}b_{ij}^2\)입니다. 대칭 \(W\), \(a=|\rho|\|W\|_2<1\)이면 10절의 로그급수와 절단 상한이 성립합니다.

증명. \(z^TBz=\operatorname{tr}B+2\sum_{i<j}b_{ij}z_iz_j\)입니다. 각 \(i\ne j\)에서 \(z_i\)의 부호만 바꾸는 짝짓기로 \(z_iz_j\)의 평균은 0입니다. 제곱편차의 서로 다른 두 쌍에 대한 곱은 어떤 첨자가 홀수 번 나타나므로 같은 짝짓기로 평균 0입니다. 같은 쌍의 곱은 1입니다. 따라서 제곱편차 평균은 \(4\sum_{i<j}b_{ij}^2=2\sum_{i\ne j}b_{ij}^2\)입니다. 비대칭 \(B\)에는 \(z^TBz=z^T(B+B^T)z/2\)를 사용해 대칭 부분의 공식을 적용합니다.

\(W=Q\operatorname{diag}(\lambda_i)Q^T\)에서 \(|\rho\lambda_i|<1\)이므로 각 스칼라 로그급수 \(\log(1-\rho\lambda_i)=-\sum_{j\ge1}(\rho\lambda_i)^j/j\)가 절대수렴합니다. 유한 고윳값 합과 급수를 교환하면 자취식입니다. 꼬리는 \(\sum_{j>K}\sum_i|\rho\lambda_i|^j/j\le n\sum_{j>K}a^j/(K+1)=na^{K+1}/((K+1)(1-a))\)입니다. ∎

Krylov 방법은 반복으로 만든 부분공간 안에서 오차나 잔차를 줄입니다. 다음 장에서는 무작위 부분공간과 스케치를 사용해 계산량을 더 줄이고, 절약한 계산에 어떤 보증이 남는지 확인합니다. N9로 이어 읽기.