보강 2. 수치선형대수 입문 - 선형대수

이 책은 같은 경고를 다섯 번 했다.

L19에서 여인수 전개를 유도해 놓고 쓰지 말라고 했다. L20에서 크래머 공식을 유도해 놓고 또 쓰지 말라고 했다. L28에서 조르당 표준형은 수치적으로 존재하지 않는다고 했다. L33에서 정규방정식을 손으로 쓰지 말라고 했다. 보강 1에서 공분산행렬을 만들지 말라고 했다.

다섯 번이면 우연이 아니다. 그런데 다섯이 정말 같은 이야기인가.

이번 글의 목표는 그것이다. 다섯을 한자리에 놓고 무엇이 같고 무엇이 다른지 가른 뒤, 그 뒤에 무엇이 있는지 본다. 새로 배우는 것은 거의 없다. 이미 가진 조각들을 다른 순서로 놓을 뿐이다.

1. 다섯 번의 경고는 두 갈래였다¶

한자리에 놓고 보자.

경고	무엇이 문제였나
L19 여인수 전개	항이 $n!$ 개다
L20 크래머 공식	행렬식을 $n+1$ 번 구해야 한다
L28 조르당 형	답이 $\sqrt{\epsilon}$ 만큼 튄다
L33 정규방정식	조건수가 제곱된다
보강 1 공분산행렬	조건수가 제곱된다

위의 둘과 아래의 셋이 다른 이야기다.

위의 둘은 비용이다. 답은 정확한데 살아 생전에 못 받는다. 셋을 나란히 적어 보자.

\text{소거} : \frac{n^3}{3}, \qquad \text{여인수 전개} : n!, \qquad \text{크래머} : (n+1)\,n!

(1)

(1)의 세 값은 $n$ 이 작을 때는 비슷하다. $n = 5$ 에서 41번, 120번, 720번이라 스무 배 차이가 날 뿐이다. 그런데 $n = 25$ 가 되면 5,208번과 $4 \times 10^{26}$ 번이 된다. 1초에 10억 번씩 계산해도 10¹⁰ 년, 곧 우주의 나이와 맞먹는다.

아래의 셋은 정확도다. 계산은 순식간에 끝나는데 나온 숫자가 틀렸다. 그리고 이쪽이 훨씬 위험하다. 오래 걸리는 것은 눈에 보이지만 틀린 것은 안 보인다.

왼쪽이 비용으로 무너지는 길, 오른쪽이 정확도로 무너지는 길이다.
왼쪽은 축이 n 이고 오른쪽은 축이 조건수다. 다른 것을 재고 있다. — Figure 1:왼쪽이 비용으로 무너지는 길, 오른쪽이 정확도로 무너지는 길이다. 왼쪽은 축이 $n$ 이고 오른쪽은 축이 조건수다. **다른 것을 재고 있다.**

2. 전방오차와 후방오차¶

정확도 이야기를 하려면 말부터 정해야 한다. $A\vv{x} = \vv{b}$ 를 풀어 $\hat{\vv{x}}$ 를 얻었다고 하자. 얼마나 틀렸는가.

가장 자연스러운 답은 참값과 비교하는 것이다.

\text{전방오차} = \frac{\lVert \hat{\vv{x}} - \vv{x} \rVert}{\lVert \vv{x} \rVert}

(2)

그런데 (2)의 값은 잴 수가 없다. 참값 $\vv{x}$ 를 모르니까 푸는 것이다.

그래서 질문을 뒤집는다¶

수치해석의 요령이 여기 있다. "내 답이 얼마나 틀렸는가"를 묻는 대신 내 답이 어떤 문제의 정확한 답인가를 묻는다.

(A + \Delta A)\,\hat{\vv{x}} = \vv{b}

(3)

(3)의 식을 만족하는 $\Delta A$ 중 가장 작은 것을 찾아 크기를 잰다.

\text{후방오차} = \frac{\lVert \Delta A \rVert}{\lVert A \rVert}

(4)

