ReflexKR

그림으로 시작하는 선형대수: NumPy와 함께

강의 개요

숫자 묶음으로 세상 적기: 벡터와 NumPy 배열

벡터 = 화살표 = 숫자 목록: 덧셈과 스칼라배

길이와 각도, 그리고 내적

일차결합과 생성(span)

행렬은 공간을 옮기는 규칙이다

행렬곱 = 변환을 이어 하기

행렬식과 역행렬

연립방정식 Ax = b

열공간·영공간·계수

기저 바꾸기

고유값과 고유벡터

대각화와 거듭제곱

직교·사영·최소제곱

SVD: 돌리고, 늘이고, 돌리기

PCA, 그리고 다음 걸음

그림으로 시작하는 선형대수: NumPy와 함께

14교시

강의 시간 10:50–11:40

SVD: 돌리고, 늘이고, 돌리기

11교시와 12교시에서 고유벡터는 아주 쓸모 있는 도구였습니다. 고유벡터 기저로 바꾸면 행렬이 "그냥 늘이기"가 되고, 거듭제곱도 쉬워졌습니다. 그런데 아쉬운 점이 있었습니다. 회전 행렬에는 실수 고유벡터가 없었고, 밀기 행렬에는 하나뿐이었습니다. 게다가 고유벡터는 정사각행렬에만 정의됩니다. 사용자 300명 × 영화 1000편 같은 데이터 표에는 쓸 수가 없습니다.

이번 교시의 특이값 분해(SVD, Singular Value Decomposition) 는 이 아쉬움을 모두 풀어 줍니다. 어떤 행렬이든, 정사각이든 아니든, 반드시 돌리기 → 축 방향으로 늘이기 → 다시 돌리기 세 단계로 쪼갤 수 있습니다. 예외가 없다는 점이 SVD를 선형대수에서 가장 많이 쓰이는 분해로 만들었습니다.

쓰이는 곳도 넓습니다. 사진을 적은 숫자로 줄이는 압축, 잡음 걷어 내기, 추천 시스템의 "취향 축" 찾기, 9교시에서 도구로만 썼던 영공간 구하기, 그리고 다음 교시의 PCA까지 모두 SVD 위에 서 있습니다.

14-1. 원은 늘 타원이 된다 (5분)

5교시에서 행렬은 격자 전체를 옮기는 규칙이었습니다. 이번에는 격자 대신 단위원(원점에서 길이 1인 점 전부)을 옮겨 봅니다.

어떤 2×2 행렬로 옮기든, 단위원은 늘 타원이 됩니다. 원이 그대로 원으로 남기도 하고(회전), 아주 납작해져 선분이 되기도 하지만(det = 0), 찌그러진 감자 모양이 되는 일은 없습니다. 격자선이 평행·등간격을 유지하는 변환이라 그렇습니다.

타원에는 가장 긴 반지름과 가장 짧은 반지름이 있고, 두 방향은 서로 수직입니다. SVD는 이 두 반지름을 찾는 이야기입니다.

타원의 긴 반지름 길이 = 첫째 특이값 σ₁
타원의 짧은 반지름 길이 = 둘째 특이값 σ₂
원 위에서 그 두 반지름으로 가는 점(서로 수직인 단위벡터) = 오른쪽 특이벡터 v₁, v₂
타원의 두 반지름 방향(역시 서로 수직) = 왼쪽 특이벡터 u₁, u₂

원 위의 수직인 두 화살표가 옮겨진 뒤에도 여전히 수직이라는 점이 핵심입니다. 보통의 행렬은 수직을 지켜 주지 않습니다(밀기 행렬로 î, ĵ를 옮겨 보면 비스듬해집니다). 그런데 원 위에는 수직이 지켜지는 특별한 한 쌍이 언제나 있습니다. 그 한 쌍이 v₁, v₂입니다.

14-2. 돌리고, 늘이고, 돌리기 (10분)

아래 실습대의 파랑(인디고) 도형은 처음에 단위원입니다. 오른쪽의 단추를 차례로 눌러 보세요.

처음 (단위원): 원과 노랑 화살표 두 개. 진한 노랑이 v₁, 옅은 노랑이 v₂입니다.
① Vᵀ: 돌리기: 원을 돌려 v₁을 x축, v₂를 y축 위에 올려놓습니다. 원은 돌려도 원이라 모양은 그대로입니다.
② Σ: 축 방향으로 늘이기: x축 방향으로 σ₁배, y축 방향으로 σ₂배 늘입니다. 이제 원이 축에 나란한 타원이 됩니다.
③ U: 다시 돌리기 = A: 타원을 통째로 돌려 제자리에 놓습니다. 이 결과가 행렬 A로 원을 옮긴 것과 정확히 같습니다.

