LOBPCG

1. 개요[편집]

LOBPCG Locally Optimal Block Preconditioned Conjugate Gradient
제안	Andrew Knyazev (2001)
푸는 문제	대칭·에르미트 $A\mathbf{x} = \lambda B\mathbf{x}$의 극단 고유쌍 $m$개
부분공간	$3m$차원 $[\,X,\ W,\ P\,]$ 위의 레일리-리츠
필요한 것	$A\mathbf{v}$ 곱, $B\mathbf{v}$ 곱, 전처리기 $T \approx A^{-1}$
구현	BLOPEX, hypre, PETSc, SLEPc, SciPy, Anasazi

크릴로프 부분공간을 쌓지 않는다. 세 덩어리만 들고 다니면서 매번 그 안에서 최선을 고른다.

LOBPCG(Locally Optimal Block Preconditioned Conjugate Gradient)는 대칭(에르미트) 일반화 고유값 문제 $A\mathbf{x} = \lambda B\mathbf{x}$의 최소(또는 최대) 쪽 고유쌍 여러 개를, 행렬-벡터 곱과 전처리기만으로, 고정 크기의 $3m$차원 부분공간 위에서 반복적으로 구하는 블록 고유값 해법기다. 2001년 안드레이 크냐제프가 정리해 발표했다.¹

란초스 알고리즘이나 아놀디 알고리즘 같은 크리로프 부분공간법은 부분공간을 한 차원씩 계속 키우다가 메모리가 감당 못 하면 재시작한다. LOBPCG는 애초에 부분공간 크기를 고정한다. 매 반복마다 현재 근사벡터 블록 $X$, 전처리된 잔차 블록 $W$, 직전 이동방향 블록 $P$ 세 개만 유지하고 그 $3m$차원 공간 위에서 레일리 몫을 최소화한다. 재시작이라는 개념 자체가 없고, 전처리기를 아무 제약 없이 끼워 넣을 수 있으며, 연산이 통째로 블록(BLAS-3)이라 GPU 컴퓨팅과 궁합이 좋다. 밀도범함수이론 코드가 콘-샴 방정식의 최저 밴드 수백 개를 뽑을 때, 또는 대형 모드 해석에서 최저 고유진동수 다발이 필요할 때 자주 등장하는 이유가 이것이다.

2. 레일리 몫 최소화에서 출발[편집]

대칭 양정부호 $B$에 대해 레일리 몫을

$$\rho(\mathbf{x}) = \frac{\mathbf{x}^{\mathsf{T}} A \mathbf{x}}{\mathbf{x}^{\mathsf{T}} B \mathbf{x}}$$

로 두면 최소 고유값 $\lambda_1$은 $\rho$의 전역 최솟값이고, 기울기는

$$\nabla \rho(\mathbf{x}) = \frac{2}{\mathbf{x}^{\mathsf{T}} B \mathbf{x}}\left(A\mathbf{x} - \rho(\mathbf{x}) B\mathbf{x}\right)$$

즉 잔차 $\mathbf{r} = A\mathbf{x} - \rho(\mathbf{x})B\mathbf{x}$가 곧 기울기 방향이다. 그러면 고유값 문제가 $\rho$에 대한 무제약 최적화가 되고, 경사하강법을 그대로 얹을 수 있다. 다만 순수 최속강하는 조건수에 비례해 기어가므로 두 가지를 덧댄다.

• 전처리. $\mathbf{w} = T\mathbf{r}$ ($T \approx A^{-1}$). 최속강하의 스텝 방향을 미리 늘여 놓는 것으로, 선형계에서 전처리 CG가 하는 일과 같다.
• 직전 방향의 재사용. $\mathbf{x}_{k-1}$을 버리지 않고 부분공간에 함께 넣는다. 이게 CG의 켤레방향에 해당하는 가속이다.

