N4 · 최소제곱의 수치계산과 정칙화#

1. 같은 관측을 설명하는 두 계수를 따로 구할 수 있을까#

두 생산활동의 효과를 세 번의 관측으로 구별한다고 합시다. 활동량과 관측값은 각각 정한 기준단위로 나누어 무차원으로 기록합니다. 첫 관측은 두 활동의 합을 재고, 나머지 두 관측은 각 활동에 약하게 반응합니다. 선형 측정모형은

\[\begin{split} A_\varepsilon=\begin{pmatrix}1&1\\\varepsilon&0\\0&\varepsilon\end{pmatrix}, \qquad b=A_\varepsilon x_*+e,\qquad0<\varepsilon\le1 \end{split}\]

입니다. \(x_*\)는 실제 효과, \(e\)는 측정오차입니다. 우선 \(e=0\), \(x_*=(1,1)^T\)로 두어 계산오차만 분리합시다. 그러면 \(b=(2,\varepsilon,\varepsilon)^T\)입니다. 두 번째와 세 번째 식은 작지만 계수를 구별하는 유일한 관측입니다.

최소제곱은 \(\|A_\varepsilon x-b\|_2^2\)를 가장 작게 만드는 \(x\)를 찾습니다. 이번에는 정확히 맞는 해가 있으므로 최소값이 0이고 \(x=x_*\)입니다. H1의 사영 정리에서 정규방정식은

\[\begin{split} \underbrace{\begin{pmatrix}1+\varepsilon^2&1\\1&1+\varepsilon^2\end{pmatrix}}_{A_\varepsilon^TA_\varepsilon}x =\begin{pmatrix}2+\varepsilon^2\\2+\varepsilon^2\end{pmatrix}. \end{split}\]

합 방향 \((1,1)\)의 고윳값은 \(2+\varepsilon^2\), 차 방향 \((1,-1)\)의 고윳값은 \(\varepsilon^2\)입니다. 따라서

\[ \sigma_1(A_\varepsilon)=\sqrt{2+\varepsilon^2},\qquad \sigma_2(A_\varepsilon)=\varepsilon,\qquad \kappa_2(A_\varepsilon)=\frac{\sqrt{2+\varepsilon^2}}{\varepsilon}. \]

작은 차이를 제곱해서 1에 더하는 계산에 주목하세요. binary64에서 \(\varepsilon=10^{-8}\)이면 \(\varepsilon^2=10^{-16}<2^{-53}\)이어서 \(1+\varepsilon^2\)가 1로 반올림됩니다. 저장한 정규행렬은 \(\begin{pmatrix}1&1\\1&1\end{pmatrix}\)가 됩니다. 원래 두 열은 독립인데 이 경로에서는 그 차이를 잃었습니다. 이후의 Cholesky가 이 행렬을 특이하다고 판정하는 것은 이미 잃은 정보를 되찾지 못한 결과입니다.

2. 세 계산 경로가 사용하는 정보#

정규방정식은 \(G=A^TA\), \(c=A^Tb\)를 만든 뒤 \(Gx=c\)를 풉니다. QR은 \(A=Q\begin{pmatrix}R\\0\end{pmatrix}\), \(Q^Tb=(c,d)^T\)로 바꾸어 \(Rx=c\)를 풉니다. SVD는 \(A=U\Sigma V^T\)에서 \(u_i^Tb\)를 각 \(\sigma_i\)로 나누고 \(v_i\) 방향을 합합니다. 정확한 산술에서 같은 최소제곱을 풀지만, 중간에 저장하는 양은 다릅니다.

경로

실수 밀집행렬의 대표 비용

완전 열계수일 때 계산의 특징

계수결핍

정규방정식·Cholesky

\(mn^2+n^3/3\)

\(A^TA\)의 형성오차와 조건수 제곱을 고려

별도 처리가 필요

Householder QR

\(2mn^2-2n^3/3\)

원래 최소제곱 자료의 작은 후진오차로 설명

피벗·완전 직교분해 등이 필요

SVD

\(O(mn^2+n^3)\)

특이값별 증폭과 절단을 직접 확인

최소노름해와 수치적 계수 처리

앞의 비용은 \(m\ge n\)에서 분해의 주항이며, SVD의 상수는 벡터를 어디까지 구하는지와 알고리즘에 따라 달라집니다. 우변 적용, 메모리 이동과 병렬화까지 포함한 실제 실행시간 표는 아닙니다.

