Fisher's linear discriminant 구현

지난 번에 알아 본 Fisher's Linear Discriminant 의 원리를 이제 C++로 구현해 보자. 또한 지난 번 글에서 matlab 으로 테스트 해 보았던 m 코드의 일부를 약간 설명한다. Matlab 에서의 matrix 를 다루는 방식이 FLD 원리 편에서 본 matrix 와 조금은 다르기 때문에 약간의 설명이 필요하다. 또한 matrix inverse 를 구하는 것을 직접 구현하기에는 무리가 따르므로 다른 사람이 구현해 놓은 것을 사용한다^[각주:1].

우리가 직접적으로 계산해야 하는 것은 다음의 수식이다.

w는 벡터이고, 우리는 이 벡터의 '방향'만이 필요하므로 위에서 r 은 1로 둔다. S_w 와 m₁, m₂ 는 다음과 같다.

우선 m₁ 과 m₂ 부터 구해 보자. 실제로 우리가 다룰 data 는 위에서 x_n 가 column vector 로 된 전체 matrix 의 형태를 띄게 된다. 즉, 다음과 같다.

각 row는 각 변수에 관한 것이 되고, 각 column 은 각 측정에 관한 것이 된다. 즉, 만약 변수가 10 였다면 X 의 행의 수가 10이 되고, 그렇게 10개의 변수를 100 번 측정했다면 X 의 열의 길이가 100 이 되는 것이다.

위와 같다. 따라서, 위에서 의미하는 mean vector는 다음과 같이 구해야 한다.

그런데 matlab에서의 mean 함수는 각각의 column 을 average 시켜서 1-by-N 개의 row vector 를 만드는 함수이다. 다음과 같다.

그래서 지난 번의 m 파일에서는 다음과 같이, transpose를 시켜서 mean 에 넣은 후 결과를 다시 transpose 를 시켰던 것이다.

