보강 1. 주성분분석과 SVD

L34의 마지막에서 우리는 "다 배운 것이 아니라 지도를 손에 넣은 것"이라고 했다. 이제 그 지도를 들고 밖으로 나가 볼 차례다. 그리고 문을 열자마자 만나는 것이 주성분분석(principal component analysis, PCA)이다.

읽는 사람 가운데 적지 않은 수가 이미 PCA를 쓰고 있을 것이다. 함수 하나를 부르면 100개짜리 좌표가 2개로 줄어든다. 그런데 그 한 줄이 무슨 일을 하는지 물으면 대개 여기까지다. “공분산행렬의 고유벡터를 구한다.”

그러면 왜 하필 공분산행렬인가. 왜 하필 고유벡터인가. 그리고 이것이 이 강의의 진짜 물음인데, 실제 라이브러리를 열어 보면 공분산행렬을 아예 만들지 않는다. 데이터 행렬에 SVD를 바로 건다. 교과서가 시키는 대로 하지 않는 것이다. 왜인가.

세 물음의 답이 전부 지난 서른네 강의 안에 있다. 새로 배울 것은 하나도 없다.

1. 문제를 정확히 적는다¶

데이터 행렬¶

$m$ 개의 표본을 각각 $n$ 개의 변수로 재었다고 하자. 행 하나가 표본 하나이고, 열 하나가 변수 하나다.

X = \begin{bmatrix} \text{---} & \vv{x}_1^{\mathsf T} & \text{---} \\ & \vdots & \\ \text{---} & \vv{x}_m^{\mathsf T} & \text{---} \end{bmatrix} \in \R^{m \times n}

(1)

(1)의 배치를 L16과 헷갈리면 안 된다. 최소제곱에서 $A$ 의 열은 "설명변수"였고 우변 $\vv{b}$ 가 따로 있었다. 여기에는 우변이 없다. $X$ 하나만 놓고 그 자체의 모양을 묻는다. 정답이 없는 문제라는 점에서 성격이 다르다.

중심화¶

먼저 평균을 뺀다. 표본평균은 열마다 하나씩 나오므로 $n$ 차원 벡터다.

\boldsymbol{\mu} = \frac{1}{m}\sum_{i=1}^{m}\vv{x}_i = \frac{1}{m}X^{\mathsf T}\vv{1}, \qquad X_c = X - \vv{1}\boldsymbol{\mu}^{\mathsf T}

(2)

(2)의 $\vv{1}$ 은 성분이 전부 1인 $m$ 차원 벡터다. $\vv{1}\boldsymbol{\mu}^{\mathsf T}$ 는 랭크 1 행렬이고(L11), 모든 행이 $\boldsymbol{\mu}^{\mathsf T}$ 로 같다. 그것을 빼면 각 열의 합이 0이 된다.

이 뺄셈을 행렬 하나로 적을 수도 있다. $\boldsymbol{\mu}^{\mathsf T} = \vv{1}^{\mathsf T}X/m$ 을 (2)에 넣고 묶으면 이렇게 된다.

X_c = X - \frac{1}{m}\vv{1}\vv{1}^{\mathsf T}X = \left(I - \frac{\vv{1}\vv{1}^{\mathsf T}}{m}\right)X

(3)

(3)의 괄호 안이 낯이 익을 것이다. $\vv{1}^{\mathsf T}\vv{1} = m$ 이므로 $\vv{1}\vv{1}^{\mathsf T}/m$ 은 L15에서 배운 한 벡터로의 투영행렬 $\vv{a}\vv{a}^{\mathsf T}/\vv{a}^{\mathsf T}\vv{a}$ 바로 그것이다. 따라서 $I$ 에서 그것을 뺀 것은 $\vv{1}$ 에 수직인 공간으로의 투영이다. 중심화는 투영이었다.

방향 하나에 자를 댄다¶

단위벡터 $\vv{v}$ 를 하나 골라 모든 점을 그 방향에 투영한다. 투영한 좌표를 모으면 벡터 하나가 나온다.

\vv{y} = X_c\vv{v}, \qquad y_i = \vv{v}^{\mathsf T}\vv{x}_{c,i}, \qquad \lVert\vv{v}\rVert = 1

(4)

(4)의 $y_i$ 가 표본 $i$ 의 그림자 위치다. 중심화를 해 두었으므로 $\vv{1}^{\mathsf T}\vv{y} = \vv{1}^{\mathsf T}X_c\vv{v} = 0$ 이고, 곧 $\vv{y}$ 의 평균이 0이다. 평균이 0이면 표본분산은 제곱합을 $m-1$ 로 나눈 것뿐이다.

\operatorname{Var}(\vv{y}) = \frac{1}{m-1}\sum_{i=1}^{m} y_i^2 = \frac{\lVert X_c\vv{v}\rVert^2}{m-1} = \vv{v}^{\mathsf T}\!\left(\frac{X_c^{\mathsf T}X_c}{m-1}\right)\!\vv{v}

(5)

(5)의 변형이 이 강의 전체의 열쇠다. $\lVert X_c\vv{v}\rVert^2 = (X_c\vv{v})^{\mathsf T}(X_c\vv{v}) = \vv{v}^{\mathsf T}X_c^{\mathsf T}X_c\vv{v}$ 로 옮긴 것이 전부인데, 이 한 번의 이동으로 데이터에 관한 물음이 행렬 하나에 관한 물음으로 바뀐다. 가운데 괄호에 이름을 붙이자.

C = \frac{X_c^{\mathsf T}X_c}{m-1} \qquad(n \times n \ \text{표본 공분산행렬})

(6)