(4)의 값은 잴 수 있다. 그리고 뜻이 분명하다. “내가 푼 것은 원래 문제가 아니라 이만큼 옆에 있는 문제였다.”

정말 잴 수 있는지 보이자. 잔차를 $\vv{r} = \vv{b} - A\hat{\vv{x}}$ 라 두고 $\Delta A$ 를 이렇게 잡아 보면 된다.

\Delta A = \frac{\vv{r}\,\hat{\vv{x}}^{\mathsf T}}{\lVert\hat{\vv{x}}\rVert^2}

(5)

(5)의 행렬을 (3)에 넣어 보자. 랭크 1 행렬을 $\hat{\vv{x}}$ 에 곱하는 것이라 계산이 한 줄로 끝난다.

(A + \Delta A)\hat{\vv{x}} = A\hat{\vv{x}} + \vv{r}\,\frac{\hat{\vv{x}}^{\mathsf T}\hat{\vv{x}}}{\lVert\hat{\vv{x}}\rVert^2} = A\hat{\vv{x}} + \vv{r} = \vv{b}

(6)

(6)의 마지막이 정확히 우리가 원한 것이다. 그리고 이 $\Delta A$ 가 가장 작은 것임도 알려져 있으므로, 후방오차를 잔차만으로 계산할 수 있다.

\frac{\lVert\Delta A\rVert}{\lVert A\rVert} = \frac{\lVert\vv{r}\rVert}{\lVert A\rVert\,\lVert\hat{\vv{x}}\rVert}

(7)

(7)의 오른쪽은 전부 손에 있는 값이다. 참해를 몰라도 된다.

둘을 잇는 부등식¶

이제 둘을 연결한다. (3)에서 $A\vv{x} = \vv{b}$ 를 빼면

A(\hat{\vv{x}} - \vv{x}) = -\Delta A\,\hat{\vv{x}} \qquad\Longrightarrow\qquad \hat{\vv{x}} - \vv{x} = -A^{-1}\Delta A\,\hat{\vv{x}}

(8)

이고, (8)의 오른쪽에 노름을 씌우면

\lVert \hat{\vv{x}} - \vv{x} \rVert \ \le\ \lVert A^{-1} \rVert\,\lVert \Delta A \rVert\,\lVert \hat{\vv{x}} \rVert

(9)

이 된다. (9)의 양변을 $\lVert\hat{\vv{x}}\rVert$ 로 나누고 $\lVert A\rVert/\lVert A\rVert$ 를 끼워 넣으면 목표한 꼴이 나온다.

\underbrace{\frac{\lVert \hat{\vv{x}} - \vv{x} \rVert}{\lVert \hat{\vv{x}} \rVert}}_{\text{전방오차}} \ \le\ \underbrace{\lVert A \rVert\,\lVert A^{-1} \rVert}_{\kappa(A)} \cdot \underbrace{\frac{\lVert \Delta A \rVert}{\lVert A \rVert}}_{\text{후방오차}}

(10)

왼쪽 화살표가 알고리즘이 책임지는 것이고 오른쪽 화살표가 우리가 원하는 것이다.
둘 사이의 환율이 조건수다. — Figure 2:왼쪽 화살표가 알고리즘이 책임지는 것이고 오른쪽 화살표가 우리가 원하는 것이다. 둘 사이의 환율이 조건수다.

다섯 경고를 다시 읽으면¶

이제 Figure 1의 세 기울기를 읽을 수 있다.

$\epsilon\kappa$ — 후방안정한 알고리즘의 한계다. 후방오차가 $\epsilon$ 이니 전방오차는 $\kappa\epsilon$ 이다. 여기까지는 어쩔 수 없다.

$\epsilon\kappa^2$ — L33의 정규방정식과 보강 1의 공분산이 여기 있다. $A^{\mathsf T}A$ 를 만드는 순간 조건수가 제곱되니, 그 뒤에 아무리 좋은 알고리즘을 써도 $\epsilon\kappa^2$ 을 넘지 못한다. 문제를 자기 손으로 나쁘게 만든 것이다.