import numpy as np
import scipy.linalg as la
def ls_qr(A, b):
    Q, R = la.qr(A, mode="economic")
    return la.solve_triangular(R, Q.T@b)

def ls_svd(A, b, rtol=0.):
    U, s, Vh = la.svd(A, full_matrices=False)
    keep = s > rtol*s[0]
    return Vh[keep].T @ ((U[:, keep].T@b)/s[keep])

eps = 1e-8
A = np.array([[1., 1.], [eps, 0.], [0., eps]])
b = np.array([2., eps, eps])
G = A.T@A
assert np.array_equal(G, np.ones((2, 2)))
for solve in (ls_qr, ls_svd):
    assert np.linalg.norm(solve(A, b)-[1., 1.]) < 1e-7
print("정규행렬은 계수를 잃지만 QR·SVD는 원래 두 열을 사용합니다.")
정규행렬은 계수를 잃지만 QR·SVD는 원래 두 열을 사용합니다.
같은 세 관측 최소제곱에서 조건수가 커질 때 정규방정식과 QR 및 SVD의 전진오차를 비교하고 Cholesky 실패를 따로 표시한 그래프

그림 97 정확한 해를 \((1,1)\)로 정한 동일 입력의 비교이다. 양축은 로그 눈금이다. 실패는 유한 오차값과 구별해 위쪽 십자로 표시하며, 오차가 0인 결과는 표시 하한에 놓았다. 최악오차 상한이 모든 실험의 기울기를 결정하지는 않는다.#

일반적인 완전 열계수 \(A\)에서도 \(\kappa_2(A^TA)=\kappa_2(A)^2\)입니다. 다만 이것이 “항상 맞는 자릿수가 절반”이라는 뜻은 아닙니다. \(u\kappa\)\(u\kappa^2\)는 작은 오차 조건 아래의 규모를 설명하는 상한에 등장합니다. 특정 입력에서는 훨씬 작은 오차도 가능하고, \(u\kappa^2\gtrsim1\)은 보장이 무너지는 영역이지 모든 입력의 실패를 판정하는 필요충분조건은 아닙니다. 예컨대 대각행렬은 매우 다른 크기의 대각도 교차곱에서 각각 저장할 수 있습니다.

3. QR을 사용해도 큰 잔차가 문제를 민감하게 만든다#

이제 알고리즘을 바꾸어도 피할 수 없는 민감도를 보겠습니다. \(A=\begin{pmatrix}1&0\\0&\varepsilon\\0&0\end{pmatrix}\), \(b=(1,0,1)^T\)이면 최소제곱해는 \(x=(1,0)^T\), 잔차는 \(r=b-Ax=(0,0,1)^T\)입니다. 셋째 관측은 어떤 계수로도 설명할 수 없습니다.

행렬의 셋째 행 둘째 성분만 \(\delta\)만큼 바꾼 \(A+E\)로 다시 풀면 목적함수는

\[ (x_1-1)^2+\varepsilon^2x_2^2+(\delta x_2-1)^2. \]

첫 성분은 \(x_1=1\)이고, 둘째 성분의 미분을 0으로 두면 \(2\varepsilon^2x_2+2\delta(\delta x_2-1)=0\)입니다. 따라서

\[ \widetilde x_2=\frac{\delta}{\varepsilon^2+\delta^2}. \]

\(|\delta|\ll\varepsilon\)에서는 약 \(\delta/\varepsilon^2\)입니다. 정규방정식을 컴퓨터에 만들지 않아도 이 정확한 해 자체\(1/\varepsilon^2\) 규모로 반응합니다. 설명하지 못한 관측 \(r\)을 작은 열 방향으로 설명하려는 변화가 들어왔기 때문입니다.

일반식도 같은 두 경로를 보여 줍니다. \(A\)가 완전 열계수이고 \(r=b-Ax\), \(A^Tr=0\)일 때 작은 자료 변화 \((E,f)\)에 대한 1차 해 변화는

\[ dx=A^+(f-Ex)+(A^TA)^{-1}E^Tr. \]

