대각화와 거듭제곱

같은 규칙을 계속 반복하는 일은 생각보다 흔합니다. 오늘 날씨로 내일 날씨를 예측하고, 그 예측으로 모레를 예측합니다. 이번 달 활성 사용자와 휴면 사용자 수로 다음 달을 계산하고, 그다음 달을 또 계산합니다. 게임 속 경제, 인구 변화, 웹 페이지의 중요도 계산도 모두 "벡터에 같은 행렬을 거듭 곱하기", 곧 Aⁿx입니다.

문제는 n이 클 때입니다. 행렬을 100번 곱하면 계산이 많을 뿐 아니라, 결과가 어디로 가는지 미리 알기도 어렵습니다. 이번 교시의 열쇠는 11교시의 고유벡터를 10교시의 기저로 쓰는 것입니다. 고유벡터 방향에서 보면 행렬은 "축마다 늘이기"일 뿐이고, 늘이기는 100번 반복해도 숫자를 100제곱하면 끝납니다.

이 생각 하나로 피보나치 수열에서 황금비가 나오는 이유, 날씨가 결국 일정한 비율에 머무는 이유가 함께 풀립니다. 마지막에는 고유값을 반복 계산만으로 찾아내는 거듭제곱법을 직접 짜 봅니다.

12-1. 고유벡터 기저에서 보면 그냥 늘이기 (8분)

11교시의 $A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$를 다시 봅시다. 고유벡터는 $\mathbf{e}_1 = (1, 1)$(고유값 3)과 $\mathbf{e}_2 = (1, -1)$(고유값 1)이었습니다. 두 방향이 서로 다르니 이 둘로 평면의 모든 벡터를 만들 수 있습니다(4교시). 곧 고유벡터 두 개가 기저가 됩니다.

$\mathbf{v} = (3, 1)$을 이 기저로 나타내면 $\mathbf{v} = 2\,\mathbf{e}_1 + 1\,\mathbf{e}_2$입니다. 이제 $A$를 곱해 봅시다. 고유벡터마다 하는 일이 "$\lambda$배" 뿐이니 계산이 아주 간단합니다.

고유벡터 좌표 $(2, 1)$이 $(6, 1)$로, 각 좌표가 자기 고유값만큼 곱해졌을 뿐입니다. 방향끼리 섞이지 않습니다. 이것을 세 단계로 나누면 다음과 같습니다.

1. 좌표 바꾸기. 표준 좌표 $(3, 1)$을 고유벡터 좌표 $(2, 1)$로 바꿉니다. 10교시에서 본 대로, 기저 행렬 $P$(열 = 고유벡터)의 역행렬 $P^{-1}$을 곱하는 일입니다.
2. 축마다 늘이기. 고유벡터 좌표를 각각 $\lambda_1$배, $\lambda_2$배 합니다. 대각선에 고유값만 있는 대각행렬 $D$를 곱하는 일입니다.
3. 좌표 되돌리기. 고유벡터 좌표를 다시 표준 좌표로 바꿉니다. $P$를 곱하는 일입니다.

6교시에서 배운 대로 행렬곱은 오른쪽부터 읽으니, 세 단계를 이어 쓰면 이렇습니다.

행렬을 이렇게 쓰는 것을 대각화(diagonalization) 라고 합니다. "적당한 좌표계에서 보면 이 행렬은 대각행렬이다"라는 뜻입니다.

12-2. 거듭제곱이 쉬워진다 (6분)

$A^2$를 대각화로 써 보면 가운데에서 $P^{-1}P = I$가 만나 사라집니다.

몇 번을 곱해도 마찬가지라 $A^n = PD^nP^{-1}$입니다. 그리고 대각행렬의 거듭제곱은 대각선 성분만 각각 거듭제곱하면 됩니다.

말로 하면 "좌표를 고유벡터로 바꾸고, 각 좌표에 $\lambda^n$을 곱하고, 되돌린다"입니다. 12-1의 $\mathbf{v} = 2\,\mathbf{e}_1 + 1\,\mathbf{e}_2$에 $A$를 10번 곱하면 이렇게 됩니다.

행렬곱 10번 대신 $3^{10}$ 한 번이면 됩니다. 더 중요한 것은 결과가 어디로 가는지 보인다는 점입니다. $3^n$은 빠르게 커지고 $1^n$은 그대로이니, $n$이 클수록 $\mathbf{e}_1$ 쪽 항이 압도합니다. 그래서 $A^n\mathbf{v}$의 방향은 점점 $(1, 1)$에 붙습니다.

