행렬함수

1. 개요

행렬함수 Matrix function
대상	정사각 $A$ 에 대해 스칼라 $f$ 를 $f(A)$ 로 확장
동치 정의	조르당 형 · 스펙트럼 위 에르미트 보간 · 코시 적분
핵심 개념	1차(primary) 행렬함수, 프레셰 미분 $L(A,E)$
범용 알고리즘	슈어-파레 (funm), 스케일링-제곱
대형 희소	$f(A)b$ — 크릴로프 · 유리 크릴로프 · 등고선 적분

$\sin(A)$ 를 물어보면 대부분 성분마다 사인을 취한 것을 떠올린다. 그건 행렬함수가 아니다.

행렬함수(matrix function)는 스칼라 함수 $f$ 와 정사각행렬 $A \in \mathbb{C}^{n\times n}$ 에 대해 정의되는 같은 크기의 행렬 $f(A)$ 로, $A$ 의 스펙트럼 위에서 $f$ 및 그 도함수와 값이 일치하는 다항식을 $A$ 에 대입한 결과다. 성분별 적용과는 완전히 다른 대상이라, MATLAB조차 이쪽에 funm·expm·sqrtm·logm 이라는 별도 이름을 준다.¹

정의가 다항식으로 환원된다는 것이 결정적이다. 케일리-해밀턴 정리에 따라 $A$ 의 거듭제곱은 전부 $I, A, \dots, A^{n-1}$ 의 선형결합이므로, 무한급수로 써 놓은 $f(A)$ 도 실은 차수 $n-1$ 이하 다항식이다. 그래서 $f(A)$ 는 $A$ 와 항상 교환하고, $A$ 의 모든 불변 부분공간을 보존하며, 고유값은 $f(\lambda_i)$ 로 사상되고 고유벡터는 그대로다.

가장 유명한 사례는 행렬 지수함수지만 $A^{1/2}$, $\log A$, $\operatorname{sign}(A)$, $A^{-1/2}$, 지수 적분기의 $\varphi$ 함수도 못지않게 굴러다닌다. 이 문서는 그 전체를 관통하는 이론과 계산 메뉴를 다루고, $e^A$ 개별 사정은 해당 문서에 맡긴다.

2. 세 가지 정의, 하나의 함수

교과서는 전혀 달라 보이는 정의 셋을 제시하는데, $f$ 가 스펙트럼 위에서 충분히 미분 가능하면 셋 다 같은 행렬을 준다.

① 조르당 형 정의. $A = ZJZ^{-1}$, $J = \operatorname{diag}(J_1,\dots,J_p)$ 로 두고 각 조르당 블록에 대해

$$f(J_k) = \begin{bmatrix} f(\lambda_k) & f'(\lambda_k) & \cdots & \dfrac{f^{(m_k-1)}(\lambda_k)}{(m_k-1)!} \\ & f(\lambda_k) & \ddots & \vdots \\ & & \ddots & f'(\lambda_k) \\ & & & f(\lambda_k) \end{bmatrix}, \qquad f(A) = Z f(J) Z^{-1}$$

로 정의한다. 도함수가 필요하다는 점이 중요하다 — 크기 $m_k$ 인 조르당 블록이 있으면 $f$ 는 그 고유값에서 $m_k-1$ 번 미분 가능해야 한다. 대각화 가능하면 $f(A) = Z\operatorname{diag}(f(\lambda_i))Z^{-1}$ 로 단순해진다.

② 에르미트 보간 정의. $A$ 의 최소다항식이 지정하는 근과 중복도를 보고, 스펙트럼 데이터 $f^{(j)}(\lambda_i)$ 를 전부 보간하는 유일한 최소 차수 다항식 $\psi$ 를 만들어 $f(A) := \psi(A)$ 로 둔다. 고유값이 모두 다를 때 이 다항식을 라그랑주 형태로 쓴 것이 실베스터 공식이다.

③ 코시 적분 정의. $f$ 가 $A$ 의 스펙트럼을 감싸는 영역에서 해석적이면

$$f(A) = \frac{1}{2\pi i}\oint_{\Gamma} f(z)\,(zI - A)^{-1}\,dz$$

로 정의한다. $\Gamma$ 는 스펙트럼을 한 번 감는 폐곡선. 섭동 해석과 노름 한계 증명에서 가장 쓸모 있고, 동시에 계산 알고리즘으로 직접 번역되는 유일한 정의다 — 아래 등고선 적분법 참조.

