쌍대각화

1. 개요[편집]

쌍대각화 Bidiagonalization
대상	임의의 $m\times n$ 행렬 ($m \ge n$)
결과	$U^{\top}AV = B$, $B$ 는 상이중대각 (특이값 보존)
표준 방법	하우스홀더 반사를 좌우 양쪽에서
비용	$4mn^2 - \tfrac{4}{3}n^3$ flops ($U,V$ 미조립)
$m \gg n$ 변형	QR 먼저 → $2mn^2 + 2n^3$ (Chan, $m \ge \tfrac{5}{3}n$)
반복법 대응물	골룹-칸-란초스 쌍대각화 (LSQR의 심장)
LAPACK	dgebrd → dbdsqr/dbdsdc

고유값 문제에서 상사변환이 하던 일을, 특이값 문제에서는 양쪽에서 서로 다른 직교행렬이 한다.

쌍대각화(bidiagonalization, 이중대각화)는 $m\times n$ 행렬 $A$ 에 좌우에서 서로 다른 직교(유니터리) 변환을 걸어 대각과 그 바로 위 한 줄만 남은 상이중대각행렬로 옮기는 축약이다.

$$U^{\top} A V = B = \begin{bmatrix} \alpha_1 & \beta_1 & & \\ & \alpha_2 & \ddots & \\ & & \ddots & \beta_{n-1} \\ & & & \alpha_n \\ & & & \\ & & \mathbf{0} & \end{bmatrix}$$

$U$ 와 $V$ 가 직교이므로 $A$ 와 $B$ 의 특이값이 완전히 같다. 특이값 분해를 계산하는 거의 모든 밀집 알고리즘이 이 형태를 1단계로 두는 이유이고, 삼중대각화가 대칭 고유값 문제의 입구인 것과 정확히 대응되는 자리를 차지한다.

두 축약의 차이는 딱 한 군데다. 고유값을 보존하려면 같은 행렬을 양쪽에 걸어야 하지만($Q^{\top}AQ$), 특이값을 보존할 때는 다른 행렬을 걸어도 된다($U^{\top}AV$). 자유도가 두 배라 한 칸 더 깎을 수 있을 것 같지만, 뒤에서 보듯 결국 부대각 한 줄에서 똑같이 멈춘다.

2. 좌우 반사를 번갈아 건다[편집]

절차는 하우스홀더 변환 반사를 왼쪽·오른쪽에서 교대로 적용하는 것뿐이다. $k$ 번째 단계에서

1. 왼쪽 반사 $U_k$ 로 $k$ 열의 대각 아래를 전부 0으로 만든다. 이건 QR 분해에서 하던 것과 완전히 같다.
2. 오른쪽 반사 $V_k$ 로 $k$ 행의 초대각(superdiagonal) 오른쪽을 전부 0으로 만든다. 여기서 한 칸을 남기는 것이 핵심 — $A(k,k{+}1)$ 까지 지우려 들면 오른쪽 곱이 방금 청소한 $k$ 열을 되살린다.

$n$ 번 반복하면 끝이다. 왼쪽 반사가 만든 0을 오른쪽 반사가 건드리지 않고, 그 반대도 성립하도록 작용 범위가 어긋나 있다는 것이 전부다. 총 비용은

$$4mn^2 - \tfrac{4}{3}n^3 \ \text{flops}$$

이고 $U, V$ 는 반사들의 곱으로만 존재한다. 실제로 조립하려면 완전 $U$($m\times m$)에 $4(m^2 n - mn^2 + n^3/3)$, $V$($n \times n$)에 약 $\tfrac{4}{3}n^3$ 이 더 든다. LAPACK dgebrd 는 관례대로 반사 벡터를 소거된 자리에 저장하고 스칼라만 따로 담아 추가 메모리를 쓰지 않는다.

2.1. 왜 여기서 멈추는가[편집]

“자유도가 두 개인데 왜 대각까지 못 가나?” 삼중대각화에서와 같은 두 겹의 답이 있다. 기계적으로는 위에서 말한 대로 왼손이 지운 자리를 오른손이 복구하기 때문이고, 원리적으로는 유한 번의 사칙연산·제곱근으로 특이값에 도달한다면 대칭행렬 $\begin{bmatrix} 0 & A \\ A^{\top} & 0\end{bmatrix}$ 의 고유값을 유한 번에 얻는 셈이라 아벨-루피니 정리에 걸린다. 특이값을 뽑는 뒷단계는 반드시 반복법이어야 한다.