일반적으로 정리하면 다음과 같습니다.

모든 행렬이 대각화되는 것은 아닙니다. 11교시의 밀기 행렬처럼 고유벡터 직선이 하나뿐이면 기저를 만들 수 없어 $P^{-1}$이 없습니다. 회전 행렬처럼 고유값이 복소수인 경우는 복소수 세계에서는 대각화되지만, 이 강의에서는 "회전 성분이 섞여 있어 돌면서 커지거나 줄어든다"는 그림으로만 읽겠습니다.

아래 실습대는 시작 벡터 x₀(청록 손잡이)에서 출발해 x₀, Ax₀, A²x₀, …, Aⁿx₀를 차례로 점으로 찍고 선으로 이어 줍니다. 옅은 점에서 진한 점으로 갈수록 나중 단계이고, 진홍 화살표가 마지막 xₙ입니다. 노랑 점선은 절댓값이 가장 큰 고유값의 고유벡터 방향입니다. 슬라이더 n을 0에서 20까지 움직이며 점들이 노랑 점선으로 끌려가는지 보세요.

'날씨 마르코프'와 '피보나치'는 숫자가 작아서, 그림을 3.5배로 키워 그렸다는 안내가 함께 나옵니다. 오른쪽과 아래 NumPy 칸의 숫자는 실제 값입니다.

12-3. 마르코프 연쇄: 내일 날씨는? (7분)

날씨가 "맑음"과 "흐림" 두 가지뿐이고, 다음 규칙을 따른다고 합시다.

오늘의 상태를 (맑을 확률, 흐릴 확률)이라는 확률 벡터 $\mathbf{x}$로 쓰면, 내일의 상태는 행렬 하나를 곱해 얻습니다.

행렬의 1열은 "오늘 맑음"에서 출발해 내일 어디로 가는가, 2열은 "오늘 흐림"에서 출발해 어디로 가는가입니다(5교시: 열 = 기저가 도착한 곳). 그래서 열마다 합이 1입니다. 오늘 확실히 맑으면 $\mathbf{x}_0 = (1, 0)$이고, 내일은 $A\mathbf{x}_0 = (0.9, 0.1)$, 모레는 $(0.86, 0.14)$입니다. 이렇게 "다음 상태가 지금 상태에만 달린" 과정을 마르코프 연쇄(Markov chain) 라고 합니다.

고유값을 구해 봅시다. 대각합은 $1.4$, 행렬식은 $0.45 - 0.05 = 0.4$입니다.

12-2의 규칙대로 읽으면, $0.4^n$ 방향은 금방 사라지고 고유값 1의 방향만 남습니다. 고유값이 정확히 1이니 그 방향의 벡터는 $A$를 곱해도 변하지 않습니다. 이것을 정상 상태(steady state) 라고 합니다. $(A - I)\mathbf{v} = \mathbf{0}$의 첫 줄은 $-0.1s + 0.5c = 0$, 곧 $s = 5c$입니다. 확률의 합이 1이 되게 맞추면 다음과 같습니다.

실습대에 나온 "오래 지나면 맑음 83.3% : 흐림 16.7% 로 수렴합니다"가 바로 이 값입니다. n = 5에서 이미 $(0.835, 0.165)$로 거의 도착했고, n = 20이면 소수 일곱째 자리까지 같아집니다.

놀라운 점은 출발점과 상관이 없다는 것입니다. 오늘 확실히 흐려도, 반반이어도 오래 지나면 같은 83.3% : 16.7%에 도착합니다. 출발점의 차이는 모두 $0.4^n$ 방향에 실려 있다가 사라지기 때문입니다. 웹 페이지의 중요도를 매기는 PageRank도 "사용자가 링크를 무작위로 따라갈 때 오래 지나면 각 페이지에 머무는 비율"이라는, 같은 구조의 정상 상태입니다.

실습대에서 x₀를 끌어 성분의 합이 1이 아닌 벡터로 만들면, 수렴하는 점도 합에 비례해 달라집니다. 예를 들어 합이 2이면 $(1.667, 0.333)$으로 갑니다. 그래도 노랑 점선 위, 맑음과 흐림의 비율 5 : 1은 그대로입니다.

12-4. 피보나치와 황금비, 그리고 소용돌이 (5분)

