행렬 거듭제곱과 거듭제곱법: 스펙트럼 반지름이 정하는 것

순환망을 100 스텝 펼쳐 학습시키면 로그에 둘 중 하나가 찍힙니다. 기울기 노름이 1.4e+08처럼 손댈 수 없이 커지거나, 3.2e-11처럼 0과 구별되지 않게 작아집니다. 학습률을 열 배 줄여도 앞의 것은 스텝만 잘게 쪼개질 뿐 여전히 커지고, 뒤의 것은 아예 움직이지 않습니다.

학습률의 문제가 아닙니다. 시간 축을 거슬러 올라가는 역전파는 스텝마다 같은 가중치 행렬을 한 번씩 곱하므로, 100 스텝이면 같은 행렬을 100번 곱한 것이 기울기에 그대로 붙습니다. 그러니 물어야 할 것은 하나입니다 — 행렬을 계속 곱하면 무슨 일이 벌어지는가.

지난 글에서 $A^k = P\Lambda^k P^{-1}$ 까지 왔습니다. 이 글은 그 식을 끝까지 밀어붙여 답을 하나의 수로 줄이고, 그 수를 구하는 방법까지 갑니다.

고유기저에서 펼친 A^k

방향마다 갈라지는 운명

행렬을 벡터에 반복해 먹이는 상황을 봅니다. $A$ 가 대각화되고 고유벡터 $\mathbf{v}_1, \dots, \mathbf{v}_n$ 이 기저를 이룬다면, 출발 벡터 $\mathbf{x}$ 를 그 기저의 좌표로 적을 수 있습니다.

$\mathbf{x} = c_1\mathbf{v}_1 + c_2\mathbf{v}_2 + \cdots + c_n\mathbf{v}_n$

$A$ 를 한 번 먹이면 각 항이 자기 고윳값만큼 늘어납니다. $A\mathbf{v}_i = \lambda_i\mathbf{v}_i$ 이므로 항끼리 섞이지 않는다는 것이 핵심입니다.

$A\mathbf{x} = c_1\lambda_1\mathbf{v}_1 + c_2\lambda_2\mathbf{v}_2 + \cdots + c_n\lambda_n\mathbf{v}_n$

$k$ 번 먹이면 그 배율이 $k$ 번 곱해질 뿐입니다.

$A^k\mathbf{x} = c_1\lambda_1^k\mathbf{v}_1 + c_2\lambda_2^k\mathbf{v}_2 + \cdots + c_n\lambda_n^k\mathbf{v}_n$

행렬의 반복이 수의 거듭제곱 $n$ 개로 갈라졌습니다. 방향마다 따로 운명이 정해지고, 그 운명은 $|\lambda_i|$ 하나가 정합니다.

고윳값에 절댓값이 붙는 이유를 짚어 둡니다. $\lambda = -2$ 이면 $\lambda^k$ 는 $-2, 4, -8, 16$ 으로 부호가 번갈아 바뀌지만 크기는 꾸준히 두 배씩 커집니다. 폭발하느냐 소실하느냐를 정하는 것은 부호가 아니라 크기입니다.

대각화가 안 되는 자리

위 전개는 첫 줄에 조건을 하나 달고 시작했습니다 — 고유벡터가 기저를 이룬다는 조건입니다. 이것이 깨지면 식 자체가 안 세워집니다. 고윳값은 있는데 그에 딸린 독립인 고유벡터가 모자란 행렬이 있고, 그런 자리를 조르당 블록이라고 부릅니다. 가장 작은 예가 이것입니다.

$J = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}, \qquad J^k = \begin{pmatrix} 1 & k \\ 0 & 1 \end{pmatrix}$

$J$ 의 고윳값은 1 하나뿐이고, 고유벡터도 $(1,0)$ 방향 하나뿐이라 2차원의 기저를 못 만듭니다. 그런데 거듭제곱을 직접 해 보면 오른쪽 위 칸이 $k$ 로 자랍니다. $k = 100$ 이면 그 칸이 100입니다.