여기서 LOBPCG가 보통의 켤레방향법과 갈라진다. 스텝 크기 $\alpha$, 관성 $\beta$를 공식으로 정해 $\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha \mathbf{w}_k + \beta \mathbf{p}_k$를 만드는 대신, $\mathrm{span}\{\mathbf{x}_k, \mathbf{w}_k, \mathbf{p}_k\}$ 전체에서 $\rho$를 최소화하는 원소를 레일리-리츠로 곧바로 찾는다. 3차원 부분공간 안에서는 이보다 더 잘할 수 없다는 뜻에서 “국소 최적(locally optimal)“이라는 이름이 붙었다. 블록으로 올리면 $X$가 $n \times m$, $W = T(AX - BX\Lambda)$, $P$도 $n \times m$이 되어 기저 $S = [\,X\ W\ P\,]$가 $n \times 3m$, 풀어야 할 축소 문제는 $3m \times 3m$짜리 작은 일반화 고유값 문제 $(S^{\mathsf{T}}AS)C = (S^{\mathsf{T}}BS)C\Theta$ 하나다.

3. 왜 $P$를 명시적으로 들고 다니는가[편집]

수학적으로는 $\mathrm{span}\{X_k, W_k, P_k\} = \mathrm{span}\{X_k, W_k, X_{k-1}\}$이다. 그러면 그냥 직전 반복 벡터 $X_{k-1}$을 넣으면 될 것 같은데, 실제로 그렇게 구현하면 알고리즘이 수렴 직전에 붕괴한다. 이유는 단순하다. 수렴할수록 $X_k \to X_{k-1}$이므로 두 블록이 거의 같은 공간을 가리키고, 기저 $[X_k\ W_k\ X_{k-1}]$이 급격히 선형종속에 가까워진다.

그래서 LOBPCG는 레일리-리츠 결과 $X_{k+1} = X_k C_X + W_k C_W + P_k C_P$에서 $X_k$ 성분을 뺀 나머지만 다음 방향으로 저장한다.

$$P_{k+1} = W_k C_W + P_k C_P$$

$P_{k+1}$은 “이번 반복에서 실제로 이동한 양”이고, 정의상 $X_k$가 이미 차지한 방향을 포함하지 않는다. 같은 공간을 훨씬 나은 조건수의 기저로 표현하는 셈이다. 이 한 줄이 LOBPCG의 실질적 핵심이며, 초보 구현이 “이론대로 짰는데 왜 안 되지”에 빠지는 지점 1순위다.²

4. 수치 불안정과 직교화[편집]

$P$를 제대로 잡아도 불안정성이 완전히 사라지지는 않는다. 근본 원인은 수렴 그 자체다. 잔차가 작아지면 $W \to 0$이고, 수렴한 고유벡터에 대해서는 $P$도 함께 작아진다. 그러면 $S = [X\ W\ P]$가 랭크 결손에 가까워지고, 그람 행렬 $S^{\mathsf{T}}BS$가 수치적으로 양정부호를 잃어 촐레스키 분해가 실패하거나, 성공하더라도 축소 고유값이 쓰레기가 된다. 블록 크기가 클수록, 그리고 고유값이 뭉쳐 있을수록 더 빨리 온다. 실전 처방은 다음과 같다.

• 블록별 $B$-정규직교화. $X$, $W$, $P$ 각각을 $B$-내적에서 정규직교화한 뒤 그람 행렬을 만든다. 그러면 $S^{\mathsf{T}}BS$의 대각 블록이 단위행렬이 되어 조건수가 크게 개선된다.
• 재직교화와 조건수 감시. 촐레스키-QR을 두 번 돌리거나, 특이값 분해 기반 정규직교화로 갈아탄다. 그람 행렬의 조건수가 임계값을 넘으면 그때만 비싼 쪽으로 전환하는 적응형이 표준이다.
• 소프트 락킹(soft locking). 수렴한 열은 $W$·$P$ 갱신에서 빼고 $X$에만 남겨 다른 벡터들의 직교화 기준으로 쓴다. 완전히 제거(hard locking)하면 부분공간이 줄어 나머지 수렴이 느려질 수 있다.
• $P$ 재시작. 그래도 그람 행렬이 무너지면 $P = 0$으로 초기화한다. 그 반복만 전처리 최속강하로 퇴화하지만 곧 회복한다. 우아하진 않아도 실무 코드에는 거의 항상 들어 있는 안전장치다.