노랑 화살표는 오른쪽 특이벡터가 가는 곳을 따라갑니다. 마지막 단계에서 노랑 화살표 끝이 타원의 긴 반지름과 짧은 반지름 끝에 정확히 닿는지 보세요.

A20.511.5

특이값: σ₁ = 2.558, σ₂ = 0.977

U0.850.53-0.530.85Vᵀ0.77-0.640.640.77

지금 보이는 것: ③ U: 다시 돌리기 = A

A의 열 끌기

NumPy로 확인하기

실습대에서 값을 정할 때마다, 서버의 NumPy가 같은 숫자로 계산한 결과입니다.

>>> A = np.array([[2, 1], [0.5, 1.5]])
>>> U, S, Vt = np.linalg.svd(A)
>>> np.round(U, 6)
array([[-0.850651, -0.525731],
       [-0.525731,  0.850651]])
>>> np.round(S, 6)
array([2.558336, 0.977198])
>>> np.round(Vt, 6)
array([[-0.767752, -0.640747],
       [-0.640747,  0.767752]])
>>> np.round(U @ np.diag(S) @ Vt, 6)
array([[2. , 1. ],
       [0.5, 1.5]])

'한 번에 재생'을 누른 뒤 ①·②·③ 단추로 한 단계씩 되짚어 보세요. ③ 단계에서 청록(1열)과 진홍(2열) 끝점을 끌어 청록을 (2, 1), 진홍을 (1, 0.5)에 놓으면 두 열이 한 직선에 놓여 σ₂가 0이 되고, 타원이 선분으로 납작해집니다.

실습대에서 해 볼 만한 것들입니다. A의 열 끝점 손잡이는 ③ 단계(= A)일 때만 보입니다.

1열을 (0, 1), 2열을 (−1, 0)에 놓으면 90° 회전입니다. 원이 원 그대로라 σ₁ = σ₂ = 1이 됩니다. 늘이기가 전혀 없는 행렬이라는 뜻입니다.
1열을 (1, 0), 2열을 (1, 1)에 놓으면 6교시의 밀기입니다. 11교시에서 고유벡터가 하나뿐이라 애먹었던 행렬인데, SVD는 문제없이 σ₁ ≈ 1.618, σ₂ ≈ 0.618을 줍니다. 12교시의 황금비가 여기서 다시 나옵니다.
1열과 2열의 자리를 바꿔 det를 음수로 만들면, ③ 단계에서 타원이 한 번 납작해졌다가 뒤집혀 나옵니다. 이때 U는 순수한 회전이 아니라 뒤집기가 섞인 회전입니다(7교시: det가 음수면 방향이 뒤집힘).

책의 기호 읽기

세 단계를 한 줄로 적은 것이 특이값 분해입니다. 행렬곱은 오른쪽부터 읽으니(6교시), $V^{\mathsf T}$ 가 먼저, $U$ 가 마지막입니다.

A = U\,\Sigma\,V^{\mathsf T} = \begin{bmatrix} | & | \\ \mathbf{u}_1 & \mathbf{u}_2 \\ | & | \end{bmatrix} \begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{bmatrix} \begin{bmatrix} - & \mathbf{v}_1^{\mathsf T} & - \\ - & \mathbf{v}_2^{\mathsf T} & - \end{bmatrix}

가장 오른쪽의 $V^{\mathsf T}$ 가 ① 돌리기, 가운데의 $\Sigma$ 가 ② 늘이기, 왼쪽의 $U$ 가 ③ 다시 돌리기입니다. U의 열은 $\mathbf{u}_1, \mathbf{u}_2$ 이고, $V^{\mathsf T}$ 의 행은 $\mathbf{v}_1, \mathbf{v}_2$ 를 눕혀 쓴 것입니다.

