2014년 8월 29일 금요일

협업 필터링(Collaborative Filtering) 예제

Matfact_newcomer
% Collaborative Filtering
% 2014. 8. 29
% refopen.blogspot.com

A = [9 2 1 8;
    8 5 6 10;
    7 9 8 9;
    2 3 4 8;
    5 9 5 7];

rank = 2;

m = mean(A, 2);
s = std(A, 1, 2);
matm = repmat(m, 1, 4);
mats = repmat(s, 1, 4);

% normalization
normA = (A - matm) ./ mats;

[u d v] = svd(normA)

vt = v';
temp = u(:, 1:rank)*d(1:rank, 1:rank)*vt(1:rank, :)

% denormalization
Ahat = matm + (temp .* mats)
u =

   -0.5699   -0.1579   -0.5892   -0.5351   -0.1297
   -0.5471   -0.3829   -0.0037    0.7322   -0.1339
    0.3887   -0.5537    0.0504   -0.1310   -0.7230
   -0.1319   -0.6232    0.4847   -0.3305    0.5001
    0.4554   -0.3654   -0.6445    0.2261    0.4388


d =

    3.0336         0         0         0
         0    2.9144         0         0
         0         0    1.5179         0
         0         0         0    0.0000
         0         0         0         0


v =

   -0.5690    0.4985   -0.4216    0.5000
    0.7368   -0.0436   -0.4530    0.5000
    0.1603    0.3410    0.7798    0.5000
   -0.3281   -0.7958    0.0948    0.5000


temp =

    0.7543   -1.2537   -0.4339    0.9334
    0.3882   -1.1743   -0.6465    1.4326
   -1.4753    0.9392   -0.3612    0.8973
   -0.6777   -0.2156   -0.6834    1.5767
   -1.3170    1.0644   -0.1417    0.3943


Ahat =

    7.6667    0.5676    3.4658    8.2999
    7.9954    4.9950    6.0085   10.0010
    7.0268    9.0288    7.9505    8.9940
    2.7065    3.7590    2.6934    7.8411
    4.3160    8.2652    6.2650    7.1538

최소자승법


최소자승법

Least square minimization

 최소자승법에 대한 글을 쓰는 것도 사실은 조심스러운 일입니다. 공학 전공자들에게 수학은 원하는 목표물을 얻기 위한 도구로 사용되지만, 수학을 전공하는 사람들에게는 존경 받는 수학자들의 이름이 거론되는 것 자체가 마음에 들지 않을 수도 있다고 생각합니다. 본 블로그에서는 공학적 응용을 위한 접근에 관점을 두고 기술합니다.

"Bendixen - Carl Friedrich Gauß, 1828" by Siegfried Detlev Bendixen - published in "Astronomische Nachrichten" 1828. Licensed under Public domain via Wikimedia Commons.
 논란의 여지가 남아 있지만, 일반적으로 최소자승법은 Carl Friedrich Gauss(이하 가우스)에 의해 1794년 발견되고, 1805년 Adrien-Marie Legendre에 의해 발간되었다고 알려져 있다. 가우스가 24살이던 1801년 소행성 ceres(2006년 국제천문연맹 회의에서 세레스의 지위가 왜행성(dwaft planet)으로 격상)의 궤도를 추정하기 위해 사용하였다고 하니 가우스의 위대함을 엿볼 수 있다.

 확률 분포에서 정규분포처럼 모수(parameter)를 추정하기 위한 방법으로 최소자승법은 폭넓은 활용을 가지고 있다. 그러면 모수를 추정하는 일은 왜 필요한가? 수 많은 데이터들을 모두 다루지 않고도 적절함 함수로 대체(fitting)할 수 있다면 그 함수의 모수를 알아내는 것만으로도 모든 데이터를 다루는 것과 별로 다르지 않은 결과를 얻을 수 있기 때문이다.

 간단한 예부터 살펴보면 다음과 같다.

 사전 정보: 데이터가 하나의 직선 위에 놓여 있을 것이다. → 직선의 방정식으로 근사화 할 수 있다.

 직선의 방정식은

 \( y = ax + b \)