2010년대 후반의 “robust LOBPCG” 구현들은 이 처리들을 체계화한 결과물이며, 초기 BLOPEX가 큰 블록에서 자주 죽던 문제가 여기서 상당 부분 해소됐다.

5. 전처리기가 성능을 지배한다[편집]

LOBPCG의 수렴 속도는 사실상 $T$가 결정한다. 단일 벡터 전처리 최속강하에 대해 크냐제프-나이마이어가 준 예리한 추정은, $\lVert I - TA \rVert_A \le \gamma < 1$일 때

$$\frac{\rho(\mathbf{x}_{k+1}) - \lambda_1}{\lambda_2 - \rho(\mathbf{x}_{k+1})} \ \le\ \sigma^2\,\frac{\rho(\mathbf{x}_k) - \lambda_1}{\lambda_2 - \rho(\mathbf{x}_k)}, \qquad \sigma = 1 - (1-\gamma)\frac{\lambda_2 - \lambda_1}{\lambda_2}$$

꼴이다. LOBPCG의 부분공간은 최속강하의 부분공간 $\mathrm{span}\{\mathbf{x}_k, \mathbf{w}_k\}$를 항상 포함하므로, 레일리-리츠의 최소성에 의해 LOBPCG의 레일리 몫은 최속강하의 것보다 결코 크지 않다. 즉 위 추정이 LOBPCG에도 그대로 상계로 적용되고, 실측 수렴은 대개 그보다 훨씬 빠르다.

읽는 법은 두 가지다. 첫째, $\gamma \to 0$(완벽한 전처리)이면 $\sigma \to \lambda_1/\lambda_2$로 스펙트럼 갭이 그대로 수렴률이 된다. 둘째, 갭 $\lambda_2 - \lambda_1$이 작으면 전처리를 아무리 잘해도 느리다 — 다만 블록으로 묶으면 실효 갭이 $\lambda_{m+1} - \lambda_m$으로 바뀌므로, 뭉친 고유값을 뚫는 정공법은 블록 크기를 원하는 개수보다 넉넉히 잡는 것이다. 실무에서 쓰는 $T$는 문제마다 다르다. 평면파 밀도범함수이론에서는 운동에너지 항 기반의 대각 전처리, 유한요소 라플라시안 계열에서는 대수 다중격자법 한 사이클, 전자기 곡선요소에서는 보조공간 전처리기가 표준이다. 전처리기가 곧 성능이라는 점에서 LOBPCG는 고유값 해법기라기보다 고유값 문제용 전처리 반복법에 가깝다.³

6. 크릴로프 계열과 무엇이 다른가[편집]

항목	란초스 / 크릴로프-슈어	LOBPCG
부분공간	계속 성장, 주기적 재시작 필요	$3m$ 고정, 재시작 개념 없음
전처리	시프트-역변환(선형계 정확 해) 필요	임의의 대칭 양정부호 $T$를 직접 투입
블록화	가능하지만 구현이 번거로움	태생부터 블록, BLAS-3
메모리	반복 횟수에 비례	$3m$개 벡터로 상수
타깃	양 끝 + 시프트로 내부 고유값	주로 극단(최소/최대) 쪽
좋은 전처리기가 없을 때	여전히 잘 작동	최속강하 수준으로 퇴화

정리하면 **“좋은 전처리기가 있으면 LOBPCG, 없으면 크릴로프”**가 가장 정직한 요약이다. 란초스는 전처리기를 자연스럽게 못 먹는다 — $T$를 끼우면 크릴로프 구조가 깨지고, 우회하려면 시프트-역변환으로 선형계를 정확히 풀어야 하는데 그건 희소행렬 인수분해 비용을 내는 것이다. 반대로 전처리기가 없거나 나쁘면 LOBPCG는 $\gamma \approx 1$이 되어 볼 것이 없어진다. 부분공간을 전처리기로 키우는 발상 자체는 데이비드슨 알고리즘과 같은 계보이며, 차이는 데이비드슨이 부분공간을 계속 키우는 반면 LOBPCG는 3블록으로 잘라 고정한다는 점이다.