(6)의 $C$ 는 대각선에 각 변수의 분산이, 그 밖에 변수 쌍의 공분산이 앉은 행렬이다. 그러나 지금 중요한 것은 통계적 이름이 아니라 $C$ 가 $X_c^{\mathsf T}X_c$ 꼴이라는 사실이다. L27에서 이 꼴을 실컷 다루었다.

그래서 문제는 이것이다¶

\max_{\lVert\vv{v}\rVert = 1} \ \vv{v}^{\mathsf T}C\vv{v}

(7)

(7)의 문제를 푸는 $\vv{v}$ 를 제1주성분이라 부른다. 제약이 붙은 최적화처럼 보이지만 라그랑주 승수를 꺼낼 필요가 없다. 답은 이미 우리 손에 있다.

앵커¶

손으로 끝까지 따라갈 수 있는 작은 자료를 하나 쓴다. 다섯 개의 점을 두 변수로 잰 것이다.

X = \begin{bmatrix} 1 & 3 \\ 3 & 1 \\ 4 & 3 \\ 5 & 3 \\ 7 & 5 \end{bmatrix}, \qquad \boldsymbol{\mu} = (4,\ 3), \qquad X_c = \begin{bmatrix} -3 & 0 \\ -1 & -2 \\ 0 & 0 \\ 1 & 0 \\ 3 & 2 \end{bmatrix}

(8)

(8)의 평균은 열의 합 $1+3+4+5+7 = 20$ 과 $3+1+3+3+5 = 15$ 를 각각 5로 나눈 것이다. 중심화한 열의 합이 $-3-1+0+1+3 = 0$ 과 $0-2+0+0+2 = 0$ 이라 제대로 되었다.

이제 $C$ 를 만든다. 제곱합은 20 과 8 이고 교차합은 $(-3)(0)+(-1)(-2)+0+0+(3)(2) = 8$ 이다.

X_c^{\mathsf T}X_c = \begin{bmatrix} 20 & 8 \\ 8 & 8 \end{bmatrix}, \qquad C = \frac{1}{4}\begin{bmatrix} 20 & 8 \\ 8 & 8 \end{bmatrix} = \begin{bmatrix} 5 & 2 \\ 2 & 2 \end{bmatrix}

(9)

(9)의 $C$ 를 보면 첫 변수의 분산이 5, 둘째가 2, 둘의 공분산이 2다. $m-1 = 4$ 로 나눈 것에 주의하자. $m$ 으로 나누는 유파도 있는데, 그러면 모든 고윳값이 같은 비율로 줄어들 뿐이라 주성분의 방향은 달라지지 않는다.

왼쪽은 평균을 빼는 단계, 가운데는 방향 하나에 그림자를 떨어뜨린 모습, 오른쪽은
자를 0도에서 180도까지 돌리며 잰 분산이다. 가장 좋은 각도와 가장 나쁜 각도가
정확히 90도 떨어져 있고, 두 값의 합이 대각합 7이다. — Figure 1:왼쪽은 평균을 빼는 단계, 가운데는 방향 하나에 그림자를 떨어뜨린 모습, 오른쪽은 자를 0도에서 180도까지 돌리며 잰 분산이다. 가장 좋은 각도와 가장 나쁜 각도가 정확히 90도 떨어져 있고, 두 값의 합이 대각합 7이다.

2. 레일리 몫 — 답은 L25와 L27에 이미 있었다¶

레일리 몫¶

$\lVert\vv{v}\rVert = 1$ 이라는 제약이 거추장스러우니 분모로 옮긴다.

R(\vv{v}) = \frac{\vv{v}^{\mathsf T}C\vv{v}}{\vv{v}^{\mathsf T}\vv{v}}

(10)

(10)의 몫을 레일리 몫이라 한다. $\vv{v}$ 를 $c$ 배 해도 분자와 분모가 똑같이 $c^2$ 배가 되므로 값이 변하지 않는다. 곧 길이는 상관없고 방향만 상관있다. 그래서 $R$ 을 최대로 하는 것과 (7)의 문제를 푸는 것이 같은 일이다.

$C$ 는 대칭이고 준정부호다¶

$C$ 에 대해 우리가 아는 것을 먼저 적어 두자. 전치를 취하면 $(X_c^{\mathsf T}X_c)^{\mathsf T} = X_c^{\mathsf T}X_c$ 이므로 대칭이다. 에너지를 계산하면 (5)의 식이 그대로 답이다.

\vv{v}^{\mathsf T}C\vv{v} = \frac{\lVert X_c\vv{v}\rVert^2}{m-1} \ \ge\ 0

(11)

(11)의 부등식이 양의 준정부호의 정의다(L27). 등호는 $X_c\vv{v} = \vv{0}$ 일 때만 성립하므로, 변수들이 서로 독립이면 양의 정부호가 되고 그렇지 않으면 고윳값 0이 생긴다. 어느 쪽이든 고윳값이 음수가 되는 일은 없다. 분산이 음수일 수 없다는 상식이 행렬의 언어로 적힌 것이다.

스펙트럼 정리를 끼워 넣는다¶

대칭행렬이므로 L25가 곧바로 쓰인다.

C = Q\Lambda Q^{\mathsf T}, \qquad Q^{\mathsf T}Q = I, \qquad \lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_n \ge 0

(12)

(12)의 분해에서 $Q$ 의 열 $\vv{q}_1, \dots, \vv{q}_n$ 은 정규직교이고 고윳값은 전부 실수다. 준정부호이므로 음수도 없다. 이 두 가지가 아래 계산을 전부 가능하게 한다.

$Q$ 의 열이 $\R^n$ 의 기저이므로 아무 $\vv{v}$ 나 $\vv{v} = Q\vv{w}$ 로 적을 수 있다. $\vv{w}$ 는 $\vv{v}$ 를 고유벡터 기저에서 본 좌표다(L31). 이것을 (10)에 넣는다.

