Main Content

비선형 시스템을 야코비 행렬을 사용하지 않고 풀거나 사용하여 풀기

이 예제에서는 비선형 연립방정식에 도함수를 제공하면 함수 실행 횟수가 줄어드는 것을 보여줍니다. 벡터 및 행렬 목적 함수 작성하기에서 설명한 것처럼, 연립방정식 F(x)의 야코비 행렬 J(x)Jij(x)=Fi(x)xj입니다. 이 도함수를 목적 함수의 두 번째 출력값으로 제공합니다.

예를 들어 multirosenbrock 함수는 n의 양의 짝수 값에 대한 로젠브록 함수(제약 조건이 있는 비선형 문제 풀기, 문제 기반 참조)의 n차원 일반화입니다.

F(1)=1-x1F(2)=10(x2-x12)F(3)=1-x3F(4)=10(x4-x32)F(n-1)=1-xn-1F(n)=10(xn-xn-12).

연립방정식 F(x)=0의 해는 점 xi=1, i=1n입니다.

이 목적 함수의 경우 모든 야코비 행렬 항 Jij(x)ij가 최대 1만큼 차이가 나는 항을 제외하고는 0입니다. i<n인 홀수 값의 경우 0이 아닌 항은 다음과 같습니다.

Jii(x)=-1J(i+1)i=-20xiJ(i+1)(i+1)=10.

이 예제의 마지막 부분에 있는 multirosenbrock 헬퍼 함수는 목적 함수 J(x)와 해당 야코비 행렬 F(x)를 생성합니다.

i<n인 홀수 값의 경우 점 xi=-1.9에서 시작하고 짝수 값 i의 경우 xi=2에서 시작하여 연립방정식을 풉니다. n=64를 지정합니다.

n = 64;  
x0(1:n,1) = -1.9; 
x0(2:2:n,1) = 2;
[x,F,exitflag,output,JAC] = fsolve(@multirosenbrock,x0);
Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

참(True) 해와 계산된 해 x 사이의 거리와 fsolve가 해 계산을 위해 행한 함수 실행 횟수를 검토합니다.

disp(norm(x-ones(size(x))))
     0
disp(output.funcCount)
        1043

fsolve가 해를 구하고, 이를 위해 함수 실행을 1000번 넘게 행합니다.

다시 연립방정식을 풉니다. 이번에는 야코비 행렬을 사용합니다. 이를 위해 'SpecifyObjectiveGradient' 옵션을 true로 설정합니다.

opts = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,F2,exitflag2,output2,JAC2] = fsolve(@multirosenbrock,x0,opts);
Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

참(True) 해와 계산된 해 x2 사이의 거리와 fsolve가 해 계산을 위해 행한 함수 실행 횟수를 다시 검토합니다.

disp(norm(x2-ones(size(x2))))
     0
disp(output2.funcCount)
    21

fsolve가 이전 해와 동일한 해를 반환하되, 이를 위해 함수 실행을 1000번 넘게 행하는 것이 아니라 약 20번 행합니다. 이 예제에서는 견고성이 향상된 것이 드러나지는 않지만 야코비 행렬을 사용하면 대체로 함수 실행 횟수는 줄어들고 견고성이 높아집니다.

헬퍼 함수

다음 코드는 multirosenbrock 헬퍼 함수를 생성합니다.

function [F,J] = multirosenbrock(x)
% Get the problem size
n = length(x);  
if n == 0, error('Input vector, x, is empty.'); end
if mod(n,2) ~= 0
   error('Input vector, x ,must have an even number of components.');
end
% Evaluate the vector function
odds  = 1:2:n;
evens = 2:2:n;
F = zeros(n,1);
F(odds,1)  = 1-x(odds);
F(evens,1) = 10.*(x(evens)-x(odds).^2); 
% Evaluate the Jacobian matrix if nargout > 1
if nargout > 1
   c = -ones(n/2,1);    C = sparse(odds,odds,c,n,n);
   d = 10*ones(n/2,1);  D = sparse(evens,evens,d,n,n);
   e = -20.*x(odds);    E = sparse(evens,odds,e,n,n);
   J = C + D + E;
end
end

참고 항목

관련 항목