피보나치 수열 $1, 1, 2, 3, 5, 8, 13, \dots$은 "앞의 두 수를 더해 다음 수를 만든다"는 규칙입니다. 연속한 두 수를 벡터로 묶으면 이 규칙도 행렬 곱셈 하나가 됩니다.

첫 줄은 $F_{n+1} = F_n + F_{n-1}$, 둘째 줄은 $F_n$을 한 칸 옮겨 적은 것입니다. 이 행렬의 고유값은 $\lambda^2 - \lambda - 1 = 0$의 해입니다.

$\lambda_1$이 바로 황금비입니다. $\lvert\lambda_2\rvert < 1$이라 그 방향은 사라지고, 벡터 $(F_{n+1}, F_n)$은 $\lambda_1$의 고유벡터 방향으로 붙습니다. 그 고유벡터는 $(\lambda_1, 1)$ 방향이라, 두 성분의 비율 $F_{n+1} / F_n$이 황금비 $1.61803\dots$로 다가갑니다. 실습대의 '피보나치'에서 n을 올리면 "비율 x_n[0] / x_n[1]"이 1, 2, 1.5, 1.667, 1.6, 1.625, …로 흔들리며 1.61803에 가까워집니다. $\lambda_2$가 음수라서 남은 오차의 부호가 매번 바뀌기 때문에 위아래로 흔들립니다.

'피보나치'에서는 '방향만 보기'가 처음부터 켜져 있습니다. 성분이 $1.618^n$배로 커져 화면 밖으로 나가 버리니, 매 단계 길이를 1로 맞춰 방향만 보여 주는 것입니다. 끄면 점들이 노랑 점선을 따라 순식간에 화면 밖으로 달아납니다.

'돌며 줄어들기' 행렬 $\begin{bmatrix} 0.8 & -0.5 \\ 0.5 & 0.8 \end{bmatrix}$은 고유값이 복소수 $0.8 \pm 0.5i$입니다. 11교시에서 본 대로 복소수 고유값은 회전 성분이 있다는 신호입니다. 복소수의 크기 $\sqrt{0.8^2 + 0.5^2} \approx 0.943$이 1보다 작으니 한 번 곱할 때마다 길이가 약 94%로 줄어듭니다. 그래서 점들이 원점을 돌며 소용돌이처럼 빨려 들어갑니다. 실수 고유벡터가 없으니 노랑 점선도 나오지 않습니다.

12-5. NumPy로 대각화와 거듭제곱 (11분)

실습대의 NumPy 칸 읽기

실습대의 처음 화면('날씨 마르코프', n = 5)에서 NumPy 칸에는 이런 코드가 나옵니다.

>>> import numpy as np
>>> A = np.array([[0.9, 0.5], [0.1, 0.5]])
>>> x0 = np.array([1, 0])
>>> np.linalg.matrix_power(A, 5)
array([[0.83504, 0.8248 ],
       [0.16496, 0.1752 ]])
>>> np.linalg.matrix_power(A, 5) @ x0
array([0.83504, 0.16496])
>>> np.round(np.linalg.eigvals(A), 6)
array([1. +0.j, 0.4+0.j])

np.linalg.matrix_power(A, n)은 A를 n번 곱한 행렬 Aⁿ입니다. A @ A @ A @ A @ A와 같은 값인데, 제곱을 거듭하는 방식으로 곱셈 횟수를 줄여 계산합니다. n = 0이면 항등행렬이 나옵니다. 거기에 @ x0을 붙이면 5일 뒤의 날씨 확률 x₅입니다. 한 번에 한 단계씩 진행해도 같은 답입니다.

>>> A @ x0
array([0.9, 0.1])
>>> A @ (A @ x0)
array([0.86, 0.14])
>>> np.linalg.matrix_power(A, 2) @ x0
array([0.86, 0.14])
>>> np.linalg.matrix_power(A, 0)
array([[1., 0.],
       [0., 1.]])

np.linalg.eigvals(A)는 11교시에서 본 대로 고유값만 돌려주고, 칸에서는 그것을 np.round(…, 6)으로 감싸 소수 여섯째 자리까지 반올림해 보여 줍니다. 그래서 '피보나치'를 고르면 array([ 1.618034+0.j, -0.618034+0.j])처럼 여섯 자리까지만 나옵니다. 1. +0.j의 빈칸은 0.4+0.j와 자리를 맞추려고 넣은 것일 뿐입니다.