$\sqrt{\epsilon}$ — L28의 조르당 형이 여기 있다. 이것은 조건수 이야기가 아니라 문제 자체가 매끄럽지 않다는 이야기인데, 5절에서 다시 본다.

3. 반복법 — 답을 향해 걸어간다¶

소거는 유한 번에 끝난다. $n^3/3$ 번이면 답이 나온다. 그런데 $n$ 이 백만이면 $3 \times 10^{17}$ 번이라 그것도 못 쓴다.

이럴 때 쓰는 것이 반복법이다. 정확한 답 대신 충분히 가까운 답을 노린다.

쪼개고 되풀이한다¶

$A$ 를 두 조각으로 나눈다. 하나는 풀기 쉬운 것, 하나는 나머지다.

A = M - N \qquad\Longrightarrow\qquad M\vv{x} = N\vv{x} + \vv{b}

(12)

(12)의 오른쪽 식은 여전히 $\vv{x}$ 를 양쪽에 갖고 있어 풀 수 없다. 그런데 오른쪽에는 옛 답을 넣고 왼쪽에서 새 답을 뽑으면 되풀이가 된다.

\vv{x}_{k+1} = M^{-1}\left(N\vv{x}_k + \vv{b}\right)

(13)

$M$ 을 무엇으로 잡느냐가 방법을 정한다. $A = D + L + U$ 로 대각·아래·위를 쪼개면

\text{야코비} : M = D, \qquad \text{가우스-자이델} : M = D + L

(14)

이다. 성분으로 풀어 쓰면 왜 싼지가 보인다. 야코비는 이렇다.

x_i^{(k+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j \neq i} a_{ij}\,x_j^{(k)} \right)

(15)

(15)의 오른쪽에는 옛 값만 들어간다. 그래서 $n$ 개를 동시에 계산할 수 있고 병렬화가 쉽다. 가우스-자이델은 한 글자만 다르다.

x_i^{(k+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j < i} a_{ij}\,x_j^{(k+1)} - \sum_{j > i} a_{ij}\,x_j^{(k)} \right)

(16)

(16)의 첫 합에 $k+1$ 이 들어 있다. 방금 갱신한 값을 즉시 쓴다. 그래서 순서를 지켜야 하지만 정보가 더 새것이라 빠르다. 둘 다 한 번 도는 데 $A$ 의 0이 아닌 성분 수만큼만 곱셈을 하므로, 성긴 행렬에서는 $n^3$ 근처에도 안 간다.

수렴 조건은 L22가 이미 준다¶

되풀이가 답으로 가는가. 오차 $\vv{e}_k = \vv{x}_k - \vv{x}$ 를 보면 바로 나온다. (13)에서 참해도 같은 식을 만족하므로 빼면

\vv{e}_{k+1} = M^{-1}N\,\vv{e}_k \qquad\Longrightarrow\qquad \vv{e}_k = \left(M^{-1}N\right)^{k}\vv{e}_0

(17)

가 된다. (17)의 오른쪽은 행렬의 거듭제곱이고, 그것이 0으로 가는 조건은 L22에서 이미 배웠다. 모든 고윳값의 절댓값이 1보다 작아야 한다.

L22의 $A^k \to 0$ 이 여기서 알고리즘의 수렴 판정이 된 셈이다.

앵커로 확인¶

이 책에서 여러 번 쓴 삼중대각행렬을 쓰자.

T = \operatorname{tridiag}(-1,\ 2,\ -1) \qquad \lambda_k(T) = 4\sin^2\!\frac{k\pi}{2(n+1)}

(18)

(18)의 고윳값이 닫힌 꼴로 알려져 있어 모든 것을 손으로 확인할 수 있다. $n = 20$ 에서 조건수가 178.06 이고, 야코비의 스펙트럼 반지름은

\rho_{\text{야코비}} = \cos\frac{\pi}{n+1} = 0.98883

(19)

이다. (19)의 값이 1에 아주 가깝다. 한 번 되풀이할 때마다 오차가 1.1%씩만 준다는 뜻이라 느리다. 가우스-자이델은 정확히 그 제곱인 0.97779 이라 두 배 빠르지만, 여전히 느리다.