2.2. $m \gg n$ 이면 순서를 바꾼다[편집]

$A$ 가 세로로 극단적으로 긴 경우 위 절차는 낭비다. 왼쪽 반사가 매번 $m$ 행 전체를 훑는데 정작 정보는 위쪽 $n\times n$ 에만 있기 때문이다. 그래서 먼저 QR 분해로 $A = QR$ 을 구하고 $n \times n$ 상삼각 $R$ 만 쌍대각화한다(R-쌍대각화, Lawson–Hanson–Chan). 비용은 $2mn^2 + 2n^3$ 이고, 두 식을 비교하면 $m \ge \tfrac{5}{3}n$ 부터 이쪽이 싸다. LAPACK dgesvd 가 내부에서 형상비를 보고 경로를 갈아타는 것이 이 지점이다.

3. 왜 $A^{\top}A$ 를 만들면 안 되는가[편집]

특이값이 $A^{\top}A$ 의 고유값의 제곱근이니 그냥 $A^{\top}A$ 를 만들어 삼중대각화하면 되지 않을까? 안 된다. $\kappa_2(A^{\top}A) = \kappa_2(A)^2$ 라서 조건수가 제곱으로 나빠지고, 작은 특이값은 그 순간 반올림 아래로 잠긴다. $\sigma_{\min}/\sigma_{\max} = 10^{-9}$ 인 행렬은 배정도에서도 곱을 만드는 순간 정보가 사라진다(정규방정식이 최소자승법에서 기피되는 것과 같은 사연이다).

재미있는 것은 쌍대각화가 그 곱을 명시적으로 만들지 않고 같은 일을 한다는 점이다. $B$ 가 상이중대각이면 $B^{\top}B$ 는 정확히 삼중대각이고, 그것이 $A^{\top}A$ 를 삼중대각화한 결과와 (부호를 빼면) 일치한다. 즉 쌍대각화는 $A^{\top}A$ 의 삼중대각화를 $A$ 위에서 제곱근을 취한 채 수행하는 절차다. 조건수를 제곱하지 않고 같은 정보를 얻는 것이 전부의 요지다.¹

4. 이중대각행렬에서 특이값 뽑기[편집]

$B$ 를 손에 쥐면 뒷단계는 삼중대각 고유값 문제와 같은 지형이 펼쳐진다.

• 암시적 영시프트 QR(dbdsqr). 골룹-라인시(1970)의 원형 알고리즘에 데멜-칸(1990)의 영시프트 변형을 얹은 것. 핵심은 작은 특이값을 상대 정확도로 계산한다는 것 — 이중대각 성분만 흔들리는 섭동에서는 모든 특이값이 상대 오차 $O(nu)$ 로 결정되며, 이 성질은 일반 행렬에는 없고 이중대각 구조에서만 나온다.
• 분할 정복(dbdsdc). $B$ 를 쪼개 재귀로 풀고 랭크-1 수정을 세큘러 방정식으로 병합한다. 특이벡터 전부가 필요한 큰 문제의 기본값이며, 자세한 것은 분할 정복 고유값 알고리즘 참고.
• dqds(Fernando–Parlett 1994, dlasq). 이중대각을 $LDL^{\top}$ 꼴로 유지한 채 시프트 있는 차분 상몫 변환을 돌린다. 뺄셈이 나오지 않도록 점화식을 짜서 상대 정확도를 유지하면서 시프트의 이득을 얻는, 이 바닥에서 손꼽히는 세공품이다.

전체 경로가 곧 dgesvd/dgesdd 다. 쌍대각화 $O(mn^2)$ 를 한 번 내고 나면 나머지가 $O(n^2)$ 급으로 떨어지는 구조라, 삼중대각화에서와 똑같이 준비운동이 본운동보다 오래 걸리는 역전이 일어난다. dgebrd 의 절반이 원리적으로 BLAS-2인 행렬-벡터 곱이라 실측 효율도 나쁘고, 그래서 현대 라이브러리는 여기서도 조밀 → 밴드 → 이중대각의 2단 축약으로 간다. 총 flops는 늘지만 시간은 줄어든다는 그 이야기가 반복된다.

