협업 필터링과 추천 시스템으로 내일 뭐 먹을지 결정해보기

협업 필터링과 추천 시스템으로 내일 뭐 먹을지 결정해보기

Collaborative Filtering and Recommendation System

 한글로 적어 놓고 나니 왠지 어색합니다. 협업 필터링(Collaborative Filtering)과 추천 시스템(Recommendation System)은 밀접한 관련이 있고, 우리가 모르는 사이에 매일 이것을 이용하고 있습니다. 실제로 Amazon에서 구매하는 물건의 약 35%가 추천으로부터 발생하고, Netflix의 영화 대여도 2/3 정도가 추천에 의해서 발생한다고 알려져 있습니다. 한국에서도 ‘와챠’라는 모바일 앱이 개인화된 추천 서비스를 제공함으로써 기존의 영화 리뷰 시장을 뒤흔들었습니다. 이 글에서는 추천 시스템 중에서 협업 필터링의 개념을 이용해서 ‘내일 뭐 먹을까?’를 결정해보는 간단한 예시를 통해 접근하고자 합니다.

 여기에서도 SVD를 적용합니다. SVD는 singular value가 큰 순서로 정렬되는 형태이므로, 중요한 정보를 압축할 수 있는 특성을 가진다는 것을 알고 있다면, Low-rank approximation으로 분류를 하는데 활용할 수 있다는 것도 생각해 볼 수 있습니다. 다시 말하면 eigenvalue가 큰 것들은 어떤 자료에서 중요한 의미를 지니는데, 그에 해당되는 eigenvector들은 모두 orthogonal합니다. orthogonal하다는 말은 직교한다는 말하고 비슷한데, 그릴 수 없는 3차원을 넘어가는 다차원에 대해서도 일반화된 개념입니다. 내적이 ‘0’이므로 서로 겹치는 부분이 없습니다. 어떤 자료를 SVD한다는 말은 결국 중요한 의미를 지니는 eigenvector공간으로 재 투영한 값들로 재구성한다는 의미로 해석하면 될 것 같습니다.

 예를 들자면, 어느 학교의 학생에 대한 정보가 몸무게, 손가락 길이, 키로 각각 주어진다고 가정해보겠습니다. 그러면 3차원 공간에 데이터를 그려볼 수 있습니다. (1,0,0)(0,1,0)(0,0,1)-x, y, z-몸무게, 손가락 길이, 키를 기저 벡터로 하는 공간에 그리는 겁니다. 그런데, 만약에 키가 모두 비슷한 학생들만 다니는 학교라고 하면 z축은 별로 의미가 없어 집니다. x-y평면에 투영해서 봐도 3차원 공간의 데이터와 별로 손실이 없을 수 있다는 의미입니다.

 추천 시스템에 활용하는 간단한 예시를 들어보겠습니다.

 어떤 모임에 4명의 사람이 있었는데, 새로운 사람이 들어온 경우에 그 사람이 어떤 성향을 가지는 사람인지 판단하고 메뉴를 추천해주는 것을 생각해 볼 수 있습니다. 아래 표와 같은 정보가 주어진다고 가정해보겠습니다.

	Korean Cuisine	Pizza	Hamburger	Chinese Cuisine
A	9	2	1	8
B	8	5	6	10
C	7	9	8	9
D	2	3	4	8

 데이터를 직관적으로 관찰해보면 중국 요리는 대부분의 사람이 좋아하고 피자를 좋아하는 사람은 대체로 햄버거도 좋아한다는 특성을 관찰할 수 있습니다.

 새로운 멤버 E와 식사를 하러 가야 하는데, 그날 마침 한국음식과 햄버거 중에 결정해야 하는 상황이라면 E는 어떤 음식을 더 좋아할까요?

	Korean Cuisine	Pizza	Hamburger	Chinese Cuisine
A	9	2	1	8
B	8	5	6	10
C	7	9	8	9
D	2	3	4	8
E	?	9	?	7

 SVD를 이용한 Low-rank approximation으로 추정해볼 수 있습니다. Rank 3으로 줄이면 아마도 한국음식, 피자, 햄버거와 비슷한 방향을 가지는 기저 벡터 공간에 재 투영될 것을 예측할 수 있습니다. 다음 예시를 참고하여 주시기 바랍니다. 

 협업 필터링 예제