책에 적힌 모양	읽는 법	뜻	파이썬에서는
$A = U\Sigma V^{\mathsf T}$	A는 U 시그마 V 전치	특이값 분해	`U, S, Vt = np.linalg.svd(A)`
$\sigma_1 \ge \sigma_2 \ge 0$	시그마 일, 시그마 이	특이값. 큰 것부터, 음수는 없음	`S[0]`, `S[1]`
$\Sigma$	대문자 시그마	대각선에 특이값을 놓은 행렬	`np.diag(S)`
$\mathbf{u}_i$	왼쪽 특이벡터	U의 i번째 열, 타원의 반지름 방향	`U[:, i]`
$\mathbf{v}_i$	오른쪽 특이벡터	V의 i번째 열, 원 위의 출발점	`Vt[i]` (Vᵀ의 행!)
$U^{\mathsf T}U = I$	U는 직교행렬	열끼리 수직이고 길이가 1	`np.allclose(U.T @ U, np.eye(2))`

세 단계를 벡터 하나로 따라가면 SVD의 가장 중요한 식이 나옵니다. $\mathbf{v}_1$ 은 ①에서 x축 위의 $(1, 0)$ 이 되고, ②에서 $(\sigma_1, 0)$ 이 되고, ③에서 $\sigma_1\mathbf{u}_1$ 이 됩니다.

A\mathbf{v}_1 = \sigma_1\mathbf{u}_1, \qquad A\mathbf{v}_2 = \sigma_2\mathbf{u}_2

11교시의 $A\mathbf{v} = \lambda\mathbf{v}$ 와 비교해 보세요. 고유벡터는 "같은 방향으로" 돌아와야 했지만, 특이벡터는 다른 방향( $\mathbf{u}$ )으로 가도 됩니다. 대신 출발하는 두 방향도 수직, 도착하는 두 방향도 수직이어야 합니다. 조건을 이렇게 바꾸었더니 모든 행렬에 답이 생겼습니다.

U와 V를 "회전"이라고 불렀지만, 정확한 이름은 직교행렬입니다. 13교시의 QR 분해에서 본 Q처럼 열끼리 수직이고 길이가 1인 행렬로, 회전이거나 뒤집기가 섞인 회전입니다. 길이와 각도를 바꾸지 않는다는 점은 같습니다.

14-3. 특이값과 고유값, 행렬식, 계수 (6분)

특이값은 앞에서 배운 개념들과 촘촘히 이어져 있습니다.

고유값과의 관계. $A^{\mathsf T}A$ 를 SVD로 풀어 쓰면 가운데의 $U^{\mathsf T}U$ 가 $I$ 가 되어 사라집니다.

A^{\mathsf T}A = (V\Sigma U^{\mathsf T})(U\Sigma V^{\mathsf T}) = V\,\Sigma^2\,V^{\mathsf T}

12교시의 대각화 $PDP^{-1}$ 와 같은 모양입니다( $V^{\mathsf T}$ 가 곧 $V^{-1}$ ). 그러니 $A^{\mathsf T}A$ 의 고유벡터가 $\mathbf{v}_i$ 이고, 고유값이 $\sigma_i^2$ 입니다. 실습대도 이 방법으로 σ와 V를 계산합니다. $A^{\mathsf T}A$ 는 늘 대칭행렬이라 11교시에서 본 대로 서로 수직인 고유벡터가 있고, 고유값도 0 이상입니다. 그래서 어떤 행렬이든 SVD가 있습니다.

행렬식과의 관계. 7교시에서 det는 넓이 배율이었습니다. 돌리기는 넓이를 바꾸지 않고, 늘이기만 넓이를 $\sigma_1\sigma_2$ 배로 만듭니다. 그래서

|\det A| = \sigma_1\,\sigma_2

입니다. 부호(뒤집혔는지)는 U와 V가 가져갑니다.

계수와의 관계. 9교시의 계수는 "0이 아닌 특이값의 개수"입니다. $\sigma_2 = 0$ 이면 ② 단계에서 y축 방향이 통째로 눌려 사라지고, 원이 선분이 됩니다. 이때 $A\mathbf{v}_2 = 0\cdot\mathbf{u}_2 = \mathbf{0}$ 이니 $\mathbf{v}_2$ 가 바로 영공간의 방향입니다. 9교시에서 Vt[-1]로 영공간을 구했던 이유가 이것입니다.

14-4. NumPy로 SVD 하기 (12분)

실습대 아래 "NumPy로 확인하기" 칸에 나오는 코드를 처음부터 따라가 봅니다. 실습대의 처음 행렬 A로 시작합니다.

