행렬 부호 함수

1. 개요[편집]

행렬 부호 함수 Matrix sign function
정의	$\operatorname{sign}(A) = A\,(A^{2})^{-1/2}$
존재 조건	$A$ 의 고유값이 허수축 위에 없을 것
기본 성질	$S^{2} = I$, $SA = AS$, 고유값은 $\pm 1$
표준 반복	뉴턴 $X_{k+1} = \tfrac12(X_k + X_k^{-1})$, 2차 수렴
도입	Roberts (1971 보고서 / 1980 출판)
주 용도	스펙트럼 사영자 · 리카티 방정식 · 스펙트럼 분할 정복

실수에서 $\operatorname{sign}(x)$ 는 정보가 1비트뿐인 시시한 함수다. 행렬로 올리면 그 1비트가 불변 부분공간 하나를 통째로 뽑아낸다.

행렬 부호 함수(matrix sign function)는 정사각행렬 $A$ 의 고유값을 실수부의 부호에 따라 $+1$ 과 $-1$ 로 보내는 행렬함수다. $A = ZJZ^{-1}$ 에서 조르당 블록을 실수부 부호로 묶어 $J = \operatorname{diag}(J_-, J_+)$ 라 쓰면

$$\operatorname{sign}(A) \;=\; Z\begin{bmatrix} -I_p & \\ & I_q\end{bmatrix}Z^{-1} \;=\; A\,(A^{2})^{-1/2}$$

이다. 고유값이 허수축 위에 있으면 정의되지 않는다 — 스칼라 $\operatorname{sign}$ 이 0에서 정의되지 않는 것과 같은 사정이다.¹

이 함수가 수치해석에서 대접받는 이유는 딱 하나, 고유분해를 하지 않고 불변 부분공간을 분리해 주기 때문이다. $\tfrac12(I - S)$ 는 좌반평면 고유값에 대응하는 불변 부분공간 위로의 사영자다. 안정/불안정 모드를 갈라야 하는 제어 문제, 스펙트럼을 반씩 잘라 병렬화해야 하는 대형 고유값 문제, 페르미온 행렬의 $\operatorname{sign}(H)b$ 를 요구하는 격자 QCD가 전부 이 한 줄에 걸려 있다.

2. 기본 성질[편집]

$S = \operatorname{sign}(A)$ 라 두면 행렬함수의 일반론에서 다음이 즉시 나온다.

• $S^{2} = I$ — 대합(involution)이다. 따라서 $S$ 는 대각화 가능하고 고유값이 $\pm1$ 뿐이며 $S^{-1} = S$.
• $SA = AS$, 그리고 $S$ 는 $A$ 의 모든 불변 부분공간을 보존한다.
• $A$ 가 실행렬이면 $S$ 도 실행렬이다(켤레쌍이 같은 값을 받으므로).
• $\operatorname{sign}(A^{-1}) = \operatorname{sign}(A)$, $\operatorname{sign}(cA) = \operatorname{sign}(c)\operatorname{sign}(A)$ ($c$ 는 0이 아닌 실수).
• $A$ 가 에르미트면 $S$ 도 에르미트이고 유니터리다. 비에르미트면 $\lVert S\rVert$ 에 상한이 없다 — 이 사실이 뒤에서 조건수 이야기로 돌아온다.

사영자 표현이 실무의 전부다.

$$P_{-} = \tfrac12\left(I - S\right), \qquad P_{+} = \tfrac12\left(I + S\right)$$

$P_-$ 의 상은 $\operatorname{Re}\lambda < 0$ 인 고유값들의 불변 부분공간이고 $P_+$ 는 그 반대다. 직교사영자가 아니라 빗각(oblique) 사영자라는 점은 기억해 둘 것 — $A$ 가 비정규 행렬이면 $\lVert P_\pm\rVert$ 가 1보다 훨씬 클 수 있고, 그 크기가 그대로 오차 증폭 인자가 된다.

3. 뉴턴 반복[편집]

