PCA, 그리고 다음 걸음

1교시에서 검색 엔진이 문장을 숫자 768개짜리 벡터로 바꾼다고 했습니다. 이런 데이터는 그림으로 그릴 수 없고, 열이 많을수록 계산도 무거워집니다. 그런데 실제 데이터의 열들은 서로 얽혀 있는 경우가 많습니다. 키가 크면 몸무게도 대체로 많이 나가는 것처럼요. 그렇다면 숫자 여러 개를 따로 들고 다니는 대신, 데이터가 가장 넓게 퍼진 방향 몇 개만 골라 그 방향의 좌표만 남겨도 정보를 거의 잃지 않습니다.

이 생각을 정리한 방법이 주성분 분석(PCA, Principal Component Analysis) 입니다. 데이터 시각화(768차원 → 2차원 지도), 머신러닝 전처리, 잡음 제거, 얼굴 사진 압축 등에서 가장 먼저 꺼내 드는 도구입니다.

PCA에는 새 개념이 거의 없습니다. 3교시의 내적과 그림자, 10교시의 기저 바꾸기, 11교시의 대칭행렬 고유벡터, 13교시의 사영, 14교시의 SVD를 한 줄로 꿰면 PCA가 됩니다. 이번 교시는 그렇게 강의 전체를 묶고, 마지막에 다음 강의인 수치해석을 예고합니다.

15-1. 데이터는 행렬이다 (5분)

데이터 표는 그대로 행렬입니다. 약속은 하나입니다. 행 하나 = 표본(사람, 사진, 문장) 하나, 열 하나 = 특성(키, 몸무게, 밝기) 하나. 다섯 사람의 키(cm)와 몸무게(kg)를 적어 봅니다.

각 행은 2차원 벡터, 즉 평면 위의 점 하나입니다. 데이터 전체는 점구름입니다.

>>> import numpy as np
>>> X = np.array([[170, 65], [160, 52], [180, 80], [175, 70], [165, 60]])
>>> X.shape
(5, 2)

(5, 2)는 "표본 5개, 특성 2개"로 읽습니다. 앞으로 데이터 행렬은 늘 이 모양(표본 수, 특성 수)으로 둡니다.

axis: 어느 방향으로 접을까

평균을 구할 때는 어느 방향으로 평균을 낼지 정해야 합니다. 이것이 axis입니다.

>>> X.mean(axis=0)
array([170. ,  65.4])
>>> X.mean(axis=1)
array([117.5, 106. , 130. , 122.5, 112.5])
>>> X.mean()
np.float64(117.7)

axis=0은 0번 축(행이 늘어선 방향, 즉 위에서 아래로)을 따라 접어서 없앱니다. 다섯 행이 하나로 접히니 열마다 평균 하나, "평균 키 170, 평균 몸무게 65.4"가 남습니다. axis=1은 1번 축(왼쪽에서 오른쪽)을 접어 사람마다 "키와 몸무게의 평균"을 줍니다. 여기서는 뜻 없는 숫자입니다. axis를 빼면 숫자 10개 전체의 평균입니다.

외우는 법: axis=k는 shape의 k번째 자리를 없앤다. (5, 2)에 axis=0을 쓰면 (2,)가, axis=1을 쓰면 (5,)가 남습니다. PCA에서는 거의 언제나 axis=0(특성마다)입니다.

15-2. 평균 빼기와 공분산 (8분)

중심화. 먼저 모든 점에서 평균을 뺍니다. 점구름 전체를 평행이동해 중심을 원점에 놓는 것입니다(2교시: 화살표는 옮겨도 같은 벡터). 이렇게 해야 "어느 방향으로 퍼졌나"를 원점에서 뻗는 벡터로 잴 수 있습니다.

분산. 한 특성이 평균에서 얼마나 퍼졌는지는 "평균과의 차이를 제곱해 평균 낸 값"인 분산으로 잽니다. 표본 $n$개로 계산할 때는 $n$ 대신 $n-1$로 나누는 것이 관례입니다(표본으로 전체를 짐작할 때 조금 작게 나오는 것을 바로잡는 약속이라고만 알아 두세요).