E가 평가를 하지 않은 빈 곳에 대해서 평균으로 채우고 Low-rank approximation을 수행해보면 이 사람은 아마도 햄버거를 더 좋아하는 사람인 것 같습니다.

 실제 시스템에서는 고려해야 할 사항이 훨씬 더 많을 것입니다. 왜냐하면 SVD의 연산 복잡도는 선형으로 증가하지 않기 때문에 100만명의 회원이 있을 경우에 간단한 방법으로 해결할 수 없을 것입니다. 의도적으로 데이터에 손상을 가하려는 사람들(평점을 일부로 높게 주거나 낮게 주는 사람들)이 있을 수도 있고 개인적인 성향에 따라 평점이 후한 사람도 있고 그렇지 않은 사람도 있을 것입니다.

협업 필터링(Collaborative Filtering) 예제

% 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

Published with MATLAB® R2013b

최소자승법

최소자승법

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으로 구현한 최소자승법의 예제를 제시합니다.

동차 선형 방정식을 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}\)의 마지막 열이 됩니다.

기계 학습 - 서론

기계 학습 - 서론

Machine learning - preliminary

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

원본 이미지: http://www.rugusavay.com/william-of-ockham-quotes/

 기계 학습을 연구하는 사람들의 공통적 철학은 '오컴의 면도날(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을 사용하였으나 최근에는 경향이 좀 바뀌었습니다. 추후에 언급할 기회가 있을 거라고 생각합니다.

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

• 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

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

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

Published with MATLAB® R2014a

최적화 문제(Optimization problem)

최적화 문제

Optimization problem

 '최적화' 라는 말이 갖는 의미가 참 모호한 것 같다. 성능을 최적화 한다는 의미로 받아 들이면 탑재된 플랫폼에 맞도록 속도나 메모리 사용량을 줄인다는 의미로 생각할 수도 있다.

 그러나 수학적 의미의 최적화 문제라 하면 주어진 제한 조건을 만족하는 최대, 최소 해를 찾는 문제로 정의된다. 제한 조건은 때로는 출력값의 범위를 제한할 수도 있고, 입력값의 범위를 제한할 수도 있으므로 입력과 출력 집합에 대한 최적의 연결(함수)을 찾는 문제라고 생각할 수도 있다.

 로봇, 컴퓨터 비전, 기계 학습과 밀접한 관련을 가지고 있으므로 자세히 알고 있으면 매우 편리하다. 아래에 든 예시 말고도, 생각할 수 있는 거의 모든 문제는 목적 함수와 제한 조건을 설정할 수 있다면, 최적화 문제로 접근이 가능하다.

• Laser Range sensor 기반의 지도 작성 및 위치 인식에 사용되는 ICP(Iterative Closest Point)
• 다수의 카메라의 특징점의 위치를 알고 있는 경우, Bundle Adjustment를 통해 최초의 값으로부터 최적의 카메라 위치, 특징점의 3차원 위치를 구하기 위해 LM(Levenberg-Marquartd) 알고리즘을 이용
• Neural Network의 Weight를 학습하기 위해 Loss function의 값을 역 전파(Back-propagation)하여 최적의 Weight값을 찾는데 이용

 Least square minimization(최소자승법) 문제를 통해서 접근하는 것이 최적화 문제에 입문하는데 효과적인 방법이 될 것이다.

SLAM(Simultaneous localization and mapping) - Simulation

SLAM 예제

 애초에는 SLAM 예제를 만들어 볼 생각이었으나, 이미 너무 잘 만들어진 시뮬레이터가 있어서 이것을 소개하는 것으로 갈음하고자 합니다.

 Austraila 의 연구자 Tim Bailey의 SLAM simulations software를 소개합니다.

 구성된 내용은

• EKF-SLAM version 1, 2
• FASTSLAM 1.0, FASTSLAM 2.0
• UKF-SLAM

 입니다.

Tim bailey의 SLAM simulations EKF-SLAM

 SLAM을 연구하고자 입문하시는 분들에게 큰 도움이 되실 거라고 생각합니다. 아래 링크를 참조하세요.

http://www-personal.acfr.usyd.edu.au/tbailey/

SLAM(Simultaneous localization and mapping)과 칼만 필터 두 번째

SLAM(Simultaneous localization and mapping) - Kalman Filter second