3. 1차 행렬함수와 가지

위 정의들은 모두 1차 행렬함수(primary matrix function)를 준다. 다가함수의 가지 하나를 골라 모든 고유값에 일관되게 적용하고, 같은 고유값에는 반드시 같은 값을 준다는 뜻이다. $f$ 가 단일가면 여기서 끝이지만, 제곱근이나 로그처럼 가지가 여럿이면 선택지가 생긴다.

고유값이 전부 서로 다르고 $A$ 가 정칙이면 각 $\lambda_i$ 에서 $\pm\sqrt{\lambda_i}$ 를 고르는 만큼, 즉 $2^n$ 개의 1차 제곱근이 존재한다. 전부 주값(principal branch)을 고른 것이 주 제곱근 $A^{1/2}$ 이고, $A$ 가 음의 실축 위에 고유값을 갖지 않으면 유일하게 존재하며 스펙트럼이 우반평면에 놓인다.

비1차 함수는 고유값이 겹칠 때, 그것도 비유도(derogatory) 행렬일 때만 생긴다. $A = I_2$ 의 1차 제곱근은 $\pm I$ 둘뿐인데, 임의의 정칙 $Z$ 에 대해 $Z\operatorname{diag}(1,-1)Z^{-1}$ 도 전부 제곱근이라 실제로는 연속체만큼 많다. 같은 고유값에 다른 가지를 배정했으니 1차가 아니고, $A$ 와 교환하지도 않는다. 반대로 제곱근이 아예 없을 수도 있다 — $\begin{bmatrix}0&1\\0&0\end{bmatrix}$ 은 어떤 행렬의 제곱도 아니다.²

로그는 절단면이 더 사납다. 주 로그 $\log A$ 는 고유값이 닫힌 음의 실축 $(-\infty,0]$ 에 없을 때만 정의되고, 고유값이 그 근처를 지나가면 조건이 급격히 나빠진다. $\log(AB) = \log A + \log B$ 는 일반적으로 거짓이며(행렬 지수함수의 비가환성과 같은 뿌리), $\log(e^A) = A$ 도 고유값 허수부가 전부 $(-\pi,\pi)$ 안에 있을 때만 보장된다.

4. 프레셰 미분과 조건수

”$A$ 를 조금 흔들면 $f(A)$ 는 얼마나 흔들리나”는 질문에 답하는 것이 프레셰 미분이다. $L_f(A,\cdot)$ 는

$$f(A+E) = f(A) + L_f(A,E) + o(\|E\|)$$

를 만족하는 $E$ 에 대한 선형 사상이다. 스칼라와 결정적으로 다른 점은 $A$ 와 $E$ 가 교환하지 않으면 $L_f(A,E) \ne f'(A)E$ 라는 것. 지수함수의 경우 $L_{\exp}(A,E) = \int_0^1 e^{A(1-s)}E\,e^{As}ds$ 라는 적분 표현이 되어 곱하기 한 번으로 끝나지 않고, 제곱근의 경우는 실베스터 방정식 $X^{1/2}L + LX^{1/2} = E$ 의 해다.

상대 조건수는 여기서 바로 정의된다.

$$\operatorname{cond}(f, A) = \frac{\|L_f(A)\|\,\|A\|}{\|f(A)\|}, \qquad \|L_f(A)\| = \max_{E \ne 0}\frac{\|L_f(A,E)\|}{\|E\|}$$

이 값이 $10^{12}$ 인 문제에서 배정도로 열 자리를 기대하는 것은 알고리즘 탓이 아니라 문제 탓이다. 후진 안정한 알고리즘조차 전진 오차는 $\operatorname{cond}(f,A)\cdot u$ 까지 허용된다는 것이 후진 오차 해석의 결론이고, 답이 이상하면 조건수부터 재봐야 하는 이유다. 덤으로 $L_f$ 자체가 쓸모 있다 — $f(A(\theta))$ 를 설계변수로 미분할 때 자동 미분 대신 닫힌 형을 쓰면 훨씬 싸다.

5. 계산 메뉴

5.1. 슈어-파레 — 범용 해법

$A = QTQ^{H}$ 로 슈어 분해하면 $f(A) = Q f(T) Q^{H}$ 이므로 상삼각행렬의 함수만 구하면 된다. 유니터리 변환은 조건수가 1이라 오차를 증폭하지 않고, 조르당 형과 달리 결함 행렬에서도 계산 가능하다.