공분산. 두 특성이 함께 움직이는 정도입니다. 평균과의 차이를 서로 곱해 평균 냅니다. 한쪽이 평균보다 클 때 다른 쪽도 크면 곱이 양수, 반대면 음수가 쌓입니다. 3교시의 내적 부호(같은 쪽 / 수직 / 반대)와 같은 그림입니다. 실제로 중심화한 두 열의 내적을 $n-1$로 나눈 것이 공분산입니다.

이것을 한 번에 모아 놓은 것이 공분산 행렬입니다. 중심화한 데이터를 $X_c$라 하면

대각선에는 각 특성의 분산이, 대각선 밖에는 공분산이 놓입니다. $C^{\mathsf T} = C$인 대칭행렬이라는 점이 중요합니다. 11교시에서 대칭행렬은 서로 수직인 고유벡터를 가진다고 했습니다.

책에 따라 공분산 행렬을 $\Sigma$로 쓰기도 합니다. 14교시의 특이값 행렬 Σ와 글자가 같으니 문맥으로 구분하세요. 이 교안에서는 $C$로 씁니다.

NumPy로 차례로 계산해 봅니다.

>>> Xc = X - X.mean(axis=0)
>>> Xc
array([[  0. ,  -0.4],
       [-10. , -13.4],
       [ 10. ,  14.6],
       [  5. ,   4.6],
       [ -5. ,  -5.4]])
>>> np.allclose(Xc.mean(axis=0), 0)
True

X는 shape이 (5, 2), X.mean(axis=0)은 (2,)입니다. 모양이 다른데 뺄셈이 되는 것은 2교시의 브로드캐스팅 덕분입니다. 길이 2짜리 평균이 다섯 행 각각에 똑같이 적용됩니다. 중심화한 뒤의 평균은 0이어야 하는데, 부동소수점 찌꺼기(-5.7e-15 같은 값)가 남을 수 있어 np.allclose로 확인했습니다.

공분산 행렬을 식대로 만든 것과 np.cov가 주는 것을 비교합니다.

>>> Xc.T @ Xc / (len(X) - 1)
array([[ 62.5,  82.5],
       [ 82.5, 110.8]])
>>> C = np.cov(X, rowvar=False)
>>> C
array([[ 62.5,  82.5],
       [ 82.5, 110.8]])
>>> np.var(X[:, 0], ddof=1)
np.float64(62.5)

키의 분산이 62.5, 몸무게의 분산이 110.8, 공분산이 82.5(양수: 키가 크면 몸무게도 큼)입니다. np.cov는 알아서 평균을 빼고 n − 1로 나눠 줍니다.

rowvar=False는 "행이 변수가 아니다(열이 변수다)"라는 뜻입니다. np.cov는 기본값이 rowvar=True라서, 아무 말 없이 넣으면 각 행을 특성으로 읽습니다. 우리 표는 행이 사람이니 꼭 rowvar=False를 줘야 합니다. np.var의 ddof=1은 n − 1로 나누라는 뜻입니다. np.var는 기본값이 n으로 나누기라 np.cov와 맞추려면 이렇게 씁니다.

15-3. 주성분 = 공분산 행렬의 고유벡터 (9분)

이제 PCA의 질문입니다. 점구름을 어느 방향에서 보면 가장 넓게 퍼져 보일까?

길이 1인 방향 $\mathbf{w}$를 하나 고르면, 중심화한 점 $\mathbf{x}$마다 $\mathbf{w}$ 위로 드리운 그림자 길이는 3교시의 내적 $\mathbf{x}\cdot\mathbf{w}$입니다. 점 전체의 그림자 길이를 한 번에 구하면 $X_c\mathbf{w}$이고, 이 그림자들의 분산을 계산하면 깔끔한 식이 나옵니다.

$\mathbf{w}^{\mathsf T}C\mathbf{w}$를 가장 크게 만드는 길이 1짜리 $\mathbf{w}$는 C의 가장 큰 고유값에 대한 고유벡터이고, 그때의 값이 바로 그 고유값 $\lambda_1$입니다. 14교시 실습대를 떠올리면 이해가 쉽습니다. 대칭행렬 C로 단위원을 옮기면 타원이 되는데, 긴 반지름 방향이 가장 크게 늘어나는 방향이고, 대칭행렬에서는 그 방향이 고유벡터입니다.