동시적 위치추정 및 지도작성과 칼만 필터 두 번째

 이제 우리는 SLAM(Simultaneous localization and mapping)과 칼만 필터를 통해서 2차원 에서 운동하는 등속 모델의 물체에 대한 추적을 할 수 있게 되었다. 그런데 궁금한 것이 있다. 관측 행렬 \(\mathbf{H}=\left[\begin{array}{cccc}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\end{array}\right]\)에서 알 수 있는 것처럼 측정할 수 있는 값은 위치 뿐이고, 속도에 대한 정보를 입력해 준 적은 없는데, 상태 벡터를 관찰해보면 속도에 대한 값이 생성되고 있다. 왜 그럴까? 공분산 행렬 \(\mathbf{P}\)을 살펴보면 예측 과정에서는 행렬의 대각 성분, 다시 말하면 각 위치와 각 속도 스스로의 항에 더해지지만, 갱신하는 과정에서 대각 성분이 아닌  곳에 연관성(correlation)이 발생하기 때문이다. 그러므로 칼만 필터는 위치에 대한 정보만 입력 받아도 상태 천이 행렬로부터 적절한 속도를 갱신하도록 한다고 생각해 볼 수 있다.

 그러면 마찬가지로 입력 받지 않은 다른 측정 값에 대해서도 갱신하는 것이 가능하지 않을까? 그렇다. SLAM의 기본 개념은 측정된 값을 이용해서 측정 되지 않은 다른 값들을 갱신하는 것이다. 현재 위치에서 측정할 수 있는 랜드 마크는 센서의 시야각과 거리의 한계로 제한이 있을 수 밖에 없다. 한정된 측정값을 이용해서 다른 상태 벡터의 값들을 좀 더 신뢰할 수 있는 값으로 갱신하는 것이다.

 아래 그림을 보면, 움직이는 로봇이 측정할 수 있는 랜드 마크의 숫자가 한정적인 경우에도 하나의 상태 벡터를 포함하고 있는 상태에서는 모든 랜드 마크와 로봇의 위치에 대한 신뢰도가 향상되는 것을 볼 수 있다. 로봇과 랜드 마크의 타원은 공분산 행렬(Covariance matrix)에 비례하여 그려진 것이므로 해당 항목의 불확실성(Uncertainty)을 의미하는 것으로 생각할 수 있다.

이미지 원본: Andrew Davison의 박사 학위 논문

 이렇게 로봇의 상태와 랜드 마크의 좌표가 하나로 합쳐진(Augmented) 상태 벡터를 사용하는 방법은 관측 범위의 제약으로 한정된 관측이 수행되는 경우에 유용하지만, 태생적으로 차원의 저주(Curse of dimensionality) 문제를 갖고 있다.

 추정하고자 하는 로봇 좌표계의 차원과 특징점의 차원이 모두 하나의 상태 벡터에 포함되므로 특징점의 개수가 증가함에 따라서, 공분산 행렬 \( \mathbf{P} \)의 크기가 \( 2^d \)에 비례하여 증가하기 때문이다.

SLAM(Simultaneous localization and mapping)과 칼만 필터

SLAM(Simultaneous localization and mapping) - Kalman Filter

동시적 위치추정 및 지도작성과 칼만 필터

 국문으로 번역된 이름이 마음에 들지 않지만, 마땅히 다르게 번역할 방법이 없어서 기존에 사용되던 것들 중에서 차용하였습니다.

 제목에서 느낄 수 있는 것처럼 이동 로봇이 미지의 세계를 방문할 때 자신의 위치추정과 지도작성을 동시에 수행하는 것을 말한다. 동시적이라고 하면 시간의 흐름상 완전한 동시성을 의미하는 것처럼 느껴지기 때문에 사실은 일관된, 연관된 정도로 번역하는 것이 자연스럽다고 생각된다.

 칼만 필터를 이용한 이동 로봇의 SLAM에 대해서 먼저 기술하고 다음으로 다른 종류의 필터(EKF, UKF, Particle Filter)에 대해 적어볼 예정입니다.

 위치 추정과 지도 작성은 동시에 진행하는 것이 타당하다. 위치를 인식하기 위해서는 지도가 정확해야 하는데, 지도의 정확성은 위치의 정확성에 의존하기 때문이다. SLAM이 어려운 이유는 그림에서 보는 것처럼 예측할 수 없는 요소들이 많기 때문이다.

왜 SLAM은 어려운 문제인가?

 세상엔 정확한 센서는 존재하지 않기 때문에 항상 노이즈를 포함하고 있다. 정확한 센서 하나만 있었더라도 이렇게 힘들게 고민할 이유가 없었을 것이다.

 위의 그림을 조건부 확률식으로 적어보면 아래와 같다.