n을 20, 50으로 올리면 정상 상태 (5/6, 1/6)에 붙습니다. 흐린 날에서 출발해도 같은 곳에 도착합니다. 합이 2인 벡터에서 출발하면 같은 비율로 두 배인 곳에 도착합니다.

>>> np.linalg.matrix_power(A, 20) @ x0
array([0.83333334, 0.16666666])
>>> np.linalg.matrix_power(A, 50) @ x0
array([0.83333333, 0.16666667])
>>> np.linalg.matrix_power(A, 50) @ np.array([0, 1])
array([0.83333333, 0.16666667])
>>> np.linalg.matrix_power(A, 30) @ np.array([2, 0])
array([1.66666667, 0.33333333])

np.diag: 대각행렬 만들기와 대각선 꺼내기

np.diag는 넘기는 것에 따라 두 가지 일을 합니다. 1차원 배열(숫자 목록)을 넘기면 그 숫자를 대각선에 놓은 대각행렬을 만들고, 2차원 배열(행렬)을 넘기면 대각선만 꺼내 줍니다.

>>> np.diag([3, 1])
array([[3, 0],
       [0, 1]])
>>> np.diag(np.array([[1, 2], [3, 4]]))
array([1, 4])

P D P⁻¹로 되돌아오는지 확인하기

12-1의 대각화를 그대로 적어 확인해 봅니다. P의 열은 고유벡터 (1, 1)과 (1, −1)입니다.

>>> S = np.array([[2, 1], [1, 2]])
>>> P = np.array([[1, 1], [1, -1]])
>>> D = np.diag([3, 1])
>>> P @ D @ np.linalg.inv(P)
array([[2., 1.],
       [1., 2.]])
>>> S @ np.array([3, 1])
array([7, 5])

np.linalg.inv가 실수 계산을 하므로 결과가 2.처럼 실수로 나왔습니다. 이번에는 NumPy가 구한 고유벡터로 같은 일을 해 봅니다. S는 대칭이니 11교시의 eigh를 씁니다.

>>> vals, vecs = np.linalg.eigh(S)
>>> P = vecs
>>> D = np.diag(vals)
>>> np.allclose(P @ D @ np.linalg.inv(P), S)
True

이제 12-2의 Aⁿ = P Dⁿ P⁻¹을 씁니다. 대각선만 거듭제곱하면 되니 vals ** 10이면 됩니다. 배열의 **는 원소마다 거듭제곱하는 연산입니다(2교시의 원소별 연산).

>>> P @ np.diag(vals ** 10) @ np.linalg.inv(P)
array([[29525., 29524.],
       [29524., 29525.]])
>>> np.linalg.matrix_power(S, 10)
array([[29525, 29524],
       [29524, 29525]])
>>> np.linalg.matrix_power(S, 10) @ np.array([3, 1])
array([118099, 118097])

12-2에서 손으로 구한 (118099, 118097)이 그대로 나왔습니다.

덤으로 하나 더 보겠습니다. 대칭행렬의 고유벡터는 서로 수직이고 길이가 1이라서, P의 역행렬이 전치 P.T와 같습니다. 이런 행렬을 직교행렬이라고 부르며, 다음 교시의 주인공입니다.

>>> np.allclose(np.linalg.inv(P), P.T)
True

정상 상태를 고유벡터로 바로 구하기

마르코프 연쇄의 정상 상태는 고유값 1의 고유벡터입니다. np.argmax는 가장 큰 값이 몇 번째 칸에 있는지 알려 줍니다. 그 번호의 열을 꺼내고(11교시: vecs[:, k]), 확률이 되도록 성분의 합으로 나눕니다.

>>> vals, vecs = np.linalg.eig(A)
>>> k = np.argmax(vals.real)
>>> k
np.int64(0)
>>> s = vecs[:, k].real
>>> s
array([0.98058068, 0.19611614])
>>> s / s.sum()
array([0.83333333, 0.16666667])
>>> A @ (s / s.sum())
array([0.83333333, 0.16666667])

NumPy가 준 고유벡터는 길이 1로 맞춘 (0.98, 0.196)이라 확률처럼 보이지 않지만, 비율은 5 : 1입니다. 합으로 나누면 (5/6, 1/6)이 되고, A를 곱해도 변하지 않습니다. 이것이 정상 상태입니다.

열마다 합이 1인지는 sum(axis=0)으로 확인합니다. axis=0은 "0번 축(행 방향)을 따라 위에서 아래로 더하라", 곧 열마다 더하라는 뜻입니다. 축은 15교시에서 다시 정리합니다.

