MATLAB 을 사용하여 원의 방정식에 대해 가우스-뉴턴 방법 적용 실험

최근에 개인적인 필요에 의해 nonlinear regression 을 접하고 있는데 아래 포스팅을 보면서 개념을 잡아가고 있습니다.

http://darkpgmr.tistory.com/58

nonlinear regression 알고리즘 중 유명한 것이 Levenberg-Marquardt 라고 하는데 Levenberg-Marquardt 알고리즘은 gradient descent 와 Gauss–Newton 방법이 합쳐진 형태라고 하더군요. 그래서 일단 gradient descent 는 아니까 Gauss–Newton 방법이 무엇인지에 대해 공부를 하기 위해 검색을 하던 중 위 블로그를 발견하게 되었고 백문이 불여일타라고 일단 코딩해 봐야 이해가 되니까 위 포스팅에서 맨 아래 있던 예제를 MATLAB 으로 돌려 봤습니다.

예제는 원 방정식에 대해 가우스-뉴턴 방법을 적용해서 파라미터를 찾는 것입니다.

원의 중심과 반지름을 a, b, r 이라고 하고 이를 벡터로 표현하면 x=[a; b; r] 입니다.

원의 방정식(F)는 아시는 바와 같이 아래 수식과 같고~ 아래 수식에서 대문자 X 와 소문자 x 가 다르다는것을 염두에 두시기 바랍니다…. 다른 문자로 쓸걸 그랬나…. ??

위 수식에 대하여 변수 a, b, r 에 대한 자코비안(J)는 아래와 같습니다.

Gauss–Newton 업데이트 수식은 아래와 같습니다. 아래 수식에서 ‘ 는 Transpose 를 의미 하고 inv() 함수는 matrix inverse 함수입니다.

x = x - inv(J(x)' * J(x))*J(x)'*F(x)

MATLAB 에서 inv 함수를 사용하는것 보다는 역슬래쉬(\) 를 사용하는게 더 추천하는 방식이기 때문에 위 수식은 아래와 같이 구현 했습니다.

x = x - (J(x)' * J(x))\J(x)'*F(x)

반지름이 3 이고 중심이 (x,y) = (1,3) 인 원의 방정식을 하나 만들었고 제가 작성한 MATLAB 코드는 아래와 같습니다.

NOP = 100;

center = [1,3];

radius = 3;

style='b-';

THETA=linspace(0,2*pi,NOP);

RHO=ones(1,NOP)*radius;

[X,Y] = pol2cart(THETA,RHO);

X=X+center(1);

Y=Y+center(2);

H=plot(X,Y,style);

grid on

axis square;

axis equal;

다음으로 해당 원의 방정식에 대해 Gauss–Newton 방정식을 적용하는 코드는 아래와 같습니다.

XX=X(:);

YY=Y(:);

x = [0;1;15];

% 원의 방정식에 대한 가우스-뉴턴 방법

for n=1:60

Y=YY(n:(n+10));

X=XX(n:(n+10));

F = sqrt((X-x(1)).^2 + (Y-x(2)).^2)-x(3);

J = [(x(1)-X)./sqrt((X-x(1)).^2 + (Y-x(2)).^2) (x(2)-Y)./sqrt((X-x(1)).^2 + (Y-x(2)).^2) -ones(size(X))];

x = x - (J' * J)\J'*F;

end

x

위 코드에서는 초기값은 0, 1, 15 로 했고 돌려보니 아래와 같이 원의 중심과 반지름을 정확히 추정함을 확인 할 수 있었습니다.

x =

1.0000

3.0000

3.0000

그런데 다들 아시겠지만 Gauss–Newton 방식은 초기값에 따라 발산의 우려가 있다고 합니다. 저도 실험을 해보니 초기값을 다른 값들로 해서 실험을 많이 하다 보니(ex x = [0;8;15];) 로 한 경우에는 추정값이 발산하더군요.