로 표현할 수 있으므로, 추정해야 하는 모수는 \(a, b\), 즉 직선의 기울기와 절편이다. 우리가 잘 알고 있는 연립 방정식으로 문제를 풀기 위해서는 두 개의 미지수가 있으므로, 적어도 2개의 식이 필요하다는 것을 알 수 있다. \(m\)이 미지수의 개수, \(n\)을 주어진 식의 수라고 하면,
  1. \(m=n\)(자명하지 않은 유일해 - nontrivial solution)
  2. \(m<n\)인 경우 (부정 - overdetermined)
  3. \(m>n\)인 경우 (불능 - underdetermined)
 최소자승법을 사용하는 경우는 위의 세 가지 예시 중 두 번째 경우에 유리하다. 다시 말하면 2개의 데이터만 있어도 직선의 방정식을 구할 수는 있으나, 그보다 더 많은 수의 식이 주어지는 경우에는 어떤 것이 최적의 식인가? 라는 질문의 답이 될 수 있다. 추정한 모수로 그린 직선과 데이터와의 차이(residual)이 가장 적은 것이 답일 것이다.

 사실 해석적 방법(Analytical method) 을 사용하면 간단한 식의 경우에는 식에서 모수와의 관계를 편미분을 통해서 바로 구해낼 수 있다. 아래 그림을 보면 residual에 대한 각 모수에 대한 편미분이 최대 또는 최소가 되는 점, 미분치가 0인 점을 찾아내면 된다는 것을 알 수 있다.



 요즘 우리의 컴퓨터는 너무나 성능이 좋고, 빅 데이터의 시대에 살아가고 있기 때문에 수치적 방법(Numerial method)를 사용하는 것이 보편화 되었다.

 수치적 방법으로 최소자승법을 푸는 방법은 다음 문제와 관련이 있다.

  • 선형 동차식(homogeneous equation)인 경우
  • 선형 비동차식(inhomogeneous equation)인 경우 
 여기서는 \(\mathbf{Ax=b}\)의 형태로 정리할 수 있는 경우 이므로 비동차식 해법을 적용한다.

 \( x_1, y_1, x_2, y_2, x_3, y_3, ... , x_n, y_n \)의 데이터가 주어진 경우

\( ax_1 + b = y_1 \)
\( ax_2 + b = y_2 \)
\( ax_3 + b = y_3 \)
\( \vdots \)
\( ax_n + b = y_n \)

식을

\(\left[\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ x_3 & 1 \\ \vdots \\ x_n & 1\end{array} \right] \left[\begin{array}{c} a \\ b \end{array} \right] = \left[\begin{array}{c} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n \end{array} \right] \)
 \( \mathbf{~~~~~A~~~~~~~~~~x~~~ = ~~~ b} \)
 \( \mathbf{~~(n\times2)~(2\times1) =  (n\times1)} \)

와 같이 정리할 수 있다.

 \(\mathbf{Ax=b}\)의 식에서 행렬 \(\mathbf{A}\)의 의사역행렬(left pseudo inverse) \(\mathbf{A^+=(A ^{T} A) ^{-1} A ^{T}}\)를 양변에 곱함으로써 \(\mathbf{x}\)를 구할 수 있다.

 여기서 SVD(Singular Value Decomposition)를 이용할 수 있는데, 행렬 \(\mathbf{A}\)의 역행렬은 \(\mathbf{A ^{-1} =V ^{T} D ^{*} U ^{T}}\)를 이므로 SVD를 통해 바로 구할 수 있기 때문이다.

 만약 MATLAB을 이용한다면 강력한 solver가 제공된다. mldivide 함수 또는 '\' 연산자를 사용하면 행렬의 크기에 따라 적절한 방법의 최소자승법에 의한 해를 구해 준다.

 MATLAB으로 구현한 최소자승법의 예제를 제시합니다.

2014년 8월 28일 목요일

동차 선형 방정식을 SVD로 풀 수 있는 이유


동차 선형 방정식을 SVD로 풀 수 있는 이유

