행렬 균형화

1. 개요

행렬 균형화 Matrix balancing
연산	대각 닮음 변환 $A \to D^{-1}AD$ (+ 치환)
목표	각 $i$ 에서 $i$ 번째 행과 열의 노름을 (거의) 같게
불변	고유값 전부 (수학적으로)
줄이는 것	$\lVert A\rVert$ → 반올림 오차의 절대 크기
원조	오스본 (1960) · 파를렛-라인슈 (1969)
구현	LAPACK xGEBAL / xGEBAK, MATLAB balance

고유값은 하나도 안 바꾸면서, 고유값 솔버가 받는 오차만 바꾼다.

행렬 균형화(matrix balancing)는 고유값 계산 전에 대각 행렬 $D$ 로 닮음 변환 $D^{-1}AD$ 를 걸어, 각 인덱스 $i$ 에 대해 $i$ 번째 행의 노름과 $i$ 번째 열의 노름이 비슷해지도록 만드는 전처리다. 닮음 변환이므로 고유값은 수학적으로 하나도 안 바뀌고, 고유벡터는 $D$ 를 곱해 되돌리면 된다.

왜 이런 걸 하나? QR 알고리즘은 후진 안정이라 계산된 고유값이 「$\lVert E\rVert \approx u\lVert A\rVert$ 인 근처 행렬 $A + E$ 의 정확한 고유값」이다. 즉 오차의 절대 크기가 $\lVert A\rVert$ 에 비례한다. 원소 크기가 $10^{-14}$ 부터 $10^{14}$ 까지 퍼진 행렬은 $\lVert A\rVert$ 가 거대해서, 크기 1 짜리 고유값이 그 거대한 노름 몫의 오차를 뒤집어쓴다. 균형화는 같은 고유값을 가진 행렬 중 노름이 작은 쪽을 골라 QR 에 넘긴다. 단위가 섞인 물리 모델(구조물의 변위와 회전, 회로의 전압과 전류, 화학 반응망의 농도)에서 이런 행렬이 흔하다.

LAPACK 의 비대칭 고유값 드라이버 xGEEV 는 기본으로 균형화를 하고, NumPy·SciPy·MATLAB 의 eig 가 모두 그 위에 있다. 사용자가 의식하지 못하는 사이 거의 모든 비대칭 고유값 계산에 끼어 있는 단계다. 동반 행렬로 다항식 근을 구할 때 이것이 결정적이라는 이야기는 그 문서에 있고, 여기서는 알고리즘·효과·역효과를 다룬다.

2. 알고리즘 — 오스본 반복과 2의 거듭제곱

치환 단계. 먼저 행과 열을 치환해, 대각 원소 외에는 0 인 행이나 열을 찾아 맨 아래·맨 위로 보낸다. 그런 행·열의 대각 원소는 그대로 고유값이므로 떼어 내고, 가운데 남은 블록만 다룬다. 비용은 거의 없고 정확도에 손해도 없다(순수한 치환이라 반올림이 없다).

스케일 단계. 남은 블록에서 인덱스 $i$ 를 하나씩 돌며 비대각 열합 $c_i = \sum_{j\ne i}\lvert a_{ji}\rvert$ 와 행합 $r_i = \sum_{j\ne i}\lvert a_{ij}\rvert$ 를 잰다. $D$ 의 $i$ 번째 성분에 $f$ 를 곱하면 열합은 $f$ 배, 행합은 $1/f$ 배가 되므로 $f \approx \sqrt{r_i/c_i}$ 로 고르면 둘이 맞는다. 이것을 모든 $i$ 에 대해 돌고, 더 바뀌는 게 없을 때까지 반복한다. 오스본(1960)이 2-노름으로 제안했고, 그 반복이 행·열 노름이 모두 같은 균형 행렬로 수렴하며 그것이 대각 닮음 중 프로베니우스 노름을 최소화하는 행렬이라는 것도 보였다.

**파를렛·라인슈(1969)**의 개량이 오늘날 표준이 됐다.

