Main Content

행렬 지수

이 예제에서는 행렬의 지수를 계산하는 19가지 방법 중 세 가지를 보여줍니다.

행렬 지수 계산에 대한 배경 정보를 알아 보려면 다음 문헌을 참조하십시오.

Moler, Cleve, and Charles Van Loan. “Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later.” SIAM Review 45, no. 1 (January 2003): 3–49. https://doi.org/10.1137/S00361445024180.

먼저 행렬 A를 만들어 보겠습니다.

A = [0 1 2; 0.5 0 1; 2 1 0]
A = 3×3

         0    1.0000    2.0000
    0.5000         0    1.0000
    2.0000    1.0000         0

Asave = A;

방법 1: 스케일링과 스퀘어링(Scaling and Squaring)

expmdemo1은 다음 문헌의 알고리즘 11.3.1을 구현한 것입니다.

Golub, Gene H. 및 Charles Van Loan. Matrix Computations, 3rd edition. Balitmore, MD: Johns Hopkins University Press, 1996.

% Scale A by power of 2 so that its norm is < 1/2 .
[f,e] = log2(norm(A,'inf'));
s = max(0,e+1);
A = A/2^s;

% Pade approximation for exp(A)
X = A;
c = 1/2;
E = eye(size(A)) + c*A;
D = eye(size(A)) - c*A;
q = 6;
p = 1;
for k = 2:q
   c = c * (q-k+1) / (k*(2*q-k+1));
   X = A*X;
   cX = c*X;
   E = E + cX;
   if p
     D = D + cX;
   else
     D = D - cX;
   end
   p = ~p;
end
E = D\E;

% Undo scaling by repeated squaring
for k = 1:s
    E = E*E;
end

E1 = E
E1 = 3×3

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

방법 2: 테일러 급수

expmdemo2는 멱급수로 지정된 행렬 지수에 대한 전형적인 정의를 사용합니다.

eA=k=01k!Ak.

A0A와 차원이 동일한 단위 행렬입니다. 실용적인 수치 계산법 측면에서 이 방법은 norm(A)가 너무 큰 경우 실행 속도가 느려지고 정확도가 떨어집니다.

A = Asave;

% Taylor series for exp(A)
E = zeros(size(A));
F = eye(size(A));
k = 1;

while norm(E+F-E,1) > 0
   E = E + F;
   F = A*F/k;
   k = k+1;
end

E2 = E
E2 = 3×3

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

방법 3: 고유값(Eigenvalue)과 고유벡터(Eigenvector)

expmdemo3은 행렬에 A=VDV-1을 충족하는 고유벡터의 완전한 집합 V가 있다고 가정합니다. 고유값의 대각 행렬을 거듭제곱하여 행렬 지수를 계산할 수 있습니다.

eA=VeDV-1.

이는 실용적인 수치 계산법 측면에서 고유벡터 행렬의 조건에 따라 정확도가 결정됩니다.

A = Asave;

[V,D] = eig(A);
E = V * diag(exp(diag(D))) / V;

E3 = E
E3 = 3×3

    5.3091    4.0012    5.5778
    2.8088    2.8845    3.1930
    5.1737    4.0012    5.7132

결과 비교하기

이 예제의 행렬에서는 세 가지 방법이 모두 효과적입니다.

E = expm(Asave);
err1 = E - E1
err1 = 3×3
10-14 ×

    0.3553    0.1776    0.0888
    0.0888    0.1332   -0.0444
         0         0   -0.2665

err2 = E - E2
err2 = 3×3
10-14 ×

         0         0   -0.1776
   -0.0444         0   -0.0888
    0.1776         0    0.0888

err3 = E - E3
err3 = 3×3
10-13 ×

   -0.0711   -0.0444   -0.0799
   -0.0622   -0.0488   -0.0933
   -0.0711   -0.0533   -0.1066

테일러 급수 실패

일부 행렬에서는 테일러 급수의 항이 0이 되기 전에 매우 커지며, 따라서 expmdemo2가 실패합니다.

A = [-147 72; -192 93];
E1 = expmdemo1(A)
E1 = 2×2

   -0.0996    0.0747
   -0.1991    0.1494

E2 = expmdemo2(A)
E2 = 2×2
106 ×

   -1.1985   -0.5908
   -2.7438   -2.0442

E3 = expmdemo3(A)
E3 = 2×2

   -0.0996    0.0747
   -0.1991    0.1494

고유값 및 고유벡터 실패

다음은 고유벡터의 완전한 집합이 없는 행렬입니다. 따라서 expmdemo3이 실패합니다.

A = [-1 1; 0 -1];
E1 = expmdemo1(A)
E1 = 2×2

    0.3679    0.3679
         0    0.3679

E2 = expmdemo2(A)
E2 = 2×2

    0.3679    0.3679
         0    0.3679

E3 = expmdemo3(A)
E3 = 2×2

    0.3679         0
         0    0.3679

참고 항목