함수최적화 기법 정리 (Levenberg–Marquardt 방법 등)

기계학습 2014. 9. 2. 16:55

※ 참고로, 아래 글 보다는 최근에 올린 [기계학습] - 최적화 기법의 직관적 이해를 먼저 읽어볼 것을 추천합니다. 그 이후에 아래 글을 읽어보시면 좋을 것 같습니다.


----------


원래는 Levenberg-Marquardt 방법에 대한 내용을 쓰고자 했는데 글을 쓰다 보니 함수 최적화 기법들을 총 정리하는 글이 되고 말았습니다.


Levenberg–Marquardt 방법은 비선형 최소 자승(nonlinear least squares) 문제를 푸는 가장 대표적인 방법입니다.


그동안 뉴턴-랩슨법(Newton–Raphson method), 가우스-뉴턴법(Gauss–Newton method), Gradient descent 방법 등 여러 함수 최적화 기법들을 소개한 바 있지만 사실 비선형 함수 최적화 문제에 있어서 가장 널리 쓰이는 방법은 Levenberg–Marquardt 방법입니다.


그런데 Levenberg–Marquardt 방법을 이해하기 위해서는 먼저 가우스-뉴턴법(Gauss–Newton method)과 Gradient descent 방법에 대한 이해가 필요합니다. 왜냐하면 Levenberg–Marquardt 방법은 가우스-뉴턴법(Gauss–Newton method)과 Gradient descent 방법이 결합된 형태로 볼 수 있기 때문입니다. 그리고 당연히 최소자승법에 대한 이해도 필요합니다.


관련하여 먼저 아래의 글들을 읽어보면 도움이 되리라 생각합니다.


[기계학습] - 최소자승법 이해와 다양한 활용예 (Least Square Method)