R(Q\vv{w}) = \frac{\vv{w}^{\mathsf T}Q^{\mathsf T}\,Q\Lambda Q^{\mathsf T}\,Q\vv{w}} {\vv{w}^{\mathsf T}Q^{\mathsf T}Q\vv{w}} = \frac{\vv{w}^{\mathsf T}\Lambda\vv{w}}{\vv{w}^{\mathsf T}\vv{w}} = \frac{\sum_{i} \lambda_i w_i^2}{\sum_{i} w_i^2}

(13)

(13)의 계산에서 한 일은 $Q^{\mathsf T}Q = I$ 를 두 번 지운 것뿐이다. 분자에서 가운데의 $Q^{\mathsf T}Q$ 가 사라져 $\Lambda$ 만 남고, 분모에서도 마찬가지다. 직교행렬은 길이를 바꾸지 않으므로(L17) 분모가 그대로 $\lVert\vv{w}\rVert^2$ 인 것이다.

최대값¶

(13)의 오른쪽을 다시 보자. $t_i = w_i^2 / \sum_j w_j^2$ 로 두면 $t_i \ge 0$ 이고 $\sum_i t_i = 1$ 이다.

R = \sum_{i=1}^{n}\lambda_i t_i, \qquad t_i \ge 0, \qquad \sum_{i=1}^{n} t_i = 1

(14)

(14)의 오른쪽은 고윳값들의 가중평균이다. 가중치를 어떻게 나누어 주든 평균이 최대값을 넘을 수는 없고 최소값 아래로 내려갈 수도 없다.

\lambda_n \ \le\ R(\vv{v}) \ \le\ \lambda_1

(15)

(15)의 두 부등식에서 위쪽 등호는 $t_1 = 1$ 일 때, 곧 $\vv{w} = \vv{e}_1$ 이고 $\vv{v} = Q\vv{e}_1 = \vv{q}_1$ 일 때만 성립한다. 라그랑주 승수도 미분도 쓰지 않았다. 정규직교 기저로 갈아 끼운 것이 증명의 전부다.

두 번째 주성분은 직교가 공짜다¶

제1주성분을 찾았으니 그다음을 찾는다. 같은 방향을 또 고르면 새로운 정보가 없으므로 $\vv{q}_1$ 에 수직인 방향 중에서 고른다.

\vv{v} \perp \vv{q}_1 \ \Longleftrightarrow\ w_1 = \vv{q}_1^{\mathsf T}\vv{v} = 0 \ \Longrightarrow\ R = \frac{\sum_{i \ge 2}\lambda_i w_i^2}{\sum_{i \ge 2} w_i^2} \ \le\ \lambda_2

(16)

(16)의 첫 동치가 요령이다. $\vv{v} = Q\vv{w}$ 이므로 $\vv{q}_1^{\mathsf T}\vv{v}$ 가 곧 $w_1$ 이고, $\vv{q}_1$ 에 수직이라는 것은 첫 좌표를 버린다는 뜻이다. 그러면 남은 합에서 (15)의 논증을 그대로 되풀이하면 되고, 최대값은 $\lambda_2$ 이며 $\vv{v} = \vv{q}_2$ 에서 달성된다.

앵커로 확인¶

(9)의 $C$ 에 특성방정식을 세운다(L21).

\det(C - \lambda I) = (5-\lambda)(2-\lambda) - 4 = \lambda^2 - 7\lambda + 6 = (\lambda - 6)(\lambda - 1)

(17)

(17)에서 $\lambda_1 = 6$ , $\lambda_2 = 1$ 이 나온다. $\lambda = 6$ 을 $C - \lambda I$ 에 넣으면 첫 행이 $-x + 2y = 0$ 이므로 $x = 2y$ 이고, $\lambda = 1$ 을 넣으면 첫 행이 $4x + 2y = 0$ 이므로 $y = -2x$ 다.

\vv{q}_1 = \frac{1}{\sqrt5}\begin{bmatrix} 2 \\ 1 \end{bmatrix}, \qquad \vv{q}_2 = \frac{1}{\sqrt5}\begin{bmatrix} -1 \\ 2 \end{bmatrix}, \qquad \vv{q}_1^{\mathsf T}\vv{q}_2 = \frac{-2+2}{5} = 0

(18)

(18)의 두 벡터가 정말 수직이다. 제1주성분의 기울기는 $1/2$ 이고 각도는 $\arctan(1/2) = 26.565^\circ$ 다. Figure 1의 오른쪽 그림에서 봉우리가 바로 이 자리에 있고, 골짜기는 $116.565^\circ$ 에 있다.

대각합이 총분산이다¶

마지막으로 하나 더. 대각합은 고윳값의 합이었다(L21).

\lambda_1 + \cdots + \lambda_n = \tr C = \sum_{j=1}^{n} C_{jj}

(19)

(19)의 오른쪽은 각 변수의 분산을 그냥 더한 것이다. 앵커에서는 $5 + 2 = 7$ 이고 왼쪽은 $6 + 1 = 7$ 이다. 좌표를 아무리 돌려도 총분산은 보존된다. 회전은 길이를 바꾸지 않으니 당연한데, 이 당연함이 다음 절에서 "설명분산 비율"의 분모가 된다.

왼쪽은 모든 방향에서 \vv{v}^{\mathsf T}C\vv{v} 를 재어 그 값만큼 뻗은 그림이다.
가장 멀리 뻗은 곳이 \vv{q}_1 이고 그 길이가 6이다. 오른쪽은 같은 정보를 타원으로
그린 것으로, 반축의 길이가 \sqrt{\lambda_i} 다. — Figure 2:왼쪽은 모든 방향에서 $\vv{v}^{\mathsf T}C\vv{v}$ 를 재어 그 값만큼 뻗은 그림이다. 가장 멀리 뻗은 곳이 $\vv{q}_1$ 이고 그 길이가 6이다. 오른쪽은 같은 정보를 타원으로 그린 것으로, 반축의 길이가 $\sqrt{\lambda_i}$ 다.

