주성분 분석(Principal Component Analysis, PCA) 은 데이터 집합 내에 존재하는 각 데이터의 차이를 가장 잘 나타내 주는 요소를 찾아 내는 방법이다. 그 구체적인 의미 및 원리, 그리고 이를 C++로 구현해 본다. 주성분 분석에 대한 개념을 위하여 그 의미를 알 수 있는 간단한 몇 가지 예제에 이어 수학적인 원리를 간략히 살펴 본 후, 그것에 대한 C++을 이용한 구현을 다룬다. 실제 C++ 코드는 후반부에 있으므로 그 부분만 필요할 경우 바로 건너 뛰어도 된다. 하지만 될 수 있으면 이 글 전체를 반복해서 읽기를 제안한다. 또한, 독자가 대학 1, 2학년 정도의 기본적인 선형 대수 내용은 알고 있다고 가정한다.
'주성분'에 대한 개념에 대한 감을 잡기 위하여 다음과 같은 상황을 생각해 보자. 도시의 특성을 주성분 분석으로 분석해 본 글도 참고하자.
강남역 8번 출구 앞을 지나는 자동차 1만대에 대하여 차의 가격, 무게, 색상, 최대탑승인원수, 배기용량, 제조회사 본사의 위도와 경도, 운전하는 사람의 키를 측정한 데이터, 자동차 문의 수, 자동차 바퀴의 수, 자동차 핸들의 모양, 차체의 높이.
위에서 측정한 특성 중 어느 변수가 1만개의 데이터를 가장 잘 나타낼 수 있을까? 즉, 각 데이터의 차이가 가장 잘 두드러지게 하는 특성은 무엇일까? 대부분의 자동차가 바퀴가 4개이기 때문에 측정한 1만개의 데이터의 대부분이 4개를 갖을 것이며, 따라서 자동차 바퀴의 수는 1만개의 데이터 각각을 구분하기에는 좋은 변수는 아니다. 유사한 이유로 자동차 핸들의 모양, 문의 수 등으로는 각 자동차를 다른 자동차와 구분하기 좋은 특성은 아닐 것이다, 동일한 값을 갖는 자동차들이 많기 때문에. 반면 가격/무게/차체의 높이 등은 각 자동차가 서로 다른 값을 갖는 경우가 많기 때문에 이 값들의 적절한 조합은 각 차를 다른 차와 구별시켜 주는 특징이 될 수 있을 것이다. 주성분 분석은 이러한 특징들을 선형 결합시킬 때 데이터들이 가장 다양한 값을 갖도록 하는 특징의 가중치를 찾아 주는 방법이다.
위의 경우 핸들의 모양이 동그라면 1, (==) 모양이면 2, 그외의 모양이면 3을 준다고 해보자. 그러면 1만개의 데이터 대부분이 이 값은 1일 것이다. (핸들 모양, 바퀴의 수)를 고려하면 대부분 (1, 4) 의 값에 대부분의 데이터가 위치할 것이다. 즉, 이 변수들로는 자동차들을 구별하기 힘들다, 대부분이 (1, 4) 에 속하기 때문에. 하지만 (차의 가격, 무게, 차체의 높이) 를 측정하면 1만개의 데이터가 다양한 값을 갖게 되며, 따라서 이 경우 각 차의 특성이 잘 드러나며 이 값으로 각각의 차를 구별하기 쉬워질 것이다.
주성분 분석은 이처럼 각 데이터의 여러 특징값들을 이용하여 데이터들이 다양한 값을 갖도록 해주는 방법이다. 주성분 분석에 의하여 각 데이터는 새로운 값을 갖게 되는데, 데이터들의 새로운 값은 가장 다양한 값을 갖게 된다. 1만개의 데이터 였다면 자동차 바퀴의 수처럼 1만개의 값 중 대부분이 4가 아니라 차체 높이 x 무게 처럼 아주 고른 값을 갖게 되는 것이다.
위에서 사용한 특징인 차체 무게, 높이 등은 데이터의 변수가 된다. 그러면 주성분 분석은 각 변수의 가중치를 찾는 방법이다. 이 때, 이 가중치는 데이터의 분산이 가장 크도록 찾아 진다. 즉, 다음과 같다.
위에서 n 은 1만개의 데이터를 측정했으므로 1만이다. p 는 변수(특성)의 개수이다. 각각의 개별 데이터는 주성분 분석에 의하여 찾아진 각 변수에 대한 가중치를 이용하여 선형 결합된다. 즉,
t1 = x11 · w1 + x12 · w2 + ... + x1p · wp
t2 = x21 · w1 + x22 · w2 + ... + x2p · wp
...
tn = xn1 · w1 + xn2 · w2 + ... + xnp · wp
처럼 각 데이터는 주성분 분석에 의하여 찾아진 가중치 w=(w1, ..., wp) 에 의하여 새로운 데이터 t=(t1, ..., tn) 으로 변형된다1. 이 때, 각 데이터는 처음에는 p-차원의 데이터였으나 w에 의하여 1차원 값인 스칼라로 변경되었다는 점을 주의 깊게 보아야 한다. 주성분 분석의 핵심은 이렇게 찾아진 가중치 w 에 의해 변환된 새로운 값인 t1, ..., tn 이 최대의 분산을 갖도록 w 를 찾는다는 점2이다. 최대의 분산을 갖는 다는 것의 의미는 데이터 값이 가장 다양하게 된다는 것으로, 1만개의 데이터가 될 수 있으면 같은 값을 갖는 데이터가 없게 된다는 것과 일맥 상통한다. 만약 자동차의 바퀴 수, 로만 생각해 보면 1만개 중 대략 80%~90%가 4를 갖게 될 것이다. 그러나 차체의 높이와 무게는 바퀴의 수보다는 다양할 것이다. 따라서 이 경우 주성분 분석을 하게 되면 바퀴에 해당하는 가중치는 차체의 높이와 무게에 해당하는 가중치보다 작게 나올 것이다.
이제 주성분 분석의 수학적 원리를 살펴 보자. 너무 깊게 들어 가지 않고 간략히 알아 보기 위하여 이 글에서는 영문 위키피디아의 주성분 분석 페이지(링크)에 있는 내용만으로 설명하고자 한다. 아래 표기 역시 그 페이지의 표기법을 따른다. 전술했듯 PCA 는 주어진 데이터 X 에 대하여
T = Xw
인 T=(t1, ..., tm) 에 대하여 각 ti 가 최대의 분산을 갖도록 w 를 찾는 방법이다. 이 때 X 는 다음과 같다.
앞에서는 w 를 마치 1개의 열벡터인 것처럼 말했으나 실제로는 자신이 찾고자 하는 개수 m 개인 열로 된 벡터이다. 그래서 T=Xw 를 행렬식으로 표현하면 다음과 같다.
이 때, 일반적으로 수학에서 괄호로 첨자를 표시할 경우 '순서'를 고려한 것을 가리키는 용례에 따라 w(1)은 제1 주성분을 찾기 위한 가중치이고 w(2)는 두 번째 주성분을 찾기 위한 가중치 등을 의미한다. 참고로 영문 위키에서 설명하는 표기법은 다음을 따른다.
살짝 번거롭게 되어 있는데, 차근차근 보자면 다음과 같다.
제 1 주성분은 다음과 같이 표현된다.
제 1 주성분 t1=(t1(1), t2(1), ..., tn(1)) 의 분산을 최대화 시키도록 w(1) = (w1(1), w2(1), ..., wp(1)) 을 찾는 것이다. 마찬가지로 그 다음의 제 2 주성분은 다음과 같이 표현된다.
이런 식으로 원하는 m 개의 주성분을 찾을 수 있다. 일반적으로 그래프로 표현하기 위해 제 1 과 제 2 주성분을 찾아 x 와 y 축으로 표현하거나 아니면 제 3 주성분까지 찾아서 z 축으로 표현해서 데이터를 시각화하곤 한다. 또한, 각 주성분이 원래의 데이터의 변화량의 몇 %를 대변하는지를 계산할 수도 있다. 또한, 각 w(i) 는 벡터의 길이가 1 인 unit vector 라는 제약 조건을 둔다.
제 1 주성분을 찾아 보자. 즉, 이제 제 1 주성분 t1=(t1(1), t2(1), ..., tn(1)) 의 분산을 최대화 시키는 w(1) = (w1(1), w2(1), ..., wp(1)) 을 찾는 것이 목표이다. 이것을 다르게 표현하면
Var(t1) = Var(t1(1), t2(1), ..., tn(1))
을 최대화 시키는 w(1) 을 찾는 것이 목표이다. 이것을 다시 좀 더 수학적으로 표현하면,
w(1) = arg max Var(t1)
||w||=1
이 된다. 즉, Var(t1) 이 최대가 되게 하는 열벡터 w 가 w(1) 이 되는 것이다. 이 때, 열벡터 w 는 그 크기가 1 인 unit vector 이다. 이것은 단순히 수학적 표기일 뿐이다. 우리는 분산의 공식이
Var(Y) = E[Y2] - (E[Y])2
임을 알고 있다. 따라서 Var(t1) = E[t12] - (E[t1])2 이다. 여기서 X 의 각 열은 평균이 0 으로 맞춰졌기 때문에 E[t1] 은 0 이 된다. 따라서
Var(t1) = E[t12]
이고, 결국
w(1) = arg max E[t12]
||w||=1
이 된다. 어차피 행 수는 n 개로 고정된 것이므로 평균을 구할 때의 1/n 은 고정된 값이므로 위의 식을 만족시키는 w(1) 은 다음과 같이 되는 것이다.
즉, 원래 arg max { 1/n sigma ...} 형식인데 1/n 은 고정된 값으로 w 에 의해 영향을 받지 않으므로 위에서 1/n 은 사라진 것이다. 그 다음부터는 표기법만 바꾼 것이므로 별도의 설명은 하지 않는다. 결론적으로
을 얻는다. 이 때, ||w||=1 이 제약 조건이있고, ||w||=wTw=wwT이므로 wTXTXw 와 wTXTXw÷wTw는 같은 값이 된다. 왜 굳이 이렇게 분모에 값이 1 인 wTw 를 표기했냐면 다음과 같이 정의되는 Rayleigh quotient 형식으로 만들기 위해서이다.
x* 는 x 의 conjugate transpose matrix 이고, M 은 complex Hermitian matrix 이다. 여기서 중요한 점은 Rayleigh quotient 의 최소값은 M의 eigenvalue의 최소값이랑 같고 그 때 x 는 그 최소값의 eigenvalue에 해당하는 eigenvector 라는 점이다. 마찬가지로 Rayleigh quotient 의 최대값은 M 의 eigenvalue의 최대값이랑 같고, 그 때 x의 값은 그 최대값의 eigenvalue에 해당하는 eigenvector 값이 된다. 이 정리를 다음의 w(1) 의 관점에서 살펴 보자.
행렬의 요소에 복소수가 없고 모든 요소가 실수인 경우 conjugate transpose 는 그냥 transpose 가 된다. 또한, Hermitian matrix 는 원래의 matrix A와 그것의 conjugate transpose 가 같은 matrix 로 정의된다. 이 때, 모든 요소가 실수로만 되어 있다면 Hermitian matrix 는 원래의 matrix A 와 그것의 transpose 인 AT가 같은 matrix 로 정의된다. 이 때, 위의 식에서 A = XTX 로 두면, AT= XTX 이므로 A 즉 XTX 는 Hermitian matrix 가 된다. 따라서 위에서 표기된 w(1) 은 Rayleigh quotient 가 되는 것이다. 따라서 위의 최대값은 XTX 의 eigenvalue의 최대값이 되고 그 때의 w 는 그 최대의 eigenvalue에 해당하는 eigenvector 가 된다.
제 2 주성분에 해당하는 가중치 w(2) 값은 원래의 matrix X 에서 제 1 주성분을 뺀 후 제 1 주성분을 구하는 것과 같이 구하면 된다고 한다. 그런데 결론적으로 w(2) 은 XTX 의 두 번째 크기의 eigenvalue에 해당하는 eigenvector 가 된다고 한다. 그 나머지도 같다고 한다. 즉, w(k)는 XTX 의 k-번째 크기의 eigenvalue 에 해당하는 eigenvector 라는 것이다.
이제 주성분 분석을 C++로 구현해 보자. 의외로 간단하다. 여기서, 위 내용을 다시 한 번 정리해 보자. 다음과 같이 데이터가 주어져 있다.
원본 데이터에서 컬럼 벡터의 평균을 0 으로 맞추어 준다. 그 후 주성분 분석은 XTX 의 eigenvector 의 크기 순으로 정렬하여 각각에 해당하는 eigenvector 를 구하면 그 순서에 해당하는 가중치가 된다. 실제 k-번째 주성분은 원본 데이터를 k-번째 가중치에 내적하면 얻어 진다.
1 #include <pca.h>
2 #include <Eigen/Dense>
3
4 #include <numeric>
5
6
7
8 // row of [_m] : variable.
9 // col of [_m] : observation.
10 // The meanining of row and col are different from those of Matlab.
11 bool pcapc12(const std::vector<std::vector<double> >& _m, std::vector<double>* const w1, std::vector<double>* const w2){
12 return pcapc12(_m, w1, w2, NULL, NULL);
13 }
14 bool pcapc12(const std::vector<std::vector<double> >& _m, std::vector<double>* const w1, std::vector<double>* const w2, double* const pc1explained, double* const pc2explained)
15 {
16 int M = (int)(_m.size());
17 const int N = (int)(_m[0].size());
18 Eigen::MatrixXd centered(N, M);
19 int m = 0; int n = 0;
20 for(m = 0; m<M; m++){
21 double sum = std::accumulate(_m[m].begin(), _m[m].end(), 0.0);
22 double mean = sum/N;
23 for(n = 0; n<N; n++){
24 double value = _m[m][n];
25 centered(n, m) = (value - mean);
26 }
27 }
28
29 Eigen::MatrixXd cov = centered.adjoint() * centered;
30
31 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(cov);
32 Eigen::VectorXd normalizedEigenValues = eig.eigenvalues() / eig.eigenvalues().sum();
33 if(pc1explained != NULL) *pc1explained = normalizedEigenValues[M-1];
34 if(pc2explained != NULL) *pc2explained = normalizedEigenValues[M-2];
35
36 Eigen::MatrixXd evecs = eig.eigenvectors();
37 Eigen::MatrixXd pcaTransform = evecs.rightCols(2);
38 for(m = 0; m<M; m++){
39 if(w1 != NULL) w1->push_back(pcaTransform(m, 1));
40 if(w2 != NULL) w2->push_back(pcaTransform(m, 0));
41 }
42
43 return true;
41 }
이 때, eigenvalue와 eigenvector 를 구하기 위해 외부 라이브러리인 Eigen (링크)을 사용했으므로 필요시 설치해 준다. 위 코드는 matlab 이나 위에서 설명한 형식과 달리 행벡터가 변수이고 열벡터가 데이터(개별 데이터 혹은 반복실험)이므로 그 점을 주의한다3. 또한 입력으로 받은 데이터는 평균이 0 으로 되어 있지 않다고 가정, 20~27 번 줄에서 평균을 0 이 되도록 해주고 있다. 나머지는 Eigen 설명서를 일어 보면 쉽게 알 수 있다. 즉. 36 번 줄에서 eigenvecor를 구해서 37번째 줄에서 그 중 두 개를 따온다. 그것이 곧 w(1) 과 w(2) 이 되는 것이다. 개인적으로 제 1, 제 2 주성분만 필요했기 때문에 위처럼 했었으며, 더 필요한 사람은 그에 맞게 코드를 수정하면 된다. 또한, 33, 34 부분은 각 주성분에 의해 실제 데이터의 변량(분산)이 어느 정도 반영되었는지를 나타내는 값이다. 위 코드를 사용해서 실제로 값을 구해 보자. 다음과 같이 tab으로 분리된 텍스트 파일을 다루어 보자. 위 코드의 주석에 있듯 각 행이 변수이고 각 열이 데이터이다. 즉, 4개의 변수로 되어 있는 데이터이다.
7 1 11 11 7 11 3 1 2 21 1 11 10
26 29 56 31 52 55 71 31 54 47 40 66 68
6 15 8 8 6 9 17 22 18 4 23 9 8
60 52 20 47 33 22 6 44 22 26 34 12 12
util.h 와 reader.h 는 개인적인 라이브러리이고, 9~16 번 줄에 의해 m 이 설정되므로 이 부분은 각자 환경에 맞게 적당히 작업한다. 즉, text file 을 한 줄씩 읽어 와서 tab 으로 자른 후 double 형으로 바꾸어서 std::vector 에 집어 넣는 것. 그 후 실행을 시키면
처럼 나온다. matlab 으로 구해 보면
이 나온다. 이 때, matlab에 사용하는 matrix 형태와 위 C++ 코드의 행/열이 바뀌었기 때문에 matlab에서는 치환을 해서 실행을 시켰다. 변수가 4개 이므로 최대 4개의 주성분이 가능하며 matlab은 이것을 전부 구해 주었다. 각 주성분은 방향이 다를 수도 있긴 하지만 그것은 크게 상관은 없으며 그것을 제외하면 같은 값이 나온 것을 확인할 수 있다. 위에서 4 번째 변수, 즉, 4번째 행의 값의 가중치가 가장 큰 것을 알 수 있다. 즉, 4번째 변수의 값이 각 데이터의 차이를 가장 잘 나타내 주는 변수라는 것을 알 수 있다. 이번에는 다음의 데이터를 분석해 보자.
두 번째 변수, 즉 두 번째 행은 모든 데이터가 같다. 따라서 두 번째 변수는 데이터들 간의 차이를 나타내기에는 전혀 필요 없는 값이다. 또한 5번째 변수 역시 거의 비슷하다. 이 변수 역시 데이터 간의 차이를 표현하기에는 그리 좋은 변수는 아니다. 따라서 두 번째 변수에 대한 가중치는 0.0 이 나올 것이고, 다섯 번째 변수에 대한 가중치는 매우 작게 나올 것이라고 예측할 수 있다. 실제 결과를 보면 다음과 같다.
예상한 대로 두 번째 변수에 대한 가중치는 10의 -16승 이 나왔다. 이것은 부동 소수점 표기에 의한 값일 뿐이며 이 경우 실제로는 0 이라고 간주할 수 있다. 또한, 예상한 것과 같이 다섯 번째 변수에 대한 가중치는 -0.0179924 나 -0.00375559 처럼 매우 작은 값이 나온 것을 알 수 있다.
위 코드는 논문 Oncopression: gene expression compendium for cancer with matched normal tissues, Bioinformatics, Volume 33, Issue 13, 1 July 2017, Pages 2068–2070 에 사용된 Fig. 1. B 인 다음 그림을 만드는데 사용하였다4. 즉, 실제로 잘 작동하는 코드이다. 아래 그림은 8개의 조직에서 총 1만 8천여 개의 유전자의 발현양을 측정한 데이터를 표현한 것이다. 원래라면 총 1만 8천여 차원에 데이터를 표현해야 하는데 그것은 시각적으로 불가능하므로, 각 데이터의 차이를 가장 잘 나타내 주는 유전자의 선형 조합에 대한 가중치를 위의 코드로 찾아서 그 선형 결합으로 찾은 제 1, 제 2 주성분을 x축과 y 축으로 표현한 것이다.
- 굵은 글씨로 되어 있는 것은 벡터 혹은 행렬이다. [본문으로]
- 무작정 분산을 제일 크게 하지는 않는다. 아래 나오겠지만 X 의 열벡터의 평균이 0 이고 가중치의 내적이 1 이 라는 제약조건이 있기 때문에 '분산을 최대화한다'라는 것이 가능한 것이다. 이러한 제약조건 없이 분산을 무작정 크게 하려면 그냥 가중치를 무작정 키우면 될 것이다. 즉, 이 경우 해가 없다. [본문으로]
- 위 코드는 인터넷의 어느 코드를 참조한 것인데 작성한 지가 오래 돼서 어느 정도 참고한 것인지 기억이 잘 안 난다, >.<"";;; [본문으로]
- 실제 논문에 들어 간 그림은 이 그림을 조금 변형한 것이다. [본문으로]
'컴퓨터 > 수학이랑' 카테고리의 다른 글
| 주성분 분석으로 알아 보는 도시 통계 (2) | 2018.05.15 |
|---|---|
| 상관계수2 (0) | 2011.08.30 |
| Fisher's linear discriminant 구현 (4) | 2010.12.03 |
| Fisher's linear discriminant 원리 (17) | 2010.11.30 |
| 통계 분석 - 평균 비교에 대하여 (0) | 2010.10.21 |