첫 항은 우변과 적합값의 직접 변화이고 둘째는 기존 잔차가 바뀐 열공간에 투영되는 효과입니다. \(r=0\)일 때에만 둘째 항이 사라집니다. 이를 단순히 “QR은 언제나 \(u\kappa\)의 전진오차”라고 요약하면 이 효과를 놓칩니다. 마지막 절에서 이 식과 유한한 섭동에 대한 상한을 증명합니다.

4. 해가 여러 개이면 무엇을 돌려줄 것인가#

\(A=\begin{pmatrix}1&1\\0&0\end{pmatrix}\), \(b=(2,1)^T\)이면 가능한 적합값은 \((x_1+x_2,0)^T\)입니다. 최소제곱은 \(x_1+x_2=2\)만 정하고 잔차 \((0,1)\)을 남깁니다. 해를 \((1+t,1-t)\)로 쓰면 제곱노름은 \(2+2t^2\)이므로 최소노름해는 \((1,1)\)입니다. 하나의 최적 적합값, 여러 계수, 유일한 최소노름 계수를 구별하세요.

일반적으로 \(A^+b=\sum_{\sigma_i>0}(u_i^Tb/\sigma_i)v_i\)가 그 최소노름해입니다. 너무 작은 \(\sigma_i\)를 0으로 취급하는 허용오차 rtol은 수치적 계수를 선택합니다. 이는 수학적으로 비영인 모든 특이값을 그대로 역수화하는 것과 다른 문제입니다. 원래 행렬의 정확한 계수를 모르는데 소프트웨어가 그것을 자동으로 알아냈다고 해석할 수 없습니다.

5. 작은 특이값을 그대로 나누지 않기#

세 독립된 측정 방향의 감도가 \(D=\operatorname{diag}(1,0.1,0.01)\)이라고 합시다. 실제 효과를 \(x_*=(1,1,1)\), 관측을 \(b=(1,0.1,0.02)\)로 두면 마지막 관측에 \(0.01\)의 오차가 있습니다. 그대로 나눈 해는 \((1,1,2)\)입니다. 마지막 관측의 작은 절대오차가 계수에서는 1이 되었습니다.

능형 또는 Tikhonov 정칙화는 \(\lambda>0\)에 대해

\[ \min_x\{\|Ax-b\|_2^2+\lambda\|x\|_2^2\} \]

를 풉니다. \(\lambda\)는 관측 제곱오차와 계수 제곱의 상대 가중치입니다. 이 예에서는 모든 좌표를 기준단위로 정규화했으므로 무차원 가중치로 쓸 수 있습니다. 서로 다른 단위의 계수에 같은 벌점을 주려면 어떤 기준척도를 사용하는지 먼저 정해야 합니다. 절편을 벌점에서 제외하는 회귀는 \(I\) 대신 다른 벌점행렬을 쓰는 별도 모형입니다.

SVD 좌표에서는 각 방향을 따로 최소화하여

\[ x_\lambda=\sum_{\sigma_i>0}\frac{\sigma_i}{\sigma_i^2+\lambda}(u_i^Tb)v_i =\sum_{\sigma_i>0}\underbrace{\frac{\sigma_i^2}{\sigma_i^2+\lambda}}_{f_i(\lambda)} \frac{u_i^Tb}{\sigma_i}v_i. \]

\(f_i\)는 정칙화하지 않은 계수에 곱하는 필터인자입니다. \(\sigma_i^2\gg\lambda\)이면 1에 가깝고, \(\sigma_i^2\ll\lambda\)이면 0에 가깝습니다. 위 \(D\)\(\lambda=10^{-4}\)를 넣으면

\[ x_\lambda=\left(\frac{10000}{10001},\frac{100}{101},1\right)^T. \]

마지막 방향의 잡음 증폭은 줄었지만 첫 두 방향에도 작은 편향을 넣었습니다. 실제 \(x_*\)를 알 수 없는 문제에서 이 \(\lambda\)가 가장 좋다는 결론은 여기서 나오지 않습니다.

일반 벌점 \(\lambda\|Lx\|^2\)의 경우 해가 유일할 필요충분조건은 \(\ker A\cap\ker L=\{0\}\)입니다. 목적함수의 이차부가 \(x^T(A^TA+\lambda L^TL)x=\|Ax\|^2+\lambda\|Lx\|^2\)이므로 공통 영방향이 없을 때 정확히 양의 정부호가 됩니다. N7의 차분 평활화가 이 형태입니다.