>>> import numpy as np
>>> A = np.array([[2, 1], [0.5, 1.5]])
>>> U, S, Vt = np.linalg.svd(A)
>>> S
array([2.55833637, 0.97719754])

np.linalg.svd는 결과를 세 개 돌려줍니다. U, S, Vt = …처럼 쉼표로 이름 셋을 늘어놓으면 차례로 하나씩 받습니다(11교시의 vals, vecs = np.linalg.eig(A)와 같은 문법입니다). 이름에 주의하세요.

S는 행렬 Σ가 아니라 특이값만 담은 1차원 배열입니다. 큰 것부터 줄 서 있습니다.
Vt는 V가 아니라 이미 전치된 Vᵀ입니다. 그래서 오른쪽 특이벡터 v₁은 Vt[0], 즉 행으로 꺼냅니다. U의 왼쪽 특이벡터는 U[:, 0], 즉 열입니다.

>>> U
array([[-0.85065081, -0.52573111],
       [-0.52573111,  0.85065081]])
>>> Vt
array([[-0.76775173, -0.64074744],
       [-0.64074744,  0.76775173]])

실습대 아래 NumPy 칸에서는 소수점 아래가 너무 길어지지 않도록 np.round(U, 6), np.round(S, 6), np.round(Vt, 6), np.round(U @ np.diag(S) @ Vt, 6)처럼 줄마다 반올림해서 보여 줍니다(그래서 칸의 S는 array([2.558336, 0.977198])입니다). np.round(배열, 자릿수)는 배열 전체를 그 자릿수로 반올림한 새 배열을 줍니다. 값을 바꾸는 것이 아니라 보기 좋게 보여 주는 용도입니다.

>>> np.round(U, 3)
array([[-0.851, -0.526],
       [-0.526,  0.851]])

세 조각을 다시 곱하면 A가 돌아와야 합니다. S가 1차원 배열이라 먼저 np.diag(S)로 대각행렬 Σ를 만듭니다. np.diag는 1차원 배열을 받으면 그 숫자들을 대각선에 놓은 정사각행렬을 만들어 줍니다.

>>> np.diag(S)
array([[2.55833637, 0.        ],
       [0.        , 0.97719754]])
>>> U @ np.diag(S) @ Vt
array([[2. , 1. ],
       [0.5, 1.5]])
>>> np.allclose(U @ np.diag(S) @ Vt, A)
True

이번에는 깔끔하게 나왔지만, 부동소수점 계산이라 1.9999999999999998 같은 찌꺼기가 섞일 때가 많습니다. 그래서 비교는 3교시부터 써 온 np.allclose로 합니다.

A @ Vt[0]과 S[0] * U[:, 0]을 비교하면 14-2의 식 Av₁ = σ₁u₁을 직접 확인할 수 있습니다.

>>> A @ Vt[0]
array([-2.1762509 , -1.34499702])
>>> S[0] * U[:, 0]
array([-2.1762509 , -1.34499702])
>>> np.allclose(A @ Vt[1], S[1] * U[:, 1])
True

U와 V가 직교행렬인지도 확인해 봅니다. 13교시에서 Q를 확인했던 방법 그대로입니다.

>>> np.allclose(U.T @ U, np.eye(2))
True
>>> np.allclose(Vt @ Vt.T, np.eye(2))
True

실습대와 부호가 다른 이유

실습대의 U는 대략 [[0.85, -0.53], [0.53, 0.85]]이고, NumPy의 U는 [[-0.85, -0.53], [-0.53, 0.85]]입니다. 첫째 열의 부호가 반대입니다. Vᵀ도 첫째 행의 부호가 반대입니다. 둘 중 하나가 틀린 것이 아닙니다.

Av₁ = σ₁u₁에서 v₁과 u₁을 함께 뒤집으면 A(−v₁) = σ₁(−u₁)이 되어 식이 그대로 성립합니다. 즉 u와 v의 짝을 함께 뒤집어도 같은 분해입니다. 실습대는 V가 순수한 회전이 되도록 부호를 골랐고, NumPy(내부의 LAPACK)는 계산 과정에서 나온 부호를 그대로 줍니다. 짝을 함께 뒤집어 실습대와 맞춰 보겠습니다.

>>> U2 = U.copy()
>>> Vt2 = Vt.copy()
>>> U2[:, 0] *= -1
>>> Vt2[0] *= -1
>>> np.round(U2, 2)
array([[ 0.85, -0.53],
       [ 0.53,  0.85]])