$F = f(T)$ 의 대각은 $f_{ii} = f(t_{ii})$ 로 즉시 나오고, 비대각은 $FT = TF$ 를 성분별로 풀어 얻는 파레 점화식으로 채운다.

$$f_{ij} = \frac{t_{ij}(f_{ii} - f_{jj}) + \sum_{k=i+1}^{j-1}\left(f_{ik}t_{kj} - t_{ik}f_{kj}\right)}{t_{ii} - t_{jj}}$$

문제는 분모다. 고유값 두 개가 가까우면 $t_{ii} - t_{jj} \approx 0$ 이라 점화식이 폭발한다. 데이비스-하이엄(2003)의 처방은 두 단계다 — 재정렬로 가까운 고유값끼리 한 대각 블록에 모으고(순서화 슈어 분해, 분리 파라미터 보통 $\delta = 0.1$), 대각 블록은 블록 중심에서의 테일러 급수로 직접 평가한 뒤 블록 사이만 실베스터 방정식으로 채운다. 비용은 스펙트럼이 잘 분리되면 $\approx 28n^3$ flops, 뭉쳐서 큰 블록이 생기면 그보다 훨씬 비싸다. 이 블록화 버전이 오늘날 슈어-파레 방법이라 부르는 것이고 MATLAB funm 이 그 구현이며, $f$ 의 고차 도함수를 사용자가 제공해야 한다는 것이 실질적 제약이다.

5.2. 스케일링-제곱과 유리근사

특정 $f$ 에는 전용 알고리즘이 훨씬 낫다. 공통 아이디어는 함수방정식으로 인수를 작게 만들고 파데 근사를 쓴 뒤 되돌리기다.

• $e^{A} = \left(e^{A/2^{s}}\right)^{2^{s}}$ — 스케일링-제곱. 하이엄(2005)의 차수 13 대각 파데가 표준.
• $\log A = 2^{s}\log\!\left(A^{1/2^{s}}\right)$ — 역스케일링-제곱. 제곱근을 반복해 $A$ 를 $I$ 쪽으로 민 뒤 $\log(1+x)$ 의 파데를 쓴다.
• 행렬 제곱근 $A^{1/2}$ — 뉴턴 반복 $X_{k+1} = \tfrac12(X_k + X_k^{-1}A)$ 는 우아하지만 수치적으로 불안정하다. 실무 표준은 슈어 형 위에서 삼각 제곱근을 재귀적으로 채우는 뵈르크-함마를링(1983) 방법.
• 행렬 부호 함수 $\operatorname{sign}(A)$ — 뉴턴 반복 $X_{k+1} = \tfrac12(X_k + X_k^{-1})$ 이 2차 수렴한다.

5.3. $f(A)b$ — 큰 문제의 유일한 길

$A$ 가 $10^{6}\times10^{6}$ 희소행렬이면 $f(A)$ 는 조밀해서 저장조차 안 된다. 다행히 실제로 필요한 것은 대개 벡터 하나에 대한 작용 $f(A)b$ 다.³

• 크릴로프 투영. 아놀디 알고리즘으로 $\mathcal{K}_m(A,b)$ 의 정규직교기저 $V_m$ 과 $H_m = V_m^{H}AV_m$ 을 만들고 $f(A)b \approx \|b\|_2 V_m f(H_m)e_1$ 로 근사한다(대칭이면 란초스 알고리즘의 3항 점화식). 다만 $\sqrt{\cdot}$ 나 $\log$ 처럼 원점 근처에서 특이한 $f$ 는 다항식 근사가 태생적으로 나빠 수렴이 느리다.
• 유리 크릴로프. 그 약점을 정면으로 치는 것이 극점 $\sigma_j$ 를 붙인 부분공간 $\operatorname{span}\{b, (A-\sigma_1 I)^{-1}b, \dots\}$ 다. 이동-역변환 한 번이 전처리기 붙은 선형계 풀이 한 번이라 반복당 비용은 비싸지만, $A^{-1/2}b$ 나 $\log(A)b$ 에서 반복 수가 한 자릿수 줄어드는 일이 흔하다.
• 등고선 적분. 코시 정의를 그대로 구적하면 $f(A)b \approx \sum_{j} w_j\,(z_j I - A)^{-1}b$ — 이동된 선형계 몇 개로 환원되고, 서로 독립이라 완벽히 병렬화된다. 직선 위 사다리꼴 구적은 대수 수렴에 그치지만, 트레페텐-바이데만-슈멜처(2006)가 정리한 대로 적분로를 왼쪽으로 열린 포물선·쌍곡선·탤벗 곡선으로 변형하면 구적점 수 $N$ 에 대해 기하급수적으로 수렴해 $N \approx 20$–$30$ 이면 기계 정밀도에 닿는다.⁴