3. 공분산을 거치지 않는 길¶

SVD 하나면 된다¶

여기까지가 교과서의 이야기다. 그런데 $C$ 를 만들지 않고도 같은 답을 얻을 수 있다. $X_c$ 에 바로 SVD를 걸면 된다(L29).

X_c = U\Sigma V^{\mathsf T}, \qquad U^{\mathsf T}U = I, \qquad V^{\mathsf T}V = I

(20)

(20)의 분해를 (6)에 그대로 넣어 보자. 가운데에서 $U^{\mathsf T}U$ 가 $I$ 로 사라진다.

C = \frac{X_c^{\mathsf T}X_c}{m-1} = \frac{V\Sigma^{\mathsf T}\underbrace{U^{\mathsf T}U}_{=\,I}\Sigma V^{\mathsf T}}{m-1} = V\!\left(\frac{\Sigma^{\mathsf T}\Sigma}{m-1}\right)\!V^{\mathsf T}

(21)

(21)의 오른쪽은 대각행렬을 직교행렬로 감싼 꼴이니 그 자체가 스펙트럼 분해다. (12)의 식과 견주면 대응이 한눈에 보인다.

Q = V, \qquad \Lambda = \frac{\Sigma^{\mathsf T}\Sigma}{m-1}, \qquad \lambda_i = \frac{\sigma_i^2}{m-1}

(22)

(22)의 대응이 이 절의 결론이다. $V$ 의 열이 주성분이고, 고윳값은 특이값의 제곱을 $m-1$ 로 나눈 것이다. 공분산행렬은 한 번도 만들지 않았다.

점수¶

각 표본이 새 좌표계에서 어디에 있는지도 공짜로 나온다. (4)의 계산을 모든 주성분에 대해 한꺼번에 하면 이렇다.

Y = X_c V = U\Sigma V^{\mathsf T}V = U\Sigma

(23)

(23)의 $Y$ 를 주성분 점수라 한다. $U$ 의 열에 $\sigma_i$ 를 곱한 것이 전부이므로 따로 곱셈을 할 필요조차 없다. 그리고 $Y^{\mathsf T}Y = \Sigma^{\mathsf T}\Sigma$ 가 대각행렬이므로 새 좌표들은 서로 상관이 0이다. 뒤엉켜 있던 변수를 풀어 놓는 것이 PCA가 하는 일이다.

앵커¶

(8)의 $X_c$ 에 특이값을 재면 $\sigma_i^2$ 가 $X_c^{\mathsf T}X_c$ 의 고윳값이므로, (9)의 행렬에서 곧바로 얻어진다. 그 행렬은 $4C$ 이고 고윳값이 24 와 4 다.

\sigma_1 = \sqrt{24} = 2\sqrt6 \approx 4.8990, \qquad \sigma_2 = 2, \qquad \frac{\sigma_1^2}{4} = 6, \qquad \frac{\sigma_2^2}{4} = 1

(24)

(24)의 마지막 두 값이 (17)에서 구한 고윳값과 같다. 두 길이 같은 곳에 닿았다. 그러면 왜 굳이 SVD를 쓰는가. 5절이 그 답이다.

4. 설명분산은 L29의 에너지다¶

비율¶

주성분을 몇 개까지 남길지 정해야 한다. 기준으로 흔히 쓰는 것이 설명분산 비율이다.

p_i = \frac{\lambda_i}{\lambda_1 + \cdots + \lambda_n} = \frac{\sigma_i^2}{\sigma_1^2 + \cdots + \sigma_n^2}

(25)

(25)의 두 번째 등식은 분자와 분모에서 $m-1$ 이 함께 지워지기 때문에 성립한다. 그래서 설명분산 비율을 계산할 때는 $m-1$ 이든 $m$ 이든 상관이 없다. 앵커에서는 $p_1 = 6/7 \approx 0.857$ , $p_2 = 1/7 \approx 0.143$ 이다.

그런데 오른쪽 표현을 L29에서 본 적이 있다. 거기서는 $\sigma_i^2$ 를 에너지라 불렀고 $\sum\sigma_i^2 = \lVert X_c\rVert_F^2$ 였다. 이름만 다르지 같은 수다.

에크하르트-영¶

같은 수인 이유가 있다. L29의 에크하르트-영 정리는 랭크 $k$ 로 근사할 때의 최소 오차를 알려주었다.

\min_{\rank B \,=\, k}\lVert X_c - B\rVert_F^2 = \sum_{i > k}\sigma_i^2, \qquad B^{\star} = \sum_{i \le k}\sigma_i\vv{u}_i\vv{v}_i^{\mathsf T}

(26)

(26)의 최적해 $B^{\star}$ 가 바로 "주성분 $k$ 개만 남긴 것"이다. 양변을 $\lVert X_c\rVert_F^2$ 로 나누면 (25)의 비율이 나온다.

\frac{\lVert X_c - B^{\star}\rVert_F^2}{\lVert X_c\rVert_F^2} = \frac{\sum_{i>k}\sigma_i^2}{\sum_{i}\sigma_i^2} = 1 - \sum_{i \le k} p_i

(27)

(27)의 등식이 두 세계를 잇는다. 왼쪽은 "복원 오차"이고 오른쪽은 "설명하지 못한 분산"인데 같은 수다. 그래서 "주성분 3개로 97%를 설명한다"는 말과 "랭크 3으로 근사하면 에너지의 97%가 남는다"는 말은 완전히 같은 말이다.

