행렬 지수함수

수치해석 계산물리 마지막 수정: 2026-08-08 04:35:52

1. 개요[편집]

행렬 지수함수 Matrix exponential
정의	$e^{A} = \sum_{k=0}^{\infty} A^k/k!$ (모든 $A$ 에서 수렴)
의미	$\dot{x}=Ax,\ x(0)=x_0$ 의 해가 $x(t)=e^{At}x_0$
표준 알고리즘	스케일링-제곱 + 대각 파데 근사
고전 문헌	Moler & Van Loan, Nineteen Dubious Ways (1978 / 2003)
구현	MATLAB `expm` · `scipy.linalg.expm`

스칼라에서 한 줄이던 것이 행렬로 오면 논문 한 편이 된다.

행렬 지수함수(matrix exponential)는 정사각행렬 $A$ 에 대해 스칼라 지수함수의 급수를 그대로 옮긴

e^{A} = \sum_{k=0}^{\infty} \frac{A^k}{k!} = I + A + \frac{A^2}{2!} + \cdots

로 정의되는 행렬함수다. 이 급수는 모든 $A$ 에 대해 절대수렴하므로 정의 자체에는 아무 문제가 없다. 존재 이유도 명확하다. 선형 상미분방정식계 $\dot x = Ax$ 의 해가 정확히 $x(t) = e^{At}x_0$ 이고, 비제차계 $\dot x = Ax + b(t)$ 의 변수변환 공식(뒤아멜 적분)도 $e^{At}$ 로 쓰인다. 구조동역학의 상태공간 형식, 제어의 이산화, 화학 반응망의 마스터 방정식, 리 군의 지수사상까지 전부 여기로 모인다.

기본 성질 몇 개는 챙겨두자. $e^A$ 는 항상 정칙이고 $(e^A)^{-1} = e^{-A}$ , $\det e^A = e^{\operatorname{tr}A}$ , $\frac{d}{dt}e^{At} = Ae^{At} = e^{At}A$ 다. 다만 결정적인 함정이 하나 있다. $AB \neq BA$ 이면 $e^{A+B} \neq e^A e^B$ 다. 이 비가환성이 연산자 분리에서 스플리팅 오차로 나타나고, 베이커-캠벨-하우스도르프 공식이 그 차이를 교환자로 정량화한다.

2. 열아홉 가지 미심쩍은 방법[편집]

Moler와 Van Loan은 1978년 SIAM Review에 “Nineteen Dubious Ways to Compute the Exponential of a Matrix”를 썼고, 25년 뒤 거의 그대로 재간행했다.¹ 제목이 이미 결론이다 — 명백해 보이는 방법 대부분이 실제로는 위험하다.

테일러 급수 직접 합산은 가장 먼저 떠오르고 가장 먼저 배신한다. 논문의 유명한 예제

A = \begin{bmatrix} -49 & 24 \\ -64 & 31 \end{bmatrix}, \qquad e^{A} \approx \begin{bmatrix} -0.7358 & 0.5518 \\ -1.4715 & 1.1036 \end{bmatrix}

를 보자. 고유값은 $-1$ 과 $-17$ 이라 답은 $O(1)$ 인데, $\|A\|_1 = 113$ 이라 급수의 최대 항은 $k \approx 113$ 부근에서 $10^{47}$ 규모까지 부푼다. 그 거대한 항들이 서로 상쇄되어 $O(1)$ 이 나와야 하니, 배정도의 16자리로는 애초에 불가능하다. 단정도로 돌리면 부호까지 틀린 결과가 나온다. 급수가 수렴한다는 사실과 그 급수를 부동소수점으로 더해도 되는지는 완전히 별개 문제라는 것.

고유분해 경로 $e^A = Ve^{\Lambda}V^{-1}$ 도 두 가지로 무너진다. $A$ 가 결함행렬(defective)이면 $V$ 자체가 없고, 존재하더라도 비정규 행렬에서는 $\kappa(V)$ 가 $10^{10}$ 을 넘길 수 있어 오차가 그 배율로 증폭된다. “그럼 조르당 표준형을 쓰면?” — 조르당 형은 항목이 불연속적으로 변하는 표준형이라 부동소수점으로는 계산할 수 없다. 반올림 하나면 $2\times2$ 조르당 블록이 서로 다른 고유값 두 개로 쪼개진다.

