바텔스-스튜어트 알고리즘

수치해석 소프트웨어 마지막 수정: 2026-08-10 04:47:52

1. 개요[편집]

바텔스-스튜어트 알고리즘 Bartels–Stewart algorithm
발표	R. H. Bartels & G. W. Stewart, 1972 (CACM, Algorithm 432)
푸는 문제	실베스터 방정식 $AX + XB = C$
전략	실 슈어 분해 → 삼각계 대입 → 되변환
비용	$O(m^3 + n^3)$ (크로네커화는 $O(n^6)$)
구현	LAPACK `dtrsyl`, MATLAB `sylvester`/`lyap`

삼각형이면 대입으로 풀린다. 그러니까 삼각형으로 만들어라. 수치 선형대수의 절반은 이 문장이다.

바텔스-스튜어트 알고리즘(Bartels–Stewart algorithm)은 실베스터 방정식 $AX + XB = C$ 를 $O(m^3+n^3)$ 에 푸는 표준 직접법이다. 1972년 리처드 바텔스와 G. W. 스튜어트가 Communications of the ACM에 Algorithm 432로 발표했고,¹ 반세기가 지난 지금도 LAPACK·MATLAB·SciPy에서 이 방정식을 풀면 결국 이 알고리즘이 돈다. 수치 선형대수에서 “발표된 그대로 살아남은” 드문 사례다.

발상은 한 줄이다. 계수행렬을 (준)삼각으로 바꿔 놓으면 미지 행렬을 열 단위로 대입해 나갈 수 있다. 여기에 필요한 삼각화가 마침 슈어 분해이고, 그것을 안정적으로 계산하는 QR 알고리즘이 1960년대에 이미 성숙해 있었기에 곧바로 실용 알고리즘이 되었다. 크로네커화가 $O(n^6)$ 을 요구하는 자리에서 $O(n^3)$ 으로 내려온 것이 이 논문의 전부이자 전부다.

2. 절차[편집]

$A$ 가 $m\times m$ , $B$ 가 $n\times n$ , $X$ 와 $C$ 가 $m\times n$ 이라 하자.

1단계 — 슈어 분해. 두 계수행렬을 각각 실 슈어 형으로 만든다.

A = U R U^{\top}, \qquad B = V S V^{\top}

$U$ , $V$ 는 직교, $R$ 과 $S$ 는 준상삼각(대각에 $1\times1$ 또는 $2\times2$ 블록). 하우스홀더 변환으로 헤센베르크 축약 후 프란시스 이중시프트 QR을 돌리는 표준 경로이며, 직교인자까지 누적하면 각각 $\approx 25n^3$ flops다.

2단계 — 변환. $AX+XB=C$ 의 좌우에 $U^{\top}$ , $V$ 를 곱하고 $Y = U^{\top}XV$ , $D = U^{\top}CV$ 로 두면

R\,Y + Y\,S = D

가 된다. 문제의 구조는 그대로인데 계수행렬만 삼각이 되었다.

3단계 — 열 단위 대입. $S$ 가 상삼각이므로 $k$ 번째 열을 보면 $YS$ 의 $k$ 열이 $\sum_{j \le k} s_{jk}\,y_j$ 이고, 따라서

\left(R + s_{kk} I\right) y_k = d_k - \sum_{j < k} s_{jk}\, y_j

가 된다. 열은 $k=1$ 부터 앞으로 훑고, 각 열 안에서는 $R + s_{kk}I$ 가 상삼각이므로 후진대입으로 성분을 푼다. 이미 구한 열들이 우변에 누적될 뿐, 새로운 미지수는 매번 한 열뿐이다.

4단계 — 되변환. $X = U Y V^{\top}$ .

$B$ 가 상삼각이 아니라 하삼각이 되도록 잡거나 $A^{\top}$ 을 슈어화하는 등 문헌마다 관례가 조금씩 다른데, 방향이 앞으로냐 뒤로냐만 바뀔 뿐 골격은 같다.

3. 실 슈어형의 2×2 블록 — 실무의 진짜 디테일[편집]

위 3단계 설명은 $R$ , $S$ 가 진짜 삼각일 때의 이야기다. 실행렬은 켤레복소 고유값을 가질 수 있고, 실수 산술만으로는 그것을 삼각화할 수 없어 $2\times2$ 대각 블록이 남는다. 그러면 대입이 성분 단위가 아니라 블록 단위가 된다.

$R$ 의 대각 블록이 $p\times p$ ( $p \in \{1,2\}$ ), $S$ 의 대각 블록이 $q \times q$ ( $q \in \{1,2\}$ )일 때, 그 교차점에서 풀어야 하는 것은 $Y$ 의 $p\times q$ 소블록에 대한 작은 실베스터 방정식이다. 이것을 크로네커화하면

\left(I_q \otimes R_{ii} + S_{jj}^{\top}\otimes I_p\right)\mathrm{vec}(Y_{ij}) = \text{(우변)}

