랜덤화 SVD

수치해석 통계 마지막 수정: 2026-08-12 04:14:26

1. 개요[편집]

랜덤화 SVD Randomized SVD
정리한 문헌	Halko–Martinsson–Tropp, SIAM Review 53(2), 2011
2단 구조	① 무작위 치역 근사 Y = AΩ → Q ② 사영 후 소형 SVD B = QᵀA
손잡이	과표본 p (보통 5~10) · 멱반복 q (보통 0~2)
기대오차	σ_k+1 의 상수배 — 최적값의 상수배 안
진짜 이점	통과 횟수 · BLAS-3 · 병렬성 (flop 수가 아니다)
구현	scikit-learn `randomized_svd`, MATLAB `svdsketch`, RSVDPACK

상위 20개 특이쌍이 필요한데 100만 차원 행렬을 통째로 분해하고 있다면, 문제는 알고리즘이 아니라 자세다.

랜덤화 SVD(randomized SVD, RSVD)는 행렬 $A$ 에 무작위 벡터 뭉치를 곱해서 그 치역(range)의 근사 정규직교기저 $Q$ 를 먼저 얻고, 원래 문제를 $Q$ 가 펼친 저차원 공간으로 눌러 내린 뒤 거기서 정직한 특이값 분해를 수행해 상위 $k$ 개 특이쌍을 구하는 알고리즘이다. 할코·마르틴손·트롭(Halko–Martinsson–Tropp)이 2011년 리뷰 논문 Finding Structure with Randomness 에서 흩어져 있던 결과를 하나의 레시피와 오차 이론으로 묶어 낸 뒤 사실상 표준이 되었다.¹

핵심 통찰은 허무할 정도로 단순하다. $A$ 가 수치적으로 랭크 $k$ 근처라면, 무작위 벡터 $\omega$ 에 대해 $A\omega$ 는 지배적인 좌특이벡터들이 펼친 공간 쪽으로 자동으로 쏠린다. 작은 특이값 방향은 곱해지면서 그만큼 작아져 사라지기 때문이다. 그러니 $k$ 개보다 조금 넉넉하게 $\ell = k + p$ 개만 던져도 치역은 거의 다 잡힌다. 나머지는 $\ell$ 차원짜리 작은 행렬을 다루는 평범한 선형대수다.

2. 알고리즘 — 두 단계가 전부다[편집]

$A \in \mathbb{R}^{m\times n}$ 에서 랭크 $k$ 근사를 원한다고 하자. $\ell = k+p$ 로 둔다.

1단계 (무작위 치역 근사). 원소가 i.i.d. 표준정규인 $\Omega \in \mathbb{R}^{n\times \ell}$ 을 뽑아

Y = A\Omega \in \mathbb{R}^{m\times \ell}, \qquad Y = QR \;\Rightarrow\; Q \in \mathbb{R}^{m\times \ell}

를 계산한다. $Q$ 는 $Y$ 의 열공간의 정규직교기저이고, 목표는 $\|A - QQ^{T}A\|$ 를 작게 만드는 것이다. QR은 안정성 때문에 하우스홀더 변환이나 재직교화를 넣은 그람-슈미트를 쓴다.²

2단계 (사영 후 소형 SVD). 이제 $\ell \times n$ 짜리 작은 행렬만 다룬다.

B = Q^{T}A, \qquad B = \tilde{U}\Sigma V^{T} \;\Rightarrow\; U = Q\tilde{U}, \qquad A \approx U\Sigma V^{T}

$\ell$ 개 중 앞의 $k$ 개만 남기면 랭크 $k$ 근사다. 비용은 $A\Omega$ 와 $Q^{T}A$ 두 번의 행렬곱 $O(mn\ell)$ (희소면 $O(\ell \cdot \mathrm{nnz})$ )에 소형 분해 $O((m+n)\ell^{2})$ 를 더한 것. 전체 SVD의 $O(mn\min(m,n))$ 과 비교하면 $\min(m,n)$ 이 $\ell$ 로 바뀐 셈이다.

$A$ 를 곱셈 연산자로만 접근하면 된다는 점이 중요하다. 원소를 명시적으로 저장하지 않고 고속 푸리에 변환이나 희소행렬 곱, 혹은 시뮬레이션 코드 한 번 돌리기로 정의된 $A$ 에도 그대로 쓴다.

3. 과표본추출 p — 왜 딱 k개로는 안 되나[편집]