• PC1(제1 주성분): C의 가장 큰 고유값 $\lambda_1$의 고유벡터. 점구름이 가장 넓게 퍼진 방향.
• PC2: 다음 고유값 $\lambda_2$의 고유벡터. C가 대칭이라 PC1과 수직입니다.
• $\lambda_i$ = 그 방향으로 잰 분산.
• 설명된 분산 비율: $\dfrac{\lambda_1}{\lambda_1 + \lambda_2}$. 전체 퍼짐 가운데 PC1 하나가 담는 몫입니다.

분모 $\lambda_1 + \lambda_2$는 C의 대각선 합 $\operatorname{Var}(\text{x}) + \operatorname{Var}(\text{y})$와 같습니다. 축을 돌려 다시 재도 전체 퍼짐은 그대로이고, 그 퍼짐을 PC1에 최대한 몰아주는 방향 잡기가 PCA입니다.

아래 실습대에는 점 60개짜리 점구름이 있습니다. 노랑 화살표가 PC1, 회색 화살표가 PC2이고, 화살표 길이는 그 방향 분산의 제곱근(표준편차)의 2배, 즉 2√λ입니다.

실습대에서 눈여겨볼 점입니다.

• 검정 손잡이(중심)를 끌어도 공분산 행렬과 PC1은 변하지 않습니다. 평균을 빼고 계산하니 점구름이 어디에 있든 퍼진 모양만 봅니다.
• 처음 화면은 기울기 30°인데 PC1은 약 (0.855, 0.518), 각도로 약 31°입니다. 점 60개는 무작위로 뽑은 표본이라 만든 방향과 딱 맞지는 않습니다. 실제 데이터에서도 PCA가 찾는 것은 "이 표본에서" 가장 넓은 방향입니다.
• '짧은 폭'을 '긴 폭'과 같게 하면 둥근 점구름이 되어 비율이 59% 안팎으로 떨어집니다. 정확히 50%가 아닌 것은 점 60개가 우연히 어느 한쪽으로 조금 더 퍼져 있기 때문이고, 이때 PC1은 '기울기'로 정한 방향과 상관없이 그 우연한 방향을 가리킵니다. 특별히 넓은 방향이 없는 데이터에서는 PC1을 믿기 어렵다는 뜻입니다. '짧은 폭'을 더 키우면 PC1이 90° 휙 돌아갑니다. 슬라이더 이름과 상관없이 PC1은 늘 실제로 더 넓은 쪽입니다.
• 'PC1 위로 사영하기'를 켜면 점마다 PC1 직선 위로 수직으로 내려앉습니다(13교시의 사영). 점마다 숫자 두 개 대신 직선 위의 위치 숫자 하나만 남는, 2차원 → 1차원 줄이기입니다.

NumPy 칸 읽기

실습대 아래 "NumPy로 확인하기" 칸은 처음 화면에서 이렇게 나옵니다.

>>> X = np.array(points)  # 실습대의 점 60개
>>> X.shape
(60, 2)
>>> np.round(X.mean(axis=0), 4)
array([0.5   , 0.5001])
>>> C = np.cov(X, rowvar=False)
>>> np.round(C, 4)
array([[1.6861, 0.8735],
       [0.8735, 0.7732]])
>>> vals, vecs = np.linalg.eigh(C)
>>> np.round(vals, 4)
array([0.244 , 2.2152])
>>> np.round(vecs[:, -1], 4)
array([-0.8553, -0.5181])
>>> np.round(vals[-1] / vals.sum(), 4)
np.float64(0.9008)
>>> z = (X - X.mean(axis=0)) @ vecs[:, -1]
>>> z.shape
(60,)

(실습대의 점 목록 points는 페이지 안에만 있어서, 이 블록은 따라 칠 수 있는 코드가 아니라 읽기용입니다.)

한 줄씩 읽어 봅니다.