가장 큰 고윳값의 크기가 정확히 1인데 크기가 자랍니다. 뒤에서 세울 "크기를 정하는 것은 고윳값뿐"이라는 말이 여기서는 성립하지 않습니다. 다만 자라는 속도가 $k$ 에 비례할 뿐이라, 지수로 자라는 쪽과는 자릿수가 전혀 다릅니다.

무작위 행렬과 가까운 고윳값

실무에서 이 예외를 얼마나 걱정해야 할까요. 답은 "거의 안 해도 된다"입니다. 대각화가 안 되려면 특성방정식이 중근을 가져야 하는데, 원소를 연속분포에서 뽑으면 그런 일이 생길 확률이 0입니다. 초기화한 가중치 행렬은 거의 언제나 대각화됩니다.

값을 치르는 자리는 따로 있습니다. 두 고윳값이 정확히 같지는 않고 아주 가까운 경우입니다. 이때 대각화는 성립하지만 $P^{-1}$ 을 만들 때 나누는 수가 작아져, 계산이 수치적으로 조르당 블록 쪽에 붙습니다. 실제로 겪는 증상도 이 글 뒤쪽의 거듭제곱법에서 그대로 나옵니다 — 두 고윳값이 붙어 있으면 수렴이 한없이 느려집니다. "대각화되느냐"는 예 아니오로 갈리지만 "얼마나 잘 대각화되느냐"는 정도 문제이고, 계산에서 걸리는 것은 언제나 뒤쪽입니다.

스펙트럼 반지름

정의와 이름

$n$ 개의 $|\lambda_i|$ 중 가장 큰 것 하나가 나머지를 압도합니다. 그 수에 이름을 붙입니다.

정의. 행렬 $A$ 의 고윳값들의 절댓값 중 가장 큰 값을 $A$ 의 스펙트럼 반지름이라 하고 $\rho(A)$ 로 적는다.

$\rho(A) = \max_i |\lambda_i|$

이름에 '반지름'이 들어간 것은 고윳값들을 복소평면에 점으로 찍었을 때 그 점들을 전부 담는 원점 중심 원의 반지름이 이 값이기 때문입니다. 회전이 섞인 행렬에서 고윳값이 복소수로 나온다는 것은 지난 글에서 봤고, 그때도 절댓값은 실수 하나로 정해집니다.

$k$ 를 키우면 $\rho(A)$ 에 해당하는 항만 살아남습니다. 가장 큰 것으로 나눠 보면 이유가 보입니다.

$\frac{A^k\mathbf{x}}{\lambda_1^k} = c_1\mathbf{v}_1 + c_2\left(\frac{\lambda_2}{\lambda_1}\right)^k\mathbf{v}_2 + \cdots$

$|\lambda_1| > |\lambda_2| \ge \cdots$ 로 번호를 매겼으므로 두 번째 항부터는 괄호 안의 크기가 1보다 작고, $k$ 제곱하면 0으로 갑니다. 남는 것은 $c_1\mathbf{v}_1$ 입니다.

$\rho(A)$	$A^k$ 가 하는 일	역전파에서
$\rho < 1$	모든 방향이 0으로 줄어든다	기울기 소실
$\rho = 1$	크기가 유지되는 방향이 있다	경계 — 안정하지만 아슬아슬하다
$\rho > 1$	최대 방향이 지수적으로 커진다	기울기 폭발

세 갈래를 가르는 문턱이 정확히 1이라는 점이 중요합니다. $0.95$ 와 $1.05$ 는 0.1밖에 차이 나지 않지만 100 스텝 뒤에는 $0.95^{100} \approx 0.0059$ 와 $1.05^{100} \approx 131.5$ 로 4만 배 넘게 벌어집니다. 지수는 작은 차이를 봐주지 않습니다.

극한으로 읽기

$\rho$ 와 행렬의 크기를 잇는 식이 하나 있습니다.

$\rho(A) = \lim_{k \to \infty} \|A^k\|^{1/k}$

