LSQR

1. 개요

LSQR Least Squares QR
푸는 문제	$\min_x \lVert Ax-b\rVert_2^2 + \lambda^2\lVert x\rVert_2^2$ — 직사각·희소·랭크 결손 허용
필요한 것	$A\mathbf v$ 와 $A^{\top}\mathbf u$ 곱만
엔진	골룹-칸 쌍대각화 + 기븐스 회전 QR
저장	길이 $m$ 벡터 1개 + 길이 $n$ 벡터 3개
정확 산술 동치	정규방정식 위의 CG (CGLS)
단조 감소	$\lVert r_k\rVert$ (반면 $\lVert A^{\top}r_k\rVert$ 는 아님)
발표	페이지·손더스, ACM TOMS 1982

정규방정식을 만들지 않고 정규방정식 위의 켤레기울기법을 돌리는 법.

LSQR는 크고 희소한 행렬 $A \in \mathbb R^{m\times n}$ 에 대한 최소자승법 문제 $\min_x\lVert Ax - b\rVert_2$ (선택적으로 감쇠항 $\lambda^2\lVert x\rVert^2$ 포함)를 푸는 크리로프 부분공간법이다. 크리스 페이지와 마이클 손더스가 1982년 ACM TOMS 에 알고리즘과 포트란 코드(Algorithm 583)를 함께 발표했다¹. $A$ 의 성분은 쓰지 않고 곱 $A\mathbf v$ 와 $A^{\top}\mathbf u$ 만 있으면 되므로, 행렬을 명시적으로 만들기 어려운 문제(투영 연산자, 합성곱, 다른 솔버의 출력)에도 그대로 붙는다.

이름은 「Least Squares + QR」 이다. 골룹-칸 쌍대각화로 작은 이중대각 행렬을 키워 가고, 그 작은 최소자승 문제를 기븐스 회전 QR 로 한 스텝에 한 번씩 갱신한다. 쌍대각화 자체의 유도와 LSQR 이 거기서 나온다는 큰 그림은 쌍대각화 문서에 있다. 이 문서는 실제로 돌릴 때 필요한 것 — 재귀식, 공짜로 나오는 추정치, 정지 규칙, 감쇠, 전처리, 그리고 역문제에서의 반수렴 — 을 다룬다.

2. 알고리즘

초기화: $\beta_1\mathbf u_1 = \mathbf b$, $\alpha_1\mathbf v_1 = A^{\top}\mathbf u_1$, $\mathbf w_1 = \mathbf v_1$, $\mathbf x_0 = 0$, $\bar\phi_1 = \beta_1$, $\bar\rho_1 = \alpha_1$. 여기서 $\alpha$, $\beta$ 는 벡터를 단위 길이로 만드는 정규화 상수다.

$k = 1, 2, \dots$ 에 대해

1. 쌍대각화 한 스텝: $\beta_{k+1}\mathbf u_{k+1} = A\mathbf v_k - \alpha_k\mathbf u_k$, $\alpha_{k+1}\mathbf v_{k+1} = A^{\top}\mathbf u_{k+1} - \beta_{k+1}\mathbf v_k$.
2. 기븐스 회전: $\rho_k = \sqrt{\bar\rho_k^2 + \beta_{k+1}^2}$, $c_k = \bar\rho_k/\rho_k$, $s_k = \beta_{k+1}/\rho_k$, $\theta_{k+1} = s_k\alpha_{k+1}$, $\bar\rho_{k+1} = -c_k\alpha_{k+1}$, $\phi_k = c_k\bar\phi_k$, $\bar\phi_{k+1} = s_k\bar\phi_k$.
3. 해 갱신: $\mathbf x_k = \mathbf x_{k-1} + (\phi_k/\rho_k)\mathbf w_k$, $\mathbf w_{k+1} = \mathbf v_{k+1} - (\theta_{k+1}/\rho_k)\mathbf w_k$.

반복당 비용은 $A$ 곱 1회, $A^{\top}$ 곱 1회, 그리고 벡터 연산 $3m + 5n$ 남짓이다. 과거 기저 벡터를 저장하지 않으므로 메모리가 반복 횟수와 무관하다. GMRES 처럼 재시작할 필요가 없는 것은 $A^{\top}A$ 가 대칭이라 3항 재귀로 끝나기 때문이다 — MINRES·SYMMLQ 가 대칭 정사각 계에서 누리는 것과 같은 특권이다.

왜 CGLS 가 아니라 LSQR 인가. 정규방정식 $A^{\top}A\mathbf x = A^{\top}\mathbf b$ 에 켤레기울기법을 $A^{\top}A$ 를 만들지 않고 돌리는 CGLS(CGNR)도 같은 크릴로프 공간 $\mathcal K_k(A^{\top}A, A^{\top}\mathbf b)$ 에서 같은 문제를 푼다. 정확 산술에서는 두 방법이 같은 반복열을 낸다. 차이는 유한 정밀도다. 페이지·손더스는 $A$ 가 악조건일 때 LSQR 이 CGLS 보다 수치적으로 믿을 만하다고 보고했고, 이것이 이 알고리즘이 표준이 된 이유다. 다만 실무 문제 상당수에서 둘의 차이는 크지 않아, 「LSQR 이 CGLS 보다 몇 배 빠르다」 같은 기대는 하지 않는 것이 맞다 — 같은 반복열이니 반복 수는 비슷하고, 이득은 악조건에서의 안정성이다.