왼쪽은 10개 변수로 잰 합성 자료의 스크리 도표다. 앞 세 개가 97%를 차지하고 나머지는
거의 평평한데, 자료를 만들 때 실제로 쓴 방향이 3개였다. 오른쪽은 랭크 k 로 실제
복원해 잰 오차와 설명분산으로 예측한 값을 겹쳐 그린 것이다. 두 곡선이 하나다. — Figure 3:왼쪽은 10개 변수로 잰 합성 자료의 스크리 도표다. 앞 세 개가 97%를 차지하고 나머지는 거의 평평한데, 자료를 만들 때 실제로 쓴 방향이 3개였다. 오른쪽은 랭크 $k$ 로 실제 복원해 잰 오차와 설명분산으로 예측한 값을 겹쳐 그린 것이다. 두 곡선이 하나다.

5. 공분산을 만들면 조건수가 제곱된다¶

이 강의를 한 문장으로 줄이면 이 절이다.

제곱된다¶

$C$ 의 고윳값이 $\sigma_i^2/(m-1)$ 이므로 조건수를 계산하면 $m-1$ 이 지워진다.

\kappa(C) = \frac{\lambda_1}{\lambda_n} = \frac{\sigma_1^2/(m-1)}{\sigma_n^2/(m-1)} = \left(\frac{\sigma_1}{\sigma_n}\right)^{\!2} = \kappa(X_c)^2

(28)

(28)의 등식은 L33에서 이미 만났다. 그때는 최소제곱의 정규방정식이 문제였고 지금은 PCA다. 문제는 다르지만 죄목은 같다. $X^{\mathsf T}X$ 라는 곱을 손으로 만든 것이다.

앵커로도 확인된다. $\kappa(X_c) = \sigma_1/\sigma_2 = 2\sqrt6/2 = \sqrt6 \approx 2.449$ 이고 $\kappa(C) = 6/1 = 6$ 이다. 정확히 제곱이다.

얼마나 손해인가¶

제곱이 왜 문제인지 자릿수로 따져 보자. 고윳값 계산에서 생기는 오차는 가장 큰 고윳값에 비례한다. 곧 $\lvert d\lambda_i\rvert \lesssim \epsilon_{\text{mach}}\lambda_1$ 이다. 우리가 원하는 것은 $\sigma_i$ 이므로 $\sigma_i = \sqrt{\lambda_i}$ 를 미분한다.

\sigma_i = \lambda_i^{1/2} \qquad\Longrightarrow\qquad \frac{d\sigma_i}{\sigma_i} = \frac{1}{2}\cdot\frac{d\lambda_i}{\lambda_i}

(29)

(29)의 관계는 $d\sigma_i = \tfrac12\lambda_i^{-1/2}d\lambda_i$ 를 $\sigma_i = \lambda_i^{1/2}$ 로 나눈 것이다. 여기에 오차의 크기를 넣는다.

\frac{\lvert d\sigma_i\rvert}{\sigma_i} \ \lesssim\ \frac{\epsilon_{\text{mach}}}{2}\cdot\frac{\lambda_1}{\lambda_i} = \frac{\epsilon_{\text{mach}}}{2}\left(\frac{\sigma_1}{\sigma_i}\right)^{\!2} \ \le\ \frac{\epsilon_{\text{mach}}\,\kappa(X_c)^2}{2}

(30)

(30)의 마지막 항이 손해의 크기다. 반면 $X_c$ 에 직접 SVD를 걸면 특이값의 오차는 $\sigma_1$ 에 비례할 뿐 제곱이 끼어들지 않는다.

\frac{\lvert d\sigma_i\rvert}{\sigma_i} \ \lesssim\ \epsilon_{\text{mach}}\cdot\frac{\sigma_1}{\sigma_i} \ \le\ \epsilon_{\text{mach}}\,\kappa(X_c)

(31)

(31)의 지수와 (30)의 지수가 1과 2로 갈린다. $\kappa = 10^6$ 이면 앞의 것은 10^-10 이고 뒤의 것은 10^-4 다. 여섯 자리가 날아간다.

갈리는 지점¶

배정밀도에서 $\epsilon_{\text{mach}} \approx 2.2\times10^{-16}$ 이다.

\kappa(X_c) = 10^{8} \qquad\Longrightarrow\qquad \kappa(C) = 10^{16} \ >\ \frac{1}{\epsilon_{\text{mach}}} \approx 4.5\times10^{15}

(32)

(32)의 지점을 넘으면 $\lambda_n$ 이 반올림 오차에 묻힌다. 그 결과가 꽤 사납다. 양의 준정부호 행렬인데 고윳값이 음수로 나온다. 수학적으로 불가능한 일이 계산에서는 일어나고, 그 뒤에 $\sqrt{\lambda_n}$ 을 하려 들면 그대로 무너진다.

실제로 재어 보면 이렇게 갈린다. 특이값을 아는 $400 \times 8$ 자료를 만들어 가장 작은 특이값을 얼마나 정확히 되찾는지 본 것이다.

$\kappa(X_c)$	$\kappa(C)$	공분산 경유	SVD 직접
10²	10⁴	$1.8\times10^{-13}$	$1.7\times10^{-16}$
10⁴	10⁸	$5.6\times10^{-10}$	$9.2\times10^{-14}$
10⁶	10¹²	$2.6\times10^{-5}$	$7.8\times10^{-13}$
10⁸	$7.6\times10^{15}$	$\mathbf{1.4\times10^{-1}}$	$1.1\times10^{-12}$
10¹⁰	$1.8\times10^{17}$	$\mathbf{1.0}$ (고윳값이 음수)	$2.7\times10^{-8}$

마지막 줄에서 공분산 경유의 상대오차가 1.0이다. 100% 틀렸다는 뜻이다. 같은 자료에서 SVD는 여덟 자리를 지키고 있다.