$p = 0$ 으로 두면 $\ell = k$ 개의 무작위 벡터가 정확히 $k$ 차원 지배 부분공간을 잡아야 하는데, 이는 $\Omega$ 의 $k\times k$ 부분블록이 병적으로 조건이 나쁘지 않기를 요구한다. 무작위 행렬은 가끔 그렇게 되고, 그러면 근사가 통째로 망가진다. $p$ 는 그 사고에 대한 보험이다. 이론이 말하는 바는 상당히 구체적이다. HMT의 평균 오차 한계는

\mathbb{E}\,\|A - QQ^{T}A\|_{F} \le \left(1 + \frac{k}{p-1}\right)^{1/2}\left(\sum_{j>k}\sigma_j^{2}\right)^{1/2}

\mathbb{E}\,\|A - QQ^{T}A\|_{2} \le \left(1 + \sqrt{\frac{k}{p-1}}\right)\sigma_{k+1} + \frac{e\sqrt{k+p}}{p}\left(\sum_{j>k}\sigma_j^{2}\right)^{1/2}

이다. 오른쪽 괄호 안의 양은 에카르트-영-미르스키 정리가 말하는 최적값 그 자체이고, 앞에 붙은 것은 순전히 상수다. 즉 랜덤화 SVD는 최적 랭크 $k$ 근사보다 상수배 나쁠 뿐이며, $p=5$ 면 $\sqrt{1+k/4}$ , $p=10$ 이면 $\sqrt{1+k/9}$ 수준이다. 게다가 실패 확률은 $p$ 에 대해 대략 $3p^{-p}$ 로 떨어져서, $p=10$ 이면 $10^{-10}$ 급이다. “확률 알고리즘이라 운이 나쁘면 틀린다”는 걱정은 $p$ 를 5에서 10으로 올리는 것으로 실질적으로 종료된다.

4. 멱반복 q — 스펙트럼이 안 떨어질 때의 생명줄[편집]

위 한계식을 다시 보면 프로베니우스 항 $\left(\sum_{j>k}\sigma_j^{2}\right)^{1/2}$ 가 스펙트럼 노름 한계에도 끼어 있다. 특이값이 뚝 떨어지는 행렬이면 이 항이 $\sigma_{k+1}$ 과 비슷해 아무 문제가 없다. 문제는 감쇠가 느릴 때다. 꼬리에 $\sigma_{k+1}$ 급 특이값이 수천 개 있으면 그 제곱합의 제곱근은 $\sigma_{k+1}$ 의 수십 배가 되고, 한계식이 쓸모없어지는 정도가 아니라 실제로 근사가 나빠진다. 무작위 벡터가 신호와 잡음을 구분할 근거가 없기 때문이다.

처방은 감쇠를 인위적으로 만드는 것이다. $A$ 대신 $(AA^{T})^{q}A$ 에 무작위 벡터를 던진다.

Y = (AA^{T})^{q}A\,\Omega

이 행렬은 좌우 특이벡터가 $A$ 와 같고 특이값만 $\sigma_j^{2q+1}$ 로 바뀐다. 특이값 비율이 $(\sigma_{k+1}/\sigma_k)^{2q+1}$ 로 지수적으로 벌어지므로, $q=1$ 이나 $2$ 만 써도 스펙트럼이 평평한 행렬이 갑자기 가파른 행렬처럼 행동한다. 실제로 오차 한계의 상수 부분이 $1/(2q+1)$ 제곱으로 눌리면서 $\sigma_{k+1}$ 로 빠르게 수렴한다.

여기에 반드시 지켜야 할 구현 디테일이 하나 있다. $(AA^{T})^{q}A\Omega$ 를 곧이곧대로 곱셈만 반복하면, 곱할 때마다 작은 특이값 방향의 정보가 $\sigma^{2q+1}$ 로 줄어들어 $q=2$ 쯤에서 이미 배정밀도 바닥( $\varepsilon_{\mathrm{mach}} \approx 10^{-16}$ )에 잠긴다. 부분공간이 지배 방향 몇 개로 붕괴해 버리는 것이다. 그래서 실전 코드는 곱할 때마다 QR로 재직교화한다 — 즉 $Y \leftarrow \mathrm{qr}(A^{T}Y)$ , $Y \leftarrow \mathrm{qr}(AY)$ 를 번갈아 돈다. 이렇게 쓰면 이 알고리즘은 사실상 무작위 시작점을 가진 부분공간 반복법이고, 실제로 랜덤화 SVD의 절반은 오래된 직교 반복법을 확률적 시작점과 오차 이론으로 새로 포장한 것이다.³

5. 진짜 이점은 flop 수가 아니다[편집]