>>> np.round(Vt2, 2)
array([[ 0.77,  0.64],
       [-0.64,  0.77]])
>>> np.allclose(U2 @ np.diag(S) @ Vt2, A)
True

.copy()는 배열을 새로 복사합니다. 복사하지 않고 U2 = U라고 쓰면 같은 배열에 이름만 하나 더 붙어서, U2를 고칠 때 U도 함께 바뀝니다. U2[:, 0] *= -1은 "0번 열 전체에 −1을 곱해 그 자리에 다시 넣기"입니다.

한 쪽만 뒤집으면 분해가 깨집니다. 두 SVD 결과를 비교할 때는 특이값 S를 비교하거나, 벡터는 부호를 무시하고(np.abs) 비교하세요.

고유값·행렬식과의 관계를 숫자로

>>> vals, vecs = np.linalg.eigh(A.T @ A)
>>> vals
array([0.95491503, 6.54508497])
>>> S ** 2
array([6.54508497, 0.95491503])
>>> np.linalg.det(A)
np.float64(2.5)
>>> S.prod()
np.float64(2.5000000000000004)

A.T @ A의 고유값이 σ²과 같습니다. eigh는 작은 것부터, svd는 큰 것부터 줄을 세운다는 점만 다릅니다. S.prod()는 배열의 숫자를 모두 곱한 값이고, det 2.5와 마지막 자리 찌꺼기만 다릅니다.

밀기 행렬과 9교시의 납작한 행렬도 넣어 봅니다. 특이값만 필요할 때는 compute_uv=False를 주면 U와 V를 계산하지 않고 S만 돌려줍니다.

>>> np.linalg.svd(np.array([[1, 1], [0, 1]]), compute_uv=False)
array([1.61803399, 0.61803399])
>>> F = np.array([[2, 1], [1, 0.5]])
>>> U, S, Vt = np.linalg.svd(F)
>>> np.round(S, 6)
array([2.5, 0. ])
>>> np.linalg.matrix_rank(F)
np.int64(1)
>>> np.allclose(F @ Vt[-1], 0)
True

납작한 행렬 F는 둘째 특이값이 0이고, 계수가 1이며, Vt[-1](v₂)을 넣으면 영벡터가 나옵니다. 둘째 특이값은 계산 환경에 따라 1e-17처럼 아주 작은 수로 나올 수 있어 반올림해서 보였습니다. np.linalg.matrix_rank도 내부에서 SVD를 하고, 아주 작은 특이값은 0으로 쳐서 셉니다.

정사각이 아닌 행렬

SVD는 정사각이 아닌 행렬에도 됩니다. 3×2 행렬을 넣어 보면 세 조각의 모양(shape)이 다릅니다.

>>> B = np.array([[1, 0], [1, 1], [0, 1]])
>>> U, S, Vt = np.linalg.svd(B)
>>> U.shape, S.shape, Vt.shape
((3, 3), (2,), (2, 2))
>>> U @ np.diag(S) @ Vt
Traceback (most recent call last):
  ...
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)

기본값(full_matrices=True)에서는 U가 3×3으로 꽉 찬 직교행렬로 나옵니다. 그런데 특이값은 2개뿐이라 U의 셋째 열은 곱할 짝이 없고, np.diag(S)(2×2)와 모양이 맞지 않아 오류가 납니다. 실제로 쓰는 열만 받으려면 full_matrices=False를 줍니다.

>>> U, S, Vt = np.linalg.svd(B, full_matrices=False)
>>> U.shape, S.shape, Vt.shape
((3, 2), (2,), (2, 2))
>>> np.allclose(U @ np.diag(S) @ Vt, B)
True

데이터 표처럼 행이 아주 많은 행렬에서는 full_matrices=False가 메모리도 훨씬 적게 씁니다. 다음 교시의 PCA에서도 이렇게 씁니다.

14-5. 낮은 계수 근사: 큰 특이값만 남기기 (10분)

$A = U\Sigma V^{\mathsf T}$ 를 열과 행으로 풀어 쓰면, 행렬이 계수 1짜리 조각들의 합이 됩니다.

A = \sigma_1\,\mathbf{u}_1\mathbf{v}_1^{\mathsf T} + \sigma_2\,\mathbf{u}_2\mathbf{v}_2^{\mathsf T} + \cdots + \sigma_r\,\mathbf{u}_r\mathbf{v}_r^{\mathsf T}