왼쪽에서 실측한 \kappa(C) 가 \kappa(X_c)^2 위에 정확히 올라앉다가, 초록 선
1/\epsilon_{\text{mach}} 를 넘는 순간부터 의미를 잃고 흩어진다. 오른쪽에서 빨간 점은
\epsilon\kappa^2 를 따라 올라가고 파란 점은 \epsilon\kappa 를 따라간다. — Figure 4:왼쪽에서 실측한 $\kappa(C)$ 가 $\kappa(X_c)^2$ 위에 정확히 올라앉다가, 초록 선 $1/\epsilon_{\text{mach}}$ 를 넘는 순간부터 의미를 잃고 흩어진다. 오른쪽에서 빨간 점은 $\epsilon\kappa^2$ 를 따라 올라가고 파란 점은 $\epsilon\kappa$ 를 따라간다.

단정밀도면 훨씬 빨리 온다¶

기계학습에서 흔히 쓰는 float32 는 $\epsilon_{\text{mach}} \approx 1.2\times10^{-7}$ 이다. (32)의 계산을 다시 하면 임계점이 이렇게 내려온다.

\frac{1}{\sqrt{\epsilon_{\text{mach}}}} = \frac{1}{\sqrt{1.19\times10^{-7}}} \approx 2900 \approx 10^{3.5}

(33)

(33)의 값은 배정밀도의 10⁸ 보다 네 자리 넘게 낮다. $200 \times 2$ 자료에서 $\sigma_1 = 1$ 로 고정하고 $\sigma_2$ 만 낮춰 가며 재 보면 경계가 정확히 그 자리에 있다.

$\kappa(X_c)$	공분산 경유	상대오차	SVD 직접
10³	$1.0033\times10^{-3}$	$0.3\%$	$1.0000\times10^{-3}$
$3\times10^{3}$	$2.9369\times10^{-4}$	$2.1\%$	$3.0000\times10^{-4}$
10⁴	$\mathbf{0}$	$\mathbf{100\%}$	$1.0000\times10^{-4}$

마지막 줄에서 $\lambda_2$ 가 반올림에 통째로 묻혀 0이 되었다. 그 뒤에 제곱근을 씌우면 $\sqrt0 = 0$ 이라 제2주성분의 크기가 0으로 보고된다. 조건수 10⁴ 는 실무에서 드물지도 않은 값이다.

6. 중심화를 빼먹으면¶

항등식 하나¶

평균을 빼지 않고 $X^{\mathsf T}X$ 를 만들면 무슨 일이 생기는지는 식으로 정확히 알 수 있다. (2)에서 $X = X_c + \vv{1}\boldsymbol{\mu}^{\mathsf T}$ 이므로 그대로 전개한다.

X^{\mathsf T}X = X_c^{\mathsf T}X_c + X_c^{\mathsf T}\vv{1}\boldsymbol{\mu}^{\mathsf T} + \boldsymbol{\mu}\vv{1}^{\mathsf T}X_c + \boldsymbol{\mu}\vv{1}^{\mathsf T}\vv{1}\boldsymbol{\mu}^{\mathsf T}

(34)

(34)의 네 항 가운데 가운데 둘이 사라진다. 중심화의 뜻이 $X_c^{\mathsf T}\vv{1} = \vv{0}$ 이기 때문이다. 그리고 $\vv{1}^{\mathsf T}\vv{1} = m$ 이므로 마지막 항이 정리된다.

X^{\mathsf T}X = X_c^{\mathsf T}X_c + m\,\boldsymbol{\mu}\boldsymbol{\mu}^{\mathsf T}

(35)

(35)의 항등식이 이 절의 전부다. 중심화를 빼먹는다는 것은 랭크 1 행렬 하나를 얹은 채로 푸는 것이다. 앵커로 검산해 보면 $X^{\mathsf T}X = \begin{pmatrix}100 & 68\\ 68 & 53\end{pmatrix}$ 이고 $X_c^{\mathsf T}X_c + 5\boldsymbol{\mu}\boldsymbol{\mu}^{\mathsf T}$ 도 같은 행렬이다.

얹힌 항이 무엇을 하는가¶

얹힌 항의 레일리 몫을 계산해 보자.

\vv{v}^{\mathsf T}X^{\mathsf T}X\vv{v} = \lVert X_c\vv{v}\rVert^2 + m\left(\boldsymbol{\mu}^{\mathsf T}\vv{v}\right)^2

(36)

(36)의 둘째 항은 $\vv{v}$ 가 $\boldsymbol{\mu}$ 방향일 때 $m\lVert\boldsymbol{\mu}\rVert^2$ 로 최대가 되고, $\boldsymbol{\mu}$ 에 수직일 때 0이 된다. 곧 평균 방향에만 $m\lVert\boldsymbol{\mu}\rVert^2$ 크기의 가짜 분산이 얹힌다.

그러면 승부는 뻔하다. 진짜 퍼짐의 크기는 $\lambda_1$ 정도인데 얹힌 것은 $m\lVert\boldsymbol{\mu}\rVert^2$ 이고, 표본이 많을수록 $m$ 이 커진다. $\lVert\boldsymbol{\mu}\rVert$ 가 퍼짐보다 크기만 하면 제1주성분은 평균 방향으로 끌려간다.

실험¶

평균이 $(8, 8)$ 이고 퍼짐은 그것과 수직인 방향인 자료 300개를 만들어 재어 보았다.

	제1주성분	특이값	설명분산 비율
중심화함	$(0.707,\ -0.707)$	$50.6,\ 8.5$	$97.3\%,\ 2.7\%$
중심화 안 함	$(0.709,\ 0.706)$	$196.8,\ 50.6$	$93.8\%,\ 6.2\%$