m1 = mean(x1')';
m2 = mean(x2')';

위 상황을 C++ 로 표현하자면 다음과 같다.

// If [v] is N-by-M matrix, [dest] is N-by-1 matrix.

// (k,1) element of [dest] is an average of elements of k-th row of [v].

bool get_mean_vector(const std::vector<std::vector<double> >& v, std::vector<double>* dest)

{

int N = (int)(v.size());

int M = (int)(v[0].size());

dest->clear();

std::vector<std::vector<double> >::const_iterator pos = v.begin();

for(; pos != v.end(); pos++){

//if((int)pos->size() != M) return false;

double sum = std::accumulate(pos->begin(), pos->end(), 0.0);

dest->push_back(sum/M);

}

return true;

}

(에러 처리는 역시 하지 않음. template 을 이용한 일반화는 하지 않는다. ) 매우 직관적임을 알 수 있다. STL 에 익숙하지 않으면 std::vector<std::vector<double> > 가 약간 혼동될 수 있으므로 간단히 설명한다. std::vector 는 template 에 기반한, 추상화된 dynamic array 이다. generic 의 편리함에 의하여, std::vector 는 요소로 그 어떤 type 도 받을 수 있으므로 std::vector<double> 을 요소를 갖는다. 따라서 위의 [v] 는 2차원 배열이라 생각할 수 있다. 또한, 위에서 설명한 바에 따라 [dest] 는 column vector 이어야 한다. 그러나 현재 std::vector 는 1차원 배열이므로 row 이거나 column vector 냐에 상관없는 형태를 갖게 된다. 뒤에서 보는 TNT::Array2D 형태는 이럴 경우, 1-by-N, 또는 N-by-1 처럼 명시적으로 표시해 주어야 한다.

이제 S_w 를 구해 보자. S_w 는 다시 보면 다음과 같다.

위에서 보면, 앞의 항과 뒤의 항은 C1 과 C2 인 것만 다를 뿐 같은 식으로 계산이 되고 있다. 또한 column vector 와 row vector 를 곱하고 있으므로 matrix 가 된다. 즉, d-by-1 x 1-by-d = d-by-d matrix 가 되는 것이다. 따라서 이 부분은 외부 라이브러리를 사용하여 계산한다. matrix 연산에 필요한 라이브러리는 약간 검색하여 TNT (Template Numerical Toolkit) 를 이용한다. TNT 는 실제로 구현 파일을 모두 제공해 주기 때문에 링크시킬 때 object 파일을 link 시키는 것이 아니라 source file 에 include 를 해서 사용하게 된다. TNT 의 구체적인 사용법은 홈페이지에 있으므로 참고하고, TNT를 이용하여 within-class-variance 를 구하는 routine 을 C++로 작성하면 다음과 같다.

// If [v] is N-by-M matrix and [m] is 1-by-N, then [dest] should have been initialized as N-by-N matrix.

// Actually [m] is N-by-1 matrix in mathematical formula.

// [dest] should be initalized by caller in advance.

bool get_within_class_variance(const std::vector<std::vector<double> >& v, const std::vector<double>& m, TNT::Array2D<double>* dest)

{

int N = (int)(v.size());

int M = (int)(v[0].size());

int row = 0; int col = 0;

for(col = 0; col < M; col++){

TNT::Array2D<double> column(N,1);

TNT::Array2D<double> column_t(1,N);

for(row = 0; row < N; row++){

column[row][0] = (v[row][col] - m[row]);

column_t[0][row] = (v[row][col] - m[row]);

}

TNT::Array2D<double> m = TNT::matmult(column, column_t);

TNT::operator+=(*dest,m);

}

return true;

}

별로 어려운 부분은 없다. 수식을 그대로 옮겼을 뿐이다.

이제 Sw 의 역행렬(inverse matrix) S_w^-1 를 구해 보자. TNT 에는 역행렬을 구해주는 함수가 직접적으로 제공되지는 않지만 그대신 TNT와 같이 관리되는 JAMA toolkit 에 solve 라는 함수가 제공이 된다. 이 함수는 Ax = B 형태의 방정식에서 x 를 구해주는 함수이다. 역행렬은 B 를 identity matrix 로 두고 문제를 풀면 되는 것이므로 다음과 같이 구현한다.

TNT::Array2D<double> I(N,N);

make_identity_matrix(N, &I);

JAMA::LU<double> Sw(Sw_);

double det = Sw.det();

if(det == 0) return false;

TNT::Array2D<double> Sw_inv = Sw.solve(I);

determinant 가 0 이면 inverse 를 구할 수 없으므로 이것만은 확인을 해주고 넘어 간다. 또한, identity matrix 를 만드는 함수는 다음과 같이 간단히 구현할 수 있다.

// [dest] should be initialized by caller in advance.

bool make_identity_matrix(int n, TNT::Array2D<double>* dest)

{

if(n < 1) return false;

int i = 0; int j = 0;

for(i = 0; i<n; i++){

for(j = 0; j<n; j++){

(*dest)[i][j] = (int)(i == j);

}

return true;

}

그 후, 실제로 [w]를 사용하여 data들을 직선으로 projection 시켰을 때의 값도 반환하면 써먹을 수 있으므로 이 값을 다음과 같이 구한다.

// [dest].T = [w].T * [m].

bool get_projected_pos(const std::vector<std::vector<double> >& m, const std::vector<double>& w, std::vector<double>* dest)

{

if(dest == NULL) return true;

int N = (int)(m.size()); // row

int M = (int)(m[0].size()); // col

int row = 0; int col = 0;

for(col = 0; col < M; col++){

std::vector<double> column;

for(row = 0; row < N; row++){

column.push_back(m[row][col]);

}

double v = std::inner_product(w.begin(), w.end(), column.begin(), 0.0);

dest->push_back(v);

}

return true;

}

[m]에서 각 column 을 뽑아 내서, [w]와의 내적을 계산하면 된다. 수학적으로는 [w]는 row vector 이지만 현재 우리는 std::vector 에 집어 넣었으므로 위처럼 inner_product 함수를 사용해서 계산하면 간단히 끝난다.

중요한 부분들의 설명이 끝났다. 실제로 위와 같은 작업들을 한 곳에 몰아서 하는 함수는 다음과 같다.

// [class1] and [class2] are matrices of which

// each row is for each metric and each column is for each data.

// [prj_i].T = [w_].T * [class_i] where .T means transpose.

// Mathematically [w_] is column vector and [prj_i] are row vectors.

bool FLD(const std::vector<std::vector<double> >& class1, const std::vector<std::vector<double> >& class2, std::vector<double>* prj1, std::vector<double>* prj2, std::vector<double>* w_)

{

std::vector<double> m1; //

std::vector<double> m2; // m_i are column vectors.

get_mean_vector(class1, &m1);

get_mean_vector(class2, &m2);

int N = (int)(m1.size());

TNT::Array2D<double> Sw1(N,N, 0.0);

TNT::Array2D<double> Sw2(N,N, 0.0);

get_within_class_variance(class1, m1, &Sw1);

get_within_class_variance(class2, m2, &Sw2);

TNT::Array2D<double> Sw_ = Sw1 + Sw2;

TNT::Array2D<double> I(N,N);

make_identity_matrix(N, &I);

JAMA::LU<double> Sw(Sw_);

double det = Sw.det();

if(det == 0) return false;

TNT::Array2D<double> Sw_inv = Sw.solve(I);

std::vector<double> m_diff_;

std::transform(m2.begin(), m2.end(), m1.begin(), std::back_inserter(m_diff_), std::minus<double>());

// [m_diff_] is a column vector.

TNT::Array2D<double> m_diff(N,1);

int i = 0;

for(i = 0; i<N; i++){

m_diff[i][0] = m_diff_[i];

}

TNT::Array2D<double> w = TNT::matmult(Sw_inv,m_diff);

std::vector<double> W;

for(i = 0; i<N; i++){

W.push_back(w[0][i]);

}

if(w_ != NULL){

*w_ = W;

}

get_projected_pos(class1, *w_, prj1);

get_projected_pos(class2, *w_, prj2);

return true;

}

많은 것들이 직관적이다. STL 의 몇 가지 알고리즘들을 여기서 설명하지는 않는다. 실제 사용법은 매우 간단한데, 다음과 같다.

// test.cpp

#include "fld.h"

#include <util.h>

#include <string>

#include <fstream>

#include <iostream>

bool load_data(const std::string& file_name, std::vector<std::vector<double> >* dest)

{

std::ifstream ifile(file_name.c_str());

if(ifile.is_open() == false) return false;

std::string line;

while(std::getline(ifile, line)){

if(line[0] == '#' || line[0] == '\n' || line[0] == '\r') continue;

std::vector<double> v = split<double,string2double>(line, '\t');

dest->push_back(v);

}

return true;

}

void show_vector(const std::vector<double>& v)

{

std::vector<double>::const_iterator pos = v.begin();

for(; pos != v.end(); pos++){

std::cout << *pos << std::endl;

}

int main(int argc, const char* argv[])

{

std::vector<std::vector<double> > x1;

std::vector<std::vector<double> > x2;

load_data("x1.txt", &x1);

load_data("x2.txt", &x2);

std::vector<double> w;

std::vector<double> p1;

std::vector<double> p2;

FLD(x1, x2, &p1, &p2, &w);

std::cout << w.size() << std::endl;

std::cout << w[0] << '\t' << w[1] << std::endl;

std::cout << "x1_p" << std::endl;

show_vector(p1);

std::cout << "x2_p" << std::endl;

show_vector(p2);

return 0;

}

위에서 x1.txt 과 x2.txt 는 다음과 같고, util.cpp, util.h, fld.h 는 다음과 같다.

위의 x1.txt 와 x2.txt 를 이용하여 위의 코드로 구한 [w] 와 Matlab 으로 구한 [w] 를 비교해 보면 다음과 같이 같은 것을 확인할 수 있다.

C++로 테스트 한 [w] 값

Matlab으로 테스트 한 [w] 값

source compile 방법

위의 예제는 TNT 와 JAMA 에서 제공해 주는 파일들을 이용한다. 따라서 이 파일들을 TNT 홈페이지에서 다운로드 받는다.

위에서 언급된 것 중 C++과 관련된 것은 모두 linux 에서 행해진 것이다 (CentOS). 특히 util.h 와 util.cpp 에는 이러저러한 잡다한 함수들이 구현되어 있는데, 그 중 몇몇은 linux 에만 있는 header 파일을 이용한다. 따라서 linux 에서와 Windows 에서의 compile 방법이 조금 달라질 수 있다.

1. Visual C++ 을 사용하는 경우 : util.h, util.cpp 파일에는 linux 에만 있는 header 를 include 한 부분이 있다. 이런 부분은 위의 예를 실행하는데 필요 없으므로 모두 주석처리 한다. 또한 Visual C++ 6.0 을 이용하여 compile 하면 STL 을 사용하는 부분에서 변수 타입의 이름이 너무 길어진다는 유명한 warning 이 나오므로 아예 warning level 을 0 으로 놓고 해본다.

2. gcc(g++) 을 사용하는 경우 : 간단하다. util.o 와 fld.o 를 만들어서 링크시키면 된다.

흐, 이젠 이 정도의 융통성을 발휘할 수는 있게 되었다. 뭐, 보다 직접적인 원인은 수치해석학 책이 지금 내 옆에 없기 때문이긴 하지만... 아직 책을 다 가져오지 않아서... [본문으로]

저작자표시

'컴퓨터 > 수학이랑' 카테고리의 다른 글

주성분 분석(PCA, Principal Component Analysis)의 개념 및 구현 (4)	2018.03.20
상관계수2 (0)	2011.08.30
Fisher's linear discriminant 원리 (17)	2010.11.30
통계 분석 - 평균 비교에 대하여 (0)	2010.10.21
엑셀에서 z-score로 p-value 계산하기 (1)	2010.10.12

진화하자 - 어디에도 소속되지 않기

Fisher's linear discriminant 구현

'컴퓨터 > 수학이랑' 카테고리의 다른 글

티스토리툴바

Fisher's linear discriminant 구현

'컴퓨터 > 수학이랑' 카테고리의 다른 글

관련글

티스토리툴바