오른쪽이 극한이라는 것이 이 식의 전부입니다. 유한한 $k$ 에서 $\|A^k\|$ 가 $\rho^k$ 보다 클 수 있고, 실제로 자주 그렇습니다. 극단적인 예가

$N = \begin{pmatrix} 0 & 10 \\ 0 & 0 \end{pmatrix}$

입니다. 고윳값이 0 하나뿐이라 $\rho(N) = 0$ 인데 $\|N\| = 10$ 이고, $N^2$ 은 그제서야 영행렬이 됩니다. 앞 절의 조르당 블록도 $\rho = 1$ 인데 $\|J^k\| \approx k$ 라 같은 일이 벌어집니다.

그래서 읽는 법을 정확히 해 둡니다. $\rho < 1$ 이 말하는 것은 "언제나 줄어든다"가 아니라 "충분히 큰 $k$ 에서는 줄어든다" 입니다. 초기 몇 스텝에서 노름이 오히려 커지는 순환망을 보고 계산이 틀렸다고 의심할 일이 아닙니다 — 그 구간은 극한이 아직 안 온 구간입니다.

행합으로 가두기

$\rho$ 를 구하려면 고윳값을 풀어야 하는데, 성분이 전부 0 이상인 행렬에는 풀지 않고 가두는 방법이 있습니다.

성분이 전부 0 이상이면 행합의 최솟값과 최댓값 사이에 $\rho$ 가 있다.

$\min_i \sum_j A_{ij} \ \le\ \rho(A)\ \le\ \max_i \sum_j A_{ij}$

$A = \begin{pmatrix} 0.5 & 0.2 \\ 0.1 & 0.4\end{pmatrix}$ 은 행합이 $0.7$ 과 $0.5$ 이므로 $\rho \le 0.7 < 1$ 이고, 고윳값을 안 풀고도 이 행렬은 반드시 줄어든다고 말할 수 있습니다. 실제 고윳값은 $0.6$ 과 $0.3$ 이라 상한 안에 들어 있습니다.

반대로 최솟값 쪽이 1을 넘으면 폭발이 확정됩니다. 다음 절에서 볼 $W$ 는 행합이 $1.5$ 와 $1.0$ 이라 둘 다 문턱을 안 가르므로, 이 방법으로는 결론이 안 나고 고윳값을 풀어야 합니다. 가두는 폭이 좁을 때만 쓸모가 있는 도구라는 뜻입니다.

순환망 기울기의 자릿수

W의 고윳값 읽기

숫자를 넣어 확인합니다. 순환 가중치가

$W = \begin{pmatrix} 0.9 & 0.6 \\ 0.2 & 0.8 \end{pmatrix}$

라고 합시다. 성분이 전부 1보다 작으니 반복해 곱하면 줄어들 것 같지만, 정하는 것은 성분이 아니라 고윳값입니다. 지난 글의 두 줄로 바로 읽습니다.

$\operatorname{tr}W = 0.9 + 0.8 = 1.7, \qquad \det W = 0.9 \times 0.8 - 0.6 \times 0.2 = 0.72 - 0.12 = 0.6$

$\lambda^2 - 1.7\lambda + 0.6 = 0 \quad \Longrightarrow \quad (\lambda - 1.2)(\lambda - 0.5) = 0 \quad \Longrightarrow \quad \lambda = 1.2,\ 0.5$

$\rho(W) = 1.2$ 입니다. 1을 넘었으므로 이 순환망은 폭발합니다. 성분을 아무리 들여다봐도 안 보이던 것이 대각합과 행렬식 두 번의 계산으로 나왔습니다.

얼마나 커지는지도 예측됩니다. $1.2^{10} \approx 6.19$, $1.2^{50} \approx 9.1\times10^3$, $1.2^{100} \approx 8.3\times10^7$ 입니다. 실제로 $W^{100}$ 의 크기를 재면 $9.5\times10^7$ 로 같은 자릿수입니다 — 어긋나는 부분은 $c_1$ 과 $\mathbf{v}_1$ 이 만드는 상수배뿐이고, 자릿수는 $\rho^k$ 가 정합니다. 작은 쪽 고윳값 $0.5$ 는 $0.5^{100} \approx 7.9\times10^{-31}$ 이라 100 스텝 뒤에는 흔적도 없습니다.