즉 크기 $pq$ 의 소연립이고, $pq \in \{1, 2, 4\}$ 다. $1\times1$ 과 $2\times2$ 뿐 아니라 양쪽이 모두 복소 켤레쌍인 자리에서는 $4\times4$ 연립이 나온다. 크기가 4 이하로 고정이라 완전 피벗팅 가우스 소거로 안전하게 풀 수 있고, 전체 비용에는 영향이 없다. 여기가 바로 “이론적으로는 자명한데 구현하면 지저분한” 대표적 구간이고, 실제로 LAPACK dtrsyl 소스의 길이 대부분을 이 분기 처리가 차지한다.²

이 소연립이 특이에 가까워지는 것은 $\lambda_i + \mu_j \approx 0$ , 즉 유일해 조건이 아슬아슬한 경우다. LAPACK은 이때 오버플로를 막으려고 우변을 스케일 인자 scale로 줄여 $R Y + Y S = \sigma D$ 를 대신 풀고 그 $\sigma$ 를 반환한다. 반환값이 1이 아니면 “문제가 나쁘다”는 신호로 읽어야 한다.

4. 비용 — $O(n^3)$ 과 $O(n^6)$ 의 거리[편집]

$m = n$ 일 때 각 단계의 비용은 대략 이렇다.

단계	비용	비고
슈어 분해 2회	$\approx 50n^3$	지배적. 직교인자 누적 포함
$C \to D$ 변환	$\approx 4n^3$	밀집 행렬곱 2회
삼각 실베스터 풀이	$\approx n^3$	열 대입
$Y \to X$ 되변환	$\approx 4n^3$	밀집 행렬곱 2회

전부 $O(n^3)$ 이고 메모리는 $O(n^2)$ 다. 크로네커화의 $O(n^6)$ · $O(n^4)$ 와 비교하면 $n=1000$ 에서 연산량이 $10^9$ 배, 메모리가 $10^6$ 배 차이다. 초 단위로 끝날 일과 수십 년 걸릴 일의 차이라, “행렬방정식은 펼쳐서 풀지 않는다”는 원칙이 여기서 완전히 굳는다.

비용을 더 깎는 개량도 있다. 헤센베르크-슈어 방법(Golub–Nash–Van Loan, 1979)은 두 행렬 중 하나를 슈어형까지 몰지 않고 헤센베르크 형까지만 축약한다. 헤센베르크 축약은 하우스홀더 변환 $n-2$ 번으로 유한 단계에 끝나 QR 반복이 아예 필요 없으므로 훨씬 싸고, 대신 각 열의 풀이가 삼각계 대입이 아니라 헤센베르크 계 풀이가 된다. 슈어 분해 한 번 값을 통째로 아끼는 셈이라 큰 쪽 행렬에 적용할수록 이득이 크다.

5. 왜 이 순서여야 하는가[편집]

세 갈래를 나란히 놓으면 이 알고리즘의 설계 논리가 선명해진다.

접근	계수행렬을 만드는가	비용	메모리
크로네커화 + 소거	만든다 ( $n^2\times n^2$ )	$O(n^6)$	$O(n^4)$
완전 대각화 $A = V\Lambda V^{-1}$	안 만든다	$O(n^3)$	$O(n^2)$
슈어 (바텔스-스튜어트)	안 만든다	$O(n^3)$	$O(n^2)$

두 번째 줄이 흥미롭다. $A$ 와 $B$ 를 대각화해 버리면 방정식이 성분별로 $(\lambda_i + \mu_j)\,\tilde{x}_{ij} = \tilde{c}_{ij}$ 로 완전히 분해되어, 나눗셈 $mn$ 번이면 끝난다. 비용도 $O(n^3)$ 이다. 그런데 아무도 이렇게 하지 않는다. 이유는 두 가지다.

존재하지 않을 수 있다. 결함(defective) 행렬은 대각화가 안 된다. 조르당 표준형은 부동소수점에서 계산 불가능한 물건이다.
존재해도 위험하다. 비정규 행렬에서 고유벡터 행렬 $V$ 의 조건수는 $10^{10}$ 을 우습게 넘길 수 있고, $\tilde{C} = V^{-1}CW$ 와 되변환 과정에서 그 조건수가 오차에 곱해진다.

반면 슈어 분해는 모든 정사각행렬에 대해 존재하고, 변환행렬이 직교라 조건수가 정확히 1이다. 대각화의 “성분별 나눗셈”이라는 편의를 조금 포기하고(삼각이므로 대입이 필요하다) 대신 무조건적 존재성과 완벽한 조건수를 얻는 교환 — 이것이 바텔스-스튜어트의 본질이다. 같은 교환이 슈어 분해를 쓰는 거의 모든 알고리즘의 밑바닥에 깔려 있다.

6. 리아푸노프 특수화와 해머링 변형[편집]

$B = A^{\top}$ 이면 리아푸노프 방정식이 되고, 슈어 분해를 한 번만 하면 되므로 비용이 절반 가까이 줄어든다. 해가 대칭이라는 사실도 이용해 상삼각 부분만 계산한다.