• np.round(…, 4): 14교시에서 본 반올림 함수입니다. 소수점 아래가 길게 늘어지지 않도록 넷째 자리까지만 보여 줍니다. 값 자체는 그대로이고, 보여 줄 때만 반올림합니다. 그래서 0.244처럼 끝자리 0이 빠진 숫자는 넷째 자리가 0이라는 뜻입니다.
• X.shape이 (60, 2): 점(표본) 60개, 특성 2개(x, y).
• X.mean(axis=0): 점구름의 중심, 즉 검정 손잡이 자리입니다. 점 좌표를 소수 셋째 자리에서 반올림해 그려서 0.5001처럼 아주 조금 어긋납니다.
• np.linalg.eigh(C): 11교시에서 본 대칭행렬 전용 고유값 함수입니다. 고유값을 작은 것부터 줄 세워 줍니다. 그래서 가장 큰 고유값은 마지막 vals[-1]이고, 그 고유벡터는 vecs의 마지막 열 vecs[:, -1]입니다. 이것이 PC1입니다.
• vals[-1] / vals.sum(): 설명된 분산 비율. 0.9008이니 실습대의 "90.1%"와 같습니다. vals.sum()은 배열의 숫자를 모두 더합니다.
• z = (X - X.mean(axis=0)) @ vecs[:, -1]: 중심화한 점들을 PC1 위로 사영한 그림자 길이. 점 60개가 숫자 60개가 되었으니 shape이 (60,)입니다.

실습대는 PC1을 (0.855, 0.518)로, NumPy는 (−0.8553, −0.5181)로 보여 줍니다. 14교시에서 본 부호 문제와 같습니다. 고유벡터는 −1을 곱해도 고유벡터이고, 같은 직선을 가리킵니다. PCA에서 중요한 것은 방향이 놓인 직선이지 화살표 끝이 어느 쪽이냐가 아닙니다.

15-4. 사영으로 차원 줄이기 (6분)

키·몸무게 데이터로 돌아와 PCA를 끝까지 해 봅니다.

>>> vals, vecs = np.linalg.eigh(C)
>>> np.round(vals, 3)
array([  0.688, 172.612])
>>> np.round(vecs, 3)
array([[-0.8,  0.6],
       [ 0.6,  0.8]])
>>> v = vecs[:, -1]
>>> round(float(vals[-1] / vals.sum()), 4)
0.996

PC1은 (0.6, 0.8) 방향, "키와 몸무게가 함께 커지는 방향"이고 전체 분산의 99.6%를 혼자 설명합니다. PC2는 (−0.8, 0.6) 방향으로, "키에 비해 몸무게가 많은지 적은지"를 잽니다. 두 방향은 수직입니다.

>>> bool(np.isclose(v @ vecs[:, 0], 0))
True

이제 사람마다 PC1 위의 그림자 길이를 구합니다. 3교시의 내적 그대로입니다.

>>> z = Xc @ v
>>> np.round(z, 2)
array([ -0.32, -16.72,  17.68,   6.68,  -7.32])
>>> np.round(np.var(z, ddof=1), 3)
np.float64(172.612)

사람마다 숫자 두 개(키, 몸무게)가 숫자 하나(PC1 좌표)가 되었습니다. 이 숫자들의 분산은 정확히 λ₁입니다. 3번 사람(17.68)이 "가장 크고 무거운" 쪽, 2번 사람(−16.72)이 반대쪽 끝입니다.

숫자 하나로 원래 점을 얼마나 되살릴 수 있을까요? 평균에서 출발해 PC1 방향으로 z만큼 가면 됩니다. 실습대에서 '사영하기'를 켰을 때 점들이 내려앉은 자리가 이것입니다.

>>> X_hat = X.mean(axis=0) + np.outer(z, v)
>>> np.round(X_hat, 1)
array([[169.8,  65.1],
       [160. ,  52. ],
       [180.6,  79.5],
       [174. ,  70.7],
       [165.6,  59.5]])

np.outer(z, v)는 14교시의 외적으로, 사람마다 z[i] * v를 한 줄씩 쌓은 (5, 2) 행렬입니다. 원래 표와 1~2 정도씩만 다릅니다. 잃어버린 것은 PC2 방향의 퍼짐, 즉 λ₂ = 0.688만큼입니다.

10교시의 눈으로 보면 PCA는 기저 바꾸기입니다. PC1, PC2를 새 기저로 삼아 좌표를 다시 재는 것이고, 그 가운데 분산이 작은 좌표를 버리는 것입니다. vecs는 열끼리 수직이고 길이가 1인 직교행렬이라 역행렬을 구할 필요 없이 Xc @ vecs가 곧 새 좌표입니다.

15-5. 전체 흐름, 그리고 SVD로 같은 답 (9분)

실습대와 비슷한 점구름 200개를 직접 만들어 PCA 전체를 한 번에 해 봅니다. 가로 폭 2, 세로 폭 0.5인 점구름을 30° 돌리고 중심을 (3, 1)로 옮깁니다.