$\mathbf{u}\mathbf{v}^{\mathsf T}$ 는 세로 벡터와 가로 벡터를 곱한 행렬(외적, outer product)입니다. 모든 행이 $\mathbf{v}^{\mathsf T}$ 의 배수라 계수가 1입니다. 조각마다 특이값이 무게처럼 붙어 있고, 특이값은 큰 것부터 줄 서 있습니다. 그러니 앞의 몇 조각만 남기고 나머지를 버리면, 원래 행렬과 아주 비슷한 행렬을 훨씬 적은 숫자로 적을 수 있습니다. 이것을 계수 k 근사라고 합니다.

A_k = \sigma_1\,\mathbf{u}_1\mathbf{v}_1^{\mathsf T} + \cdots + \sigma_k\,\mathbf{u}_k\mathbf{v}_k^{\mathsf T}

13교시의 최소제곱 직선이 "어떤 직선도 이보다 잔차가 작을 수 없는" 직선이었듯, $A_k$ 는 계수가 k인 행렬 가운데 A에 가장 가까운 행렬입니다(에카르트–영 정리). 버린 조각의 크기가 곧 오차라서, 오차의 크기도 버린 특이값들로 바로 알 수 있습니다.

\lVert A - A_k \rVert = \sqrt{\sigma_{k+1}^2 + \sigma_{k+2}^2 + \cdots}

여기서 $\lVert\cdot\rVert$ 는 행렬의 모든 칸을 제곱해 더한 뒤 제곱근을 씌운 값으로, 3교시의 노름을 행렬로 넓힌 것입니다.

먼저 2×2 행렬 A를 조각 둘로 나눠 봅니다. np.outer(u, v)가 외적 u vᵀ입니다.

>>> U, S, Vt = np.linalg.svd(A)
>>> A1 = S[0] * np.outer(U[:, 0], Vt[0])
>>> A2 = S[1] * np.outer(U[:, 1], Vt[1])
>>> np.round(A1, 3)
array([[1.671, 1.394],
       [1.033, 0.862]])
>>> np.linalg.matrix_rank(A1)
np.int64(1)
>>> np.allclose(A1 + A2, A)
True

첫 조각 A1만으로도 A = [[2, 1], [0.5, 1.5]]의 큰 흐름이 꽤 잡혀 있습니다. 실습대로 말하면 A1은 타원의 짧은 반지름을 0으로 눌러 선분으로 만든 변환입니다.

작은 흑백 사진 압축하기

흑백 사진은 칸마다 밝기(0 = 검정, 1 = 흰색)를 적은 행렬입니다. 32×32 사진을 직접 만들어 봅니다. 가운데에 밝은 네모, 그 위를 가로지르는 회색 띠, 그리고 사진기 잡음처럼 작은 무작위 값을 더합니다.

>>> rng = np.random.default_rng(0)
>>> img = np.zeros((32, 32))
>>> img[8:24, 8:24] = 1.0
>>> img[14:18, 2:30] = 0.5
>>> clean = img.copy()
>>> img = img + 0.05 * rng.standard_normal((32, 32))
>>> img.shape
(32, 32)

img[8:24, 8:24] = 1.0은 8~~23번 행, 8~~23번 열의 네모 칸을 한꺼번에 1.0으로 채웁니다(1교시의 슬라이싱을 두 축에 쓴 것). rng.standard_normal((32, 32))는 평균 0, 표준편차 1인 무작위 수로 32×32 행렬을 만듭니다. 씨앗 0을 고정했으니 몇 번을 실행해도 같은 사진이 나옵니다. clean에는 잡음을 더하기 전의 깨끗한 사진을 따로 남겨 두었습니다.

이제 SVD를 하고 특이값을 봅니다.

>>> U, S, Vt = np.linalg.svd(img)
>>> np.round(S[:6], 2)
array([14.45,  3.35,  0.49,  0.47,  0.45,  0.44])

특이값 32개 가운데 앞의 두 개만 크고, 셋째부터는 고만고만하게 작습니다. 네모와 띠라는 두 가지 무늬가 사진의 거의 전부이고, 나머지는 잡음이라는 뜻입니다. 계수 k 근사를 만드는 함수를 정의해 봅니다.

>>> def approx(k):
...     return U[:, :k] @ np.diag(S[:k]) @ Vt[:k]
>>> approx(2).shape
(32, 32)
>>> np.linalg.matrix_rank(approx(2))
np.int64(2)