6. 절단과 조기종료도 필터가 된다#

TSVD는 선택한 문턱 \(\tau\) 아래의 특이값을 버립니다. 이때 \(f_i=1\) 또는 0입니다. Landweber 반복은 \(x_0=0\)에서

\[ x_{k+1}=x_k+\omega A^T(b-Ax_k) \]

로 갱신합니다. \(0<\omega<2/\sigma_1^2\)이면 비영 특이방향에서 수렴하고, \(k\)번 뒤의 필터는

\[ f_i^{(k)}=1-(1-\omega\sigma_i^2)^k. \]

잡음이 있는 작은 특이방향을 완전히 역수화하기 전에 멈추면 정칙화 효과가 있습니다. 필터가 0과 1 사이에서 단조 증가한다고 말하려면 더 강한 \(0<\omega\le1/\sigma_1^2\)를 사용해야 합니다. 그보다 큰 안정적인 보폭에서는 일부 필터가 진동할 수 있습니다.

특잇값별로 능형의 매끄러운 감쇠, TSVD의 문턱, Landweber 조기종료의 필터를 비교하는 그래프

그림 98 가로축은 로그 특잇값, 세로축은 정칙화하지 않은 계수에 곱하는 값이다. Landweber에는 \(\omega=1\), \(0<\sigma\le1\)을 사용해 필터가 \([0,1]\)에 머무른다.#

적합벡터를 \(\widehat b=H_\lambda b\)라 쓰면 \(H_\lambda=U\operatorname{diag}(f_i)U^T\)이고 \(\operatorname{tr}H_\lambda=\sum_i f_i\)입니다. 이를 선형 평활의 유효자유도라고 부릅니다. 여기서는 자료에 대한 선형 반응의 자취라는 대수적 정의이며, 통계적 자유도와의 관계는 잡음모형을 도입한 뒤 다룹니다. \(\lambda\downarrow0\)에서 계수 \(r\)로, \(\lambda\to\infty\)에서 0으로 갑니다. 열 수 \(n\)으로 가는 것은 완전 열계수일 때입니다.

7. 정칙화 모수는 무엇을 보고 선택할까#

한 번 계산한 SVD로 여러 \(\lambda\)의 계수·적합·잔차를 구할 수 있습니다. 또

\[\begin{split} \left\|\begin{pmatrix}A\\\sqrt\lambda I\end{pmatrix}x- \begin{pmatrix}b\\0\end{pmatrix}\right\|^2 =\|Ax-b\|^2+\lambda\|x\|^2 \end{split}\]

이므로 확대행렬의 QR로 풀어도 됩니다. 확대행렬의 특이값은 \(\sqrt{\sigma_i^2+\lambda}\)이며, 영 특이값까지 \(n\)개로 채우면 조건수는 정확히 \(\sqrt{(\sigma_1^2+\lambda)/(\sigma_n^2+\lambda)}\)입니다.

L-곡선은 \((\|Ax_\lambda-b\|,\|x_\lambda\|)\)를 로그축에 그려 맞춤과 계수 크기의 교환관계를 보여 줍니다. 곡선이 반드시 뚜렷한 L 모양이거나 유일한 모서리를 갖는 것은 아닙니다. GCV는

\[ \operatorname{GCV}(\lambda)= \frac{\|b-H_\lambda b\|^2/m}{(1-\operatorname{tr}H_\lambda/m)^2} \]

로 잔차를 보정합니다. 분모가 0인 경우는 제외합니다. 실제로 한 관측씩 뺀 능형 적합의 예측오차는 \(r_i/(1-h_{ii})\)이며, GCV는 이 개별 분모를 평균 \(1-\operatorname{tr}H/m\)으로 바꾼 것입니다. \(h_{ii}\)가 매우 다르면 이 근사가 좋지 않을 수 있습니다. 오차의 분포와 독립성에 관한 통계적 정당화는 아직 가정하지 않았습니다.

Picard 그림은 \(|u_i^Tb|\), \(\sigma_i\), \(|u_i^Tb|/\sigma_i\)를 나란히 봅니다. 관측 계수가 어느 지점부터 감소하지 않는다면 작은 특이값으로 나누는 단계가 그 성분을 키울 수 있습니다. 유한차원에서 이 그림이 특정 감소율을 보이지 않는다는 이유만으로 “해가 의미 없다”고 단정하지는 않습니다. 유용한 문턱과 모수는 관측오차 크기, 요구 정확도와 계수의 의미를 함께 사용해 정해야 합니다.