>>> rng = np.random.default_rng(0)
>>> base = rng.standard_normal((200, 2)) * [2.0, 0.5]
>>> th = np.radians(30)
>>> R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
>>> X = base @ R.T + [3, 1]
>>> X.shape
(200, 2)

* [2.0, 0.5]는 브로드캐스팅으로 첫째 열을 2배, 둘째 열을 0.5배 합니다. R은 6교시의 회전 행렬입니다. 그런데 R @ base가 아니라 base @ R.T라고 쓴 데 주의하세요. 점이 행으로 들어 있어서, 점 하나하나에 R을 곱하려면 오른쪽에서 전치를 곱해야 합니다. (행 벡터 p에 대해 p Rᵀ = (R p)ᵀ입니다.)

>>> mean = X.mean(axis=0)
>>> Xc = X - mean
>>> C = np.cov(X, rowvar=False)
>>> vals, vecs = np.linalg.eigh(C)
>>> np.round(vals, 3)
array([0.254, 3.856])
>>> pc1 = vecs[:, -1]
>>> np.round(pc1, 3)
array([-0.859, -0.513])
>>> if pc1[0] < 0:
...     pc1 = -pc1
>>> round(float(np.degrees(np.arctan2(pc1[1], pc1[0]))), 1)
30.8
>>> round(float(vals[-1] / vals.sum()), 3)
0.938

부호가 음수로 나와 if로 뒤집었습니다. 비교하기 편하게 "x 성분이 양수" 쪽으로 맞추는 흔한 습관입니다. np.arctan2(y, x)는 벡터 (x, y)의 각도를 라디안으로 주고, np.degrees로 도(°)로 바꿨습니다(3교시). 만들 때 넣은 30°에 가까운 30.8°가 나왔고, 분산(가로 폭² = 4, 세로 폭² = 0.25)도 대략 맞습니다.

SVD로도 같은 답 (14교시와 연결)

14교시에서 AᵀA의 고유벡터가 V, 고유값이 σ²이라고 했습니다. 공분산 행렬은 Xcᵀ Xc를 n − 1로 나눈 것이니, 중심화한 데이터의 SVD에서 곧바로 PCA가 나옵니다.

>>> U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
>>> U.shape, S.shape, Vt.shape
((200, 2), (2,), (2, 2))
>>> np.round(Vt[0], 3)
array([-0.859, -0.513])
>>> np.round(S ** 2 / (len(X) - 1), 3)
array([3.856, 0.254])
>>> np.allclose(S ** 2 / (len(X) - 1), vals[::-1])
True

• Vt[0](첫째 오른쪽 특이벡터)이 PC1입니다(부호까지 같게 나왔지만, 늘 그렇다는 보장은 없습니다).
• S ** 2 / (n - 1)이 고유값입니다. SVD는 큰 것부터, eigh는 작은 것부터라서 vals[::-1]로 순서를 뒤집어 비교했습니다. [::-1]은 "끝에서부터 하나씩 거꾸로" 꺼내는 슬라이싱입니다.
• 행이 200개라 full_matrices=False를 줬습니다. 안 주면 U가 200×200이 됩니다.

실무의 PCA 라이브러리(예: scikit-learn)는 공분산 행렬을 만들지 않고 이렇게 SVD로 계산합니다. Xcᵀ Xc를 곱하는 과정에서 생기는 부동소수점 오차를 피할 수 있기 때문인데, 이것이 바로 다음 강의 수치해석의 주제입니다.

특성이 많을 때

특성이 3개 이상이어도 코드는 그대로입니다. 셋째 특성을 "첫째 + 둘째 + 약간의 잡음"으로 만들어 보면, 셋째 열은 사실상 새 정보가 없습니다.

>>> t = rng.standard_normal(200)
>>> X3 = np.column_stack([X, X[:, 0] + X[:, 1] + 0.1 * t])
>>> X3.shape
(200, 3)
>>> vals3 = np.linalg.eigvalsh(np.cov(X3, rowvar=False))
>>> np.round(vals3[::-1], 2)
array([11.11,  0.27,  0.  ])
>>> np.round(np.cumsum(vals3[::-1]) / vals3.sum(), 4)
array([0.9764, 0.9997, 1.    ])