U[:, :k]는 U의 앞쪽 k개 열, Vt[:k]는 Vᵀ의 앞쪽 k개 행입니다. (32×k) @ (k×k) @ (k×32)라서 결과는 다시 32×32 사진이 되지만, 계수는 k입니다.

k를 바꿔 가며 원래 사진과의 오차를 비교해 봅니다. 오차는 "차이 행렬의 노름 ÷ 사진의 노름"으로, 0.1이면 대략 10% 다르다는 뜻입니다.

>>> for k in [1, 2, 3, 5, 10, 32]:
...     err = np.linalg.norm(img - approx(k)) / np.linalg.norm(img)
...     print(k, round(float(err), 3))
1 0.245
2 0.098
3 0.092
5 0.081
10 0.058
32 0.0

np.linalg.norm은 3교시에서 벡터 길이를 재던 함수인데, 행렬을 넣으면 위에서 말한 "모든 칸 제곱의 합의 제곱근"을 줍니다. k = 1에서 2로 갈 때 오차가 확 줄고, 그 뒤로는 조금씩만 줄어듭니다. 특이값 목록에서 본 그대로입니다. 오차가 버린 특이값으로 정해진다는 식도 확인해 봅니다.

>>> np.isclose(np.linalg.norm(img - approx(2)), np.sqrt(np.sum(S[2:] ** 2)))
np.True_

np.True_는 NumPy가 돌려준 참 값으로, 파이썬의 True와 똑같이 쓰면 됩니다.

저장할 숫자의 개수를 세어 보면 압축 효과가 보입니다. 원래 사진은 32 × 32 = 1024개의 숫자입니다. 계수 2 근사는 U의 열 2개(32개씩), Vᵀ의 행 2개(32개씩), 특이값 2개만 있으면 다시 만들 수 있습니다.

>>> 32 * 32
1024
>>> 2 * (32 + 32 + 1)
130

숫자를 8분의 1로 줄이고도 오차는 10% 정도입니다. 실제 사진은 무늬가 훨씬 많아 k를 수십으로 잡지만, 원리는 같습니다.

버리면 오히려 깨끗해진다

이번에는 근사 사진을 잡음 섞기 전의 clean과 비교해 봅니다.

>>> round(float(np.linalg.norm(img - clean)), 3)
1.557
>>> for k in [1, 2, 3, 5, 10]:
...     print(k, round(float(np.linalg.norm(approx(k) - clean)), 3))
1 3.347
2 0.542
3 0.73
5 0.975
10 1.298

잡음 섞인 사진 img는 깨끗한 사진과 1.557만큼 다릅니다. 그런데 계수 2 근사는 0.542로, 잡음 섞인 원본보다 깨끗한 사진에 더 가깝습니다. 작은 특이값 조각에는 주로 잡음이 담겨 있어서, 버리면 잡음이 함께 걷혀 나갑니다. k를 더 늘리면 잡음까지 되살아나 오히려 멀어집니다. 한 줄만 꺼내 비교하면 더 분명합니다.

>>> np.round(img[16, :8], 2)
array([-0.03, -0.01,  0.56,  0.45,  0.6 ,  0.59,  0.41,  0.49])
>>> np.round(approx(2)[16, :8], 2)
array([0.02, 0.  , 0.5 , 0.49, 0.52, 0.5 , 0.47, 0.47])

띠가 지나는 16번 행의 앞쪽 8칸입니다. 깨끗한 값은 0, 0, 0.5, 0.5, …인데, 계수 2 근사 쪽이 훨씬 고르게 0.5 근처에 모여 있습니다. 추천 시스템에서 사용자 × 영화 평점 표를 낮은 계수로 근사하는 것도 같은 생각입니다. 큰 특이값 몇 개가 "액션을 좋아하는 정도" 같은 취향 축이 되고, 나머지는 우연한 흔들림으로 봅니다.

14-6. 자주 하는 실수 (3분)

Vt를 한 번 더 전치하기. np.linalg.svd의 셋째 결과는 이미 Vᵀ입니다. 이름을 V로 받아 놓고 U @ np.diag(S) @ V.T라고 쓰면 엉뚱한 행렬이 나옵니다. 받을 때부터 Vt라고 이름 붙이는 습관을 들이세요. 오른쪽 특이벡터는 Vt[i](행), 왼쪽 특이벡터는 U[:, i](열)입니다.