켤레기울기 — 조건수에 제곱근이 붙는다¶

$A$ 가 대칭 양의 정부호이면(L27) 훨씬 좋은 방법이 있다.

L27에서 최급강하법을 보았다. 그때 수렴 인자가 $\dfrac{\kappa-1}{\kappa+1}$ 이었고, 조건수가 크면 1에 붙어 지그재그로 기어갔다. 켤레기울기법은 같은 문제에서

\text{최급강하} : \frac{\kappa - 1}{\kappa + 1} \qquad\text{vs}\qquad \text{켤레기울기} : \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}

(20)

의 인자를 갖는다. (20)의 차이는 제곱근 하나뿐인데 그것이 전부를 바꾼다. $\kappa = 178$ 에서 왼쪽은 0.9888 이고 오른쪽은 0.8606 이다. 자릿수 하나를 더 얻는 데 왼쪽은 205번, 오른쪽은 15번이 든다.

그리고 하나 더 있다. 켤레기울기는 $n$ 번 안에 정확히 끝난다. 매 단계마다 서로 $A$ -직교인 방향을 하나씩 잡아 나가는데, $n$ 개를 다 잡으면 공간을 다 덮기 때문이다.

왼쪽에서 켤레기울기가 k = 20 에서 절벽처럼 떨어진다. 야코비와 가우스-자이델은
250번을 돌아도 10^{-3} 근처다. 오른쪽은 세 감소 인자를 n 에 따라 그린 것이다. — Figure 3:왼쪽에서 켤레기울기가 $k = 20$ 에서 절벽처럼 떨어진다. 야코비와 가우스-자이델은 250번을 돌아도 10^-3 근처다. 오른쪽은 세 감소 인자를 $n$ 에 따라 그린 것이다.

4. 고윳값을 실제로 구하는 법¶

이제 고윳값이다. L21에서 우리는 $\det(A - \lambda I) = 0$ 을 풀어 고윳값을 구했다. 컴퓨터는 절대 그렇게 하지 않는다. 왜인지부터 보자.

특성다항식이라는 함정¶

근이 $1, 2, \dots, 20$ 인 다항식을 만들어 계수 하나를 아주 조금 흔들어 보자. $x^{19}$ 의 계수는 -210 인데, 이것을 상대적으로 10^-10 만큼 건드린다.

근 스무 개 중 열 개가 복소수가 된다. 가장 많이 움직인 근은 2.18 이나 옮겨 간다.

같은 크기로 행렬 쪽을 흔들면 어떻게 되는가. $\diag(1,\dots,20)$ 의 성분을 같은 비율로 흔들면 고윳값은 $4 \times 10^{-9}$ 만 움직인다.

\text{다항식 계수를 흔들 때} \ :\ 2.18 \qquad\text{vs}\qquad \text{행렬을 흔들 때} \ :\ 4 \times 10^{-9}

(22)

(22)의 두 값은 $5 \times 10^8$ 배 차이가 난다.

같은 크기의 상대 섭동이다. 왼쪽에서는 근이 복소평면으로 날아가고 오른쪽에서는
고윳값이 제자리를 지킨다. — Figure 4:같은 크기의 상대 섭동이다. 왼쪽에서는 근이 복소평면으로 날아가고 오른쪽에서는 고윳값이 제자리를 지킨다.

행렬을 다항식으로 바꾸는 순간 문제가 나빠진다.

이것은 L33·보강 1의 $A^{\mathsf T}A$ 와 같은 종류의 실수다. 푸는 방법이 나쁜 것이 아니라 푸는 대상을 자기 손으로 나쁘게 만든 것이다.

그래서 방향이 반대다. numpy.roots 는 다항식의 근을 구할 때 동반행렬을 만들어 그 고윳값을 구한다.