np.linalg.eigvalsh는 eigh에서 고유값만 주는 함수입니다. np.cumsum은 앞에서부터 차례로 더한 누적합입니다. 주성분 하나로 97.6%, 둘로 99.97%가 설명되니, 3차원 데이터를 2차원으로 줄여도 거의 잃는 것이 없습니다. 실제 데이터에서는 이 누적 비율이 90%나 95%를 넘는 지점까지 주성분을 남기는 방식을 많이 씁니다.

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

rowvar=False를 빼먹기. 행이 표본인 데이터에 그냥 np.cov(X)를 쓰면 표본끼리의 공분산이 나옵니다. shape을 보면 바로 알 수 있습니다.

>>> X = np.array([[170, 65], [160, 52], [180, 80], [175, 70], [165, 60]])
>>> np.cov(X).shape
(5, 5)
>>> np.cov(X, rowvar=False).shape
(2, 2)

eigh의 순서를 잊기. np.linalg.eigh는 작은 고유값부터 줍니다. vecs[:, 0]은 PC1이 아니라 가장 좁은 방향입니다. PC1은 vecs[:, -1]입니다.

단위가 다른 특성을 그대로 섞기. 분산은 단위에 따라 크기가 달라집니다. 같은 데이터에서 키만 미터로 바꿔 적으면 키의 분산이 10000분의 1이 되어, PC1이 몸무게만 가리키게 됩니다.

>>> Xm = X * [0.01, 1]
>>> vals_m, vecs_m = np.linalg.eigh(np.cov(Xm, rowvar=False))
>>> np.round(vecs_m[:, -1], 3)
array([0.007, 1.   ])
>>> Xs = (X - X.mean(axis=0)) / X.std(axis=0, ddof=1)
>>> vals_s, vecs_s = np.linalg.eigh(np.cov(Xs, rowvar=False))
>>> np.round(vecs_s[:, -1], 3)
array([0.707, 0.707])

단위가 다른 특성을 섞을 때는 각 열을 평균을 빼고 표준편차로 나눠(표준화) 모두 "분산 1"로 맞춘 뒤 PCA를 합니다. X.std(axis=0, ddof=1)은 열마다의 표준편차입니다.

평균을 빼지 않고 SVD하기. SVD로 PCA를 할 때 중심화를 빼먹으면, 첫째 특이벡터가 "가장 넓게 퍼진 방향"이 아니라 원점에서 점구름 중심을 향하는 방향이 되어 버립니다.

>>> np.round(np.linalg.svd(X)[2][0], 3)
array([-0.933, -0.361])
>>> np.round(X.mean(axis=0) / np.linalg.norm(X.mean(axis=0)), 3)
array([0.933, 0.359])

PC1인 (0.6, 0.8)이 아니라, 평균 (170, 65.4)를 정규화한 방향과 거의 같아졌습니다. 공분산은 np.cov가 알아서 평균을 빼 주지만, SVD는 그렇지 않습니다.

15-7. 15교시를 돌아보며 (4분)

숫자 묶음 하나에서 시작해 여기까지 왔습니다. 각 교시의 핵심과, 그것을 NumPy로 쓰는 한 줄입니다.

돌아보면 두 가지 그림이 강의 전체를 꿰고 있습니다. 하나는 행렬 = 공간을 옮기는 변환(5교시)이라는 그림입니다. 행렬식은 그 변환의 넓이 배율이고, 역행렬은 되돌리기, 고유벡터는 방향이 안 바뀌는 축, 특이값은 원이 타원이 될 때의 반지름이었습니다. 다른 하나는 Av = 열들의 일차결합(4교시)이라는 그림입니다. 연립방정식, 열공간, 최소제곱이 모두 "b를 A의 열로 만들 수 있나, 없으면 가장 가깝게는?"이라는 질문이었습니다.

이 두 그림이 머릿속에 있으면, 처음에 덮었던 선형대수 책을 다시 펼쳐도 기호 뒤의 그림이 보일 것입니다.

15-8. 다음 강의 예고: 수치해석 (3분)

이 강의 내내 몇 번이고 마주친 장면이 있습니다. 0.30000000000000004, 1.9999999999999996, np.allclose, "inv 대신 solve를 쓰세요"(8교시), "특이값이 1e-17처럼 나올 수 있습니다"(14교시). 수학으로는 깔끔한 계산이 컴퓨터 안에서는 조금씩 어긋납니다. 다음 강의 수치해석은 바로 이 어긋남을 정면으로 다룹니다.