[기계학습] - 뉴턴법/뉴턴-랩슨법의 이해와 활용(Newton's method)

[기계학습] - Gradient Descent 탐색 방법


이 글에서는 먼저 위 내용들을 포함하여 관련 내용들을 전반적으로 정리한 후 Levenberg-Marquardt 방법에 대해 살펴보도록 하겠습니다.


1. 최소자승법 & 가중최소자승법

2. Gradient Descent 탐색 방법

3. 뉴턴법 (뉴턴-랩슨법)

4. 가우스-뉴턴법

5. Levenberg-Marquardt 방법

6. 마무리하며



1. 최소자승법 & 가중최소자승법


관측값을 (xi, yi), i=1,...,n, 모델 파라미터를 p = ( p1, p2, ..., pm) , 모델을 y = f(x, p), 관측값과 모델과의 오차(residual)를 ri라 할때, 오차 제곱합이 최소화 되도록 모델 파라미터 p를 정하는 방법을 최소자승법(least square method)이라 한다.


 --- (1)


만일 관측값의 신뢰도(중요도)가 서로 다를 경우에는 각각의 오차 제곱에 가중치 wi를 곱해서 최소화시키는 방법을 사용하는데 이러한 문제를 가중최소자승(weighted least squares, WLS)문제라 부른다 (통계 쪽에서는 보통 wi = 1/σi2로 잡음. 단, σi2은 yi의 분산).


 --- (2)


이 때 식 (1), 식 (2)의 에러합 부분을 행렬 형태로 표현하면 다음과 같다.


 --- (3)


 --- (4)


단, r = [r1, ..., rn]T, W는 wi를 대각원소로 갖는 대각행렬 (W = diag(w1, ..., wn)).


최소자승법은 결국 식 (3) 또는 식 (4)에 의해 주어지는 에러 함수를 최소화시키는 p를 구하는 문제로 볼 수 있으며 이는 p에 대한 미분이 E'(p) = 0, Ew'(p) = 0인 p를 구함으로써 얻어진다. 그런데 선형(linear) 최소자승문제의 경우에는 벡터미분을 이용하여 이러한 해를 직접적으로 구할 수 있다 (closed-form solution이 존재).


☞ 최소자승 문제에 있어서 모델 f(x,p)가 모델 파라미터에 대해 선형인 경우 선형 최소자승 문제(linear least squares problem)라 부르고 그렇지 않은 경우 비선형 최소자승문제(non-linear least squares problem)라 부른다. 예를 들어 f(x,p) = p1*sinx + p2*cosx라면 f(x,p) 자체는 비선형 함수이지만 파라미터 p = (p1,p2)에 대해서는 선형이므로 f(x,p)에 대한 최소자승 문제는 선형(linear) 문제가 된다.


만일 f(x,p)가 p에 대해 선형인 경우에는 f(x,p)를 아래와 같이 p에 대한 일차 결합 형태로 표현할 수 있다.


 --- (5)


따라서, 선형 최소자승문제의 경우에는 아래와 같은 행렬 표현이 가능해진다.


 --- (3)'


 ---(4)'


단,


 --- (6)


벡터미분을 이용하여 E(p), Ew(p)를 p로 미분한 후 0으로 놓으면 


 --- (7)


 --- (8)


가 된다(식 (7), (8)의 양변에 transpose를 취한 후 p에 대해 정리, W는 대각행렬이므로 WT=W, 벡터미분에 대해서는 벡터미분과 행렬미분 글 참조).


따라서, 선형 최소자승문제와 선형 가중최소자승문제의 해는 다음과 같이 closed-form 형태로 구해진다.


 --- (9)


--- (10)


[적용예제] 예를 하나 들어보면, 세점 (x1,y1), (x2,y2), (x3.y3)를 삼각함수 f(x) = p1*sin(x) + p2로 근사시키고자 할 때 최소자승 해는 다음과 같이 구해진다.


--- (11)


[음함수 최소자승법] 만일 모델이 양함수(explicit function) 형태가 아니라 f(x,p) = 0와 같은 음함수(implicit function) 형태일 경우에는 최소자승 문제의 형태 및 풀이방법이 조금 달라진다. 이 경우 ri = f(xi,p)이며, f(x,p)가 p에 대해 선형이면 최소자승문제를 ∑ri2 = ∥Ap∥2와 같이 표현할 수 있다. 그리고, ∑ri2 = ∥Ap∥2를 최소화시키는 p는 A의 SVD(singular value decomposition)을 A = U∑VT라 할 때 최소 특이값(singular value)에 대응하는 V의 열벡터가 된다. 자세한 내용은 [선형대수학 #5] 선형연립방정식 풀이 글 참조.


이상과 같이 선형(linear) 최소자승 문제의 경우에는 pseudo inverse나 SVD(singular value decomposition)을 이용하면 해를 손쉽게 구할 수 있다. 하지만 비선형(non-linear) 최소자승 문제의 경우에는 closed-form solution이 없기 때문에 점진적으로 해를 찾는 방법(iterative minimization)을 사용해야 한다. 즉, 어떤 초기값 p0에서 시작해서 점차적으로 E(p)를 최소화시키는 방향으로 p를 갱신하는 방법을 사용해야 한다. 그리고 이러한 iterative minimization 기법으로는 아래에 설명할 gradient descent 방법, 가우스-뉴턴법, Levenberg-Marquardt 방법 등이 있다.



2. Gradient Descent 탐색 방법


어떤 다변수 함수 E(x1,x2,...,xn)의 그레디언트(gradient)는


 --- (12)


로 정의되며 함수값이 가장 가파르게 증가하는 방향을 나타낸다. 그리고 그레디언트의 크기는 그 증가의 가파른 정도를 나타낸다.


☞ 여기서 다변수 함수라는 말을 어렵게 생각할 필요는 없으며 예를 들어 사람의 외적인 아름다움이 키와 팔의 길이, 다리의 길이로 결정된다고 했을 때, '아름다움 = f(키,팔길이,다리길이)'와 같이 함수값이 여러 요인에 의해 결정된다는 의미이다. 키, 팔다리 길이와 아름다움의 관계를 실제 구하는 것은 어렵겠지만 만일 그 관계가 예를 들어 'f(키,팔길이,다리길이) = 키^3 - 팔길이^2 + 다리길이^2'라 한다면 f의 그레디언트는 ∇f = (3*키^2, -2*팔길이, 2*다리길이)가 된다. 따라서 어떤 사람 A의 키, 팔길이, 다리길이가 키 = 1, 팔길이 = 0.5, 다리길이 = 0.7과 같다면 A 위치에서의 그레디언트 값은 (3, -1, 1.4)이므로 현 시점에서 이 사람이 가장 아름다워지는 방향은 3:1:1.4의 비율로 키는 증가하고 팔길이는 감소, 다리길이는 증가하는 방향이다.


Gradient descent 방법은 아래 식과 같이 어떤 초기값 x0 = (x10,...,xn0)부터 시작하여 그레디언트의 반대 방향으로 계속 내려가다 보면 함수의 최소값을 찾을 수 있다는 방법이다 (자세한 내용은 Gradient Descent 탐색 방법 글 참조).


 --- (13)


만일 함수의 최대값을 찾고 싶을 경우에는 xk+1 = xk + λk∇E(xk) 형태의 gradient ascent 방법을 사용한다.


[최소자승 문제에의 적용]

이제 gradient descent 방법을 앞서 설명한 최소자승 문제에 적용해 보자. 선형 최소자승 문제의 경우에는 gradient descent 방법과 같은 근사적 방법을 사용할 필요가 없이 직접 해를 구할 수 있으므로 여기서는 비선형 최소자승 문제에 한해 생각한다. 관측값을 (xi, yi), i=1,...,n, 모델 파라미터를 p = ( p1, p2, ..., pm) , 모델을 y = f(x, p), 오차(residual)를 r(x,p) = y - f(x,p)라 하면 최소화시키고자 하는 목적함수는 다음과 같다.


 --- (14)



 --- (15)


따라서, 비선형 최소자승 문제에 대한 gradient descent 해는 모델 파라미터 p에 대한 초기 추정값 p0=(p10,...,pm0)에서 시작하여 p를 다음 수식에 의해 반복적으로 갱신함으로써 얻어진다.


 --- (16)


단, Jr은 r의 야코비언(Jacobian) 행렬.


참고로 비선형 가중최소자승(weighted least squares) 문제의 경우에 gradient descent 해는 다음 수식에 의해 구해진다.


 --- (17)



3. 뉴턴법 (뉴턴-랩슨법)


뉴턴법(뉴턴-랩슨법)은 어떤 함수 E(x) = 0의 해를 찾기 위해 임의의 초기값 x0로부터 시작하여 다음 수식에 따라 x를 반복적으로(iterative) 갱신함으로써 해를 찾아가는 방법이다.


 --- (18)


그 원리는 미분이 접선의 기울기임을 이용하여 현재의 x에서 E(x)의 부호와 기울기 E'(x)의 부호에 따라 x를 증가시킬지 감소시킬지를 결정하고 그 기울기의 크기에 따라서 x를 얼마나 많이 증가(감소)시킬지를 결정하는 방식이다.


만일 뉴턴법을 함수 E(x)의 최소(최대)값을 찾는 최적화(optimization) 문제에 적용할 경우에는 다음과 같이 E'(x) = 0의 해를 찾는 문제로 식을 변형하여 적용할 수 있다.


 --- (19)


또한 뉴턴법을 다변수 함수의 최적화 문제로 확장하면 식 (19)를 아래와 같이 수정하여 적용할 수 있다.


 --- (20)


또는 아래와 같이 수렴속도를 조절하기 위한 step size를 포함하는 형태도 가능하다 (ν>0).


 --- (21)


단, xk=(x1k,...,xnk), ∇E(xk)는 x = xk에서의 E의 그레디언트(gradient) 값, HE(xk)는 x = xk에서의 E의 헤시안(Hessian) 값을 나타낸다.


 --- (22)


[최소자승 문제에의 적용]

식 (20) 또는 식 (21)을 이용하면 뉴턴법을 비선형 최소자승 문제에 직접 적용하는 것도 가능은 하다. 하지만 E(p) = ∑ri(p)2가 두 번 미분가능해야 하며 E(p)의 2차 미분까지 계산해야 하는 부담이 존재한다. 따라서 비선형 최소자승 문제의 경우에는 뉴턴법 대신에 일차 미분만으로 해의 계산이 가능한 가우스-뉴턴법이 사용된다.



4. 가우스-뉴턴법


가우스-뉴턴법은 일종의 뉴턴법의 변형된 형태로서 비선형 최소자승 문제에 대한 대표적인 최적화 방법 중 하나이다. 앞서 뉴턴법을 최적화 문제에 적용할 경우에는 2차 미분이 필요하지만 가우스-뉴턴법을 이용하면 1차 미분만으로 해를 찾을 수 있다.


앞서와 같이 관측값을 (xi, yi), i=1,...,n, 모델 파라미터를 p = ( p1, p2, ..., pm) , 모델을 y = f(x, p), 오차(residual)를 ri(p) = yi - f(xi,p)라 하고 최소화시키고자 하는 목적함수를 E(p)라 하자.


 --- (23)


이 때, E(p)를 최소화시키는 즉, E'(p) = 0으로 만드는 가우스-뉴턴 해는 모델 파라미터 p에 대한 초기 추정값 p0=(p10,...,pm0)에서 시작하여 p를 다음 수식에 의해 반복적으로 갱신함으로써 얻어진다.


 --- (24)


단, Jr은 Jr(pk)를 간략히 표현한 것으로서 pk에서의 r의 야코비언(Jacobian) 행렬값을 나타낸다.


 --- (25)


식 (24)에서 (JrTJr)-1JrT는 사실 Jr의 pseudo inverse이다. 따라서 식 (24)를 좀더 간략히 표현하면 다음과 같다.


 --- (26)


여기서 가우스-뉴턴법의 식 (26)를 잘 보면 뉴턴법의 식 (18)을 다변수 벡터함수의 경우로 일반화시킨 형태임을 알 수 있다.


참고로, 비선형 가중최소자승(weighted least squares) 문제의 경우에 가우스-뉴턴 해는 다음 수식에 의해 구해진다.


 --- (27)


[가우스-뉴턴법의 원리] 가우스-뉴턴법의 핵심 원리는 비선형 함수를 지역적으로 선형함수로 근사하여 해를 구하는데 있다. 오차벡터 r(p) = [r1(p) ... rn(p)]T를 테일러 전개를 이용하여 pk 근처에서 선형함수로 근사하면 r(p) ~ r(pk) + Jr(pk)(p - pk)가 된다(2차 이상의 항은 무시). 그런데, 선형근사된 오차에 대해 오차 제곱합 ∥r(pk)+Jr(pk)(p-pk)∥2을 최소화시키는 p를 구해보면 식 (26)의 p = pk - pinv(Jr(pk))r(pk)가 나옴을 알 수 있다. 즉, 가우스-뉴턴법은 현재의 파라미터 추정값(pk) 근처에서 에러함수를 선형근사하여 최소자승 해(pk+1)를 구하고, 구해진 해 근처에서 다시 에러 함수를 선형근사하여 최소자승 해를 구하는 방식으로 점진적으로 해를 찾아간다.


[적용예제] 예를 들어 세 점 (x1,y1), (x2,y2), (x3,y3)를 삼각함수 y = cos(ax + b)로 근사시키는 문제를 가우스-뉴턴법으로 풀어보자. 이 경우 모델 파라미터는 p = (a, b)이며 먼저 r(p)과 Jr(p)를 구해보면 다음과 같다.


 --- (28)


이제 적절한 초기 추정값 p0 = (a0, b0)부터 시작하여 다음 식에 따라 p를 수렴할 때까지 갱신함으로써 해를 찾는다.


 --- (29)



5. Levenberg-Marquardt 방법


앞서 언급한 바와 같이 Levenberg–Marquardt 방법은 가우스-뉴턴법(Gauss–Newton method)과 gradient descent 방법이 결합된 형태로서 해로부터 멀리 떨어져 있을 때는 gradient descent 방식으로 동작하고 해 근처에서는 가우스-뉴턴 방식으로 해를 찾는다. 그런데 Levenberg–Marquardt 방법은 가우스-뉴턴법보다 안정적으로 해를 찾을 수 있으며(초기값이 해로부터 멀리 떨어진 경우에도 해를 찾을 확률이 높음) 비교적 빠르게 해에 수렴하기 때문에 비선형 최소자승문제에 있어서는 대부분 Levenberg–Marquardt 방법이 사용된다.


Levenberg-Marquardt 방법은 Levenberg 알고리즘(1944년)을 1963년에 Marquardt가 보완한 것이다. 빠른 비교 및 참조를 위해 gradient descent 방법, 가우스-뉴턴법, Levenberg 방법, Levenberg-Marquardt 방법의 업데이트 식을 함께 적어보면 다음과 같다.


[Gradient descent 방법]

 --- (16)'


[가우스-뉴턴법]

 --- (24)'


[Levenberg 방법]

 --- (30)


[Levenberg-Marquardt 방법]

 --- (31)


단,

 --- (32)


먼저 gradient descent 방법은 그레디언트의 반대방향으로 이동하되 그레디언트의 크기에 비례한 step size 만큼씩 이동하면서 해(에러함수를 최소화시키는 극소점)를 찾아가는 방법이다.


반면에 가우스-뉴턴법은 함수의 그레디언트와 곡률(curvature)을 같이 고려하면서 해를 찾아가는 방식이다(식에서 JrTJr는 이차미분인 Hessian의 의미를 가지며(Hessian에 대한 근사행렬) 함수의 곡률을 나타낸다). 즉, 이동할 step size를 (그레디언트의 크기)/(곡률의 크기)로 결정하는 방식인데, 기울기(그레디언트)가 크더라도 곡률이 크면(기울기의 변화가 급격하면) 조금만 이동하고, 곡률이 작다면(기울기의 변화가 거의 없으면) 좀더 크게 이동함으로써 극소점을 찾아가는 방식이다. 따라서, 가우스-뉴턴법은 gradient descent 방법보다 훨씬 정확하고 빠르게 해를 찾을 수 있다는 장점을 갖는다. 그런데, 가우스-뉴턴법을 살펴보면 계산 과정에 (JrTJr)의 역행렬 계산을 필요로 한다. 따라서 (JrTJr)가 singular 행렬(역행렬이 존재하지 않는 행렬)에 근접한 경우에는 계산된 역행렬이 수치적으로 불안정하여 해가 발산할 수 있는 문제점이 있다.


Levenberg 방법은 가우스-뉴턴법을 개선하여 JrTJr에 항등행렬(identity matrix)의 상수배 μI를 더함으로써 발산의 위험성을 낮추고 보다 안정적으로 해를 찾을 수 있도록 한 방법이다. 상수 μ(μ>0)를 damping factor라 부르는데 μ값이 작으면 Levenberg 방법은 가우스-뉴턴법과 유사해지고 μ값이 크면 gradient descent 방법과 유사해진다. 수식에서 μ→∞ 면 (JrTJr+μI)-1→1/μI이므로 μ가 커지면 Levenberg 방법은 스텝의 크기가 1/μ인 gradient descent 방법과 유사해짐을 알 수 있다. 그런데, Levenberg 방법에서 damping factor μ는 고정된 값이 아니라 매 iteration마다 바뀌는 값으로서 현재 안정적으로 해에 수렴하고 있을 경우에는 작은 값을 주고 해를 잘 찾지 못하고 있을 경우에는 큰 값을 주는 방법을 사용한다. 좀더 구체적으로 설명하면, 현재 단계에서 계산된 에러 E(pk)가 이전의 에러 E(pk-1)에 비해 잘 감소하고 있으면 μk에 작은 값을 사용하여 가우스-뉴턴 방식으로 해를 찾고 오히려 에러가 증가하거나 에러 감소가 충분치 않은 경우에는 μk를 증가시켜서 gradient descent 방식으로 해를 찾는 방식이다.


Levenberg-Marquardt 방법은 기존의 Levenberg의 방법을 1963년에 Marquardt가 좀더 개선시킨 방법을 일컫는다. 원래의 Levenberg 방법은 μ가 큰 경우에 step size가 1/μ인 gradient descent 방법과 유사해짐은 앞서 설명한 바 있다. 하지만 이 경우 수렴 속도가 느린 문제점이 있다. Marquardt는 이러한 문제를 보완하기 위해 항등행렬 I 대신에 diag(JrTJr) 더해주는 방식을 제안하였다(diag(A)는 A의 대각원소는 유지하고 나머지 원소들의 값을 0으로 만든 대각행렬을 나타냄). JrTJr은 원래 헤시안(Hessian)에 대한 근사행렬의 의미를 갖기 때문에 JrTJr의 대각 원소들은 각 파라미터 성분(pi)에 대한 곡률(curvature)를 나타낸다. 즉, Levenberg-Marquardt 방법은 가우스-뉴턴법의 singular 문제를 피하면서도 μ가 큰 경우에도 곡률(curvature)을 반영하여 효과적으로 해를 찾을 수 있도록 한 방법이다.


☞ Levenberg-Marquardt 방법은 가우스-뉴턴법의 damped version으로서 다른 말로 damped least-squares (DLS) 방법이라고도 불린다.


마지막으로 다시 한번 Levenberg-Marquardt 방법을 정리해 보면 아래와 같다.


[Levenberg-Marquardt 알고리즘]

관측값을 (xi, yi), i=1,...,n, 모델 파라미터를 p = ( p1, p2, ..., pm) , 모델을 y = f(x, p), 오차(residual)를 ri(p) = yi - f(xi,p)라 하고 최소화시키고자 하는 목적함수를 E(p) = ∑ri(p)2라 하자. E(p)를 최소화시키는 해는 모델 파라미터 p에 대한 초기 추정값 p0=(p10,...,pm0)에서 시작하여 p를 다음 수식에 의해 반복적으로 갱신함으로써 얻어진다.


 --- (33)


만일 목적함수가 E(p) = ∑wi*ri(p)2 형태인 경우 즉, 가중최소자승(weighted least squares) 문제의 경우에는 업데이트 식이 다음과 같다.


 --- (34)


단, W = diag(w1,...,wn), r=[r1,...,rn]T, p=[p1,...,pm]T, Jr은 rp로 미분한 야코비언(Jacobian) 행렬


 --- (35)


식 (33), (34)에서 μ는 매 iteration마다 변하는 damping factor로서 μ를 어떻게 계산하느냐에 따라서 다양한 방식의 Levenbarg-Marquardt 알고리즘 구현이 가능하다. 구체적인 방법은 조금씩 다르겠지만 그 기본적인 원리는 현재의 μ값을 가지고 p를 갱신했을 때 E(p)가 증가했다면 E(p)가 감소할 때까지 μ를 계속 증가시키고, E(p)가 감소하면 μ를 원래의 기본적으로 설정된 값(최소값)으로 낮추는 원리이다.


실제 구체적인 구현 알고리즘, Matlab 코드 및 C/C++ 구현 코드에 대해서는 아래 자료들을 참조한다 (제가 직접 돌려보거나 그런 것은 아니니 저한테 사용법은 묻지 마시길..).



6. 마무리하며


마지막으로 초기값에 대한 얘기를 하고 글을 마무리합니다. 앞서 나온 gradient descent 방법, Gauss-Newton 방법, Levenberg-Marquardt 방법 등은 공통적으로 파라미터(p)에 대한 초기 추정값을 필요로 합니다. 그리고 기본적으로 이러한 방법들은 local 해를 찾는 방법들이기 때문에 초기값을 어떻게 주느냐에 따라서 그 결과가 달라질 수 있습니다. 따라서 초기값은 가능한 한 실제 해에 가까운 값을 줄수록 좋은 결과를 얻을 수 있습니다. 즉, 정확도는 떨어지더라도 대략적이나마 해를 구할 수 있는 방법이 있다면 먼저 그 방법으로 해를 찾은 후에 그 해를 초기값으로 해서 앞서의 최적화 기법들을 적용하는 것이 효과적입니다. 정 초기값에 대한 정보가 없다면 파라미터 공간을 일정한 간격으로 샘플링한 후 그 중 목적함수를 최소화시키는 파라미터 값을 초기값으로 설정할 수도 있을 것입니다.


☞ 위에 설명한 여러 함수 최적화 기법들의 실제 적용예에 대해서는 [기계학습] - 대수적 에러와 기하학적 에러(Algebraic & Geometric Error) 글을 참조하시기 바랍니다.


☞ 막상 글을 써 놓고 보니 수식만 잔뜩 있는 글이 되고 말았습니다 ㅠㅠ. 그림도 넣고 실제 예제도 돌려보면서 결과를 비교해 보면 더 좋았을 텐데요.. 하지만 글이 너무 길어지고 또 기력이 다한 관계로 이쯤에서 마무리합니다. 너무 이론적인 내용이라 읽기가 쉽지는 않겠지만 공부하는 학생이나 알고리즘 개발자, 연구자 등에게는 알아두면 도움이 되리라 생각합니다.



by 다크 프로그래머