Ref. Open: Least Square Minimization

2014년 8월 29일 금요일

최소자승법

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\)을 주어진 식의 수라고 하면,

\(m=n\)(자명하지 않은 유일해 - nontrivial solution)
\(m<n\)인 경우 (부정 - overdetermined)
\(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월 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