>>> 0.1 + 0.2
0.30000000000000004
>>> A = np.array([[2, 1], [0.5, 1.5]])
>>> S = np.linalg.svd(A, compute_uv=False)
>>> round(float(S[0] / S[1]), 3)
2.618
>>> round(float(np.linalg.cond(A)), 3)
2.618
>>> round(float(np.linalg.cond(np.array([[1, 1], [1, 1.0001]]))))
40002

마지막 두 줄은 조건수입니다. 14교시의 특이값으로 정의되는데(σ₁ ÷ σ₂, 타원이 얼마나 납작한가), 이 값이 크면 입력이 아주 조금만 흔들려도 답이 크게 흔들립니다. 거의 납작한 두 번째 행렬은 조건수가 4만이 넘습니다.

다음 강의에서 다룰 내용을 미리 적어 둡니다.

• 부동소수점: 컴퓨터가 실수를 적는 방식과, 오차가 생기고 쌓이는 이유.
• 조건수와 안정성: 문제 자체가 예민한 경우와, 계산 방법이 오차를 키우는 경우를 구분하기. inv(A) @ b보다 solve가 나은 진짜 이유.
• LU 분해: np.linalg.solve가 속에서 하는 일(가우스 소거법을 행렬로).
• 반복법: 행렬이 수백만 × 수백만이라 한 번에 풀 수 없을 때, 답에 조금씩 다가가는 방법(야코비, 가우스–자이델, 켤레기울기법).
• 거듭제곱법과 QR 알고리즘: 12교시에 맛본 거듭제곱법을 다듬어, np.linalg.eig가 실제로 고유값을 구하는 방법까지.

이번 강의에서 "그림으로 무엇을 하는가"를 익혔다면, 다음 강의에서는 "컴퓨터가 그것을 어떻게, 얼마나 정확하게 하는가"를 배웁니다.

연습 문제 (5분)

1. 실습대에서 '기울기'를 0°로 놓으면 공분산 행렬의 대각선 밖(공분산)은 얼마에 가까워질까요? 그때 PC1은 어느 방향일까요?
2. 키·몸무게 데이터 X에서 PC2 좌표(Xc @ vecs[:, 0])의 분산을 구하고, λ₂와 같은지 확인하세요.
3. 15-5의 X3에서 누적 비율이 99%를 넘으려면 주성분이 몇 개 필요한가요?

풀이

>>> Xc = X - X.mean(axis=0)
>>> vals, vecs = np.linalg.eigh(np.cov(X, rowvar=False))
>>> z2 = Xc @ vecs[:, 0]
>>> bool(np.isclose(np.var(z2, ddof=1), vals[0]))
True
>>> np.round(vals[0], 3)
np.float64(0.688)
>>> int(np.argmax(np.cumsum(vals3[::-1]) / vals3.sum() > 0.99)) + 1
2

1번: 기울기가 0°이면 점구름이 x축에 나란히 퍼져 x와 y가 함께 움직이지 않으니 공분산이 0 가까이 되고, PC1은 (1, 0) 근처(또는 (−1, 0))를 가리킵니다. 공분산 행렬이 거의 대각행렬이 되어, 대각선의 분산이 그대로 고유값입니다.

2번: PC2 좌표의 분산은 λ₂ = 0.688과 같습니다. 각 주성분 좌표의 분산이 그 고유값이라는 15-4의 결과가 PC2에서도 성립합니다.

3번: np.cumsum(…) > 0.99는 누적 비율이 99%를 넘는 자리마다 True인 배열이고, np.argmax는 처음 True가 나오는 위치(0부터 셈)를 줍니다. 1을 더해 개수로 바꾸면 2개입니다.

이것으로 「그림으로 시작하는 선형대수」 15교시를 모두 마칩니다. 다음 강의 드론 관제실의 수치해석 기초에서는 컴퓨터가 실수를 어떻게 저장하는지, solve가 속에서 어떻게 계산되는지(가우스 소거와 LU), 그 답을 얼마나 믿을 수 있는지(조건수), 그리고 공식으로 풀리지 않는 방정식·적분·미분방정식을 어떻게 근사해 푸는지를 다룹니다.