같은 정칙화 문제에서 잔차노름과 계수노름의 L 곡선 및 정칙화 계수별 GCV를 보여 주는 두 패널

그림 99 \(A=[D;0]\), \(b=(1,0.1,0.02,0.003)\)의 결정론적 예이다. 두 패널 모두 같은 \(\lambda\) 격자를 사용한다. L-곡선은 두 노름의 로그축, 오른쪽은 양축 모두 로그축이며 GCV 최솟값을 관측모형의 보편적 최적값이라고 해석하지 않는다.#

8. 관측 하나를 추가하거나 빼면#

기존 \(G=A^TA\)가 양의 정부호이고 새 관측행이 \(a^T\), 새 값이 \(y\)라면 \(G\)\(G+aa^T\), \(A^Tb\)\(A^Tb+ay\)로 바뀝니다. \(x=G^{-1}A^Tb\)에서

\[ x_{\rm new}=x+\frac{G^{-1}a}{1+a^TG^{-1}a}(y-a^Tx) \]

입니다. 새 관측의 예측오차를 방향 \(G^{-1}a\)로 보정합니다. 역행렬을 저장해야 한다는 뜻은 아니며, N3의 Givens로 \(R\)에 행을 추가하는 제곱근 방식이 같은 문제를 풉니다.

관측을 빼는 경우에는 \(G-aa^T\)이고 양의 정부호를 유지할 조건은 \(a^TG^{-1}a<1\)입니다. 분모 \(1-a^TG^{-1}a\)가 작아지면 다운데이팅은 민감해집니다. 이 조건은 \(G^{-1/2}(G-aa^T)G^{-1/2}=I-vv^T\), \(v=G^{-1/2}a\)의 한 고윳값 \(1-\|v\|^2\)에서 나옵니다. 예컨대 \(G=I\), \(a=e_1\)을 제거하면 첫 방향을 결정할 정보가 없어져 특이해집니다.

9. 연습과 전체 풀이#

1. \(A_\varepsilon\)\(b=(2,\varepsilon,\varepsilon+\eta)\)를 넣고 정확한 최소제곱해를 구하세요.

풀이. 원래 해 \((1,1)\)의 변화 \(d\)에 대해 정규방정식 우변 변화는 \((0,\varepsilon\eta)\)입니다. 두 식을 더하면 \((2+\varepsilon^2)(d_1+d_2)=\varepsilon\eta\), 빼면 \(\varepsilon^2(d_2-d_1)=\varepsilon\eta\)입니다. 따라서 \(d_1=\frac12(\varepsilon\eta/(2+\varepsilon^2)-\eta/\varepsilon)\), \(d_2=\frac12(\varepsilon\eta/(2+\varepsilon^2)+\eta/\varepsilon)\)입니다. 차 방향에 \(1/\varepsilon\)의 증폭이 있습니다.

2. \(\kappa_2(A)=10^6\)일 때 \(u\kappa\)\(u\kappa^2\)의 크기 및 해석을 쓰세요.

풀이. \(u\approx1.11\times10^{-16}\)이므로 각각 약 \(1.11\times10^{-10}\), \(1.11\times10^{-4}\)입니다. 상수와 잔차 효과를 제외한 규모로는 10자리와 4자리 정도지만 보장된 실제 자릿수의 등식은 아닙니다. 실제 상대오차를 알아야 \(-\log_{10}(\text{상대오차})\)를 보고할 수 있습니다.

3. \(A=[1\ 1]\), \(b=2\)의 능형해와 \(\lambda\downarrow0\) 극한을 구하세요.

풀이. 목적함수는 \((x_1+x_2-2)^2+\lambda(x_1^2+x_2^2)\)입니다. 두 정상조건을 빼면 \(\lambda(x_1-x_2)=0\), 따라서 \(x_1=x_2=t\)입니다. \((2+\lambda)t=2\)여서 \(t=2/(2+\lambda)\to1\)입니다. 극한은 여러 정확해 중 최소노름해입니다.