p(x) = x^n + c_{n-1}x^{n-1} + \cdots + c_0 \quad\Longrightarrow\quad C = \begin{bmatrix} -c_{n-1} & -c_{n-2} & \cdots & -c_0 \\ 1 & 0 & \cdots & 0 \\ & \ddots & & \vdots \\ & & 1 & 0 \end{bmatrix}

(23)

(23)의 $C$ 는 특성다항식이 정확히 $p$ 인 행렬이다. 다항식에서 행렬로 가지, 그 반대로 가지 않는다.

거듭제곱법 — L22가 이미 증명해 두었다¶

그러면 어떻게 구하는가. 가장 단순한 방법은 이미 배웠다.

L22의 3절에서 $A^k\vv{x}$ 를 고유기저로 펼치고 $\lambda_1^k$ 로 나눴다.

\frac{A^k\vv{x}}{\lambda_1^k} = c_1\vv{s}_1 + c_2\!\left(\frac{\lambda_2}{\lambda_1}\right)^{k}\!\vv{s}_2 + \cdots \ \longrightarrow\ c_1\vv{s}_1

(24)

(24)의 식이 곧 알고리즘이다. $A$ 를 자꾸 곱하고 길이를 1로 맞추면 지배 고유벡터가 나온다. L22는 이것을 "가장 큰 것만 살아남는다"고만 부르고 알고리즘이라 부르지 않았을 뿐이다.

수렴 속도도 그때 이미 나왔다. 가장 느리게 죽는 항이 $\lambda_2$ 쪽이므로 인자가 $\lvert\lambda_2/\lambda_1\rvert$ 다.

고윳값은 레일리 몫으로 뽑는다(보강 1에서 쓴 그것이다).

\lambda_1 \approx \vv{v}_k^{\mathsf T}A\vv{v}_k \qquad (\lVert\vv{v}_k\rVert = 1)

(25)

(25)의 값은 고유벡터보다 두 배 빨리 수렴한다. 왜 그런지 적어 두자. $A$ 가 대칭이라 하고 $\vv{v} = \vv{q}_1 + \delta\vv{w}$ 로 두자. 여기서 $\vv{w}$ 는 $\vv{q}_1$ 에 수직인 단위벡터다. 분자를 펼치면

\vv{v}^{\mathsf T}A\vv{v} = \lambda_1 + 2\delta\,\vv{q}_1^{\mathsf T}A\vv{w} + \delta^2\,\vv{w}^{\mathsf T}A\vv{w}

(26)

가 되는데, (26)의 가운데 항이 사라진다. $A\vv{q}_1 = \lambda_1\vv{q}_1$ 이고 $\vv{q}_1 \perp \vv{w}$ 이기 때문이다.

\vv{q}_1^{\mathsf T}A\vv{w} = (A\vv{q}_1)^{\mathsf T}\vv{w} = \lambda_1\,\vv{q}_1^{\mathsf T}\vv{w} = 0

(27)

(27)에서 $A$ 가 대칭이라는 것을 썼다. 분모도 $1 + \delta^2$ 이므로 결국 1차 항이 양쪽에서 다 사라진다.

\vv{v}^{\mathsf T}A\vv{v} / \vv{v}^{\mathsf T}\vv{v} = \lambda_1 + O(\delta^2)

(28)

(28)의 결과가 두 배 빠른 이유의 전부다. 고유벡터가 $\delta$ 만큼 틀렸어도 고윳값은 $\delta^2$ 만큼만 틀린다. 대칭이 아니면 이 이득이 없다.

왼쪽에서 파란 점이 회색 이론선 위에 정확히 앉는다. 오른쪽에서 주황(고윳값)의 기울기가
파랑(고유벡터)의 정확히 두 배다. — Figure 5:왼쪽에서 파란 점이 회색 이론선 위에 정확히 앉는다. 오른쪽에서 주황(고윳값)의 기울기가 파랑(고유벡터)의 정확히 두 배다.

QR 알고리즘 — 유사변환을 되풀이한다¶

거듭제곱법은 고윳값을 하나만 준다. 전부 구하려면 다른 것이 필요하다.

QR 알고리즘은 놀랍도록 간단하다. 분해하고, 순서를 바꿔 다시 곱한다.

