논문

하루에도 수만개의 글자를 읽고 있습니다. 하루에도 수백장의 종이를 들춰 읽습니다.
이것은 그 읽기에 대한 일기입니다.

The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems

H. P. Gavin, The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems (http://people.duke.edu/~hpgavin/ce281/lm.pdf)

1 Introduction

주어진 함수 $ \hat{y} (t; p) $의 파라미터 $ p $ 를 결정하고자 할 때, 관측된 데이터 $ (t_i, y_i) $ 가 어느 한 파라미터로 결정된 함수와 얼마나 잘 들어맞는지를 측정하기 위한 방법으로 관측된 데이터 $y_i$와 함수 $ \hat{y}(t_i) $와의 오차를 가중치 제곱의 합을 정의이용하는 것은 아주 일반적인 이야기이다. 이 스칼라 값을 chi-squared error 라고 부르며 수식으로 표현하면 아래와 같이 표현된다.

$$\chi ^ 2(p) = \sum^m_{i=1}\left[ \frac{y(t_i) – \hat{y}(t_i ; p)}{w_i}\right] \\
= \mathbf{(y – \hat{y}(p))^TW(y- \hat{y}(p))} \\
= \mathbf{y^TWy – 2y^Tw\hat{y} + \hat{y}^T w \hat{y}}$$

여기서 $ w_i $는 관측 시의 오차 수준을 나타낸다. 수식에 의해 행렬 W는 $ W_{ii} = w_i $를 성분으로 갖는 대각행렬이다. 우리가 원하는 적절한 파라미터 p는 이 chi-squared error를 최소화하는 p로, $ \chi^2 $를 작아지게끔 현재 p값을 h만큼 반복적으로 변화시켜가며 찾고자 한다.

문제는 현재의 p값에서 $ \chi^2 $를 작아지게 하는 h를 어떻게 찾는지다.

2 The Gradient Descent Method

단순하면서도 어느 정도 작동하는 이 방법은, $ \chi^2 $의 그래디언트를 계산하여 이를 직접 이용한다.

$$\frac{\partial}{\partial \mathbf{p}} \chi^2 = \mathbf{ \mathbf{(y – \hat{y}(p))^T W \frac{\partial}{\partial \mathbf{p}} (y- \hat{y}(p))}} \\
= – \mathbf{ \mathbf{(y – \hat{y}(p))^T W \frac{\partial \hat{y}(p)}{\partial \mathbf{p}}}} \\
= – \mathbf{ (y-\hat{y})^T W J}$$

값이 낮아져야 하기 때문에 그래디언트의 반대 방향과 변화시킬 정도를 나타내는 $ \alpha $를 이용하여 h를 결정한다.

$$\mathbf{ h_{gd} = \alpha J^T W (y-\hat{y}) }$$

3 The Gauss-Newton Method

Gauss-Newton 방법은 함수를 quadratic으로 간주하고 접근하는 방식이다. 먼저 $ \hat{y} $의 Taylor series를 구해보면,

$$\mathbf{ \hat{y}(p+h) \approx \hat{y}(p) + \left[ \frac{\partial \hat{y}}{\partial p} \right]h = \hat{y} + J h }$$

와 같다. 이 식을 이용하여 $ \chi^2(\mathbf{p + h}) $를 계산한다. 이는 다시 말하면 현재 p에 기반하여 p+h의 $ \chi^2 $의 추정 값을 계산한 것이다. $ \chi^2 $가 최소화가 되는 값을 가지려면 $ \partial \chi^2 / \partial \mathbf{h} = 0 $ 이어야 한다.

이를 계산하면,

$$\frac{\partial}{\partial \mathbf{h}} \chi^2 \mathbf{(p+h) \approx -2(y – \hat{y})^T W J + 2h^T J^T W J} = 0$$

이고, 이를 정리하면 변화시켜야 할 h는 아래와 같이 구할 수 있다.

$$\mathbf{[J^T W J] h_{gn} = J^T W (y – \hat{y}) }$$

4 The Levenberg-Marquardt Method

Levenber-Marquardt 방법은 변화시켜야 할 h를 구하기 위해서 위의 gradient descent 방법과 Gauss-Newton 방법을 혼합하여 사용한다.

$$\mathbf{ [J^TWJ + \lambda I] h_{lm} = J^T W (y – \hat{y})}$$

여기서 $ \lambda $가 작은 값이면 Gauss-Newton이 큰 가중치를, 큰 값이면 gradient descent에 더 영향을 많이 받는 값으로 h가 계산된다. 계산을 반복하는 가운데, 계산이 잘못되고 있다고 판단되면 $ \lambda$ 를 증가시기고 잘 되고 있으면 이를 감소시킬 것이다.

4.1 Numerical Implementation

사실 Levenberg Marquardt의 많은 변종들이 있지만, 이 논문에서는 그 중 하나만 다루고 있다. 먼저 h를 구하는 방법으로는 아래의 식을 사용한다.

$$\mathbf{ [J^TWJ + \lambda \text{diag}(J^T W J)] h_{lm} = J^T W (y – \hat{y})}$$

구해진 h를 평가하는 방법으로는,

$$\rho_i(\mathbf{h}) = (\chi^2 (\mathbf{p}) – \chi^2(\mathbf{p+h})) / (2\mathbf{h^T(\lambda_i h + J^TW(y-\hat{y}(p)))})$$

를 사용한다. 계산 전에 미리 정해놓은 값 $ \epsilon $ 보다 $ \rho(\mathbf{h}) $가 크면 p + h 는 이전보다 좋은 위치라고 판단하고, 그보다 작으면 계산 결과가 좋지 않다고 판단할 것이다. 그리고 그에 따라 $ \lambda $를 변화시킬 것이다.

4.1.1 Initialization and update of L-M parameter, $ \lambda $, and parameters p

h를 계산하는 방법은 아래의 3가지 방법 중 하나를 선택하여 사용한다.

  1. 미리 $ \lambda_0 $ 를 초기화한다.
    h를 $ \mathbf{ [J^TWJ + \lambda \text{diag}(J^T W J)] h_{lm} = J^T W (y – \hat{y})} $으로 계산 한 뒤,
    만약 $ \rho_i(\mathbf{h}) \gt \epsilon $ 이면, $ \mathbf{p \leftarrow p + h} $, $ \lambda_{i+1} = \max[\lambda_i / L_{\downarrow}, 10^{-7}]$
    그렇지 않으면 $ \lambda_{i+1} = \min[\lambda_i L_{\uparrow}, 10^7] $
  2. 미리 $ \lambda_0 = \lambda_0 \max[\text{diag}[\mathbf{J^T W J }]] $, $ \lambda_0$를 초기화한다.
    h를 $ \mathbf{ [J^TWJ + \lambda I] h_{lm} = J^T W (y – \hat{y})} $으로 계산한 뒤,
    추가적으로 $ \mathbf{ \alpha = ((J^T W (y – \hat{y}(p)))^T h) /
    (( \chi^2(p+h) – \chi^2(p)) / 2 + 2(J^T W (y – \hat{y}(p))))^Th} $
    를 계산한다.
    만약 $ \rho_i(\alpha\mathbf{h}) \gt \epsilon $ 이면, $ \mathbf{p \leftarrow p + \alpha h} $, $ \lambda_{i+1} = \max[\lambda_i / (1+\alpha), 10^{-7}]$
    그렇지 않으면 $ \lambda_{i+1} = \lambda_i + |\chi^2(p + \alpha h) – \chi^2(p) | / (2 \alpha) $
  3. 미리 $ \lambda_0 = \lambda_0 \max[\text{diag}[\mathbf{J^T W J }]] $, $ \lambda_0$를 초기화한다.
    h를 $ \mathbf{ [J^TWJ + \lambda I] h_{lm} = J^T W (y – \hat{y})} $으로 계산한 뒤,
    만약 $ \rho_i(\alpha\mathbf{h}) \gt \epsilon $ 이면, $ \mathbf{p \leftarrow p + h} $, $ \lambda_{i+1} = \max[1/3, 1- (2 \rho_i – 1)^3], \nu = 2 $
    그렇지 않으면 $ \lambda_{i+1} = \lambda_i \nu_i, \nu_{i+1} = 2 \nu_i $

4.1.2 Computation and rank-1 update of the Jacobian, $ [\partial \mathbf{y} / \partial \mathbf{p} ] $

우리가 함수의 모양을 알지 못하기 때문에 위치가 이동할 때마다 이전에 계산했던 Jacobian 행렬은 올바르지 못하게될 가능성이 높다. 따라서, 처음 계산할 때, 짝수번째 반복에, 그리고 $ \chi^2(p+h) \gt \chi^2(p) $ 일 때, Jacobian 행렬을 재계산해준다.

위의 조건에 해당하지 않을 경우에는, Broyden의 rank-1 방법을 이용하여 Jacobian 행렬의 추정 값으로 갱신하여준다.

$$\mathbf{ J = J + ((\hat{y}(p+h) – \hat{y}(p) – Jh)h^T) / (h^T h)}$$

4.1.3 Convergence criteria

앞에서 이야기한 방법대로 반복적으로 p를 변화시켜가며 최소지점을 찾다가 다음의 조건을 하나라도 만족하면 최소지점을 찾았다고 판단하고 계산을 중단한다.

  • 그래디언트에서 수렴하였을 경우 : $ \max | \mathbf{J^T W(y – \hat{y})| \lt \epsilon_1 } $
  • 파라미터가 수령하였을 경우 : $ \max | h_i / p_i | \lt \epsilon_2 $
  • $ \chi^2 $가 수렴하였을 경우 : $ \chi^2 / (m -n + 1) \lt \epsilon_3 $

4.4 Numerical Examples

추가적으로 3개의 예제를 설명하고 있으니 논문을 참고하도록 하자.


Add a Comment Trackback

2 Comments

  1. 3. Gauss-Newton method에서 p값을 최소화 시키는게 아니라 chi^2를 최소화 되는 조건인데
    오타가 있는 것 같습니다.

  2. 지적 감사합니다!
    급하게 작성하는 경우가 많다보니 전반적으로 꼼꼼하지 못한 부분이 많네요.
    바로 수정하였습니다. :)

댓글 남기기

This site uses Akismet to reduce spam. Learn how your comment data is processed.