4. \(D=\operatorname{diag}(1,0.01)\), \(b=D(1,1)^T\), \(\lambda=0.01\)에서 둘째 열을 100배 스케일하고 같은 계수 벌점을 주면 원래 좌표의 답도 같은가요?

풀이. 원래 둘째 계수는 \(0.0001/(0.0001+0.01)=1/101\)입니다. \(x=S z\), \(S=\operatorname{diag}(1,100)\)로 바꾸면 \(DS=I\)이고 같은 \(\lambda\|z\|^2\) 벌점의 둘째 \(z\)\(0.01/1.01\)입니다. 되돌린 \(x_2=100z_2=100/101\)로 다릅니다. 같은 문제를 유지하려면 새 벌점은 \(\lambda\|Sz\|^2\)이어야 합니다.

5. \(\omega\sigma^2=3/2\)인 Landweber 방향의 처음 세 필터를 구하세요.

풀이. \(1-\omega\sigma^2=-1/2\)이므로 \(f_1=3/2\), \(f_2=3/4\), \(f_3=9/8\)입니다. 1로 수렴하지만 1의 위아래를 번갈아 지납니다. 안정적인 보폭과 단조 필터 조건이 다릅니다.

6. \(G=\operatorname{diag}(2,1)\), 기존 해 \(x=(1,1)\)\(a=(1,0)\), \(y=2\)를 추가하세요.

풀이. \(G^{-1}a=(1/2,0)\), \(a^TG^{-1}a=1/2\), 새 관측오차는 1입니다. 변화는 \((1/3,0)\)이므로 새 해는 \((4/3,1)\)입니다. 직접 \(G+aa^T=\operatorname{diag}(3,1)\), 새 우변 \((4,1)\)을 풀어도 같습니다.

7. \(A=\operatorname{diag}(1,0.1)\)의 능형 유효자유도를 쓰고 단조성을 확인하세요.

풀이. \(\mathrm{df}(\lambda)=1/(1+\lambda)+0.01/(0.01+\lambda)\)입니다. 도함수는 \(-1/(1+\lambda)^2-0.01/(0.01+\lambda)^2<0\)입니다. \(\lambda\downarrow0\)에서 2, 무한대로 가면 0입니다.

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

최소제곱과 정칙화의 해는 H4의 SVD로 좌표별로 설명할 수 있습니다. 계산의 안정성에는 N1의 반올림·후진오차N3의 Householder 누적오차를 적용합니다. 아래에서는 실수행렬을 쓰며, 복소수에서는 전치를 수반으로 바꾸면 됩니다.

정리 1. 최소노름해와 조건수 제곱#

계수 \(r\)\(A\)의 모든 최소제곱해는 \(A^+b+\ker A\)입니다. 그 중 최소노름해는 유일하게 \(A^+b\)입니다. 완전 열계수이면 \(\kappa_2(A^TA)=\kappa_2(A)^2\)입니다.

증명. 완전한 직교기저의 SVD 좌표 \(z=V^Tx\), \(c=U^Tb\)를 씁니다. 목적함수는 \(\sum_{i=1}^r(\sigma_i z_i-c_i)^2+\sum_{i>r}c_i^2\)입니다. 첫 합은 각 \(z_i=c_i/\sigma_i\)에서 유일하게 0이고, 영 특이값 방향의 \(z_i\)는 목적함수에 나타나지 않습니다.

이 자유 방향들의 집합이 \(\ker A\)입니다. \(\|x\|^2=\sum_i z_i^2\)이므로 자유 성분을 모두 0으로 해야 최소노름이며 유일합니다. 완전 열계수에서는 \(A^TA=V\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2)V^T\)가 양의 정부호이므로 최대·최소 특이값의 비가 \(\sigma_1^2/\sigma_n^2\)입니다. ∎

정리 2. 잔차를 포함하는 최소제곱 섭동#

\(A\)가 완전 열계수이고 \(s=\sigma_n(A)>0\), \(\|E\|_2=e<s\)라 합시다. \(x=A^+b\), \(\widetilde x=(A+E)^+(b+f)\), \(r=b-Ax\)이면

\[ \|\widetilde x-x\|_2\le \frac{\|f\|_2+e\|x\|_2}{s-e} +\frac{e\|r\|_2}{(s-e)^2}. \]

그 1차 도함수는 3절의 식입니다.