$S^2 = I$ 는 곧 $S$ 가 $X^2 - I = 0$ 의 해라는 뜻이다. 여기에 뉴턴-랩슨법을 그대로 적용하면

$$X_{0} = A,\qquad X_{k+1} = \tfrac12\left(X_k + X_k^{-1}\right)$$

이 나오고, 이것이 사인 함수 계산의 표준 도구다. 모든 반복이 $A$ 의 행렬함수이므로 서로 교환하고, 스펙트럼 위에서 스칼라 반복 $x_{k+1} = \tfrac12(x_k + 1/x_k)$ 로 완전히 분리된다. 즉 행렬 문제의 수렴이 스칼라 문제의 수렴으로 환원된다.

오차 항등식도 닫힌 형이다.

$$X_{k+1} - S \;=\; \tfrac12\,X_k^{-1}\left(X_k - S\right)^{2}$$

$S^2 = I$ 와 교환성만 쓰면 두 줄로 나온다. 2차 수렴이 여기서 바로 읽힌다.

전역 거동은 케일리 변환이 깔끔하게 설명한다. $y = \dfrac{x-1}{x+1}$ 로 두면 뉴턴 사상이 정확히

$$y_{k+1} = y_k^{2}$$

가 된다.² 허수축이 단위원 $\lvert y\rvert = 1$ 로 옮겨지므로, 허수축을 벗어난 모든 초기값이 수렴하고 수렴 속도는 $\lvert y_0\rvert$ 가 1에서 얼마나 떨어져 있느냐로 결정된다. 고유값이 허수축에 붙어 있으면 $\lvert y_0\rvert \approx 1$ 이라 제곱을 아무리 해도 움직이지 않는다 — 이 알고리즘의 어려움이 어디서 오는지가 한 눈에 보인다.

반복 자체는 안정하다. 고정점 $S$ 에서 반복 사상의 프레셰 미분이 $L(E) = \tfrac12(E - SES)$ 인데, 계산해 보면 $L \circ L = L$ 로 멱등이라 오차가 증폭되지 않고 그 자리에 머문다. 뒤에 나올 행렬 제곱근의 뉴턴 반복이 불안정한 것과 정면으로 대비되는 지점이다.

4. 스케일링 — 초반 정체를 없애는 법[편집]

이론상 2차 수렴인데 실제로 돌려 보면 초반 10~20회가 아무 일도 하지 않는 것처럼 흘러가는 경우가 흔하다. 고유값 크기가 넓게 퍼져 있으면 큰 쪽은 $x \mapsto x/2$ 로, 작은 쪽은 $x \mapsto 1/(2x)$ 로 사실상 선형 수축만 하기 때문이다. 처방은 매 단계 스칼라 $\mu_k$ 를 곱해 스펙트럼을 1 근처로 밀어 넣는 것이다.

$$X_{k+1} = \tfrac12\left(\mu_k X_k + \mu_k^{-1}X_k^{-1}\right)$$

$\operatorname{sign}(\mu X) = \operatorname{sign}(X)$ 이므로 답은 바뀌지 않는다. 실무에서 쓰이는 세 가지:

스케일링	$\mu_k$	비고
행렬식 (Byers)	$\lvert\det X_k\rvert^{-1/n}$	고유값 절댓값의 기하평균을 1로. LU 가 이미 있어 사실상 공짜
노름	$\left(\lVert X_k^{-1}\rVert / \lVert X_k\rVert\right)^{1/2}$	스펙트럼 반지름의 상·하한을 균형
스펙트럼	$\left(\rho(X_k)\,\rho(X_k^{-1})\right)^{-1/2}$	최적이지만 $\rho$ 추정 비용이 붙음

행렬식 스케일링이 가장 널리 쓰인다. 역행렬을 구하느라 어차피 LU 를 만들었으니 $\det$ 은 대각 곱 한 번이고, 잘 조건화된 문제에서 반복 횟수를 수십 회에서 대여섯 회로 줄인다. 수렴이 가까워지면(예: $\lVert X_{k+1}-X_k\rVert$ 가 임계 이하) 스케일링을 꺼야 한다 — 2차 수렴 구간에서 $\mu_k$ 를 계속 곱하면 반올림 잡음만 얹힌다.