Solving homogeneous linear equation using SVD

 최소 자승법을 이야기 하면서 의사 역행렬을 이용해 over-constrained eq.을 푸는 방법을 소개 했는데 정작 왜 그것이 정답이 되는지는 얼버무리고 넘어간 것 같아서 이 글을 게시합니다.

 SVD(Singular Value Decomposition)은 선형 대수학에서 매우 중요한 의미를 가집니다. eigenvalue, eigenvector가 중요한 것처럼 singular value들은 eigenvalue와 관계가 있기 때문입니다. 컴퓨터 비전이나 기계 학습에서 선형 대수학이 빠질 수 없는 이유는 바로 선형 변환에 불변하는 특성을 가지는 eigenspace 해석을 통해서 중요한 성질을 구할 수 있기 때문입니다. PCA(Principal Component Analysis)도 결국 SVD를 이용해서 주요 성분을 얻어내는 것이고, 기계 학습에서 sparse coding으로 중요한 정보 만을 이용해서 영상을 분류하는 방법에도 적용될 수 있습니다.

 행렬 \(\mathbf{A}\)의 singular value를 \(\sigma\)라고 하면, \(\mathbf{A^*A}\)와 \(\mathbf{AA^*}\)의 0이 아닌 eigenvalue, \(\lambda\)와 \(\sigma\) = \(\sqrt{\lambda}\)의 관계를 가집니다. 여기서 \(\mathbf{A^*}\)는 \(\mathbf{A}\)의 켤레전치행렬(conjugate transpose) 입니다. 간단하게 상수라고 생각하면 제곱 한 것을 분해한 것과 그냥 분해한 것의 차이 정도로 대략 감을 잡으시면 될 것 같습니다.

 그럼 SVD로 어떻게 동차 방정식의 해를 구할 수 있는지 보겠습니다.

 \(\mathbf{Ax} = 0\)

와 같은 동차 방정식이 주어지는 경우, 자명하지 않은 해(non trivial solution)을 구한는 것이 목적입니다. (\(\mathbf{x}=0\)이면 항상 참이 되니 자명한 해이므로 관심이 없습니다.)

 그래서 자명한 해가 아닌 해를 구하기 위해서 제약 조건을 줍니다. 보통은 \(||\mathbf{x}||=1\) 과 같이 벡터의 크기를 정해줍니다. 우변이 0이기 때문에 벡터의 크기는 아무래도 관계 없기 때문입니다.

 이제 행렬 \(\mathbf{A=UDV^T}\)로 분해하면,

\(\mathbf{||UDV^Tx||=||DV^Tx||}\) - ①
\(\mathbf{||x||=||V^Tx||}\) - ②

로 둘 수 있습니다. 왜냐하면, SVD로 구한 \(\mathbf{U, V^T}\)행렬은 각 열이 모두 orthogonal 하기 때문입니다. 이 특성은 회전 행렬과 같습니다. 벡터 공간에서 임의의 벡터를 아무리 회전 시켜도 그 크기가 변하지 않는다는 것을 생각해보면 식 ①, ②가 성립하는 것을 알 수 있습니다.

\(\mathbf{y=V^Tx}\) - ③

을 도입해 보면, 최초의 문제는 결국 '\(\mathbf{||Dy||}\)를 최소화하되 \(\mathbf{||y||}=1\)을 만족 시켜라.'는 문제로 바뀌게 됩니다.

 여기서 SVD의 특성을 다시 한 번 생각해보면, 분해된 행렬 \(\mathbf{UDV^T}\)에서 \(\mathbf{D}\)는 대각 행렬이고 원소들은 singular value가 큰 순서대로 정렬되기 때문에 \(\mathbf{||Dy||}\)를 가장 작게 만들면서 \(\mathbf{||y||}=1\)을 만족 시킬 수 있는 방법은 \(\mathbf{y}=(0, 0, ..., 0, 1)^T\)로 두는 것입니다.

 그러므로, \(\mathbf{x=V^Ty}\)는 결국 \(\mathbf{V^T}\)의 마지막 열이 됩니다.

2014년 8월 27일 수요일