3. 공짜로 나오는 추정치와 정지 규칙

LSQR 의 실무적 매력은 해 말고도 추가 비용 없이 진단량이 나온다는 점이다.

• 잔차 노름: $\lVert\mathbf r_k\rVert = \bar\phi_{k+1}$ — 단조 감소한다.
• 정규방정식 잔차: $\lVert A^{\top}\mathbf r_k\rVert = \bar\phi_{k+1}\alpha_{k+1}\lvert c_k\rvert$ — 최소자승 해의 최적성 조건이 $A^{\top}\mathbf r = 0$ 이므로 이게 진짜 수렴 척도다. 단조가 아니다.
• $\lVert A\rVert_F$ 추정: 지금까지 나온 $\alpha_i$, $\beta_i$ 의 제곱합의 제곱근(하한이며 점점 커진다).
• $\mathrm{cond}(A)$ 추정과 $\lVert\mathbf x_k\rVert$.

이것들로 정지 규칙 세 개를 판정한다. 사용자가 atol, btol, conlim 을 준다.

규칙	조건	의미
S1	$\lVert r\rVert \le \texttt{btol}\,\lVert b\rVert + \texttt{atol}\,\lVert A\rVert\lVert x\rVert$	해가 있는 계($Ax = b$)를 데이터 오차 수준까지 풀었다
S2	$\lVert A^{\top}r\rVert \le \texttt{atol}\,\lVert A\rVert\lVert r\rVert$	최소자승 해를 $A$ 의 오차 수준까지 풀었다
S3	$\mathrm{cond}(A) \ge \texttt{conlim}$	추정 조건수가 한계를 넘었다 — 더 가면 해가 잡음을 증폭한다

atol·btol 은 「$A$ 와 $b$ 가 몇 자리까지 믿을 만한 데이터인가」로 정하라는 것이 원 논문의 권고다. 측정값 $b$ 가 세 자리 정확하면 btol = 1e-3 보다 엄격하게 풀어 봐야 잡음을 맞추는 일이다. 기본값 1e-6/1e-8 을 무심코 쓰기보다 데이터 정확도와 연결하는 게 V&V 정신에 맞다.

실측: $2000\times500$ 희소 행렬(밀도 1% 무작위 + 대각 0.1, 2-노름 조건수 6.6), 잡음 $10^{-3}$ 섞인 $b$. atol = btol = 1e-12 에서 SciPy LSQR 이 62번 만에 S2 로 정지했고, 밀집 SVD 최소자승 해와의 상대 차이가 $1.0\times10^{-14}$ 였다. LSMR 도 62번, $1.3\times10^{-14}$. 이때 LSQR 이 돌려준 조건수 추정은 92 였는데, 프로베니우스 노름 기반 추정이라 2-노름 조건수 6.6 보다 크게 나온 것이다. 절댓값이 아니라 추세로 읽어야 한다².

4. 감쇠와 전처리

감쇠(damping). $\lambda > 0$ 이면 LSQR 은

$$\min_x \left\lVert \begin{pmatrix} A \\ \lambda I \end{pmatrix}\mathbf x - \begin{pmatrix} \mathbf b \\ 0 \end{pmatrix}\right\rVert_2$$

를 푼다. 곧 티호노프 정규화다. 확장 행렬을 만들 필요 없이 반복마다 회전을 하나 더 넣는 것으로 처리되어 비용이 거의 늘지 않는다. 여러 $\lambda$ 를 시험하는 L-곡선 탐색에서 그래서 편하다. $\lVert L\mathbf x\rVert$ 같은 일반형 정규화는 $L$ 이 가역이면 변수 치환 $\mathbf y = L\mathbf x$ 로 표준형으로 바꿔 쓴다.

전처리. 최소자승에서는 오른쪽 전처리만 쓴다. $\min_y\lVert AM^{-1}\mathbf y - \mathbf b\rVert$ 를 풀고 $\mathbf x = M^{-1}\mathbf y$ 로 되돌리면 같은 문제다. 왼쪽에 곱하면 잔차 노름의 정의가 바뀌어 다른 최소자승 문제가 된다. 흔한 선택은 열 노름으로 나누는 대각 스케일링(값싸고 의외로 효과가 크다), 불완전 QR, 그리고 무작위 사영을 이용한 스케치 전처리(Blendenpik·LSRN 계열 — $SA$ 의 QR 에서 $R$ 을 얻어 $M = R$ 로 쓴다)다. 전처리기 문서의 일반론이 그대로 적용되지만, 목표가 $A^{\top}A$ 가 아니라 $AM^{-1}$ 의 특이값 분포를 모으는 것임을 기억할 것.

5. 역문제에서의 반수렴 — 반복 횟수가 정규화다