A_k = Q_k R_k \qquad\Longrightarrow\qquad A_{k+1} = R_k Q_k

(29)

(29)의 두 줄이 전부다. 그런데 왜 이것이 통하는가. $R_k = Q_k^{\mathsf T}A_k$ 를 대입해 보면 정체가 드러난다.

A_{k+1} = R_k Q_k = Q_k^{\mathsf T}A_k Q_k = Q_k^{-1}A_k Q_k

(30)

(30)의 마지막이 유사변환이다(L28). 유사변환은 고윳값을 보존하므로 $A_k$ 는 아무리 되풀이해도 처음과 같은 고윳값을 갖는다. 변하는 것은 비대각 성분뿐이고, 그것이 0으로 간다.

수렴 속도는 이웃한 고윳값의 비로 정해진다. 정확히는 $(i+1, i)$ 성분이 $\lvert\lambda_{i+1}/\lambda_i\rvert^k$ 로 죽으므로, 전체가 0에 가는 속도는

\max_i \left\lvert \frac{\lambda_{i+1}}{\lambda_i} \right\rvert

(31)

이 정한다. (31)의 값이 1에 가까우면, 곧 가까운 고윳값이 어딘가에 하나라도 있으면 전체가 느려진다.

시프트가 판을 바꾼다¶

느린 것은 고칠 수 있다. $\mu$ 를 하나 빼고 하고 다시 더하면 된다.

A_k - \mu I = Q_k R_k, \qquad A_{k+1} = R_k Q_k + \mu I

(32)

(32)의 변형도 여전히 유사변환이라 고윳값은 그대로다. 그런데 비가 $\left\lvert\dfrac{\lambda_{i+1}-\mu}{\lambda_i-\mu}\right\rvert$ 로 바뀌므로, $\mu$ 를 어느 고윳값 가까이 두면 그 비가 아주 작아진다.

$\mu = A_k[n,n]$ 으로 두는 것이 가장 흔한 선택이다. 그 성분이 이미 고윳값에 가까워지고 있기 때문이다. 결과가 극적이다.

되풀이	시프트 없음	시프트 있음
1	1.2	1.2
3	0.9	$3.4\times10^{-1}$
4		$3.7\times10^{-3}$
5		$4.4\times10^{-9}$
6		$6.5\times10^{-27}$

자릿수가 매번 세 배로 늘어난다. 세제곱 수렴이다. 시프트 없이는 40번을 돌아도 10^-12 에 못 간다.

왼쪽에서 행렬이 대각으로 굳어 간다. 대각선의 세 수가 곧 고윳값이다.
오른쪽에서 시프트가 있는 쪽은 여섯 걸음 만에 바닥에 닿는다. — Figure 6:왼쪽에서 행렬이 대각으로 굳어 간다. 대각선의 세 수가 곧 고윳값이다. 오른쪽에서 시프트가 있는 쪽은 여섯 걸음 만에 바닥에 닿는다.

5. 조르당 형이 왜 위험한지, 정확한 언어로¶

L28에서 물었다. 조르당 형은 수학적으로 분명히 존재하는데 왜 컴퓨터가 못 잡는가. 그때는 $\sqrt{\epsilon}$ 법칙을 관찰로만 보였다. 이제 이름을 붙일 수 있다.

고윳값에도 조건수가 있다¶

$A$ 의 고윳값 $\lambda$ 에 대해 오른쪽 고유벡터 $\vv{x}$ 와 왼쪽 고유벡터 $\vv{y}$ 가 있다. 왼쪽 고유벡터란 $\vv{y}^{\mathsf T}A = \lambda\vv{y}^{\mathsf T}$ , 곧 $A^{\mathsf T}\vv{y} = \lambda\vv{y}$ 를 만족하는 것이다.

어디서 나오는지 유도해 보자. $A$ 를 $\varepsilon E$ 만큼 흔들면 고윳값과 고유벡터가 함께 움직인다.