기계 학습 - 서론

기계 학습 - 서론

Machine learning - preliminary


 기계 학습을 한 마디로 정의하기는 어려울지도 모르겠습니다. 그러나 학습, 추론, 결정을 요구하는 시스템을 수학적 기반으로 연구하는 분야라는 데에는 이견이 없을 듯 합니다. 최근에 하드웨어 성능의 향상(특히 GPGPU-General-purpose computing on graphics processing units), 빅 데이터의 경향에 따라 산업계, 학계의 지대한 관심을 받고 있는 분야입니다. 그에 따라 연구 그룹도 우후죽순 생겨나고 있습니다. 본 블로그에서는 이 서론을 시작으로 몇 가지 활용에 대해서 기술하고자 합니다.

 기계 학습을 연구하는 사람들의 공통적 철학은 '오컴의 면도날(Occam’s Razor)'이라는 철학자의 말을 빌어보는 것이 좋을 것 같습니다. 사실 오컴은 사람 이름이 아닌 지명입니다. 레오나르도 다 빈치 라는 이름이 빈치에서 태어난 레오나르도라는 의미인 것과 마찬가지로 '윌리엄'은 영국 오컴 지방의 프란체스코회 수사였습니다.

 오컴의 면도날의 의미는 문제를 풀 수 있는 방법이 여러 가지가 있을 경우에 가장 간단한 방법이 최선의 답이라는 철학적 원리를 바탕으로 합니다. 경제성의 원리라고도 불립니다. 기계 학습을 전공하는 사람들 중에서 최근 각광 받는 Deep learning 분야를 연구 하는 사람들은 아마도 복잡한 수식이나 모델링을 통한 해법보다는 자극(Stimulation)에 의해서 활성화(Activation)되는 인간의 뉴런과 시냅스를 모사하고 학습하는 것이 오컴의 면도날을 충족하는 가장 바람직한 정답이라고 생각했을 수 있습니다.

원본 이미지: http://vashonsd.org/teacherweb/floyd/index.php/floydtemp/pages/C1191/P90
뉴런과 시냅스를 모사한 Nueral Network
 신경망 회로(Nueral network)의 개념이 최초로 도입되던 시절에는 지금에 비해서 컴퓨터의 성능이 보잘 것 없어서 실용화 하기에는 불가능할 것처럼 보였습니다. 그러나 최근에 제시되는 결과는 간단한 분류에서 문자 인식, 음성 인식, 얼굴 인식처럼 점점 고도화 되고 상용화 가능성을 보이고 있습니다.

  신경망 회로(Neural network)에 관해 생소하다면 역전파(Back-propagation)을 통해 작은 수의 레이어가 있는 경우에 한해 접근해 보는 것이 좋을 것 같습니다.

신경망 회로에서 역전파를 통한 학습


신경망 회로에서 역전파를 통한 학습

Learning through back-propagation in 3-layers neural network

 layer가 3개인 신경망 회로(이하 NN)는 input, hidden, output layer로 구성됩니다. 학습 단계는 output layer에서 원하는 출력을 얻기 위해 hidden layer가 활성화 되도록 학습하는 것이 목적입니다.

 예를 들면, 어떤 물체가 연필인지 아닌지를 구분하는 NN을 만들고자 할 경우에 loss function은

\(loss = 연필 - output\)