5. 역행렬을 피하는 변형[편집]

반복마다 $X_k^{-1}$ 을 구하는 것은 $O(n^3)$ 이고, GPU·분산 환경에서는 역행렬보다 곱셈이 훨씬 싸다. 그래서 곱셈만 쓰는 변형이 있다.

• 뉴턴–슐츠: $X_{k+1} = \tfrac12 X_k\left(3I - X_k^{2}\right)$. 역행렬이 없고 2차 수렴이지만 $\lVert I - A^{2}\rVert < 1$ 일 때만 수렴한다. 그래서 보통 스케일링 붙인 뉴턴으로 근처까지 간 뒤 갈아탄다.
• 핼리(파데 계열): $X_{k+1} = X_k\left(3I + X_k^{2}\right)\left(I + 3X_k^{2}\right)^{-1}$. 3차 수렴. 케니–라웁(1991)이 정리한 대로, $\operatorname{sign}$ 의 반복법 전체가 $\operatorname{sign}(x) = x/(1-(1-x^2))^{1/2}$ 에 대한 파데 근사 족으로 통일된다. 뉴턴과 뉴턴–슐츠는 그중 두 점일 뿐이다.
• QDWH 계열: 역행렬 대신 QR 분해 두 번으로 같은 유리반복을 구현한다. 나카츠카사–하이엄(2013)이 극분해와 부호 함수를 이 방식으로 묶어, 통신량이 적고 GPU 에서 잘 도는 스펙트럼 분할 정복을 만들었다. 역행렬을 QR 로 바꾸는 것이 현대 통신 회피 알고리즘의 상투적 수법이다.

6. 무엇에 쓰나[편집]

① 스펙트럼 사영과 불변 부분공간. $P_- = \tfrac12(I - S)$ 의 열공간에서 QR 분해로 정규직교기저 $U$ 를 뽑으면 $U^{H}AU$ 가 안정 부분의 축약이다. 고유값을 하나도 계산하지 않고 스펙트럼을 반으로 가른 셈. 이것을 재귀로 반복하면 스펙트럼 분할 정복이 되고, 원소별 통신이 적어 대규모 병렬 고유값 문제의 대안 경로가 된다(바이–데믈–구 1997).

② 대수 리카티 방정식. $A^{\top}X + XA - XGX + Q = 0$ 을 풀려면 해밀턴 행렬

$$H = \begin{bmatrix} A & -G \\ -Q & -A^{\top}\end{bmatrix}$$

의 안정 불변 부분공간이 필요하다. $W = \operatorname{sign}(H)$ 를 구하면 그 부분공간이 $W + I$ 의 영공간이므로, 블록으로 쪼갠 최소제곱 문제 하나로 $X$ 가 나온다. 로버츠(1971)의 원 논문이 정확히 이 용도였고, 리카티 방정식과 최적 제어에서 슈어 벡터 기반 방법과 나란히 쓰인다.

③ 실베스터·리아푸노프 방정식. $A$ 와 $B$ 의 스펙트럼이 각각 우·좌반평면에 있으면

$$\operatorname{sign}\!\begin{bmatrix} A & C \\ 0 & -B\end{bmatrix} = \begin{bmatrix} I & Y \\ 0 & -I\end{bmatrix}, \qquad AY + YB = 2C$$

가 성립한다. 블록 삼각행렬 하나에 뉴턴 반복을 돌리면 실베스터 방정식의 해가 우변 위쪽에 떨어진다는 뜻이고, 리아푸노프 방정식은 그 특수한 경우다.

④ 격자 QCD. 오버랩 페르미온의 디랙 연산자에 $\operatorname{sign}(H)$ 가 들어 있는데, $H$ 가 $10^{7}$ 차원 희소행렬이라 $S$ 를 만들 수가 없다. 필요한 건 $\operatorname{sign}(H)b$ 뿐이므로 란초스 알고리즘이나 졸로타레프 유리근사로 처리하며, 시뮬레이션 전체 계산량의 대부분을 여기가 먹는다.