조건부 확률식으로 전개한 SLAM 문제

 상태 벡터 \( \mathbf{x}_{k-1} \)에서 \( \mathbf{x}_{k} \)로 이동하는 로봇이 \( \mathbf{m}_{i}, \mathbf{m}_{j}\) 랜드 마크를 관측한 값이 측정치 \( \mathbf{z}_{k,j}\) 로 입력된다. \( \mathbf{u}_k \)는 현재 로봇의 조종 명령이다.

 조건부 확률 식을 말로 풀이하면 " \( \mathbf{z}_{1:t}, \mathbf{u}_{1:t} \)가 만족 되는 경우에 즉, 최초 시점 1에서 현재 시점 \( t \)까지의 관측 값과 조종 명령이 주어질 때, 현재 위치를 의미하는 상태 벡터 \( \mathbf{x}_t \)와 지도를 의미하는 \( \mathbf{m} \)이 어떤 확률 분포를 가지는가?" 로 서술할 수 있다.

 확률은 그 합이 '1'이 되어야 하기 때문에 조건이 주어지지 않은 경우의 현재 상태와 지도에 대한 확률과 합하면 조건부 확률 값은 '1'이 된다. 위와 같은 조건부 확률 식으로 전개해 두고 나면 앞서 기술한 베이지안 추정 방법이나 마르코프 위치 인식 방법을 사용할 수 있게 되므로 매우 유리한 점이 생긴다.

 상태 천이 행렬이 선형(Linearity: Homogeneity와 superposition을 만족)인 경우에 적용할 수 있다. 선형 상태 천이를 한다는 것은 상태 천이 함수가 선형이라는 말이며 등속, 등가속 운동처럼 이후의 운동의 예측할 수 있는 경우를 말한다.

 그림에서 상태 천이 함수는 \( A \) 행렬이다. 등속 모델인 경우엔 \( \mathbf{A}=\left[\begin{array}{cccc}1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\end{array}\right] \), 상태 벡터는 \(\mathbf{x}=\left[\begin{array}{c}x \\ y \\ \dot{x} \\ \dot{y} \end{array}\right] \) 로 설정하면 2차원 좌표에서 등속으로 운동하는 물체에 대한 칼만 필터의 추정식이 완성된다. 실제 시스템 모델이 맞도록 공정 잡음과 측정 잡음을 설정하면 된다. 

 간단하게 몇 줄의 코드 만으로도 칼만 필터는 훌륭하게 동작합니다. 아래에 2차원에서 움직이는 물체에 대한 칼만 필터 예제를 게시합니다. 

칼만 필터 예제

% 2D Kalman filter example
% 2014. 8. 21
% refopen.blogspot.kr

2차원에서 움직이는 참 값 생성

true = [0:0.03:pi/2; sin(0:0.03:pi/2)];
close all;
plot(true(1, :), true(2, :), 'b*-');

파라미터 초기화

x = [0; 0; 0; 0];
A = [1 0 1 0;
    0 1 0 1;
    0 0 1 0;
    0 0 0 1];
sigmax = 0.01;
sigmay = 0.01;
sigmaxdot = 0.01;
sigmaydot = 0.01;
Q = [sigmax.^2 0 0 0;
    0 sigmay.^2 0 0;
    0 0 sigmaxdot.^2 0;
    0 0 0 sigmaydot.^2];
H = [1 0 0 0; 0 1 0 0];
R = [0.05 0;
    0 0.05];

xhatk = x;
Phatk = Q;

kfest = zeros(size(true));
measure = zeros(size(true));

예측과 갱신 반복

for n=1:size(true, 2)   %% 참 값으로 만든 횟수 만큼

    % 측정값 생성
    % 표준 편차가 0.03인 가상의 노이즈를 더한 가상의 측정값
    zk = true(:, n) + randn(2, 1)*0.03;
    measure(:, n) = zk;

    % 예측 과정
    xbark = A*xhatk;
    Pbark = A*Phatk*A' + Q;

    % 갱신 과정
    K = (Pbark*H'*(H*Pbark*H'+R)^-1);
    xhatk = xbark + K*(zk - H*xbark);
    Phatk = (eye(4) - K*H)*Pbark;

    kfest(:, n) = xhatk(1:2);

end

hold on;
plot(measure(1, :), measure(2, :), 'k*-');
plot(kfest(1, :), kfest(2, :), 'r*-');