7. 구현과 실무 메모[편집]

원 저자의 참조 구현이 BLOPEX(Block Locally Optimal Preconditioned Eigenvalue Xolvers)이고, 이것이 hypre와 PETSc에 통합되면서 대규모 병렬 코드로 퍼졌다. SLEPc에는 EPSLOBPCG로, SciPy에는 scipy.sparse.linalg.lobpcg로, Trilinos의 Anasazi에도 들어 있다. ARPACK이 크릴로프 쪽의 사실상 표준이듯, 전처리 블록 쪽의 기본 선택지가 LOBPCG다.

• 블록 크기는 원하는 개수보다 크게. 뭉친 고유값을 놓치지 않으려면 $m$을 목표 개수의 1.2~1.5배로 잡는 것이 흔한 관행이다.⁴
• 초기 $X$는 난수 + $B$-정규직교화. 특정 방향으로 치우친 초기값은 그 성분이 없는 고유벡터를 영영 못 찾게 만들 수 있다.
• 잔차 판정은 상대 기준으로. $\lVert A\mathbf{x} - \theta B\mathbf{x}\rVert / (\lVert A\mathbf{x}\rVert + \lvert\theta\rvert\lVert B\mathbf{x}\rVert)$ 같은 형태를 쓴다. 절대 잔차로 판정하면 스케일이 다른 문제에서 기준이 무의미해진다. 이 잔차가 실제로 무엇을 보장하는지는 후진 오차 해석 참고.
• 내부 고유값은 그냥 포기하거나 시프트-역변환. LOBPCG는 스펙트럼 양 끝을 위한 도구다. 내부 대역이 필요하면 크릴로프-슈어 방법에 시프트-역변환을 얹는 쪽이 낫다.

8. 관련 문서[편집]

• 고유값 문제 · 레일리 몫 · 거듭제곱법 · 역반복법
• 란초스 알고리즘 · 아놀디 알고리즘 · 크릴로프-슈어 방법 · 데이비드슨 알고리즘
• 전처리기 · 다중격자법 · 크리로프 부분공간법 · 조건수
• 촐레스키 분해 · 특이값 분해 · QR 분해 · 희소행렬
• 밀도범함수이론 · 콘-샴 방정식 · 밴드 구조 계산 · 평면파 기저
• 모드 해석 · 축소차수모델 · 경사하강법
• ARPACK · LAPACK · GPU 컴퓨팅 · 후진 오차 해석

9. Footnotes[편집]

1. 크냐제프의 원 논문 제목이 “Toward the optimal preconditioned eigensolver”인데, 제목에 “toward”가 붙은 논문치고 이렇게 오래 살아남은 경우가 드물다. 20년 넘게 지난 지금도 대규모 전처리 고유값 계산의 기본 선택지다. ↩

2. 논문 수식만 보고 $[X_k\ W_k\ X_{k-1}]$로 짜면 초반 수십 반복은 멀쩡히 돌다가 잔차가 $10^{-6}$쯤에서 갑자기 튀어오른다. “수렴할수록 망가지는 알고리즘”이라는 인상을 남기고 사람들을 란초스로 돌려보내는 원흉. 사실은 기저가 붕괴한 것뿐이다. ↩

3. 전처리기의 품질 $\gamma$는 실전에서 잴 수 없는 양이다. 그래서 “AMG 한 사이클이면 $\gamma \approx 0.2$쯤 되겠지” 같은 감으로 돌리고, 안 되면 사이클을 늘린다. 이론이 있는데도 결국 튜닝인 것은 전처리기 세계의 오랜 전통. ↩

4. 블록 크기를 키우면 반복 횟수는 줄지만 반복당 비용이 $m$에 선형, 축소 문제 풀이가 $O(m^3)$으로 늘어난다. 그래서 “$m$을 키우면 빨라진다”는 어느 지점부터 반드시 뒤집힌다. 그 지점은 문제·하드웨어마다 다르고, 그래서 결국 돌려 보는 수밖에 없다. ↩