• $f$ 를 기수(radix)의 거듭제곱으로 제한한다. 이진 부동소수점이면 $f = 2^m$. 2의 거듭제곱 곱셈은 지수부만 바꾸므로 반올림 오차가 전혀 생기지 않는다(언더플로·오버플로 범위를 넘지 않는 한).¹ 균형화 단계 자체는 정확하다는 뜻이고, 그래서 $f$ 를 정확히 $\sqrt{r/c}$ 로 맞추지 않고 「2 배 범위 안」이면 만족한다.
• 1-노름을 쓴다. 계산이 싸고, 2의 거듭제곱 제한 때문에 어차피 거친 균형이라 노름 선택이 결과에 큰 차이를 안 만든다.
• 멈춤 조건. 한 번의 스케일이 $(c_i + r_i)$ 를 최소 5% 줄이지 못하면((c+r)/f < 0.95 s 판정) 그 인덱스는 건드리지 않는다. 몇 번의 스윕이면 끝난다.

고유벡터는 균형화된 행렬 $A' = D^{-1}P^TAPD$ 의 고유벡터 $v'$ 에서 $v = PDv'$ 로 되돌린다(LAPACK xGEBAK).

3. 효과 — 같은 고유값, 다른 정확도

$8\times8$ 무작위 행렬 $M$ 을 $D = \operatorname{diag}(1, 10^s, 10^{2s}, \dots, 10^{7s})$ 로 $A = D^{-1}MD$ 처럼 일부러 망가뜨렸다. 고유값은 $M$ 과 똑같다. 균형화를 하는 경로(LAPACK geev)와 치환만 하는 경로(gees, 슈어 분해)로 고유값을 구해 $M$ 의 고유값과의 최대 상대오차를 쟀다:

$s$	$\lVert A\rVert_F$	균형화 뒤 $\lVert A'\rVert_F$	오차 (균형화)	오차 (균형화 없음)
1	$9.6\times10^{6}$	7.9	$2.7\times10^{-15}$	$2.3\times10^{-13}$
2	$9.5\times10^{13}$	7.8	$5.6\times10^{-15}$	$2.2\times10^{-13}$
3	$9.5\times10^{20}$	7.5	$1.6\times10^{-15}$	$8.5\times10^{-10}$
4	$9.5\times10^{27}$	7.3	$5.8\times10^{-15}$	$1.5\times10^{-8}$

균형화는 노름을 원래 $M$ 수준($\lVert M\rVert_F = 7.3$)까지 되돌리고, 오차는 $s$ 와 무관하게 기계 정밀도에 머문다. 균형화 없는 쪽은 $s = 4$ 에서 7자리를 잃는다. 흥미로운 건 균형화가 찾은 $D$ 가 원래의 $10^{s}$ 거듭제곱 구조를 2의 거듭제곱으로 근사해 거의 되찾는다는 점이다 — 행렬 성분만 보고 「누가 단위를 바꿨는지」 알아낸 셈이다.

다만 이 예제는 균형화에 유리하게 만든 것이다. 비정규성 자체가 큰 행렬(비정규 행렬)은 대각 닮음으로 고칠 수 없고, 균형화는 노름을 줄일 뿐 행렬의 본질적인 비정규성을 없애지는 못한다. 균형화 뒤에도 남는 고유값 조건수가 정확도의 진짜 한계다.

4. 역효과 — 작은 원소를 키우면 생기는 일

균형화는 작은 원소가 「진짜 작다」는 정보를 존중하지 않는다. 반올림 수준의 원소($10^{-16}$ 크기)를 정밀한 정보로 취급해 그것과 균형을 맞추려 들면, 그 작은 원소를 거대하게 키우는 스케일이 나올 수 있다. 그러면 고유값은 괜찮아도 고유벡터가 망가진다.

MATLAB 문서가 eig(A, 'nobalance') 옵션을 설명하며 드는 행렬이 대표 예다($\varepsilon$ 은 기계 엡실론):

$$A = \begin{pmatrix} 3 & -2 & -0.9 & 2\varepsilon \\ -2 & 4 & 1 & -\varepsilon \\ -\varepsilon/4 & \varepsilon/2 & -1 & 0 \\ -0.5 & -0.5 & 0.1 & 1 \end{pmatrix}$$

고전적인 파를렛-라인슈 반복(EISPACK balanc 와 같은 1-노름·0.95 판정)을 그대로 구현해 돌리면 스케일 인자가 $(1,\ 1,\ 7.5\times10^{-9},\ 3.4\times10^{7})$ — 최대·최소 비가 약 $4.5\times10^{15}$ — 로 나온다. 균형화된 행렬로 고유쌍을 구해 되돌리면 정규화된 고유벡터의 잔차 $\max_i\lVert Av_i - \lambda_iv_i\rVert/\lVert A\rVert_2$ 가 $6.9\times10^{-2}$ 다. 균형화 없이 슈어 분해에서 직접 구하면 $2.0\times10^{-16}$ 이다. 고유값 넷(5.5616, 1.4384, 1, −1)은 양쪽 다 맞는데 고유벡터는 한쪽이 두 자리도 안 맞는다.