중심화하지 않은 제1주성분과 평균 방향의 각이 $0.009^\circ$ 다. 평균 방향 그 자체다. 그리고 진짜 제1주성분과의 각은 $89.87^\circ$ 로, 있을 수 있는 최악의 오답이다.

더 고약한 것은 아무 경고도 나지 않는다는 점이다. 설명분산 93.8%는 오히려 그럴듯해 보인다. 첫 특이값 196.8을 제곱하면 38747인데, $m\lVert\boldsymbol{\mu}\rVert^2 = 300 \times 128.9 = 38674$ 다. 거의 전부가 평균에서 온 것이다.

왼쪽에서 빨간 직선(중심화하지 않은 제1주성분)이 구름을 가로지르지 않고 원점에서
평균으로 뻗는다. 오른쪽 그림은 평균이 커질수록 각도 오차가 어떻게 변하는지 잰 것으로,
\lVert\boldsymbol{\mu}\rVert 가 퍼짐 3을 넘는 지점에서 곧장 90도로 튄다. — Figure 5:왼쪽에서 빨간 직선(중심화하지 않은 제1주성분)이 구름을 가로지르지 않고 **원점에서 평균으로 뻗는다.** 오른쪽 그림은 평균이 커질수록 각도 오차가 어떻게 변하는지 잰 것으로, $\lVert\boldsymbol{\mu}\rVert$ 가 퍼짐 3을 넘는 지점에서 곧장 90도로 튄다.

7. 무엇을 최소로 하는가¶

PCA는 회귀가 아니다¶

앵커에 세 개의 직선을 그어 보자. $S_{xx} = 20$ , $S_{yy} = 8$ , $S_{xy} = 8$ 이다.

\text{(y on x)}\ \frac{S_{xy}}{S_{xx}} = \frac{8}{20} = 0.4, \qquad \text{(PC1)}\ \frac{1}{2}, \qquad \text{(x on y)}\ \frac{S_{yy}}{S_{xy}} = \frac{8}{8} = 1

(37)

(37)의 세 기울기가 전부 다르다. 왼쪽은 L16의 최소제곱으로 $y$ 를 $x$ 로 설명한 것이고, 오른쪽은 역할을 바꿔 $x$ 를 $y$ 로 설명한 뒤 그림에서 읽은 기울기다. 가운데가 제1주성분이고 언제나 두 회귀선 사이에 있다.

무엇을 최소로 하느냐가 다르기 때문이다. 직선 방향을 $\vv{v}$ , 그 법선을 $\vv{n}$ 이라 하면

\sum_{i}\text{(수직거리)}^2 = \lVert X_c\vv{n}\rVert^2 = \sigma_2^2, \qquad \sum_{i}\text{(세로거리)}^2 = \lVert X_c(1,\,-k)^{\mathsf T}\rVert^2

(38)

가 된다. (38)의 왼쪽에서 수직거리 제곱합을 최소로 하는 것은 $\lVert X_c\vv{n}\rVert$ 을 최소로 하는 것이고, 그것은 $\vv{n} = \vv{q}_2$ 일 때 $\lambda_2$ 를 $m-1$ 배 한 값, 곧 $\sigma_2^2 = 4$ 다. PCA가 최소로 하는 것은 수직거리다. 반면 최소제곱은 세로거리를 최소로 하고 그 값은 4.8이다.

\text{PC1}:\ \sum d_\perp^2 = 4.00, \qquad \text{(y on x)}:\ \sum d_\perp^2 = 4.14

(39)

(39)의 두 값을 견주면 수직거리로 재었을 때 PC1이 이긴다. 세로거리로 재면 4.80 대 5.00으로 반대가 된다. 둘 다 옳고, 다른 문제를 푼 것이다.

같은 다섯 점에 그은 세 직선과, 각각이 짧게 만들려 애쓰는 거리다. 가운데는 수직으로
떨어뜨린 선분이고 오른쪽은 세로로 떨어뜨린 선분이다. — Figure 6:같은 다섯 점에 그은 세 직선과, 각각이 짧게 만들려 애쓰는 거리다. 가운데는 수직으로 떨어뜨린 선분이고 오른쪽은 세로로 떨어뜨린 선분이다.

눈금을 바꾸면 답이 바뀐다¶

변수 $j$ 의 단위를 바꾸는 것은 열에 상수를 곱하는 일이다.

X_c \longrightarrow X_c D \qquad\Longrightarrow\qquad C \longrightarrow DCD, \qquad D = \diag(d_1, \dots, d_n)

(40)

(40)의 변환은 $C$ 를 닮음이 아니라 합동으로 바꾼다(L32). 닮음이라면 고윳값이 보존될 텐데 합동은 부호만 보존한다. 그러니 주성분이 바뀌는 것이 당연하다.

앵커에서 둘째 변수를 10배 해 보자.

D = \diag(1, 10) \ \Longrightarrow\ DCD = \begin{bmatrix} 5 & 20 \\ 20 & 200 \end{bmatrix}, \qquad \lambda = 202.03,\ 2.97

(41)

(41)의 새 주성분은 $(0.101,\ 0.995)$ 로 거의 $y$ 축이고, 설명분산은 98.6%로 치솟는다. 원래 $26.57^\circ$ 였던 것이 $84.20^\circ$ 가 되었다. 자료는 하나도 달라지지 않았는데 답이 달라졌다.

그래서 단위가 제각각인 변수를 다룰 때는 각 열을 표준편차로 나눈 뒤 PCA를 건다. 그러면 $C$ 대신 상관행렬이 들어간다.

R = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}, \qquad \rho = \frac{S_{xy}}{\sqrt{S_{xx}S_{yy}}} = \frac{8}{\sqrt{160}} \approx 0.6325, \qquad \lambda = 1 \pm \rho

(42)