특성다항식 경로도 마찬가지다. 케일리-해밀턴 정리에 따라 $e^A$ 는 $I, A, \dots, A^{n-1}$ 의 선형결합으로 쓸 수 있지만, 그 계수를 얻으려면 특성다항식 계수를 계산해야 하고(파데예프-르베리에 등), 다항식 계수는 근에 대해 극단적으로 나쁜 조건을 갖는다. 이론적으로 우아한 길이 수치적으로 최악인 대표 사례.

3. 스케일링-제곱 + 파데 근사[편집]

실무 표준은 지수 법칙 하나를 지렛대로 쓴다. $e^{A} = \left(e^{A/2^s}\right)^{2^s}$ 이므로,

$\|A/2^s\|$ 가 충분히 작아지도록 $s$ 를 고른다.
작아진 행렬에 대각 파데 근사 $r_m(X) = q_m(X)^{-1}p_m(X)$ 를 적용한다.
결과를 $s$ 번 제곱한다.

파데를 쓰는 이유는, 같은 행렬곱 예산에서 테일러보다 정확도가 훨씬 높고 오차가 $O(X^{2m+1})$ 로 떨어지기 때문이다. 게다가 $\|X\|$ 가 작은 영역에서만 쓰이므로 상쇄 문제가 발생하지 않는다. Higham(2005)의 표준 구현은 차수 13 대각 파데를 쓰며, $\|A/2^s\| \le \theta_{13} \approx 5.37$ 이 되도록 $s = \max(0, \lceil \log_2(\|A\|_1/\theta_{13})\rceil)$ 로 잡는다. 비용은 행렬곱 6회 + 선형계 풀이 1회 + 제곱 $s$ 회, 즉 $O(n^3)$ 상수배 몇 개다. MATLAB expm 과 scipy.linalg.expm 이 이 계열이다.

제곱 단계는 공짜가 아니다. 각 제곱에서 상대오차가 최대 두 배로 늘 수 있으므로 $s$ 번 제곱하면 최악의 경우 $2^s$ 배까지 증폭될 여지가 있다. 그래서 $s$ 를 무작정 키워 파데 차수를 낮추는 전략은 나쁘고, 실제 구현은 ” $s$ 를 최소로 유지하면서 파데 차수는 감당 가능한 선에서 최대로”라는 균형점을 찾는다. $\theta_{13}$ 같은 문턱값이 이론 오차 한계에서 역산되어 표로 박혀 있는 이유다.

4. 험프와 조건수[편집]

$A$ 의 모든 고유값이 좌반평면에 있으면 $\|e^{At}\| \to 0$ 이다. 그런데 가는 길이 단조롭지 않다. $t\to\infty$ 의 감쇠율은 스펙트럼 횡좌표 $\alpha(A) = \max\operatorname{Re}\lambda_i$ 가 결정하지만, $t \to 0^+$ 에서의 초기 증가율은 수치적 횡좌표 $\mu(A) = \lambda_{\max}\!\left(\tfrac{1}{2}(A + A^{\top})\right)$ 가 결정한다. 비정규 행렬에서는 $\alpha < 0 \ll \mu$ 가 흔하고, 그 결과 $\|e^{At}\|$ 가 과도구간에 크게 부풀었다가 뒤늦게 감쇠한다. 이것이 험프(hump)다. 유동 안정성에서 “고유값은 전부 안정인데 유한 시간 동안 교란이 수천 배로 자라는” 과도 성장이 정확히 이 현상이다.

험프는 알고리즘도 괴롭힌다. 스케일링으로 $A/2^s$ 를 만들면 제곱 단계가 험프의 꼭대기를 지나가게 되고, 중간 결과의 노름이 최종 답보다 훨씬 커지면서 상대오차가 그 비율만큼 증폭된다(overscaling). Al-Mohy와 Higham(2009)은 스케일 판정에 $\|A\|$ 대신 $\alpha_p(A) = \max\left(\|A^p\|^{1/p}, \|A^{p+1}\|^{1/(p+1)}\right)$ 를 쓰는 개선을 제안했다. 비정규 행렬에서 이 양이 $\|A\|$ 보다 훨씬 작을 수 있어 불필요한 스케일링을 피한다. 현재 SciPy의 expm 이 이 알고리즘이다.²

조건수 자체도 짚고 넘어가야 한다. $e^A$ 의 상대 조건수는 프레셰 미분 $L(A,E)$ 의 노름으로 정의되며, 정규행렬이면 $\|A\|$ 수준이지만 비정규행렬에서는 그보다 몇 자릿수 클 수 있다. 답이 이상하면 알고리즘을 의심하기 전에 조건수를 먼저 재보는 것이 순서다.