>>> A.sum(axis=0)
array([1., 1.])

피보나치

>>> F = np.array([[1, 1], [1, 0]])
>>> np.linalg.matrix_power(F, 10)
array([[89, 55],
       [55, 34]])
>>> x = np.array([1, 0])
>>> for n in range(1, 8):
...     x = F @ x
...     print(n, x, x[0] / x[1])
1 [1 1] 1.0
2 [2 1] 2.0
3 [3 2] 1.5
4 [5 3] 1.6666666666666667
5 [8 5] 1.6
6 [13  8] 1.625
7 [21 13] 1.6153846153846154
>>> np.linalg.eigvals(F)
array([ 1.61803399+0.j, -0.61803399+0.j])
>>> (1 + np.sqrt(5)) / 2
np.float64(1.618033988749895)

F¹⁰의 성분이 피보나치 수 89, 55, 34입니다. for문의 비율이 위아래로 흔들리며 황금비로 다가가는 것도 12-4에서 말한 그대로입니다.

12-6. 거듭제곱법: 반복만으로 고유값 찾기 (5분)

12-2의 "가장 큰 고유값의 방향이 결국 이긴다"를 거꾸로 이용하면, 특성방정식을 풀지 않고도 가장 큰 고유값을 찾을 수 있습니다. 아무 벡터에서 출발해 A를 계속 곱하면 된다는 것이 거듭제곱법(power iteration) 입니다.

다만 그냥 곱하기만 하면 피보나치처럼 숫자가 끝없이 커지거나 0으로 줄어듭니다. 그래서 매번 길이를 1로 맞춥니다(실습대의 '방향만 보기'와 같은 일). 한 번 곱할 때 길이가 몇 배가 되었는지가 고유값의 근사입니다.

>>> x = np.array([1.0, 0.0])
>>> for k in range(1, 11):
...     y = F @ x
...     lam = np.linalg.norm(y) / np.linalg.norm(x)
...     x = y / np.linalg.norm(y)
...     print(k, round(float(lam), 6))
1 1.414214
2 1.581139
3 1.612452
4 1.617215
5 1.617914
6 1.618017
7 1.618031
8 1.618034
9 1.618034
10 1.618034

round(float(lam), 6)은 줄마다 np.float64(…)가 붙지 않게 파이썬 숫자로 바꾸고 여섯째 자리까지 보여 주려고 쓴 것입니다. 여덟 번 만에 황금비가 여섯 자리까지 맞았습니다. 이때 x는 고유벡터의 근사입니다.

>>> x
array([0.85067229, 0.52569636])
>>> np.round(x[0] / x[1], 4)
np.float64(1.6182)

얼마나 빨리 맞아 가는지는 두 고유값의 크기 비율 |λ₂| / |λ₁|이 정합니다. 피보나치는 0.618 / 1.618 ≈ 0.38이라 한 번에 오차가 0.38배쯤으로 줄어 빠르게 맞습니다. 두 고유값의 크기가 비슷하면 아주 느려집니다.

실제로 np.linalg.eig는 이보다 훨씬 정교한 방법(QR 알고리즘)으로 모든 고유값을 한꺼번에 구합니다. 그래도 거듭제곱법은 행렬이 수백만×수백만처럼 너무 커서 eig를 쓸 수 없을 때 지금도 쓰이는 방법이고, PageRank도 원래 이렇게 계산했습니다. 이런 "반복해서 답에 다가가는 계산"은 이 강의 다음의 수치해석에서 본격적으로 다룹니다.

12-7. 자주 하는 실수 (3분)

A ** n은 행렬 거듭제곱이 아닙니다. 배열의 **는 원소마다 거듭제곱입니다. 6교시에서 A * B가 행렬곱이 아니었던 것과 같은 이야기입니다.

>>> S ** 3
array([[8, 1],
       [1, 8]])
>>> np.linalg.matrix_power(S, 3)
array([[14, 13],
       [13, 14]])

정수가 넘쳐도 알려 주지 않습니다. NumPy의 int64는 약 9.2 × 10¹⁸까지만 담습니다. 피보나치 수 F₁₀₁은 그보다 훨씬 큰데, 오류도 경고도 없이 엉뚱한 숫자가 나옵니다.

>>> np.linalg.matrix_power(F, 90) @ np.array([1, 0])
array([4660046610375530309, 2880067194370816120])
>>> np.linalg.matrix_power(F, 100) @ np.array([1, 0])
array([1298777728820984005, 3736710778780434371])