이중대각 구조가 왜 그렇게 대접받는지는 한 문장으로 정리된다. $B$ 의 성분은 특이값을 상대 정확도로 결정한다. 성분에 $|\delta b_{ij}| \le \epsilon |b_{ij}|$ 규모의 상대 섭동만 가하면 모든 특이값이 $\sigma_i(1+O(n\epsilon))$ 안에서 움직인다 — 가장 작은 특이값도 자기 크기에 비례해서만 흔들린다는 뜻이다. 일반 행렬에는 이런 성질이 없다(절대 섭동 $\|E\|_2$ 만큼 움직이므로 작은 특이값은 통째로 사라질 수 있다). 그래서 이중대각까지 오는 것이 축약이 아니라 정보의 보존으로 취급되고, 뒷단계 알고리즘들이 뺄셈 하나까지 따져 가며 상대 정확도를 지키려 애쓰는 것이다.

5. 야코비 SVD — 정확도를 사는 다른 길[편집]

쌍대각화가 유일한 경로는 아니다. 축약을 아예 하지 않고 열 쌍을 골라 $2\times2$ 회전으로 직교화하는 것을 반복하는 단측 야코비(one-sided Jacobi, dgesvj)가 있다.

	쌍대각화 + QR/분할정복	단측 야코비
비용	$O(mn^2)$ 한 방	스윕 반복, 보통 몇 배 느림
작은 특이값	이중대각 성분 기준 상대 정확도	열 스케일링 기준 상대 정확도
병렬화	축약이 병목	회전을 짝지어 잘 흩어짐
실무 위치	기본값	정확도가 계약 조건일 때

차이가 나는 지점은 $A = DX$ 처럼 열마다 규모가 수십 자릿수씩 다른 행렬이다. 데멜-베셀리치(1992)는 이런 행렬에서 야코비가 $\kappa(X)$ 에만 의존하는 상대 오차를 내는 반면, 쌍대각화는 축약 단계에서 이미 규모 정보를 섞어 버려 그 보장이 깨진다는 것을 보였다. 유한요소 강성행렬이나 병렬 실험 데이터처럼 열 단위 규모 차이가 물리적으로 의미 있는 경우 이 차이가 실재한다. 물론 대부분의 문제에서는 dgesdd 가 압도적으로 빠르므로, 야코비는 “느려도 좋으니 작은 특이값을 믿을 수 있게”라는 요구가 붙을 때 꺼내는 카드다.

6. 란초스 쌍대각화 — LSQR의 심장[편집]

$A$ 가 거대한 희소행렬이면 $O(mn^2)$ 축약은 시작조차 못 한다. 반복법 대응물이 골룹-칸-란초스 쌍대각화다. $\beta_1 u_1 = b$, $\alpha_1 v_1 = A^{\top}u_1$ 로 시작해

$$\beta_{k+1} u_{k+1} = A v_k - \alpha_k u_k, \qquad \alpha_{k+1} v_{k+1} = A^{\top} u_{k+1} - \beta_{k+1} v_k$$

를 돌리면, 행렬-벡터 곱 $Av$ 와 $A^{\top}u$ 만으로 정규직교 기저 $\{u_k\}, \{v_k\}$ 와 하이중대각 $B_k$ 가 자라난다. $A$ 의 성분에는 손도 대지 않는다 — 란초스 알고리즘이 삼중대각을 만드는 것과 정확히 같은 구조이고, 실제로 이것은 $A^{\top}A$ 에 대한 란초스와 수학적으로 동치다.

여기서 LSQR(Paige–Saunders 1982)이 나온다. $\min_x \|Ax - b\|_2$ 를 $k$ 차 크릴로프 부분공간으로 제한하면 $\min_y \|\beta_1 e_1 - B_k y\|_2$ 라는 작은 이중대각 최소자승 문제가 되고, 이건 기븐스 회전으로 한 스텝에 갱신된다. 정규방정식에 켤레기울기법을 그냥 돌린 것(CGNR)과 정확 산술에서는 같은 반복열을 내지만, $A^{\top}A$ 를 만들지 않으므로 유한 정밀도에서 훨씬 낫다. 잔차 대신 $\|A^{\top}r\|$ 를 단조 감소시키는 변형이 LSMR(2011)이고, 정칙화 항 $\lambda\|x\|^2$ 를 넣는 것도 $B_k$ 에 행 몇 개를 덧붙이는 것으로 끝나 티호노프 정규화와 궁합이 좋다.