7. 조건수와 한계[편집]

부호 함수는 허수축 근처에서 태생적으로 위험하다. 스칼라 $\operatorname{sign}$ 이 원점에서 불연속이라는 사실이 그대로 올라온 것이라, 고유값 하나가 허수축을 가로지르는 순간 답이 $\pm1$ 사이를 점프한다. 정성적으로 조건수는 스펙트럼과 허수축 사이의 거리에 반비례하며, 거리가 0이면 함수 자체가 정의되지 않는다.

더 고약한 것은 비정규 행렬이다. 고유값이 허수축에서 충분히 떨어져 있어도 $\lVert S\rVert$ 가 $10^{8}$ 씩 나오는 일이 있고, 이 경우 실제로 문제의 난이도를 재는 척도는 고유값이 아니라 의사스펙트럼이 허수축에서 얼마나 떨어져 있느냐다. 뉴턴 반복이 아무리 안정해도 문제 자체가 나쁘면 후진 오차 해석이 허용하는 전진 오차가 $\operatorname{cond}\cdot u$ 이므로 답은 그만큼 나쁘게 나온다. “알고리즘이 이상하다”고 말하기 전에 조건수를 재는 것이 여기서도 첫 번째 할 일이다.³

마지막으로 비용. 뉴턴 한 번이 역행렬 하나($\approx 2n^3$ flops)이므로 69회면 $12n^3$$18n^3$ 수준이고, 슈어 분해로 같은 일을 하는 것($\approx 25n^3$)보다 반드시 싸다고 말하기 어렵다. 부호 함수의 진짜 강점은 flops 가 아니라 구조다 — 역행렬과 행렬곱만으로 이루어져 있어 블록화·병렬화·GPU 이식이 쉽고, 슈어 분해는 그렇지 않다.

8. 관련 문서[편집]

• 행렬함수 · 행렬 제곱근 · 행렬 지수함수 · 파데 근사
• 고유값 문제 · 슈어 분해 · 비정규 행렬 · 의사스펙트럼
• 리카티 방정식 · 리아푸노프 방정식 · 실베스터 방정식 · 최적 제어
• 극분해 · QR 분해 · 통신 회피 알고리즘 · 분할 정복 고유값 알고리즘
• 조건수 · 후진 오차 해석 · 격자 QCD · 뉴턴-랩슨법

9. Footnotes[편집]

1. 그래서 ”$\operatorname{sign}(A)$ 를 구해 달라”는 요청에 대해 프로그램이 할 수 있는 가장 정직한 응답은 종종 “당신 행렬은 허수축 위에 고유값이 있습니다”이다. 제어 쪽에서는 이게 버그가 아니라 시스템이 안정도 불안정도 아닌 임계 상태라는 물리적 진단이라, 오히려 알고리즘이 문제를 대신 찾아 준 셈이 된다. ↩

2. 확인은 대입 한 번이다. $x' = \frac{x^2+1}{2x}$ 를 넣으면 $\frac{x'-1}{x'+1} = \frac{(x-1)^2}{(x+1)^2}$ 이다. 뉴턴 반복이 “복소평면에서 제곱하기”와 켤레라는 사실은, 같은 반복을 실수축이 아닌 곳에서 돌렸을 때 나타나는 뉴턴 프랙탈이 왜 그렇게 생겼는지도 동시에 설명한다. ↩

3. 조건수가 나쁘다는 말은 “정확한 답을 못 준다”가 아니라 “입력의 마지막 자리를 흔들면 답이 그만큼 흔들린다”는 뜻이다. 부호 함수의 경우 입력을 흔들었더니 안정 부분공간의 차원이 바뀌었다면, 그건 계산의 실패가 아니라 모형이 그 정도로 아슬아슬하다는 보고다. 이 구분을 못 하면 엉뚱한 곳에서 정밀도를 올리며 시간을 태우게 된다. ↩