5. 정작 필요한 건 $e^A b$ 다[편집]

큰 문제에서 진실은 이렇다. $e^A$ 전체를 원하는 사람은 거의 없다. 강성 방정식을 시간적분하든 반응망을 풀든, 실제로 필요한 것은 특정 벡터에 대한 작용 $e^{A}b$ 다. $A$ 가 $10^6 \times 10^6$ 희소행렬이면 $e^A$ 는 조밀해서 저장조차 불가능하지만 $e^Ab$ 는 벡터 하나다.

표준 접근은 크리로프 부분공간법이다. 아놀디 알고리즘으로 $\mathcal{K}_m(A,b) = \operatorname{span}\{b, Ab, \dots, A^{m-1}b\}$ 의 정규직교기저 $V_m$ 과 $m\times m$ 상헤센베르크 $H_m = V_m^{\top}AV_m$ 을 만든 뒤

e^{A}b \approx \|b\|_2\, V_m\, e^{H_m} e_1

로 근사한다. $m \approx 20$ – $50$ 이면 충분한 경우가 많고, 이 작은 $H_m$ 에는 스케일링-제곱 파데를 마음껏 쓰면 된다. 큰 문제를 작은 문제로 투영하는 크리로프 부분공간법의 정신 그대로다.³ 파데조차 쓰지 않고 스케일링 + 절단 테일러만으로 $e^Ab$ 를 직접 계산하는 Al-Mohy–Higham(2011)의 expmv 계열도 널리 쓰인다(SciPy의 expm_multiply).

$e^{A}$ 전체가 정말로 필요한 대표적 상황은 제어계의 이산화다. 연속시간계 $\dot x = Ax + Bu$ 를 샘플링 주기 $\Delta t$ 로 영차 유지(ZOH) 이산화하면 $x_{k+1} = A_d x_k + B_d u_k$ 이고, $A_d = e^{A\Delta t}$ , $B_d = \int_0^{\Delta t} e^{A\tau}d\tau\, B$ 다. 두 행렬을 따로 구하는 대신 확대행렬 $\begin{bmatrix} A & B \\ O & O\end{bmatrix}\Delta t$ 의 지수 하나를 계산해 블록으로 읽어내는 것이 표준 요령이다. 상태 차원이 수십 정도라 $O(n^3)$ 이 전혀 부담되지 않는 세계이고, 그래서 여기서는 expm 을 그냥 부르면 된다.

이 기계 위에 얹힌 것이 지수 적분기다. 선형 강성항을 $\varphi$ 함수(지수함수의 사촌들)로 정확히 처리하고 비선형 항만 명시적으로 다루므로, 강성 때문에 시간 스텝이 묶이는 문제에서 안정성 제약을 크게 풀 수 있다. 물론 스텝마다 행렬함수 작용을 계산해야 하니 스텝당 비용은 훨씬 비싸다 — 스텝 수를 줄여 얻는 이득이 그 비용을 넘을 때만 남는 장사다.

6. 관련 문서[편집]

후속판 제목이 “Twenty-Five Years Later”인데, 25년 동안 방법 개수는 19개 그대로였다. 저자들이 서문에서 밝히길 “여전히 미심쩍다”는 판정도 대부분 유지됐다. 수치해석 논문 중 제목만으로 결론이 전달되는 드문 사례. ↩
“노름을 재는 방법을 바꿨더니 정확도가 올라갔다”는 문장은 수치해석 바깥 사람에게는 사기처럼 들린다. 하지만 스케일 $s$ 는 정수라 한 칸 잘못 잡으면 제곱 횟수가 통째로 하나 달라지고, 비정규행렬에서 $\|A\|$ 는 실제 필요보다 몇 배 과대평가된 스케일을 요구한다. 재는 자를 바꾸는 것이 곧 알고리즘을 바꾸는 것이다. ↩
다만 $A$ 가 비정규이고 스펙트럼이 넓게 퍼져 있으면 크릴로프 수렴이 눈에 띄게 느려진다. 이럴 때는 시간 구간을 잘게 쪼개 $e^{A\Delta t}$ 를 여러 번 적용하는 스텝 분할이 오히려 총 행렬-벡터 곱 수를 줄인다. “한 방에 끝내는 것이 항상 싸지는 않다”는, 이 바닥에서 반복되는 교훈. ↩