랜덤화 SVD를 “빠른 SVD”로만 소개하면 절반은 오해다. 밀집 행렬을 메모리에 올려놓고 상위 몇 개만 뽑는 상황이라면 란초스 알고리즘 계열(svds, PROPACK)이 flop 수로는 대개 더 적게 쓰고 정확도도 더 높다. 랜덤화가 이기는 지점은 다른 데 있다.

통과 횟수(pass). 행렬이 메모리에 안 들어가고 디스크·네트워크·스트림에서 흘러올 때, 비용을 지배하는 것은 산술이 아니라 데이터를 몇 번 훑느냐다. 기본형은 $A$ 를 2번( $Y=A\Omega$ , $B=Q^{T}A$ ), 멱반복을 쓰면 $2q+2$ 번 통과한다. 란초스는 반복 횟수만큼 통과한다.
한 번 통과(single-pass)와 스트리밍. $Y = A\Omega$ 와 $Z = A^{T}\Psi$ 를 동시에 누적해 두면 $A$ 를 딱 한 번만 보고도 근사를 만들 수 있다. 스케치는 갱신 가능하므로 $A$ 가 $A + \Delta A$ 로 조금씩 흘러들어오는 스트리밍 환경에서도 그대로 굴러간다. 정확도는 2패스 판본보다 떨어지지만, 행렬을 두 번 볼 방법이 애초에 없는 문제에서는 비교 대상이 없다.
병렬성과 BLAS-3. $A\Omega$ 는 그냥 밀집 행렬곱 한 방이라 캐시·GPU·분산 환경 어디서든 최고 효율로 돈다. 반면 란초스/아놀디는 매 반복이 앞 반복에 의존하고 반복마다 전역 동기화(내적)가 필요해 구조적으로 순차적이다. 랜덤화의 승부처는 산술 복잡도가 아니라 동기화 횟수다.

이 관점을 잡으면 변형들도 자연스럽게 읽힌다. $\Omega$ 를 가우시안 대신 SRFT/SRHT(부분표본 푸리에·아다마르 변환)로 잡으면 $A\Omega$ 를 $O(mn\log \ell)$ 에 계산할 수 있고, 희소 부호 행렬(CountSketch류)로 잡으면 희소 $A$ 에서 비용이 nnz에 비례한다. 대칭 준정부호 행렬에는 $Q$ 대신 니스트룀(Nyström) 근사가 같은 스케치에서 더 좋은 오차를 준다.

6. 랭크를 모를 때 — 적응적 치역 탐색[편집]

지금까지는 $k$ 를 안다고 가정했다. 실제로는 “오차가 $\varepsilon$ 이하가 되는 최소 랭크”를 원하는 경우가 더 많다. 이때 쓰는 것이 적응적 치역 탐색(adaptive randomized range finder)이다.

핵심은 무작위성을 오차 추정에도 재활용하는 것이다. 무작위 벡터 $\omega$ 를 새로 뽑아 잔차 $\|(I - QQ^{T})A\omega\|$ 를 재면, 그 값이 $\|(I-QQ^{T})A\|$ 의 확률적 추정치가 된다. 벡터를 $r$ 개 쓰면

\|(I-QQ^{T})A\| \le 10\sqrt{\tfrac{2}{\pi}}\ \max_{i=1,\dots,r}\big\|(I-QQ^{T})A\omega^{(i)}\big\|

가 확률 $1 - 10^{-r}$ 로 성립한다. 그래서 실전 알고리즘은 무작위 벡터를 블록 단위로 계속 추가하며 $Q$ 를 확장하다가, 최근 $r$ 개의 잔차가 전부 임계값 아래로 떨어지면 멈춘다. “랭크를 정하는 일”과 “오차를 재는 일”이 같은 연산으로 처리된다는 점이 이 방법의 우아한 구석이다. 반대로 란초스 알고리즘에서 이런 사후 오차 추정은 훨씬 성가시다.

주의할 점은 이 추정이 이미 뽑은 $\Omega$ 와 독립인 새 벡터로 이루어져야 한다는 것이다. $Y$ 를 만든 벡터를 재사용하면 잔차가 구조적으로 0에 가깝게 나와 오차를 심각하게 과소평가한다. 코드에서 난수 스트림을 아끼려다 이 함정에 빠지는 일이 실제로 있다.

7. 주의사항 — 이건 정밀 도구가 아니다[편집]