증명. 단위벡터 \(v\)에 대해 \(\|(A+E)v\|\ge\|Av\|-\|Ev\|\ge s-e>0\)이므로 \(\widetilde A=A+E\)도 완전 열계수입니다. \(\widetilde A^T\widetilde A\widetilde x=\widetilde A^T(b+f)\)에서 \(b=Ax+r=\widetilde Ax-Ex+r\)를 대입하면

\[ \widetilde A^T\widetilde A(\widetilde x-x) =\widetilde A^T(f-Ex)+E^Tr. \]

마지막 항은 \(A^Tr=0\)을 사용했습니다. 역행렬을 적용하면 첫 항의 연산자는 \(\widetilde A^+\)이고 노름은 \(1/\sigma_n(\widetilde A)\le1/(s-e)\)입니다. 둘째의 연산자 노름은 그 제곱입니다. 삼각부등식으로 상한을 얻습니다. \((E,f)\)\(t(E,f)\)로 바꾸어 \(t\)로 나누고 \(t\to0\)으로 보내면 \(dx=A^+(f-Ex)+(A^TA)^{-1}E^Tr\)입니다. 역행렬의 연속성은 \(e<s\) 영역에서 N1의 Neumann 논법으로 보장됩니다. ∎

\(x\ne0\), \(\cos\theta=\|Ax\|/\|b\|\), \(\nu=\|A\|\|x\|/\|Ax\|\)라 두면 \(1\le\nu\le\kappa\)입니다. \(\|E\|\le\epsilon_A\|A\|\), \(\|f\|\le\epsilon_b\|b\|\)인 1차 상한은

\[ \frac{\|dx\|}{\|x\|}\le \frac{\kappa}{\nu\cos\theta}\epsilon_b +\left(\kappa+\frac{\kappa^2\tan\theta}{\nu}\right)\epsilon_A. \]

이는 바로 앞 식의 각 항을 나눈 것입니다. \(\|r\|/\|Ax\|=\tan\theta\)는 적합값과 잔차의 직교성에서 나옵니다. \(x=0\)에서는 상대오차 대신 위 정리의 절대오차식을 씁니다.

정리 3. QR 최소제곱 계산의 후진 해석#

정상 범위에서 Householder를 \([A\ b]\)에 적용했다고 하자. N3의 국소오차 분석으로 정확한 직교 \(Q_*\)와 작은 \(E,f\)가 있어 \(Q_*^T(A+E)=\begin{pmatrix}\widehat R\\0\end{pmatrix}\), \(Q_*^T(b+f)=\begin{pmatrix}\widehat c\\\widehat d\end{pmatrix}\)라고 쓸 수 있습니다. 이어 정상 범위의 삼각대입이 완료되면 계산해는 작은 추가 섭동을 받은 원래 최소제곱 문제의 정확한 해입니다.

증명. N2의 삼각대입 후진오차에서 \((\widehat R+F)\widehat x=\widehat c\), \(|F|\le\gamma_{n+1}|\widehat R|\)입니다. \(\widetilde A=A+E+Q_*\begin{pmatrix}F\\0\end{pmatrix}\), \(\widetilde b=b+f\)로 두면

\[\begin{split} Q_*^T(\widetilde b-\widetilde A\widehat x)=\begin{pmatrix}0\\\widehat d\end{pmatrix},\qquad Q_*^T\widetilde A=\begin{pmatrix}\widehat R+F\\0\end{pmatrix}. \end{split}\]

따라서 잔차가 열공간에 직교하고 \(\widehat x\)는 정확한 최소제곱해입니다. 추가 행렬변화의 Frobenius 노름은 직교 불변성으로 \(\|F\|_F\le\gamma_{n+1}\|\widehat R\|_F\)입니다. 원래 \(A\)의 변환에는 N3의 단계별 상한을 \(A\) 열블록에, 우변에는 \(b\)에 적용할 수 있어 각각 차원 상수와 \(u\)에 비례하는 상대 후진오차를 얻습니다. 이 결론이 전진오차까지 작게 만들려면 정리 2의 민감도도 작아야 합니다. ∎

정리 4. 필터, 확대행렬, 극한#

\(\lambda>0\)인 능형해는 유일하며 5절의 SVD 식을 만족합니다. \(\lambda\downarrow0\)에서 \(A^+b\)로 수렴합니다. 확대행렬 조건수와 유효자유도는 6·7절의 식입니다. Landweber 필터는 \(0<\omega<2/\sigma_1^2\)에서 비영 특이방향을 모두 복원합니다.