90제곱까지는 맞지만 100제곱의 결과는 둘째 수가 첫째 수보다 크게 나와 피보나치 수가 아닙니다. 아주 큰 정수가 필요하면 dtype=object로 파이썬 정수를 담게 할 수 있습니다(느리지만 정확).

>>> Fo = np.array([[1, 1], [1, 0]], dtype=object)
>>> np.linalg.matrix_power(Fo, 100) @ np.array([1, 0])
array([573147844013817084101, 354224848179261915075], dtype=object)

대각화되지 않는 행렬. 밀기 행렬에서 eig가 준 P는 두 열이 같은 직선이라(11교시) 행렬식이 0에 가깝습니다. 그 역행렬을 구하면 오류 대신 10¹⁵ 같은 거대한 숫자가 나옵니다. 이런 결과를 믿으면 안 됩니다.

>>> vals, vecs = np.linalg.eig(np.array([[1, 1], [0, 1]]))
>>> np.linalg.det(vecs)
np.complex128(2.2204460492503185e-16+0j)

연습 문제 (5분)

1. 12-1의 대각화를 이용해 $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}^5$를 손으로 구하세요. (힌트: $3^5 = 243$. $A^n$의 성분은 $\tfrac{3^n + 1}{2}$과 $\tfrac{3^n - 1}{2}$입니다.)
2. 앱 사용자가 "활성"과 "휴면" 두 상태라고 합시다. 활성 사용자는 다음 달에 80%가 활성으로 남고, 휴면 사용자는 30%가 돌아옵니다. 행렬 $T$를 세우고, 오래 지나면 활성 사용자 비율이 얼마가 되는지 구하세요.
3. 12-6의 거듭제곱법을 $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$에 30번 돌리면 어떤 고유값과 고유벡터에 다가갈지 예상하고 확인하세요.

풀이

1번: $A^5 = P\,\mathrm{diag}(243, 1)\,P^{-1}$이고, 성분은 $(243 + 1)/2 = 122$, $(243 - 1)/2 = 121$입니다.

2번: 열이 "이번 달 상태에서 출발해 다음 달에 어디로 가는가"이니 $T = \begin{bmatrix} 0.8 & 0.3 \\ 0.2 & 0.7 \end{bmatrix}$입니다. 고유값은 대각합 1.5, 행렬식 $0.56 - 0.06 = 0.5$에서 1과 0.5입니다. 고유값 1의 고유벡터는 $-0.2a + 0.3b = 0$에서 $a = 1.5b$, 합이 1이 되게 맞추면 $(0.6, 0.4)$입니다. 오래 지나면 활성 60%, 휴면 40%입니다.

3번: 가장 큰 고유값 3과 그 고유벡터 $(1, 1)$ 방향, 곧 길이 1로 맞춘 $(0.707, 0.707)$에 다가갑니다.

>>> np.linalg.matrix_power(S, 5)
array([[122, 121],
       [121, 122]])
>>> T = np.array([[0.8, 0.3], [0.2, 0.7]])
>>> np.linalg.eigvals(T)
array([1. +0.j, 0.5+0.j])
>>> vals, vecs = np.linalg.eig(T)
>>> s = vecs[:, np.argmax(vals.real)].real
>>> s / s.sum()
array([0.6, 0.4])
>>> np.linalg.matrix_power(T, 30) @ np.array([0, 1])
array([0.6, 0.4])
>>> x = np.array([1.0, 0.0])
>>> for _ in range(30):
...     x = S @ x
...     x = x / np.linalg.norm(x)
>>> x
array([0.70710678, 0.70710678])
>>> x @ S @ x
np.float64(3.0000000000000004)

2번의 마지막 줄은 모든 사용자가 휴면인 상태에서 출발해도 30달 뒤에는 (0.6, 0.4)에 도착한다는 확인입니다. 3번의 x @ S @ x는 길이 1인 x에 대해 "S가 x를 몇 배 하는가"를 재는 식으로, 고유값 3에 부동소수점 찌꺼기가 조금 붙어 나왔습니다. 정확히 비교하려면 np.isclose(x @ S @ x, 3)을 씁니다.

다음 교시에는 서로 직교하는 벡터가 왜 그렇게 편한지 보고, 데이터 점들에 가장 잘 맞는 직선을 찾는 최소제곱법을 배웁니다. 정확한 해가 없는 연립방정식에서 "가장 가까운 답"을 찾는 방법입니다.