같은 재귀가 대규모 부분 SVD에도 쓰인다. 상위 몇 개 특이삼중항만 필요한 주성분 분석이나 잠재의미분석에서 svds 계열이 도는 것이 이 알고리즘이고, 암시적 재시작을 얹은 IRLBA(Baglama–Reichel 2005)가 대표 구현이다. 대가도 란초스와 동일하다 — 유한 정밀도에서 $u_k, v_k$ 의 직교성이 무너지고 유령 특이값이 튀어나온다. 완전 재직교화나 부분 재직교화(그람-슈미트의 그 재직교화)를 얹어야 하고, 그 비용이 반복당 $O(k(m+n))$ 로 누적된다. 밀집 쪽의 랜덤화 SVD가 “정확도를 조금 포기하고 BLAS-3 한 방”으로 가는 것과 대비되는 선택지다.²

7. 실무에서 밟는 지뢰[편집]

한 번 정리하고 가자. 아래가 밀집 SVD 한 번에 실제로 도는 것들이다.

단계	하는 일	LAPACK
(선택) QR 선처리	$m \ge \tfrac{5}{3}n$ 일 때 $R$ 만 남김	dgeqrf
쌍대각화	$U^{\top}AV = B$	dgebrd
반사 조립	$U, V$ 를 필요한 만큼만	dorgbr / dormbr
이중대각 SVD	$B = U_B \Sigma V_B^{\top}$	dbdsqr 또는 dbdsdc
드라이버	위를 묶은 것	dgesvd / dgesdd

밴드행렬이면 축약을 처음부터 밴드 구조에 맞춰 도는 dgbbrd 가 따로 있고, 열 규모 차이가 심하면 앞 절의 dgesvj 로 간다. 여기서 자주 터지는 것들이 있다.

• 부대각을 지우려는 유혹. 이중대각의 $\beta_i$ 를 반사로 더 지우려는 시도는 반드시 이미 만든 0을 되살린다. 코드가 무한 루프를 도는 대신 조용히 틀린 답을 뱉는 형태로 나타나서 더 나쁘다.
• 디플레이션 문턱값. $|\beta_i|$ 가 작으면 그 자리에서 $B$ 를 쪼개는데, 문턱을 절대값으로 잡으면 상대 정확도라는 이중대각 최대의 장점이 날아간다. dbdsqr 이 문턱을 상대 규모로 잡는 이유다.
• $U$ 를 다 만들지 마라. 최소자승법이나 저랭크 근사에는 $U$ 의 앞 $n$ 열이면 충분한데, 전체 $m \times m$ 을 조립하면 $m = 10^6$ 에서 메모리가 통째로 날아간다. dgesvd 의 jobu 인자가 존재하는 이유가 이것이고, 초심자가 가장 자주 밟는다.
• 란초스 쌍대각화의 중복 특이값. 재직교화 없이 돌리면 같은 특이값이 여러 번 수렴한 것처럼 보인다. 진짜 중복도인지 유령인지는 반복 기록을 봐야 알 수 있다.

8. 관련 문서[편집]

• 특이값 분해 · 삼중대각화 · 하우스홀더 변환 · QR 분해
• 그람-슈미트 · 기븐스 회전 · 란초스 알고리즘 · 크리로프 부분공간법
• 최소자승법 · 티호노프 정규화 · 주성분 분석 · 랜덤화 SVD
• 분할 정복 고유값 알고리즘 · 조건수 · 후진 오차 해석 · LAPACK

9. Footnotes[편집]

1. 골룹과 칸의 1965년 논문 제목이 “Calculating the singular values and pseudo-inverse of a matrix”인데, 요지가 정확히 ”$A^{\top}A$ 를 만들지 마라”였다. 5년 뒤 골룹-라인시가 실제 구현까지 붙여 Numerische Mathematik 에 알골 코드로 발표했고, 그게 지금 dgesvd 안에 그대로 살아 있다. 두 논문을 뭉뚱그려 “골룹-라인시 이중대각화”라고 부르는 문헌이 많은데, 축약 자체의 저작권은 골룹-칸 쪽이다. ↩

2. 축약은 후진 안정하다. 계산된 $\hat B$ 는 어떤 정확한 직교 $\tilde U, \tilde V$ 에 대해 $\hat B = \tilde U^{\top}(A + E)\tilde V$ 를 정확히 만족하고 $\|E\|_F \le c(m,n)u\|A\|_F$ 다. 바일 부등식을 얹으면 특이값이 최대 $\|E\|_2$ 만큼만 움직이므로, 큰 특이값은 걱정할 필요가 없다. 걱정거리는 언제나 $\sigma_{\max}$ 대비 극도로 작은 특이값이고, 그래서 데멜-칸의 상대 정확도 결과가 그토록 귀하게 취급된다. ↩