(42)의 고윳값이 $1+\rho$ 와 $1-\rho$ 인 것은 $R$ 의 대각합이 2이고 행렬식이 $1-\rho^2$ 이기 때문이다. 그리고 $2 \times 2$ 상관행렬의 제1주성분은 $\rho$ 가 무엇이든 언제나 $(1,1)/\sqrt2$ 다. 표준화가 선택이 아니라 판단임을 보여 주는 예다.

8. 자주 하는 오해¶

① 공분산행렬을 만들어 고유분해하는 것이 정석이다. 아니다. 종이 위에서는 그것이 정의지만 계산에서는 조건수를 제곱해 버린다. 5절이 답이다. $\kappa(X_c)$ 가 10⁸ 을 넘으면 작은 고윳값이 음수로 나오기까지 한다. 데이터 행렬에 SVD를 거는 것이 정석이다.

② 주성분은 원래 변수 중에서 고른 것이다. 아니다. 주성분은 모든 변수의 선형결합이다. $\vv{q}_1$ 의 성분이 각 변수의 기여도이고, 성분 대부분이 0이 아니므로 "몇 번 변수를 골랐다"고 말할 수 없다. 변수를 실제로 골라야 한다면 그것은 다른 문제이고 다른 도구가 필요하다.

③ 설명분산이 높으면 좋은 주성분이다. 아니다. 6절의 실험이 반례다. 중심화를 빼먹으면 설명분산이 오히려 올라간다. 평균이라는 아무 정보도 없는 방향이 분산을 독차지하기 때문이다. 설명분산은 "잘 맞았다"가 아니라 "이 방향이 길다"는 뜻일 뿐이다.

④ 주성분의 부호에 뜻이 있다. 아니다. $\vv{q}$ 가 고유벡터면 $-\vv{q}$ 도 고유벡터이고, 두 답 사이에 우열이 없다. 같은 자료를 다른 라이브러리에 넣으면 부호가 뒤집혀 나오는 일이 흔하다. 부호가 아니라 축과 점수의 상대적 배치를 읽어야 한다.

⑤ PCA로 줄인 좌표가 회귀에 가장 좋은 변수다. 아니다. PCA는 목표변수를 본 적이 없다. 분산이 큰 방향과 목표를 잘 설명하는 방향은 아무 관계가 없고, 목표에 중요한 정보가 마침 분산이 작은 방향에 들어 있으면 그것을 버리게 된다. 목표를 함께 보고 싶다면 부분최소제곱 같은 다른 방법이 있다.

⑥ 주성분이 직교하니 원래 변수들도 독립이 된다. 아니다. 새 좌표의 상관이 0인 것은 맞다. (23)에서 $Y^{\mathsf T}Y$ 가 대각행렬이기 때문이다. 그러나 상관이 0인 것과 통계적으로 독립인 것은 다르다. 공분산은 이차 정보만 보므로, 비선형으로 묶여 있는 관계는 PCA가 보지 못한다.

마치며...¶

이번 강의에서는 주성분분석이 무엇을 하는 계산인지 처음부터 끝까지 따라가 보았다.

분산이 최대인 방향을 찾는 문제는 레일리 몫의 최대화이고, 답은 최대 고윳값의 고유벡터다. 증명은 정규직교 기저로 갈아 끼우는 것이 전부였다(L25).
$C = X_c^{\mathsf T}X_c/(m-1)$ 은 대칭이고 양의 준정부호이므로 고윳값이 실수이고 음수가 아니다(L27). 두 번째 주성분이 첫 번째와 직교하는 것은 공짜로 따라온다.
$X_c = U\Sigma V^{\mathsf T}$ 를 걸면 $V$ 가 주성분이고 $\lambda_i = \sigma_i^2/(m-1)$ 이다. 점수는 $U\Sigma$ 라서 곱셈조차 필요 없다(L29).
설명분산 비율 $\sigma_i^2/\sum\sigma_j^2$ 은 L29의 에너지와 같은 수이고, 주성분 $k$ 개를 남기는 것은 에크하르트-영의 랭크 $k$ 근사 그 자체다.
공분산행렬을 만들면 조건수가 제곱된다. 그래서 라이브러리는 그것을 만들지 않는다.
중심화를 빼먹으면 제1주성분이 평균 방향이 된다. 그러고도 설명분산은 더 좋아 보인다.

돌아보면 이 강의에는 새 정리가 하나도 없었다. 스펙트럼 정리는 L25에 있었고, 준정부호는 L27에 있었고, SVD와 에너지는 L29에 있었고, 조건수의 제곱은 L33에 있었고, 투영은 L15에 있었고, 합동과 닮음의 구별은 L32에 있었다. 여러분이 한 줄로 호출하던 함수 안에서 벌어지는 일이, 지난 서른네 강의의 내용 그대로였다.

다음은 보강 2, 수치선형대수 입문이다. 이번 강의에서 "옳은 공식과 쓰는 공식은 다르다"는 교훈이 다섯 번째로 나왔다. 다섯 번이면 우연이 아니다. 그 다섯을 조건수라는 하나의 개념으로 묶고, 큰 행렬의 고윳값을 실제로 어떻게 구하는지 살펴보자. 그리고 L22에서 벡터에 행렬을 자꾸 곱하던 그 애니메이션이 사실은 알고리즘이었다는 것도 그때 밝힌다.

이번 강의의 내용을 파이썬으로 확인해 보려면 L35 실습 노트북으로 넘어가면 된다. 주성분 함수를 스무 줄로 직접 만들어 앵커의 손계산과 대조하고, 자를 돌려 가며 분산이 6과 1 사이를 오가는 것을 보고, 조건수를 키워 가며 공분산 경유와 SVD 직접이 어디서 갈라지는지 직접 잴 수 있다.