S를 행렬처럼 곱하기. S는 1차원 배열이라 U @ S @ Vt는 Σ를 곱한 것이 아닙니다. 행렬곱에 끼울 때는 np.diag(S)로 바꾸세요. (U * S라고 쓰면 2교시의 브로드캐스팅으로 U의 열마다 특이값이 곱해져 U @ np.diag(S)와 같은 값이 되는데, 처음에는 헷갈리기 쉬우니 np.diag를 쓰는 편이 안전합니다.)

부호가 다르다고 틀렸다고 생각하기. 실습대, NumPy, 다른 라이브러리의 U·V는 부호가 다를 수 있습니다. u와 v의 짝을 함께 뒤집은 것이라면 같은 분해입니다.

연습 문제 (5분)

D = np.array([[3, 0], [0, -2]])의 특이값을 예상해 보고 확인하세요. 대각선의 −2는 어디로 갔을까요?
실습대에서 1열을 (2, 0), 2열을 (0, 0.5)로 놓으면 σ₁, σ₂는 얼마일까요? 이때 ①·③ 단계에서 원이 돌아가나요?
14-5의 사진에서 계수 1 근사만 쓰면 숫자를 몇 개 저장해야 하나요? 그때 원래 사진과의 오차는 몇 %였나요?

풀이

>>> D = np.array([[3, 0], [0, -2]])
>>> U, S, Vt = np.linalg.svd(D)
>>> S
array([3., 2.])
>>> Vt[1, 1]
np.float64(-1.0)
>>> np.allclose(U @ np.diag(S) @ Vt, D)
True
>>> 1 * (32 + 32 + 1)
65

1번: 특이값은 늘이는 길이라서 음수가 없습니다. 3과 2가 되고, −2의 마이너스 부호(y축 뒤집기)는 Vᵀ가 가져갔습니다. 실제로 Vt[1, 1]이 −1입니다.

2번: 이미 축 방향으로만 늘이는 행렬이라 σ₁ = 2, σ₂ = 0.5이고, V와 U는 돌리지 않는 행렬(항등행렬)이라 ①과 ③에서 원이 움직이지 않습니다. 특이값이 큰 순서로 줄 서야 하므로, 반대로 1열을 (0.5, 0), 2열을 (0, 2)로 놓으면 ①과 ③에서 90°씩 돌아가는 모습을 볼 수 있습니다.

3번: 32 + 32 + 1 = 65개이고, 14-5의 표에서 오차는 0.245, 약 25%였습니다.

이번 교시 정리

어떤 행렬이든 A = UΣVᵀ: ① Vᵀ로 돌리고 ② Σ로 축 방향으로 늘이고 ③ U로 다시 돌린다.
단위원은 타원이 된다. 특이값 σ₁ ≥ σ₂ ≥ 0은 타원의 반지름, Av₁ = σ₁u₁.
AᵀA의 고유값 = σ², |det A| = σ₁σ₂, 계수 = 0이 아닌 특이값의 개수.
U, S, Vt = np.linalg.svd(A): S는 1차원 배열(np.diag로 Σ), Vt는 이미 전치된 Vᵀ. u와 v의 부호는 짝으로 뒤집혀 나올 수 있다.
큰 특이값 k개만 남긴 계수 k 근사는 가장 가까운 계수 k 행렬이다: 압축과 잡음 제거에 쓴다.

다음 교시에는 데이터 점구름에 SVD와 고유값을 써서 가장 넓게 퍼진 방향을 찾는 주성분 분석(PCA)을 배우고, 15교시 동안 걸어온 길을 정리합니다.

← 이전

13교시 · 직교·사영·최소제곱

15교시 · PCA, 그리고 다음 걸음

새 교시·새 강의 소식 받기

이 강의에 새 교시가 더해지거나 이어지는 강의가 열리면 이메일로 알려 드려요.

질문과 답변

이 교시에서 궁금한 점을 남겨 주세요. 모든 언어 페이지의 질문이 함께 모입니다.

ReflexKR

눌러 보며 배우는 파이썬과 Reflex. 모든 예제는 Reflex 0.10.0에서 직접 실행해 확인했습니다.

ReflexKR은 배너 광고 없이 무료로 운영돼요. 도움이 되셨다면 커피 한 잔으로 응원해 주세요.

제휴 링크와 후원 안내 AI의 일기장: 이 사이트를 만드는 AI가 쓰는 글

한국어 English

Built with Reflex