(A + \varepsilon E)\,\vv{x}(\varepsilon) = \lambda(\varepsilon)\,\vv{x}(\varepsilon)

(33)

(33)의 양변을 $\varepsilon$ 으로 미분하고 $\varepsilon = 0$ 을 넣으면

E\vv{x} + A\vv{x}' = \lambda'\vv{x} + \lambda\vv{x}'

(34)

이 된다. (34)에는 모르는 것이 $\vv{x}'$ 과 $\lambda'$ 둘인데, 왼쪽 고유벡터를 곱하면 $\vv{x}'$ 이 통째로 사라진다. $\vv{y}^{\mathsf T}A = \lambda\vv{y}^{\mathsf T}$ 이기 때문이다.

\vv{y}^{\mathsf T}E\vv{x} + \underbrace{\vv{y}^{\mathsf T}A\vv{x}'}_{=\ \lambda\vv{y}^{\mathsf T}\vv{x}'} = \lambda'\,\vv{y}^{\mathsf T}\vv{x} + \lambda\,\vv{y}^{\mathsf T}\vv{x}'

(35)

(35)의 양쪽에서 $\lambda\vv{y}^{\mathsf T}\vv{x}'$ 이 지워지므로 $\lambda'$ 이 곧바로 나온다.

\lambda' = \frac{\vv{y}^{\mathsf T}E\vv{x}}{\vv{y}^{\mathsf T}\vv{x}} \qquad\Longrightarrow\qquad \lvert \delta\lambda \rvert \ \le\ \frac{\lVert E\rVert}{\lvert \vv{y}^{\mathsf T}\vv{x} \rvert}\,\varepsilon

(36)

(36)의 분모가 전부를 결정한다. 그래서 이렇게 이름을 붙인다.

\kappa(\lambda) = \frac{1}{\lvert \vv{y}^{\mathsf T}\vv{x} \rvert} \qquad (\lVert\vv{x}\rVert = \lVert\vv{y}\rVert = 1)

(37)

(37)의 값이 그 고윳값의 조건수다. 좌우 고유벡터가 나란할수록 작고, 수직에 가까울수록 크다.

세 경우를 재 보자.

$A$	$\kappa(\lambda_1)$
대칭 $\begin{pmatrix}4&1\\1&3\end{pmatrix}$	1
비대칭 $\begin{pmatrix}4&1\\0&3\end{pmatrix}$	1.414
거의 결함 $\begin{pmatrix}3&1\\10^{-6}&3\end{pmatrix}$	500

대칭행렬에서 1이다. 좌우 고유벡터가 같으니 $\vv{y}^{\mathsf T}\vv{x} = 1$ 이다. 이것이 L25의 스펙트럼 정리가 주는 또 하나의 선물이다. 대칭행렬의 고윳값은 언제나 완벽하게 조건이 좋다.

결함이면 조건수가 무한대다¶

결함 행렬로 가면 좌우 고유벡터가 수직에 가까워진다. 완전히 결함이면 정확히 수직이 되어 $\vv{y}^{\mathsf T}\vv{x} = 0$ 이고, (37)의 값이 무한대다.

조건수가 무한대라는 것은 (11)의 부등식이 아무 말도 못 한다는 뜻이다. 그때는 1차 근사가 통하지 않고 다른 규칙이 나타난다. 조르당 블록 크기가 $m$ 이면

\lvert \delta\lambda \rvert \ \sim\ \delta^{1/m}

(38)

이다. (38)에서 $m = 2$ 이면 $\sqrt{\delta}$ 이고, $\delta = \epsilon_{\text{mach}} \approx 10^{-16}$ 이면 10^-8 이다. L28에서 측정한 기울기 0.5가 바로 이것이다.

6. 자주 하는 오해¶

① 잔차가 작으면 답이 정확하다. 아니다. 잔차 $\lVert A\hat{\vv{x}} - \vv{b}\rVert$ 가 작다는 것은 후방오차가 작다는 뜻이지 전방오차가 작다는 뜻이 아니다. (11)의 부등식에 조건수가 곱해져 있다. 조건수가 10¹² 이면 잔차가 10^-16 이어도 답은 10^-4 만큼 틀릴 수 있다.