여기에 얹히는 것이 해머링(Hammarling) 변형(1982)이다. 안정한 $A$ 와 $Q = C^{\top}C$ 에 대해 $A^{\top}P + PA + C^{\top}C = 0$ 을 풀 때, $P$ 를 거치지 않고 촐레스키 분해 인수 $L$ ( $P = L^{\top}L$ )을 직접 계산한다. 이유는 정확도다. 그라미안처럼 특이값이 $10^{16}$ 에 걸쳐 퍼져 있는 해는 $P$ 를 먼저 만들면 반올림 때문에 음의 고유값을 갖고 나오고, 그러면 뒤이은 촐레스키가 실패하거나 평형 절단의 한켈 특이값이 NaN이 된다. 인수를 직접 전파하면 정부호성이 구조적으로 보장된다.³

발상 자체는 칼만 필터의 제곱근 필터와 완전히 같다 — 정부호 행렬을 다루는 코드는 가능하면 행렬이 아니라 그 인수를 들고 다닌다. 축소차수모델의 평형 절단이 한켈 특이값을 $\sqrt{\lambda_i(W_cW_o)}$ 가 아니라 두 촐레스키 인수의 곱 $L_o L_c^{\top}$ 의 특이값으로 계산하는 것도 같은 이유다.

이산시간판 $A^{\top}XA - X + Q = 0$ 에도 대응물이 있다. 슈어 형으로 바꾼 뒤 대입하는 골격은 같지만, 각 단계에서 나오는 소연립의 계수가 합 $r_{ii} + s_{jj}$ 가 아니라 곱 $r_{ii}s_{jj} - 1$ 꼴이라 세부가 다르다. 바로(Barraud, 1977)의 알고리즘이 표준이며 MATLAB dlyap이 이 경로다.

7. 수치적 성질과 구현[편집]

계산된 해 $\hat X$ 는 잔차 관점에서 만족스럽다. 대략

\|A\hat{X} + \hat{X}B - C\|_F \lesssim u\left(\|A\|_F + \|B\|_F\right)\|\hat{X}\|_F

수준의 한계가 성립한다( $u$ 는 단위 반올림). 다만 이것이 전진 오차가 작다는 뜻은 아니다. $\hat X$ 가 참 해에서 얼마나 벗어나는지는 $\mathrm{sep}(A,-B)$ 로 나눈 값에 지배되고, sep이 작으면 잔차가 아무리 작아도 해는 얼마든지 틀릴 수 있다. 알고리즘은 무죄, 문제가 유죄인 전형적 구도다. 그래서 라이브러리들이 scale이나 sep 추정치를 함께 반환한다.

실제 구현 지도는 이렇다.

LAPACK — dtrsyl(실수)·ztrsyl(복소)이 3단계(삼각 실베스터 풀이)만 담당한다. 슈어 분해는 dgees로 따로 부르고, 변환·되변환은 dgemm으로 사용자가 조립한다. 최근 버전에는 재귀 블록화된 dtrsyl3이 추가되어 레벨-3 BLAS 비중을 높였다.
재귀 블록화 — 욘손·코그스트룀의 RECSY가 대표적. 삼각 실베스터 풀이를 절반씩 쪼개 재귀하면 대부분의 작업이 행렬곱으로 바뀌어, 캐시 계층에서 몇 배가 빨라진다. 알고리즘의 flop 수는 그대로인데 실측 속도만 오르는 전형적 현대화.
SLICOT — 제어 쪽 표준 라이브러리. SB04MD(실베스터), SB03MD(리아푸노프), 해머링 변형까지 갖추고 있다.
MATLAB — sylvester(A,B,C)가 $AX+XB=C$ 를, lyap/dlyap(Control System Toolbox)이 리아푸노프계를 푼다. SciPy는 solve_sylvester·solve_lyapunov가 같은 경로다.

8. 관련 문서[편집]

1972년 논문의 제목은 Solution of the matrix equation AX + XB = C [F4]이다. 대괄호 안의 F4는 CACM 알고리즘 분류 코드다. 요즘 기준으로는 제목이 아니라 함수 시그니처에 가깝지만, 실제로 그 시절 Collected Algorithms of the ACM은 논문이라기보다 지금의 표준 라이브러리에 해당하는 물건이었다. ↩
” $1\times1$ 아니면 $2\times2$ 겠지”라고 얕보고 짜기 시작하면 $4\times4$ 분기에서 반드시 한 번 데인다. 게다가 $2\times2$ 블록의 경계를 잘못 판정해 켤레쌍을 쪼개면 그 자리에서 실수 산술이 붕괴한다. 이 방정식을 직접 구현해 보는 것은 좋은 연습이지만, 프로덕션에는 dtrsyl을 부르는 게 정답이다. ↩
반대로 말하면, $P$ 를 다 구해 놓고 나서 촐레스키를 돌리는 순간 이미 늦었다는 뜻이다. 부동소수점에서 한 번 잃은 정보는 사후에 복구되지 않는다. “결과가 이론상 정부호니까 괜찮겠지”는 수치해석에서 가장 위험한 문장 중 하나다. ↩