legend('true', 'measure', 'kf est.');

Published with MATLAB® R2014a

이동 로봇을 위한 칼만 필터

이동 로봇을 위한 칼만 필터

Kalman Filter for a mobile robot

 EBS 다큐멘터리 '자본주의란 무엇인가?'에서 세계적인 석학들에게 '자본주의란?' 명제에 대해 답변해주기를 요청하는 장면에서 모두 머리를 절래절래 흔드는 장면이 나온다.

 마찬가지의 관점에서 공학자로서 칼만 필터에 대해 기술한다는 것은 굉장히 괴로운 일이다. 깊은 통찰이 담겨있는 근본적인 원리이기 때문이다. 이미 인터넷을 검색하면 무수히 많은 튜토리얼과 활용예가 있음에도 불구하고 굳이 다시 한 번 거론하는 이유는 본 블로그에 게시된 글들과의 연관성을 위해서 SLAM과 관련된 내용을 게시하기 이전에 필요하기 때문이다.

 베이지안 추정 필터와 관련된 글을 먼저 확인 바랍니다.

 칼만 필터를 이동 로봇에 적용하는 것이 적합하다고 생각되는 이유는 물리적인 특성을 유지하는 물체이기 때문이다. 다시 말하면 뉴턴의 운동 법칙을 위배하는 물체는 이 세상에 없기 때문에 (빛보다 충분히 느린 물체의 경우에) 이동 로봇은 물론이고, 잠수정, 비행체를 포함한 모든 물체는 일정한 패턴을 가지고 움직인다. 갑자기 서울에서 있던 물체가 부산에 나타날 수는 없다는 것이다. 혹자는 금융에 칼만 필터를 적용하려는 시도를 하는 경우가 있는데, 본인이 생각하기엔 불가능하다. 왜냐하면 금융은 모델 자체가 갖는 특성이 물리적인 의미가 없는 경우가 많기 때문이다. (수 많은 파생상품은 투자 대상 조차도 형태가 없는 경우가 있다.)

 아래의 글은 칼만 필터의 예측과 갱신과정에 대해서 이해하고 있는 독자를 대상으로 합니다. 칼만 필터의 기본적인 개념은 위키 백과의 글을 참고하시기 바랍니다.

 칼만 필터가 실제로 좋은 방향으로 상태를 예측할 수 있는지는 다음 칼만 게인이 의미하는 바를 해석함으로써 알 수 있다.

Kalman gain

 그림에서 보인 수식은 칼만 게인이 갱신 과정에서 어떤 역할을 수행하는지 해석한 것이다. 칼만 게인은 공정 잡음으로부터 발생된 오차 공분산 행렬 값 \(P_k^-\)과 측정 잡음 행렬 \(R_k\)와 관련이 있다. 측정 잡음이 '0'으로 수렴하면 칼만 게인은 \(H^{-1}\)가 되므로 상태 벡터의 갱신 값은 측정치를 따르게 된다. 반대로 오차 공분산 행렬 \(P_k^-\) 값이 '0'으로 수렴하게 되면 칼만 게인은 '0'이 되어 상태 벡터의 예측치를 따르게 된다.

 이처럼 칼만 게인은 알고 있는(또는 알고 있다고 믿는) 공정 잡음과 측정 잡음 예측치와 측정치 중에서 현재 상태에 가장 최적의 값을 추론한다.

 이를 이용한 위치 추정과 지도 작성에 관한 문제는 SLAM이라는 제목으로 다시 게재할 예정이므로 칼만 필터에 관해서는 여기에서 줄인다.

우분투 14.04 LTS에서 Nvidia cuda 6.0 설치

우분투 14.04 LTS에서 Nvidia CUDA 6.0 설치

 우분투 14.04에서 CUDA 6.0 설치를 위해서는 몇 가지 조작이 필요하다.

• 새로 설치하는 경우라면 다음 명령을 확인해본다.

1. apt-get install build-essential

• cuda 6.0 run 패키지를 다운로드 한다. 현재는 다음 버전이 최신 버전으로 확인됨.http://developer.download.nvidia.com/compute/cuda/6_0/rel/installers/cuda_6.0.37_linux_64.run
• 터미널을 열고 run 패키지의 압축을 해제한다. 압축을 해제하지 않고 run 패키지를 바로 실행 시켜도 설치가 되어야 하지만 아직은 호환성에 문제가 있는 것으로 보임.