증명. \(z=V^Tx\) 좌표에서 목적함수는 \(\sum_{i\le r}((\sigma_i z_i-c_i)^2+\lambda z_i^2)+\sum_{i>r}\lambda z_i^2+\sum_{i>r}c_i^2\)입니다. 각 이차식의 양의 이차계수는 \(\sigma_i^2+\lambda\) 또는 \(\lambda\)이며, 최소점은 \(z_i=\sigma_i c_i/(\sigma_i^2+\lambda)\) 또는 0입니다. 유일성과 식을 얻습니다. 유한 개 항에서 \(\lambda\downarrow0\) 극한을 취하면 비영 방향은 \(c_i/\sigma_i\), 영방향은 0으로 남습니다.

확대행렬의 Gram 행렬은 \(A^TA+\lambda I\)이므로 특이값과 조건수 식이 따릅니다. \(Ax_\lambda=\sum f_i c_i u_i\)에서 \(H_\lambda\)를 얻고 직교기저에서 자취를 계산하면 \(\sum f_i\)입니다.

Landweber의 각 비영 좌표는 \(z_{i,k+1}=(1-\omega\sigma_i^2)z_{i,k}+\omega\sigma_i c_i\)입니다. \(z_{i,0}=0\)에서 유한 등비합을 계산하면 \(z_{i,k}=[1-(1-\omega\sigma_i^2)^k]c_i/\sigma_i\)입니다. 영 특이값 방향은 0에서 변하지 않습니다. 모든 비영 방향에서 \(|1-\omega\sigma_i^2|<1\)이 정확히 위 보폭 조건으로 보장됩니다. ∎

정리 5. rank-one 갱신과 한 관측 제거#

가역 \(G\)에서 \(1+v^TG^{-1}u\ne0\)이면

\[ (G+uv^T)^{-1}=G^{-1}-\frac{G^{-1}uv^TG^{-1}}{1+v^TG^{-1}u}. \]

증명. 오른쪽에 \(G+uv^T\)를 곱합니다. \(I\) 외에 남는 항은 \(G^{-1}uv^T\)에 곱해지는 계수 \(1-(1+v^TG^{-1}u)/(1+v^TG^{-1}u)=0\)이므로 양쪽 역행렬입니다. \(u=v=a\)와 새 우변을 대입하고 기존 \(Gx=A^Tb\)를 쓰면 8절의 계수 갱신식을 얻습니다.

능형의 \(G=A^TA+\lambda I\succ0\)에서 \(i\)번째 행 \(a_i^T\)를 제거해도 \(G-a_ia_i^T\succ0\)입니다. 나머지 행의 제곱합에 \(\lambda I\)가 더해져 있기 때문입니다. 따라서 \(h_{ii}=a_i^TG^{-1}a_i<1\)입니다. 제거된 해 \(x^{(-i)}\)에 앞 역행렬식을 마이너스 부호로 적용하면

\[ x^{(-i)}=x-\frac{G^{-1}a_i}{1-h_{ii}}(b_i-a_i^Tx). \]

왼쪽에 \(a_i^T\)를 곱하여 예측값을 빼면 \(b_i-a_i^Tx^{(-i)}=r_i/(1-h_{ii})\)입니다. 이로써 7절에서 GCV와 비교한 정확한 관측 제거 오차식이 증명됩니다. ∎

여러 행을 갱신할 때의 Woodbury 공식도 같은 곱 검산으로 \(G^{-1}-G^{-1}U(C^{-1}+VG^{-1}U)^{-1}VG^{-1}\)을 얻습니다. 여기에는 \(G,C,C^{-1}+VG^{-1}U\)의 가역성이 필요합니다. 이 식은 작은 보조계가 민감하거나 차감이 큰 경우의 정확도를 자동으로 보장하지 않으며, QR·Cholesky 인자 갱신의 잔차와 함께 확인해야 합니다.

최소제곱의 수치해와 정칙화로 바꾼 문제의 해를 구별했습니다. 다음 장에서는 고윳값을 계산하는 반복을 살펴보며, 후보값의 잔차로 얻을 수 있는 정보를 확인합니다. N5로 이어 읽기.