Gauss-Newton Method

1. Introduction

Least Squares Method(최소제곱법 또는 최소자승법)는 주어진 데이터에 가장 잘 맞는 모델의 parameter를 구하는 방법 중 하나로 이름에서 알 수 있듯이 잔차의 제곱합을 목적함수로 한다. Linear regression에서 다루었던 것 처럼, 선형 모델의 경우에는 잔차의 제곱합이 기울기가 0이되는 지점에서 최소값을 갖기 때문에 least squares를 통해 쉽게 해를 구할 수 있다. 하지만 모델이 더이상 선형이 아닐경우에는 기울기가 0이되는 점이 유일하지 않을 수 있다. 

그림 1. 중심이 (0, 0)이고 반지름이 1인 원(red line)으로 부터 sampling된 반원 데이터(blue points). 평균이 0이고 표준편차가 $\sigma=0.01$인 noise를 포함한다.

그림 1과 같이 중심이 (0, 0)이고 반지름의 길이가 1인 원(red line)으로 부터 sampling된 반원 데이터(blue points)가 주어졌다고 해 보자. 원을 a, b, r의 세 개의 parameter로 정의할 때 이 데이터에 가장 잘 맞는 원은 아래와 같은 목적함수를 최소화 하는 원이다.

$$ model : \sqrt{(x_i-a)^2+(y_i-b)^2} = r $$

$$ optimal circle = [a^*, \ b^*, \ r^*] $$

$$ a^*, b^*, r^* =  \underset{a,b,r}{argmin} \sum_i^N{\Big{|} \Big{|} r - \sqrt{(x_i-a)^2+(y_i-b)^2} \Big{|} \Big{|}^2} $$

위 데이터에 대하여 목적 함수를 그려보면 그림 2와 같다. 그림 2에서는 반지름의 길이를 1로 고정하고 원의 중심만을 변수로 두어 목적함수를 plot하였다.

그림 2. 반원 데이터에 대한 목적함수의 가시화. 반지름을 1로 고정한 상태에서 원의 중심값 변화에 따른 목적함수 변화.

그림 2에서 확인할 수 있듯이, 이 문제에서는 (0, 0) 부근에서 global optimal이 존재하고, (0, 1) 부근에서 추가적인 local optimal이 존재한다. 따라서 더이상 목적함수의 gradient가 0이되는점이 하나임이 보장되지 않기 때문에 일반적으로 non-linear 모델의 해를 least squares 방식으로 구할때는 gradient descent와 비슷하게 iteration을 통해 수렴시키는 방법을 사용한다.

2. Newton method

Newton method는 다음과 같은 1D 함수에서 해($x \ where \ f(x)=0$)를 찾는 방법 중 하나이다.

그림 3. Newton method의 예. 초기 시작 점(red point)에서의 함수의 접선이 x축(black line)과 만나는 점(green point)을 다음 시작점의 x좌표로 한다. 그런 다음 해당 x좌표에 해당하는 함수위의 점(orange point)에서 같은 연산을 수행하고, 이 점이 수렴할 때 까지 반복한다.

그림 3에서와 같이 어떤 시작 점(red point)에서의 접선이 x축(black line)과 만나는 점(green point)의 x좌표를 새로운 시작점(orange point)의 x좌표로하여 같은 연산을 수렴할 때 까지 수행한다. 수식은 아래와 같다.