② 조건수가 크면 알고리즘이 나쁜 것이다. 아니다. 조건수는 문제의 성질이다. 세상에서 가장 좋은 알고리즘을 써도 줄지 않는다. 알고리즘이 책임지는 것은 후방오차 하나뿐이다.

③ 반복법이 소거보다 언제나 빠르다. 아니다. 반복법은 행렬이 크고 성기고 조건수가 나쁘지 않을 때 이긴다. $n$ 이 작거나 조건수가 크면 소거가 낫다. 그리고 Theorem 1의 조건이 깨지면 반복법은 아예 발산한다. 소거는 적어도 답은 준다.

④ 켤레기울기는 $n$ 번이면 끝나니 유한 알고리즘이다. 이론적으로는 그렇지만 실무에서는 아니다. 반올림 때문에 $A$ -직교성이 무너져 $n$ 번에 정확히 끝나지 않는다. 그리고 애초에 $n$ 이 백만이면 $n$ 번도 못 돈다. $\sqrt{\kappa}$ 로 정해지는 근사 수렴이 실제로 쓰는 성질이다.

⑤ 고윳값을 구하려면 특성다항식을 풀어야 한다. 정반대다. 4절이 답이다. 실무는 QR 알고리즘을 쓰고, 다항식의 근을 구할 때조차 동반행렬을 만들어 고윳값 쪽으로 간다.

⑥ 대칭행렬은 계산이 쉬우니 별 이득이 없다. 이득이 크다. 대칭이면 고윳값의 조건수가 정확히 1이고, eigh 가 eig 보다 빠르며, 결함이 있을 수 없어 5절의 함정이 아예 없다. 대칭인지 먼저 확인하는 것이 L32의 흐름도가 그 자리에 있는 이유다.

마치며...¶

이번 글에서는 이 책이 다섯 번 한 경고를 하나로 묶었다.

실패하는 방식은 둘이다. 비용과 정확도. 여인수와 크래머는 앞엣것이고 조르당·정규방정식·공분산은 뒤엣것이다.
정확도 이야기는 전방오차 $\le$ 조건수 $\times$ 후방오차 한 줄로 정리된다. 후방오차는 알고리즘의 몫이고 조건수는 문제의 몫이다.
반복법의 수렴은 스펙트럼 반지름 $< 1$ 이 정한다. L22의 $A^k \to 0$ 이 그대로 알고리즘의 판정 기준이 되었다.
켤레기울기가 특별한 것은 조건수에 제곱근이 붙기 때문이다.
고윳값은 특성다항식으로 구하지 않는다. 거듭제곱법은 L22에서 이미 증명해 두었고, QR 알고리즘은 유사변환을 되풀이하는 것(L28)이며, 시프트를 넣으면 세제곱으로 수렴한다.
그리고 조르당 형이 왜 위험한지 정확한 언어로 말할 수 있다. 결함 고윳값은 조건수가 무한대다. 알고리즘을 바꿔서 될 일이 아니라 문제를 바꿔야 한다.

돌아보면 이 책의 후반부가 계속 같은 이야기를 하고 있었다. 존재하는 것과 손에 쥘 수 있는 것은 다르다. 그 차이를 재는 자가 조건수다.

다음은 보강 3 — 그래프 라플라시안과 스펙트럴 클러스터링이다. 고윳값이 시간의 흐름을 예측하고(L23) 확률의 정착점을 알려주었는데(L24), 이번에는 그래프를 자르는 선까지 알려준다.

이번 글의 내용을 파이썬으로 확인해 보려면 보강 2 실습 노트북으로 넘어가면 된다. 세 반복법을 같은 문제에 걸어 감소 곡선을 겹쳐 보고, 거듭제곱법이 정말 $\lvert\lambda_2/\lambda_1\rvert$ 로 수렴하는지 재고, 특성다항식과 행렬을 같은 크기로 흔들어 어느 쪽이 무너지는지 직접 볼 수 있다.