활성함수의 도함수가 끼는 자리

실제 역전파에 곱해지는 것은 $W$ 혼자가 아닙니다. 한 스텝의 상태가 $h_t = \sigma(W h_{t-1} + U x_t)$ 이면, 한 스텝을 거슬러 갈 때 곱해지는 행렬은

$\frac{\partial h_t}{\partial h_{t-1}} = \operatorname{diag}\bigl(\sigma'(z_t)\bigr)\,W$

입니다. 여기서 $\operatorname{diag}(\sigma'(z_t))$ 는 그 스텝의 활성함수 도함수를 대각선에 늘어놓은 행렬이고, 한 스텝을 통과할 때마다 이득을 한 번 더 깎는 역할을 합니다.

깎이는 정도가 활성함수에 따라 다릅니다. tanh의 도함수는 $1 - \tanh^2$ 이라 값이 0과 1 사이이고, 입력이 조금만 커져도 빠르게 0에 붙습니다. 그러면 곱해지는 행렬의 크기가 $W$ 자신보다 작아지므로, tanh를 쓰는 순환망에서는 폭발보다 소실이 훨씬 흔합니다. 위의 $W$ 처럼 $\rho$ 가 1.2인 행렬을 써도 도함수가 평균 0.5쯤 곱해지면 실효 이득이 0.6이 되어 오히려 줄어듭니다.

ρ 하나로는 못 읽는 것

그런데 이 대각행렬은 스텝마다 다릅니다. $z_t$ 가 그 시점의 입력과 상태로 정해지기 때문입니다. 그래서 100 스텝을 거슬러 갈 때 곱해지는 것은 같은 행렬의 100제곱이 아니라

$\prod_{t} \operatorname{diag}\bigl(\sigma'(z_t)\bigr)\,W$

처럼 서로 다른 행렬 100개의 곱입니다. 고유기저에서 갈라 쓰는 이 글의 방법은 같은 행렬을 반복할 때만 통하므로, 여기서는 $\rho$ 하나로 값을 예측할 수 없습니다.

남는 것은 상한뿐입니다. 곱의 크기는 크기의 곱을 못 넘으므로, 도함수의 최댓값을 $m$ 이라 하면 전체 이득이 $(m\,\|W\|)^T$ 이하입니다. 이 상한이 1보다 작으면 소실은 확정이지만, 상한이 1보다 크다고 폭발이 확정되지는 않습니다. 기울기 폭발·소실을 다루는 방법 — 클리핑, 게이트 구조, 초기화 — 이 이 글의 몫이 아닌 이유도 여기에 있습니다. 여기서 답한 것은 하나입니다. 왜 하필 지수적인가. 같은 것을 반복해 곱하기 때문이고, 반복 곱셈은 고유기저에서 수의 거듭제곱이기 때문입니다.

거듭제곱법

절차

$k$ 가 커지면 $A^k\mathbf{x}$ 가 $\mathbf{v}_1$ 방향으로 쏠린다는 사실은 지금까지 골칫거리였습니다. 그런데 최대 고유벡터를 찾고 싶다면 그것이 그대로 알고리즘입니다. 아무 벡터나 잡고 계속 곱하면 저절로 그 방향이 됩니다.

문제는 곱할수록 길이가 $\rho^k$ 로 커지거나 작아져 수치가 넘친다는 것뿐이니, 매 스텝 길이를 다시 맞춰 주면 됩니다. 이 절차를 거듭제곱법이라고 합니다 — 행렬을 반복해 곱하는 것만으로 최대 고윳값과 그 고유벡터를 얻는 방법입니다.

거듭제곱법. $\mathbf{x}_0$ 를 아무렇게나 잡고 다음을 반복한다.

$\mathbf{y}_{k} = A\mathbf{x}_{k-1}, \qquad \mathbf{x}_{k} = \frac{\mathbf{y}_{k}}{\|\mathbf{y}_{k}\|}$

$\mathbf{x}_k$ 는 $\mathbf{v}_1$ 방향으로 수렴하고, $\|\mathbf{y}_k\|$ 는 $|\lambda_1|$ 로 수렴한다.

$A = \begin{pmatrix} 2 & 1 \\ 1 & 2\end{pmatrix}$ 로 손으로 돌려 봅니다. 정규화는 첫 성분으로 나누는 방식을 쓰겠습니다 — 길이를 1로 맞추는 것보다 손계산이 깨끗하고, 방향만 보는 데는 차이가 없습니다.

$\mathbf{x}_0 = (1, 0)$ 에서 출발합니다.

$A(1, 0) = (2, 1) \ \longrightarrow\ \text{2로 나눠} \ (1,\ 0.5)$

$A(1, 0.5) = (2.5,\ 2) \ \longrightarrow\ \text{2.5로 나눠} \ (1,\ 0.8)$

$A(1, 0.8) = (2.8,\ 2.6) \ \longrightarrow\ \text{2.8로 나눠} \ (1,\ 0.9286)$

$A(1, 0.9286) = (2.9286,\ 2.8571) \ \longrightarrow\ \text{2.9286으로 나눠} \ (1,\ 0.9756)$

지난 글에서 이 행렬의 고윳값이 $3$ 과 $1$, 최대 고유벡터가 $(1,1)$ 이라는 것을 이미 구해 두었습니다. 표가 정확히 그리로 갑니다. 고윳값을 푸는 방정식을 세운 적이 없는데도 곱셈과 나눗셈만으로 같은 답이 나왔습니다.

레일리 몫

나눈 값을 고윳값 추정으로 쓰는 것보다 나은 방법이 있습니다. 지금 벡터가 $A$ 방향으로 얼마나 늘어났는지를 직접 재는 것입니다.

$R(\mathbf{x}) = \frac{\mathbf{x}^\top A\mathbf{x}}{\mathbf{x}^\top \mathbf{x}}$

이 값을 레일리 몫이라고 합니다. $\mathbf{x}$ 가 정확히 고유벡터이면 분자가 $\lambda\,\mathbf{x}^\top\mathbf{x}$ 가 되어 값이 $\lambda$ 와 같아집니다. 고유벡터가 아니면 그 근처의 값을 냅니다.

같은 반복에 이 열을 붙이면 차이가 바로 보입니다.

$k$	$\mathbf{x}_k$	나눈 값	오차	레일리 몫	오차
1	$(1,\ 0.5)$	2	1	2.8	0.2
2	$(1,\ 0.8)$	2.5	0.5	2.9756	0.0244
3	$(1,\ 0.9286)$	2.8	0.2	2.99709	0.0029
4	$(1,\ 0.9756)$	2.9286	0.0714	2.99971	0.00029

나눈 값의 오차는 한 스텝에 3분의 1씩 줄고, 레일리 몫의 오차는 9분의 1씩 줍니다. 대칭 행렬에서 레일리 몫의 오차는 방향의 오차를 제곱한 만큼이기 때문입니다. 방향이 소수 둘째 자리까지 맞으면 고윳값은 넷째 자리까지 맞는다는 뜻이고, 같은 반복 횟수로 자릿수를 두 배 얻는 셈입니다.

정규화와 멈출 때

손계산에서 쓴 "첫 성분으로 나누기"는 실제 코드에 그대로 쓰면 안 됩니다. 최대 고유벡터의 첫 성분이 0에 가까우면 그 나눗셈이 터지기 때문입니다. 512차원 무작위 행렬에서 첫 성분은 흔히 최대 성분의 10분의 1쯤이고, 하필 아주 작게 뽑히면 나눈 값에 반올림 오차가 크게 실립니다. 길이로 나누는 쪽은 그런 자리가 없습니다 — 벡터가 영벡터가 아닌 한 분모가 작아지지 않습니다.

멈출 때는 스텝 수를 정해 두는 대신 잔차를 봅니다. 지금 추정값 $\hat\lambda$ 에 대해

$r = \|A\mathbf{x} - \hat\lambda \mathbf{x}\|$

를 재는 것입니다. $\mathbf{x}$ 의 길이가 1이고 $A$ 가 대칭이면 $r$ 안에 진짜 고윳값이 하나 있다는 것이 보장되므로, 이 값이 곧 오차의 상한입니다. 고윳값을 모르는 채로 오차를 아는 유일한 길이라 실무 구현은 거의 다 이 값으로 멈춥니다.

수렴 속도와 스펙트럼 갭

1등과 2등의 비

위 표의 오차가 대략 3분의 1씩 줄었습니다. 우연이 아닙니다. 앞에서 $\lambda_1^k$ 로 나눈 식을 다시 봅니다.

$\frac{A^k\mathbf{x}}{\lambda_1^k} = c_1\mathbf{v}_1 + c_2\left(\frac{\lambda_2}{\lambda_1}\right)^k\mathbf{v}_2 + \cdots$

목표에서 벗어난 양은 두 번째 항이 지배하고, 그 크기는 $|\lambda_2/\lambda_1|^k$ 입니다. 이 예에서는 $|1/3| = 0.333$ 이므로 한 스텝마다 오차가 3분의 1로 줄어듭니다.

수렴 속도. 거듭제곱법의 오차는 $k$ 스텝 뒤 $\left|\dfrac{\lambda_2}{\lambda_1}\right|^k$ 에 비례해 줄어든다.

이 비를 스펙트럼 갭이라고 부르기도 합니다 — 1등과 2등이 벌어진 정도라는 뜻입니다. 갭이 크면 금방 붙고, 두 고윳값이 비슷하면 한없이 느려집니다. $B = \begin{pmatrix} 3 & 1 \\ 1 & 3\end{pmatrix}$ 은 고윳값이 $4$ 와 $2$ 라 비가 $0.5$ 인데, 같은 출발점에서 돌리면 8 스텝에도 $(1, 0.9922)$ 로 앞의 예가 6 스텝에 도달한 자리에 겨우 옵니다.

시프트로 갭을 벌리기

갭이 데이터에 주어진 값이라고 포기할 일은 아닙니다. $A$ 대신 $A - \mu I$ 를 돌리면 고윳값이 전부 $\mu$ 만큼 평행이동하고 고유벡터는 그대로입니다. 이것을 시프트라고 합니다.

$\lambda_1 = 5.75, \lambda_2 = 4.48$ 인 행렬은 비가 $0.78$ 이라 느립니다. 여기에 $\mu = 2.5$ 를 빼면 두 값이 $3.25$ 와 $1.98$ 이 되어 비가 $0.61$ 로 내려갑니다. 차이는 그대로인데 비가 줄어든 것이고, 나누는 쪽이 작아졌으니 당연합니다.

$\mu$ 를 고르는 감각은 한 줄입니다 — 노리지 않는 고윳값 무리를 원점 근처로 모으는 값을 고릅니다. 고윳값이 전부 양수인 행렬이면 그 무리가 0부터 퍼져 있으므로 아래쪽을 잘라 내는 양수 $\mu$ 가 통하고, 양끝에 걸쳐 있는 행렬이면 이미 가운데가 0 근처라 시프트로 얻을 것이 별로 없습니다. 잘못 고르면 비가 오히려 커지므로 공짜가 아니라 값을 아는 만큼 버는 방법입니다.

역거듭제곱법과 켤레쌍

같은 절차를 $A^{-1}$ 에 돌리면 무엇이 나올까요. $A\mathbf{v} = \lambda\mathbf{v}$ 의 양변에 $A^{-1}/\lambda$ 를 곱하면 $A^{-1}\mathbf{v} = (1/\lambda)\mathbf{v}$ 이므로, 고유벡터는 그대로이고 고윳값만 역수가 됩니다. 그러면 원래 가장 작았던 고윳값이 역수를 취해 가장 커지므로, 거듭제곱법이 그쪽으로 수렴합니다. 이것이 역거듭제곱법입니다. 실제로는 역행렬을 만들지 않고 매 스텝 $A\mathbf{y} = \mathbf{x}$ 를 풀어 씁니다.

수렴하지 않는 자리도 분명히 해 둡니다. $|\lambda_1| = |\lambda_2|$ 이면 비가 1이라 오차가 줄지 않습니다. 실수 행렬에서 이 일이 벌어지는 흔한 이유가 고윳값이 켤레복소수 한 쌍으로 나오는 경우입니다 — 켤레는 절댓값이 같기 때문입니다. 지난 글의 90도 회전 행렬이 정확히 그런 예이고, 벡터를 계속 곱하면 붙는 대신 계속 돕니다. 이때는 한 스텝이 아니라 두 스텝을 묶어 보아야 무엇이 도는지가 보입니다.

한 가지 더 있습니다. 출발 벡터가 하필 $c_1 = 0$ 이면, 즉 $\mathbf{v}_1$ 성분이 아예 없으면 그 방향은 영원히 생기지 않습니다. 실무에서 무작위 벡터로 시작하는 이유가 이것입니다 — 무작위로 뽑은 벡터의 특정 방향 성분이 정확히 0일 확률은 0이고, 설령 0이더라도 부동소수점 반올림이 곧 아주 작은 성분을 만들어 냅니다.

코드로 확인하기

손계산을 코드로

import numpy as np

A = np.array([[2., 1], [1, 2]])

x = np.array([1., 0])
for k in range(6):
    y = A @ x
    x = y / y[0]                  # 첫 성분으로 정규화
    print(k + 1, np.round(x, 6), round(y[0], 6))
# 6 [1.       0.99726] 2.991803  —  (1,1) 방향, 고윳값 3으로 간다

스펙트럼 반지름과 폭발 여부는 두 줄이면 됩니다.

W = np.array([[0.9, 0.6], [0.2, 0.8]])
rho = max(abs(np.linalg.eigvals(W)))
print(rho)                                          # 1.2  — 1을 넘는다

for k in (10, 50, 100):
    print(k, np.linalg.norm(np.linalg.matrix_power(W, k)), rho ** k)
# 100 95385693.8...   82817974.5...   자릿수가 같다

512×512에서의 시간과 정확도

이제 손으로는 못 하는 크기로 갑니다. 무작위 대칭 행렬에 방향 하나를 심어 $\lambda_1 = 5.75$, $\lambda_2 = 4.48$ 이 되게 만들고, 거듭제곱법과 전체 고윳값 분해를 나란히 돌렸습니다.

방법                    시간        결과
np.linalg.eigh        0.027초     고윳값 512개 전부
거듭제곱법 30 스텝      0.002초     최대 고윳값 하나

 k    첫 성분 정규화 오차    레일리 몫 오차
10        9.84e-01           3.86e-01
20        6.83e-02           1.39e-03
30        1.08e-04           6.21e-06
40        2.09e-04           3.33e-08

읽을 것이 셋입니다. 첫째, 레일리 몫이 앞 절의 예측대로 한 자릿수 이상 앞섭니다. 20 스텝에서 두 열이 이미 50배 벌어져 있습니다.

둘째, $k = 40$ 에서 왼쪽 열의 오차가 오히려 커졌습니다. 앞에서 말한 첫 성분 정규화의 위험이 그대로 나타난 것입니다 — 이 행렬의 최대 고유벡터에서 첫 성분은 $-0.0134$ 로 최대 성분의 9분의 1이라, 그 수로 나누면서 반올림 오차가 실려 그 아래로는 못 내려갑니다.

셋째, 시간을 보면 512×512에서는 전체 분해가 느리지 않습니다. 0.027초면 충분히 빠릅니다. 거듭제곱법의 값어치는 이 크기에서 나오지 않고, 행렬을 통째로 메모리에 올릴 수 없거나 곱셈만 가능한 크기에서 나옵니다.

시프트가 줄이는 스텝 수

같은 행렬에 $\mu$ 만 바꿔 넣고 레일리 몫의 오차가 $10^{-8}$ 아래로 내려갈 때까지의 스텝 수를 셌습니다.

$\mu$	$\|\lambda_2/\lambda_1\|$	스텝 수
0	0.779	43
1.0	0.733	34
2.0	0.661	26
2.5	0.609	23

비가 0.78에서 0.61로 줄자 스텝이 절반 가까이 빠졌습니다. 이 행렬은 고윳값이 전부 양수라 아래쪽을 잘라 내는 시프트가 그대로 통하는 경우이고, 그래서 $\mu$ 를 키울수록 좋아집니다. 양끝에 고윳값이 걸쳐 있는 행렬에서 같은 표를 만들면 반대로 스텝이 늘어납니다 — 시프트는 스펙트럼의 생김새를 알 때만 버는 방법입니다.

정리

• $A^k\mathbf{x}$ 를 고유기저에서 펼치면 $\sum c_i\lambda_i^k\mathbf{v}_i$ 이고, 방향마다 자기 고윳값의 거듭제곱만큼 늘어납니다. 고유벡터가 기저를 못 이루면 이 전개 자체가 성립하지 않고, $\rho = 1$ 인데 크기가 $k$ 에 비례해 자라는 조르당 블록이 그 예입니다.
• 스펙트럼 반지름 $\rho(A) = \max|\lambda_i|$ 하나가 반복 곱셈의 운명을 정합니다. 다만 $\rho = \lim\|A^k\|^{1/k}$ 은 극한이라, 유한한 $k$ 에서는 $\|A^k\|$ 가 $\rho^k$ 를 넘을 수 있습니다.
• 성분이 전부 0 이상이면 행합의 최솟값과 최댓값 사이에 $\rho$ 가 갇힙니다. 가두는 폭이 문턱 1을 안 가를 때만 고윳값을 풀면 됩니다.
• 순환망의 기울기가 지수적인 이유가 이것입니다. 다만 실제로 곱해지는 것은 $\operatorname{diag}(\sigma'(z_t))W$ 이고 스텝마다 다르므로, $\rho$ 하나로는 못 읽고 상한만 남습니다. tanh의 도함수가 1을 안 넘어 소실 쪽이 더 흔합니다.
• 거듭제곱법은 그 지배 현상을 거꾸로 쓴 것입니다. 곱하고 정규화하기를 반복하면 최대 고유벡터가 나오고, 고윳값은 나눈 값보다 레일리 몫으로 읽는 편이 오차가 제곱으로 작아집니다. 멈추는 기준은 잔차 $\|A\mathbf{x} - \hat\lambda\mathbf{x}\|$ 입니다.
• 수렴 속도는 $|\lambda_2/\lambda_1|^k$ 입니다. $A - \mu I$ 로 시프트하면 이 비를 줄일 수 있고, 512×512 실험에서 0.779를 0.609로 내리자 43 스텝이 23 스텝이 됐습니다. 절댓값이 같은 켤레쌍에서는 수렴하지 않습니다.
• 고윳값 방정식을 풀지 않고 행렬-벡터 곱만으로 최대 고윳값을 얻는다는 점이 실전에서 중요합니다. 512×512에서는 전체 분해가 더 빠르고, 이 방법이 값어치를 갖는 것은 행렬을 통째로 못 들고 있는 크기에서부터입니다.

곱셈 몇 번으로 최대 고윳값을 얻는 이 절차는 뒤에서 다시 쓰입니다. 층의 가중치를 그 최댓값으로 나눠 크기를 묶어 두는 스펙트럼 정규화가 매 스텝 필요로 하는 수가 바로 이것이고, 그 이야기는 「중급 17번 · 행렬 노름」의 몫입니다. 다음 글은 방향을 조금 틀어, 고윳값이 애초에 얌전하게 나오는 행렬 — 대칭 행렬을 봅니다. AI에서 만나는 행렬이 왜 대부분 대칭인지, 그리고 대칭이면 무엇이 공짜로 따라오는지입니다.

---

읽어주셔서 감사합니다. 😊