이 되어야 합니다.

 convex optimization 문제로 만들기 위해 \(loss = (연필 - output)^2\)으로 놓고 미분하면 \(\frac{\delta loss}{\delta output} = 2(연필 - output)\)이 됩니다.

 그러면 NN이 원하는 출력 값과 얼마나 다른 값을 출력 하는지를 알 수 있습니다. 출력값의 차이가 크다면 input->hidden->output layer로 이르는 뉴런들이 올바르게 활성화될 수 있도록, 오차만큼 역전파(Back-propagation)을 통해서 weight값을 고쳐줍니다. 이 과정의 반복을 학습이라 부릅니다. 그림을 보면 이해가 쉽습니다.

  • Feed-forward 과정입니다. 입력 뉴런과 weight의 값들의 합이 다음 layer의 입력으로 들어가고 activation function의 임계치에 따라서 활성화되는지 그렇지 않은지가 결정됩니다. activation function은 과거에는 모두 sigmoid function을 사용하였으나 최근에는 경향이 좀 바뀌었습니다. 추후에 언급할 기회가 있을 거라고 생각합니다.
  • Backward propagation의 과정입니다. loss function의 값은 곧바로 weight를 학습하는데 사용할 수가 없기 때문에, 다시 말하면 출력값이 원하는 값과 다르게 나오는데 영향을 미친 뉴런이 어떤 것이고 얼마 만큼의 가중치를 가졌는지 알 수 없기 때문에 미분치를 역으로 전파시킵니다. 그림에서 \(\delta\)는 이 미분치를 의미하고 chain rule로 입력 layer까지 전파됩니다. 이제 weight를 갱신하고 다시 이 과정을 반복합니다. 이때, \(\eta\)는 learning rate으로 부르는데 학습의 시간을 결정합니다. 크게 설정하면 빠르게 학습되지만 over-shoot이 생겨서 오차 범위에 정착하기가 힘들기 때문에 수렴하는 속도를 보고 변경하는 것이 일반적입니다.











원본 이미지: http://galaxy.agh.edu.pl/~vlsi/AI/backp_t_en/backprop.html

 신경망의 기본 개념은 인간의 뉴런과 시냅스를 모사한다는 거창한 말에 비해서 이론적으로 어렵거나 구현이 까다롭지 않습니다.(물론 layer의 수가 적을 때의 이야기 입니다. 영상에 적용하는 convolutional neural network은 이보다 엄청나게 많은 수의 뉴런과 weight가 존재하기 때문에 구현 자체가 어려운 경우가 많습니다.)

 학습은 방법보다는 어떤 응응 분야에 적용할 것 인지를 염두에 두어야 합니다. 왜냐하면 학습 데이터에 너무나 잘 들어맞는(over-fitting) NN을 가지고 일반적인 데이터에 적용하면 분류 혹은 검출률이 낮을 수 있기 때문입니다. 반대로 적당히 학습되어있는 경우(under-fitting)엔 반대의 경향이 나타날 수 있습니다. 그래서 학습할 목적에 따라 오검출율(FAR-False Alarm Rate)이 얼만큼 허용 되는지 판단해야 할 수 있습니다.

 C# 예제를 게시한 MSDN Magazine이 있어서 이 곳을 소개합니다. 참고하시면 좋을 것 같습니다.

 http://msdn.microsoft.com/en-us/magazine/jj658979.aspx

구글 캠퍼스 서울 세계에서 3번째로 문 연다.

 새 정부 들어서 창조 경제와 관련된 이슈가 많은데, 정치적인 목적이 아닌 실제로 창업자들을 위한 프로그램이 지원되었으면 좋겠습니다.

 http://www.zdnet.co.kr/news/news_view.asp?artice_id=20140827104637

2014년 8월 26일 화요일

MATLAB으로 구현한 최소자승법 예제

lsq1

MATLAB으로 구현한 최소자승법 예제

% Least square minization demo
% 2014. 8. 26
% refopen.blogspot.com

추정하고자 하는 모수 참 값

a = 3.5;
b = 105.3;

모수의 참값으로 그려본 직선

x = 0:0.1:100;
y = a*x + b;

plot(x, y,'LineWidth', 3);
hold on;

측정되는 가상 데이터 생성

x_perturbed = x + randn(1, 1001)*10;
y_perturbed = y + randn(1, 1001)*10;

plot(x_perturbed, y_perturbed, 'MarkerSize', 2, 'Marker', '*', 'Color', 'r', 'LineStyle', 'none');
hold on;

최소자승법으로 모수를 추정

[param_est] = [x_perturbed' ones(1001, 1)] \ y_perturbed';

y_est = param_est(1)*x + param_est(2);
plot(x, y_est, 'Color','g' , 'LineWidth', 3);
legend('true', 'measured data', 'Least Square Estimation');
param_est =

    3.1268
  123.7713