6. 현장에서 만나는 $f$ 들

• $e^{At}$, $\varphi_k(hA)$ — 선형 시스템의 정확한 시간 전파. 지수 적분기와 제어계 이산화의 심장이다.
• $A^{1/2}$, $A^{-1/2}$ — 공분산 $\Sigma$ 에서 표본을 뽑으려면 $\Sigma = BB^{H}$ 인 인수가 필요한데, 대칭 제곱근 $\Sigma^{1/2}$ 는 촐레스키 분해와 달리 유일하고 대칭이다. 앙상블 칼만 필터의 제곱근 필터(ETKF), 가우시안 프로세스 표본 생성, 백색화 $\Sigma^{-1/2}x$, 밀도범함수이론의 뢰딘 직교화 $S^{-1/2}$ 가 전부 여기다. 희소 $\Sigma$ 에서는 촐레스키의 충전(fill-in)을 피해 $\Sigma^{1/2}\xi$ 를 크릴로프로 직접 계산한다.
• $\operatorname{sign}(A)$ — $\tfrac12(I - \operatorname{sign}(A))$ 는 좌반평면 고유값의 스펙트럼 사영자다. 안정 불변 부분공간을 고유분해 없이 뽑아내므로 리카티 방정식과 최적 제어의 해법이 되고, 격자 QCD 오버랩 페르미온에서는 $\operatorname{sign}(H)b$ 가 계산량의 대부분이다.
• $\log\det A$ — 가우시안 우도와 최대우도추정이 요구한다. $\operatorname{tr}\log A$ 로 바꿔 확률적 대각합 추정(허친슨)과 크릴로프를 엮는 것이 대형 문제의 표준.

7. 관련 문서

• 행렬 지수함수 · 실베스터 공식 · 케일리-해밀턴 정리 · 최소다항식
• 슈어 분해 · 조르당 표준형 · 비정규 행렬 · 의사스펙트럼
• 파데 근사 · 보간과 근사 · 직교다항식
• 크리로프 부분공간법 · 아놀디 알고리즘 · 란초스 알고리즘 · 전처리기
• 조건수 · 후진 오차 해석 · 지수 적분기 · 리카티 방정식 · 수반행렬

8. Footnotes

1. 실제 버그의 단골이다. NumPy에서 np.exp(A) 는 성분별 지수, scipy.linalg.expm(A) 가 행렬 지수함수다. 둘은 $A$ 가 대각행렬일 때만 일치하는데, 하필 테스트를 대각행렬로 짜는 바람에 몇 달 뒤에 터지는 시나리오가 이 바닥의 국룰. ↩

2. 존재 조건은 영고유값에 붙은 조르당 블록 크기들이 서로 짝지어질 수 있는지가 결정한다. $J_2(0)$ 하나만 덩그러니 있으면 짝이 없어 실패한다. “제곱근이 없는 행렬”이 존재한다는 사실 자체가 스칼라 직관이 처음 배신당하는 지점이다. ↩

3. 성분별 오해의 반대편에는 ”$f(A)$ 만 구하면 뭐든 된다”는 과신이 있다. $f(A)$ 는 대개 조밀행렬이라 $n = 10^5$ 만 돼도 저장에 80GB다. 필요한 건 거의 언제나 $f(A)b$ 이고, 그걸 깨닫는 순간 알고리즘 선택지가 통째로 바뀐다. ↩

4. 곡선을 휘면 왜 빨라지는가는 사다리꼴 구적의 오차 이론이 답한다. 사다리꼴 공식은 피적분함수가 해석적으로 확장되는 띠의 폭에 비례하는 지수 속도로 수렴하는데, 적분로를 왼쪽으로 열면 그 띠를 넓게 확보할 수 있다. “경로를 바꿨더니 수렴이 대수에서 지수로 올라갔다”는, 복소해석이 수치해석에 주는 가장 짭짤한 선물 중 하나. ↩