$$ x_{new} = x - \frac{f(x)}{f'(x)} $$

이를 반복하면 아래의 그림 4와 같이 수렴하게 된다.

그림 4. Newton methd를 이용한 $f(x)=0$의 해 찾기

같은 논리로 이 newton method를 다음과 같이 수행하면 $f'(x)=0$인 점을 찾을 수 있다.

$$ x_{new} = x - \frac{f'(x)}{f''(x)}  $$

그림 5. Newton method를 이용한 $f'(x)=0$의 해 찾기. (left) $x_0=-3$에서 시작한 수렴 결과와 (right) $x_0=3$에서 시작한 수렴 결과. 각각 가장 가까운 minimum 값에 수렴한다.

Newton method는 간단한 수식을 통해 최적화를 수행할 수 있으나 이는 그림 5에서 확인할 수 있듯이 초기값에 따라 수렴 결과가 달라질 수 있다는 단점이 있다.

추가로, 2차 미분을 통해 극소값을 구할 수도 있지만 함수의 값이 SSE 에러와 같이 항상 양의 값을 갖는 경우에는 1차 미분으로도 아래의 그림과 같이 해를 구할 수 있다.

그림 6. 함수값이 항상 양의 값을 가질 때 1차 미분을 통한 최소해 찾기의 예

물론, 해에 도달하였을 때 기울기가 0에 가까우므로 그 다음 update 값이 overshooting 될 수 있다. 따라서 적절한 조취가 필요하며 이는 이후에 다루도록 한다.

3. Gauss-Newton method

Gauss-Newton method는 Newton method를 연립방정식으로 확장한 것으로 생각할 수 있다. 

$$ f_1(x_1, x_2, x_3) = 0 \\ f_2(x_1, x_2, x_3) = 0 \\ \vdots \\ f_n(x_1, x_2, x_3) = 0  $$

위와 같은 다변수 연립 방정식이 있을 때, 이를 $ F \textbf{x} = 0$으로 표현할 수 있으며, 여기서 $ F \textbf{x}$가 항상 양의 값을 갖는다고 가정하면, 그림 6에서와 같이 1차 미분을 통해 해를 구할 수 있다.

Matrix $F$의 $\textbf{x}$에 대한 1차 미분은 다음과 같이 jacobian으로 표현할 수 있다.

$$ J = \begin{bmatrix} \frac{\partial{f_1}}{\partial{x_1}} \ \frac{\partial{f_1}}{\partial{x_2}} \ \frac{\partial{f_1}}{\partial{x_3}} \\ \vdots \\ \frac{\partial{f_n}}{\partial{x_1}} \ \frac{\partial{f_n}}{\partial{x_2}} \ \frac{\partial{f_n}}{\partial{x_3}} \end{bmatrix} $$ 

이를 이용하여 Newton method에 적용하면 paramter의 update는 다음과 같이 진행된다.

$$ \textbf{x}_{new} = \textbf{x} - J^{-1}F $$

이 때, 일반적으로 jacobian이 square matrix가 아니므로 pseudo inverse를 이용하여 다음과 같이 update 수식을 나타낼 수 있다.

$$ \textbf{x}_{new} = \textbf{x} - (J^T J)^{-1} J^TF $$

 

4. Circle fitting using Gauss-Newton method

다시 circle fitting 문제로 돌아와서, 이를 앞서 설명한 Gauss-Newton method로 fitting 해보자. 문제를 다시 정의하자면 다음과 같다.

$$ model : \sqrt{(x_i-a)^2+(y_i-b)^2} = r $$

$$ optimal circle = [a^*, \ b^*, \ r^*] $$

$$ a^*, b^*, r^* =  \underset{a,b,r}{argmin} \sum_i^N{\Big{|} \Big{|} r - \sqrt{(x_i-a)^2+(y_i-b)^2} \Big{|} \Big{|}^2} $$

따라서 F를 다음과 같이 설정할 수 있다.

$$ F = \begin{bmatrix} \sqrt{(x_1-a)^2+(y_1-b)^2} -r \\ \vdots \\ \sqrt{(x_n-a)^2+(y_n-b)^2} - r  \end{bmatrix} = \begin{bmatrix} F_1 \\ \vdots \\ F_n \end{bmatrix} $$

그러면, jacobian은 다음과 같이 계산된다.

$$ J = \begin{bmatrix} \frac{\partial{F_1}}{\partial{a}} \ \frac{\partial{F_1}}{\partial{b}} \ \frac{\partial{F_1}}{\partial{r}} \\ \vdots \\ \frac{\partial{F_n}}{\partial{a}} \ \frac{\partial{F_n}}{\partial{b}} \ \frac{\partial{F_n}}{\partial{r}} \end{bmatrix} $$ 

$$ J = \begin{bmatrix} \frac{a-x_1}{\sqrt{(x_1-a)^2+(y_1-b)^2}} \ \frac{b-y_1}{\sqrt{(x_1-a)^2+(y_1-b)^2}} \ -1 \\ \vdots \\ \frac{a-x_n}{\sqrt{(x_n-a)^2+(y_n-b)^2}} \ \frac{b-y_n}{\sqrt{(x_n-a)^2+(y_n-b)^2}} \ -1  \end{bmatrix} $$

따라서, 다음과 같은 반복 update를 통해 특정 parameter에 수렴할 수 있으며 일종의 근사해를 구할 수 있다.

$$ \begin{bmatrix} a_{new} \\ b_{new} \\ r_{new} \end{bmatrix} =  \begin{bmatrix} a \\ b \\ r \end{bmatrix} - (J^T J)^{-1} J^TF$$

위 식을 통해 초기 paramter 값이 각각 [0, -0.5, 0.5], [0, 1, 0.5]일 때 circle fitting을 수행한 결과는 아래와 같다.

 

그림 7. Gauss-Newton 방식을 통한 circle fitting의 예

Reference

[1] Gavin, H.. “The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems.” (2013).

[2] N. Chernov, C. Lesort, "Least Squares Fitting of Circles", Journal of Mathematical Imagging and Vision., Vol. 23, pp. 239-252, 2005.