특이값은 항상 과소평가된다. 사영 $QQ^{T}$ 는 노름을 늘리지 않으므로 $\sigma_j(B) \le \sigma_j(A)$ 가 모든 $j$ 에서 성립한다. 즉 RSVD가 내놓는 특이값은 참값의 하한이다. “특이값 몇 자리까지 맞아야 한다”는 요구가 있으면 이 방법이 아니다.
꼬리 쪽 특이쌍은 못 믿는다. $\ell = k+p$ 개를 계산했어도 신뢰할 수 있는 것은 앞의 $k$ 개다. 뒤쪽 $p$ 개는 애초에 오차를 흡수하라고 넣은 여유분이다.
간격이 좁으면 벡터가 흔들린다. $\sigma_k \approx \sigma_{k+1}$ 이면 개별 특이벡터는 랜덤화와 무관하게 원래 불안정하다. 이건 알고리즘 탓이 아니라 문제의 조건수 탓이다.
결정론적 재현성. 난수 시드를 고정하지 않으면 같은 입력에서 매번 다른 출력이 나온다. 검증 및 확인 관점에서 시드는 반드시 로그에 남겨야 하는 입력값이다.

8. 어디서 쓰이나[편집]

주성분 분석과 절단 특이값 분해의 대규모 판본이 첫 번째 고객이다. scikit-learn의 TruncatedSVD(algorithm='randomized')가 기본값으로 이걸 쓰고 있어서, 사용자 상당수는 자기가 랜덤화 알고리즘을 돌리고 있다는 사실도 모른다. 두 번째는 적합직교분해와 축소차수모델이다. 스냅숏 행렬이 수십 GB인 대형 전산유체역학 해석에서 POD 모드를 뽑을 때, 스냅숏을 디스크에서 한 번 흘려보내며 스케치만 누적하는 방식이 자연스럽다. 그 밖에 커널 행렬 근사, 계층 행렬(H-행렬) 압축, 압축센싱 및 최소자승법의 스케치 기반 가속이 같은 도구상자에서 나온다.

이론적 뿌리는 존슨-린덴슈트라우스 보조정리 계열의 차원축소 결과지만, 둘을 같은 것으로 말하면 곤란하다. JL은 점들 사이의 거리를 보존하는 사영을 말하고, 랜덤화 SVD가 필요로 하는 것은 특정 부분공간 위에서 노름을 보존하는 부분공간 매장(subspace embedding)이다. 후자가 더 강한 요구이며, 가우시안·SRFT가 그 요구를 만족한다는 것이 별도로 증명되어야 했다.⁴

9. 관련 문서[편집]

아이디어 자체는 그보다 앞선다. Frieze–Kannan–Vempala(1998)의 표본추출 기반 저계수 근사, Sarlós(2006)의 스케치, Rokhlin–Szlam–Tygert(2009)의 무작위 PCA가 각각 조각을 들고 있었다. HMT의 공은 새 알고리즘을 발명한 것이라기보다 “무작위로 던져도 되는 이유”를 한 권 분량의 오차 이론으로 정리해 실무자가 안심하고 쓸 수 있게 만든 것이다. 실제로 논문 본문의 절반이 확률 부등식이다. ↩
여기서 그람-슈미트 문서의 교훈이 그대로 재생된다. $Y = A\Omega$ 의 열은 이미 지배 방향으로 쏠려 있어 서로 거의 평행하다. 즉 조건수가 나쁜 행렬을 직교화하는 상황이고, 고전 그람-슈미트를 그냥 쓰면 직교성이 $O(u\kappa^{2})$ 로 무너진다. “무작위라 잘 흩어져 있겠지”라는 직관은 정확히 반대다. ↩
그래서 나이 지긋한 수치해석자들이 “이거 그냥 직교 반복법 아니냐”고 말하는 것도 틀린 말은 아니다. 다만 무작위 시작점 덕에 확률적 오차 보증이 붙었고, 반복 횟수를 정확도가 아니라 스펙트럼 감쇠에 맞춰 고정할 수 있게 되었다는 게 실무적으로는 큰 차이다. 옛 알고리즘 + 새 이론 = 새 알고리즘, 은 수치해석에서 꽤 흔한 방정식이다. ↩
데이터 과학 발표에서 “JL 보조정리에 의해 랜덤 사영은 안전합니다”라는 문장이 만능 주문처럼 쓰이는데, 원 논문의 상수를 대입해 보면 요구되는 차원이 실제 쓰는 것보다 훨씬 크게 나오는 경우가 많다. 실전에서 잘 되는 이유는 JL 상수가 느슨해서지, JL이 정확히 그걸 보장해서가 아니다. ↩