1. mkdir ~/Downloads/nvidia_installers;
   cd ~/Downloads
   ./cuda_6.0.37_linux_64.run -extract=~/Downloads/nvidia_installers;

• 기존에 설치된 nvidia 드라이버를 삭제한다.

• sudo apt-get --purge remove nvidia-*

• 이제 X window를 ctrl+alt+F1으로 탈출하고 터미널로 로그인한 뒤 다음을 통해 설치를 진행한다.

• cd ~/Downloads/nvidia_installers;
  sudo service lightdm stop
  sudo killall Xorg
  sudo ./NVIDIA-Linux-x86_64-331.62.run 

• nvidia 드라이버가 설치되면 이제 X window를 ctrl+alt+F1으로 탈출하고 터미널로 로그인한 뒤 다음을 통해 설치를 진행한다.

• sudo modprobe nvidia
  sudo ./cuda-linux64-rel-6.0.37-18176142.run
  sudo ./cuda-samples-linux-6.0.37-18176142.run

• 설치가 제대로 되었는지 확인

• cd /usr/local/cuda/samples
  sudo chown -R <username>:<usergroup> .
  cd 1_Utilities/deviceQuery
  make .
  ./deviceQuery    

• X window로 복귀

• sudo service lightdm start

• 설치가 제대로 되었는지 확인할 수 있다.

• lsmod | grep nv

Dijkstra, Best First Search, A* search

Dijkstra, Best First Search, A* search

전역 경로 계획

 거창하게 전역 경로 계획(Global path planning)이라고 이름을 붙이기는 했지만, 지역 경로 계획에 대비되는 용어일 뿐이다.

 지역 경로 계획은 가까운 곳의 장애물을 회피하기 위한 방법이라면 전역 경로 계획은 출발점에서 도착점까지 가장 가까운 거리를 찾는 것을 목적으로 한다. 그러므로 전역 경로 계획으로 생성된 경로에서 지역적으로 발생하는 장애물을 회피하는 방법을 복합적으로 사용하는 것이 바람직하다. 전역 경로 계획은 컴퓨터 공학 그래프 이론에서 출발한 것들이기 때문에 지역 경로 계획과는 근본적으로 조금 다른 성향을 가진다.
 실상이 없는 노드와 에지의 연결체에서 최단 경로를 찾는 것은 로봇의 모션 플래닝의 관점과는 조금 맞지 않을 수도 있다.

 이외에도 다른 방법들이 많이 있지만 다음 세가지에 국한해서 예시를 든다.

• Dijkstra
• Best First Search
• A* Search

 워낙 오래되고 많이 사용되는 방법들이기 때문에 특별히 많은 설명이 필요하지는 않을 것 같다. 다음그림으로 세 가지 방법의 차이점을 명확히 이해할 수 있다.

Dijkstra, BFS, A* 비교(왼쪽부터)

 Occupancy grid와 같은 평면 지도를 작성했다면 위의 그림처럼 경로 계획을 하는 것이 적합한 것임을 이해할 수 있다. 각 셀을 노드로 셀간의 연결(4개 또는 8개)을 에지로 생각하면 그래프와 같기 때문이다.

 Dijkstra에서는 인접 셀에 거리 측정 방법(Manhattan distance, Euclidean distance)으로 측정된 비용에 따라 도착점 까지의 비용이 산출된다. 각 인접셀에 도달하는 비용은 동일하기 때문에 도착점까지의 가능한 모든 경로에 대해서 비용을 산출한 뒤에야 최단 경로가 생성되지만 최단 거리를 보장해준다는 장점이 있다.

 BFS에서는 도착점까지의 거리가 가장 가까운 지점을 먼저 방문한다. Dijkstra보다 빠르지만 그림에서 보는 것처럼 최단경로를 찾지 못할 수 도 있다는 단점이 있다.

 A* search algorithm은 상기 두 가지 방법을 적절히 혼용한 것으로 볼 수 있다. 출발점으로 부터 인접셀의 비용을 계산하되 도착점까지의 거리가 적을 것으로 생각되는 지점을 먼저 계산한다. 아래 그림에서 경로를 계산하는 순서를 확인할 수 있다.

 게시물 작성중에 최근에 새로운 것이 있는지 확인해본 결과 흥미로운 것을 발견했다. 기존의 Dijkstra's algorithm을 GPGPU를 이용해서 비약적으로 성능을 확장시킨 것이다. 가능한 인접셀을 모두 병렬로 방문함으로써 반복의 횟수를 절감했다.