악조건 역문제에서 LSQR 은 독특하게 쓰인다. 크릴로프 기저가 큰 특이값 방향부터 차례로 잡기 때문에, 초반 반복은 신호를, 후반 반복은 잡음을 맞춘다. 그래서 오차가 줄다가 다시 커지는 반수렴(semi-convergence)이 나타나고, 반복 횟수 $k$ 자체가 정규화 매개변수 노릇을 한다.

실측: 폭 0.03 가우시안 블러 커널을 $200\times200$ 으로 이산화(조건수 $1.2\times10^{19}$ — 사실상 랭크 결손), 두 봉우리 참해, 데이터에 1% 상대 잡음. 감쇠 없이 반복 $k$ 에서 멈춘 해의 상대오차:

$k$	1	2	5	10	20	40	80	160	400
상대오차	0.195	0.090	0.027	0.014	0.028	0.128	0.812	7.39	39.5

10번째에서 1.4% 로 최저점을 찍고, 400번이면 참해의 40배 크기 잡음이 된다. 같은 데이터로 밀집 최소자승 해를 그냥 구하면 상대오차가 $6\times10^{9}$ 이다. 멈추는 시점이 곧 답의 품질이므로, 이런 문제에서는 위의 S1–S3 대신 잡음 수준을 아는 경우 불일치 원리($\lVert r_k\rVert \approx$ 잡음 노름에서 정지), 모르는 경우 L-곡선·일반화 교차검증으로 $k$ 를 고른다. 절단 특이값 분해가 「특이값 몇 개를 쓸까」를 고르는 것과 같은 결정을 SVD 없이 하는 셈이다³.

이것이 LSQR 이 지진 토모그래피·의료 영상 재구성의 표준 솔버가 된 이유다. 행렬은 광선 경로라 극도로 희소하고 수백만 열이며, 데이터에는 잡음이 있고, 정규화가 필수다. 대수적 재구성법(ART·SIRT)과 같은 자리를 놓고 경쟁하는데, 반복 몇십 번에 의미 있는 해가 나오는 쪽은 대개 LSQR 이다.

6. LSMR 과 다른 친척

• LSMR(퐁·손더스 2011): 같은 쌍대각화 위에서 $\lVert A^{\top}\mathbf r_k\rVert$ 를 최소화한다. 정확 산술에서 $A^{\top}A$ 위의 MINRES 와 동치이고, $\lVert A^{\top}\mathbf r_k\rVert$ 와 $\lVert\mathbf r_k\rVert$ 가 둘 다 단조 감소한다. LSQR 의 $\lVert A^{\top}\mathbf r_k\rVert$ 는 들쭉날쭉해서 S2 판정이 늦게 걸리는 경우가 있는데, LSMR 은 일찍 멈춰도 안전하다. 그래서 「일찍 끊을 것 같으면 LSMR」이 요즘 권고다.
• CRAIG: 해가 있는 계에서 최소 노름 해 $\min\lVert\mathbf x\rVert$ s.t. $A\mathbf x = \mathbf b$ 를 푼다. 정확 산술에서 $AA^{\top}$ 위의 CG(CGNE)와 동치.
• LSQR 을 써서는 안 되는 경우: $A$ 가 정사각 대칭이면 MINRES·켤레기울기법이 $A$ 곱 하나로 같은 일을 한다(LSQR 은 곱 둘에 조건수 제곱 공간을 탐색한다). 정사각 비대칭 가역 계에서도 GMRES·BiCGSTAB 가 대개 빠르다. LSQR 의 자리는 직사각·랭크 결손·정규화가 필요한 문제다.
• 밀집 문제라면 반복법을 쓸 이유가 없다. 수천 열 이하면 QR 분해·특이값 분해 한 방이 더 빠르고 정확하다.

7. 관련 문서

• 쌍대각화 · 기븐스 회전 · 크리로프 부분공간법 · 란초스 알고리즘
• 최소자승법 · 켤레기울기법 · MINRES · SYMMLQ · GMRES · 무어-펜로즈 유사역행렬
• 역문제 · 티호노프 정규화 · 절단 특이값 분해 · 지진 토모그래피 · 대수적 재구성법
• 전처리기 · 희소행렬 · LSMR · L-곡선 · 불일치 원리

8. Footnotes

1. 손더스는 MINRES·SYMMLQ(페이지와 공저, 1975)와 LSMR·MINRES-QLP 까지 크릴로프 최소자승 계열을 공저자들과 함께 반세기에 걸쳐 정리해 왔다. 스탠퍼드 SOL 페이지에 포트란·C·MATLAB·파이썬·줄리아 구현이 전부 올라와 있다. ↩

2. 「조건수 92 래, 괜찮네」와 「조건수 6.6 이래, 괜찮네」는 결론이 같아서 다행이지, $10^{8}$ 근처에서는 한 자리 차이가 「믿을 만함」과 「믿으면 안 됨」을 가른다. 추정치는 추정치다. ↩

3. 반복을 더 돌렸더니 답이 나빠지는 경험은 처음에 꽤 킹받는다. 수렴은 신에게 맡긴다는 국룰이 통하지 않는 몇 안 되는 동네다 — 여기선 「언제 멈출지」가 사람 몫이다. ↩