제임스·랑구·로워리(2014)는 이 현상을 분석해, 고전적 정지 조건이 「이미 거의 균형인 행렬」을 알아보지 못하고 반올림 수준 원소에 끌려 다닌다는 것을 보이고 정지 조건 변경을 제안했다. SciPy 에 딸린 현재 LAPACK 의 dgebal 에 같은 행렬을 넣어 보면 스케일 인자가 전부 1 — 즉 아예 건드리지 않는다. 그래서 오늘날 NumPy 의 eig 는 이 행렬에서 잔차 $6\times10^{-16}$ 을 낸다. 「균형화가 고유벡터를 망친다」는 말은 구현과 버전에 따라 참이기도 거짓이기도 하다.²

실무 요령:

• 행렬에 반올림 수준의 「사실상 0」이 섞여 있고 고유벡터가 필요하면, 균형화 없는 결과와 한 번 비교해 본다(MATLAB 'nobalance', LAPACK xGEEVX 의 BALANC='N').
• 고유값만 필요하면 균형화는 거의 항상 이득이다.
• 단위계가 섞인 모델이라면 물리적으로 의미 있는 무차원화를 먼저 하는 것이 어떤 자동 균형화보다 낫다.³

5. 다른 「균형화」들과의 구분

같은 이름을 쓰는 다른 문제가 둘 있어 혼동이 잦다.

• 행렬 스케일링 / 평형화(equilibration). 선형 연립방정식 $Ax = b$ 를 풀 때 $D_1AD_2$ 로 행·열을 따로 스케일하는 것(LAPACK xGEEQU). 닮음 변환이 아니므로 고유값을 바꾸지만 해 $x$ 는 되돌릴 수 있다. 조건수를 줄이는 전처리기의 가장 단순한 형태다.
• 이중 확률 스케일링(싱크혼 문제). 양의 행렬을 $D_1AD_2$ 로 행합·열합이 모두 1 인 행렬로 만드는 것. 최적수송의 싱크혼 알고리즘이 이것이다. 교대로 행·열을 맞춘다는 점에서 오스본 반복과 형제지간이지만, 대상(행합 vs 행·열 노름의 일치)과 변환(양쪽 독립 vs 닮음)이 다르다.

일반화 고유값 문제 $Ax = \lambda Bx$ 에는 왼쪽·오른쪽을 따로 스케일하는 워드(1981)의 펜슬 균형화가 있다(QZ 알고리즘 참조). 2차 고유값 문제에서는 그보다 계수 행렬 자체의 노름을 맞추는 스케일링이 더 중요하다는 것이 그 문서의 결론이다.

6. 관련 문서

• QR 알고리즘 · 동반 행렬 · QZ 알고리즘 · 2차 고유값 문제 · 고유값 문제
• 후진 오차 해석 · 고유값 조건수 · 비정규 행렬 · 조건수 · 전처리기
• LAPACK · 슈어 분해 · 최적수송 · 무차원수

7. Footnotes

1. 16진 부동소수점을 쓰던 IBM 360 시절에는 기수가 16 이라 $f$ 도 16의 거듭제곱이었다. EISPACK 코드에 RADIX 라는 변수가 굳이 따로 있는 이유다. 하드웨어가 바뀌어도 「지수부만 건드린다」는 원리는 그대로 살아남았다. ↩

2. 같은 코드가 머신마다 다른 고유벡터를 낸다고 「재현성 위기」를 외치기 전에 링크된 LAPACK 구현과 버전부터 확인하자. 이런 조용한 차이가 범인일 때가 생각보다 있다. ↩

3. 「일단 SI 로 다 넣고 돌려」가 국룰이지만, 강성 $10^{9}$ N/m 과 질량 $10^{-3}$ kg 이 한 행렬에 들어가는 순간 고유값 솔버는 단위계의 뒷수습을 하게 된다. 균형화는 그 뒷수습을 자동으로 해 주는 장치이고, 무차원화는 애초에 어지르지 않는 습관이다. ↩