데이비드슨 알고리즘

1. 개요[편집]

데이비드슨 알고리즘 Davidson Algorithm
제안	Ernest R. Davidson, 1975
적용 행렬	대각 우세한 대칭·에르미트 희소행렬
구하는 것	최저(또는 특정) 고유쌍 몇 개
핵심 장치	대각 전처리 보정 벡터
주 무대	CI, CASSCF, TDDFT, 콘-샴 대각화
일반화	야코비-데이비드슨(JD)

란초스는 부분공간을 “자라는 대로” 받아들이지만, 데이비드슨은 다음에 어느 방향을 넣을지 골라서 넣는다.

데이비드슨 알고리즘(Davidson algorithm)은 대각 성분이 비대각 성분보다 뚜렷하게 우세한 대형 대칭(에르미트) 희소행렬에서 최저 몇 개의 고유쌍을 구하는 부분공간 반복 알고리즘으로, 잔차 벡터를 행렬의 대각 성분만으로 전처리해 다음 탐색 방향을 만든다는 점이 정체성이다. 1975년 어니스트 데이비드슨이 배치상호작용(CI) 행렬의 바닥상태를 겨냥해 발표했고, 이후 반세기 동안 전자구조 계산 코드의 사실상 기본 대각화기 자리를 지켰다.

배경에는 양자화학 행렬의 특수한 사정이 있다. CI 행렬의 대각 원소는 각 슬레이터 행렬식의 에너지라서 수십 하트리에 걸쳐 넓게 퍼져 있는 반면, 비대각 원소(2전자 적분)는 상대적으로 작다. 즉 $A \approx D$가 꽤 쓸 만한 근사다. 이 사정을 그대로 알고리즘에 박아 넣은 것이 데이비드슨이고, 그래서 범용 고유값 해법이 아니라 대각 우세 행렬 전용 특화 도구다.¹

2. 부분공간 전개와 레일리-리츠[편집]

알고리즘의 뼈대는 고유값 문제의 표준 도구인 레일리-리츠 사영이다. 정규직교 시행 벡터 $V_k = [\mathbf{b}_1 \cdots \mathbf{b}_k]$가 있을 때, 원래의 $N \times N$ 문제를 부분공간으로 눌러 작은 문제를 만든다.

$$\tilde{A} = V_k^{\mathsf{T}} A V_k, \qquad \tilde{A}\mathbf{y} = \theta \mathbf{y}$$

여기서 $\theta$가 리츠값(고유값 근사), $\mathbf{x} = V_k \mathbf{y}$가 리츠벡터다. 이 근사가 얼마나 틀렸는지는 잔차가 말해 준다.

$$\mathbf{r} = (A - \theta I)\mathbf{x}$$

$\|\mathbf{r}\|$이 문턱값 아래로 내려가면 종료, 아니면 $\mathbf{r}$에서 만든 새 방향을 부분공간에 추가하고 다시 돈다. 여기까지는 란초스 알고리즘이나 아놀디와 똑같다. 갈림길은 “$\mathbf{r}$로부터 무엇을 넣느냐”에서 생긴다.

3. 대각 전처리 — 싸구려 뉴턴 보정[편집]

이상적인 보정 벡터 $\boldsymbol{\delta}$는 $A(\mathbf{x}+\boldsymbol{\delta}) = \theta(\mathbf{x}+\boldsymbol{\delta})$를 만족해야 하니, 2차 항을 버리면

$$(A - \theta I)\boldsymbol{\delta} = -\mathbf{r}$$

이다. 즉 고유값 문제를 뉴턴-랩슨법으로 푸는 한 스텝에 해당한다. 그런데 이 방정식을 제대로 풀려면 결국 $N \times N$ 계를 풀어야 하니 배보다 배꼽이 크다. 데이비드슨의 처방은 무식할 정도로 간단하다 — $A$를 그 대각 $D = \mathrm{diag}(A)$로 바꿔치기한다.

$$\delta_i = \frac{r_i}{\theta - A_{ii}}$$

나눗셈 $N$번이 전부다. 이 $\boldsymbol{\delta}$를 기존 기저에 대해 그람-슈미트로 직교화·정규화해 $\mathbf{b}_{k+1}$로 삼는다.² 성능의 근원은 분모에 있다. 잔차 성분 중 대각 에너지가 현재 리츠값 $\theta$에 가까운 성분은 작은 분모로 나뉘어 크게 증폭되고, 멀리 떨어진 성분은 억눌린다. 다시 말해 “지금 상태와 강하게 섞일 만한 배치”를 골라 부분공간에 밀어 넣는다. 전처리기가 선형계에서 하는 일과 정확히 같은 짓을 고유값 문제에서 하는 셈이다.

그래서 수렴 속도가 대각 우세도에 통째로 의존한다. 대각 산포가 넓고 비대각이 작은 CI 행렬에서는 10~30회 반복으로 $10^{-8}$ 수렴이 흔하지만, 대각이 거의 상수인 행렬(예: 잘못 표현된 평면파 기저 해밀토니안)에서는 $(\theta I - D)^{-1}$이 거의 스칼라 배가 되어 보정이 잔차와 방향이 같아진다. 이때 데이비드슨은 그냥 크리로프 부분공간을 쌓는 것과 다를 바 없어지고, 오히려 재직교화 비용만 더 무는 란초스 알고리즘의 열화판이 된다.

4. 란초스와 무엇이 다른가[편집]

항목	란초스	데이비드슨
부분공간	크리로프 $\mathcal{K}_m(A,\mathbf{v})$ 고정	전처리로 매번 새로 고른 방향
직교화	3항 점화식(이론상 $O(N)$ 저장)	항상 전 기저에 대해 완전 직교화
필요 정보	$A\mathbf{v}$만	$A\mathbf{v}$ + 대각 성분 $A_{ii}$
수렴 지배 요인	스펙트럼 상대 간격	대각 우세도
유령 고유값	재직교화 안 하면 발생	구조상 완전 직교화라 없음

핵심 대비는 이렇다. 란초스는 부분공간이 알고리즘에 의해 결정되므로 저장을 극단적으로 아낄 수 있는 대신 어느 방향이 유용한지 고를 자유가 없다. 데이비드슨은 매 스텝 방향을 골라서 넣는 대신, 그 방향이 크리로프 구조를 갖지 않아 3항 점화식의 축복을 포기하고 $A\mathbf{b}_j$를 전부 저장·직교화해야 한다. 문제 구조가 좋은 힌트($D$)를 공짜로 준다면 후자의 거래가 압도적으로 유리하다는 것이 반세기의 결론이다.

5. 블록·재시작·올슨 보정[편집]

실전 구현은 몇 겹의 살이 더 붙는다.

• 블록(다중 근): 여러 개의 최저 고유쌍이 동시에 필요하면(들뜬 상태 계산이 그렇다) 리츠벡터 $p$개에 대한 잔차 $p$개로 보정 벡터를 한꺼번에 만든다. 축퇴 상태를 한 번에 잡을 수 있고 행렬-벡터 곱이 행렬-행렬 곱이 되어 캐시 효율도 오른다.
• 재시작(collapse): 부분공간 차원 $k$가 커지면 $kN$ 크기의 벡터 두 세트($\mathbf{b}_j$와 $A\mathbf{b}_j$)를 들고 있어야 해서 메모리가 먼저 죽는다. 보통 $k$가 상한(예: 원하는 근 수의 10배)에 닿으면 현재 리츠벡터들만 남기고 부분공간을 접은 뒤 다시 키운다. 여기에 올슨 보정(Olsen, 1990)이 더해진다. 순진한 대각 보정은 $\boldsymbol{\delta}$가 현재 리츠벡터 $\mathbf{x}$와 상당히 겹쳐 새 정보가 거의 없는 경우가 생긴다. 올슨은 그 겹침을 미리 빼도록

$$\boldsymbol{\delta} = -(D-\theta I)^{-1}\mathbf{r} + \epsilon\,(D-\theta I)^{-1}\mathbf{x}, \qquad \epsilon = \frac{\mathbf{x}^{\mathsf{T}}(D-\theta I)^{-1}\mathbf{r}}{\mathbf{x}^{\mathsf{T}}(D-\theta I)^{-1}\mathbf{x}}$$

로 고쳐, 보정이 $\mathbf{x}$에 대해 직교하도록 만들었다. 대형 CI 구현이 사실상 표준으로 채택한 형태이며, 다음 절의 야코비-데이비드슨으로 가는 징검다리이기도 하다.

6. 야코비-데이비드슨과 정확한 보정방정식[편집]

여기서 유명한 역설이 하나 나온다. 대각 근사가 나쁘다면 $(A-\theta I)$를 정확히 쓰면 되지 않나? 그런데 그러면 $\boldsymbol{\delta} = -(A-\theta I)^{-1}\mathbf{r} = -\mathbf{x}$가 되어, 새 방향이 이미 부분공간 안에 있는 벡터라 아무 진전이 없다.³ 전처리를 “더 정확하게” 만들었더니 알고리즘이 멈추는 것이다.

슬라이데레흐트·판데르포르스트(Sleijpen & van der Vorst, 1996)의 야코비-데이비드슨(JD)은 이 함정을 정면으로 피한다. 보정을 $\mathbf{x}$의 직교 여공간 안에서만 찾도록 사영을 양쪽에 걸어 준다.

$$(I - \mathbf{x}\mathbf{x}^{\mathsf{T}})(A - \theta I)(I - \mathbf{x}\mathbf{x}^{\mathsf{T}})\,\boldsymbol{\delta} = -\mathbf{r}, \qquad \boldsymbol{\delta} \perp \mathbf{x}$$

이 보정방정식을 정확히 풀면 레일리 몫 반복 수준의 3차 수렴이 나오고, 안쪽을 크리로프 부분공간법으로 대충만 풀어도 실용적으로 충분하다. 데이비드슨의 대각 근사는 JD 보정방정식을 가장 싸구려로 근사한 특수 케이스로 재해석된다. 안쪽 반복의 정확도와 바깥 반복 횟수를 맞바꾸는 이 구조 덕에, JD는 대각 우세가 없는 문제나 역반복법이 필요한 내부 고유값 문제까지 커버한다.

7. 실무에서의 자리[편집]

전자구조 코드를 열어보면 데이비드슨은 거의 항상 거기 있다. Quantum ESPRESSO의 diagonalization='david'는 기본값이고, VASP의 ALGO = Normal은 블록 데이비드슨이다. 밀도범함수이론 계산의 콘-샴 방정식은 자체무결장(SCF) 사이클마다 대각화를 반복하므로, 이전 사이클의 궤도를 시행 벡터로 재활용하면 초기 추측이 거의 정답이라 반복이 몇 번 만에 끝난다. 배치상호작용·CASSCF·TDDFT에서도 행렬을 절대 저장하지 않고 $\sigma$-벡터($A\mathbf{b}$) 생성 루틴만 직접 구현하는 것이 관례다.⁴ 요컨대 데이비드슨은 “행렬을 만들지 않는 것”과 “대각 성분은 공짜로 안다”는 두 조건이 동시에 성립하는 계산화학 생태계에 완벽히 맞춰진 도구다.

8. 관련 문서[편집]

• 고유값 문제 · 란초스 알고리즘
• 크리로프 부분공간법 · 전처리기
• 배치상호작용 · 하트리-폭 방법
• 밀도범함수이론 · 콘-샴 방정식
• VASP · Quantum ESPRESSO
• 희소행렬 · 뉴턴-랩슨법 · 준-뉴턴법
• 평면파 기저 · 완전기저 극한

9. Footnotes[편집]

1. 그래서 수치해석 전공자와 양자화학 전공자가 이 알고리즘을 두고 미묘하게 다른 온도를 보인다. 전자는 “대각 근사에 의존하는 휴리스틱”이라고 하고, 후자는 “우리 행렬에서 이거 없으면 계산이 안 돌아간다”고 한다. 둘 다 맞다. ↩

2. 부호 규약이 문헌마다 갈리는데, 어차피 새 벡터는 직교화 후 정규화되므로 $\boldsymbol{\delta}$의 전체 부호는 결과에 영향이 없다. 다만 코드 리뷰에서 이걸로 30분씩 싸우는 건 국룰. ↩

3. 실제로 데이비드슨 코드를 개선한답시고 전처리기를 완전 LU로 바꿨다가 “수렴이 멈췄다”며 버그를 찾아 헤매는 사례가 주기적으로 재발한다. 버그가 아니라 수학이 그렇게 생겼다. 이 현상이 JD 논문의 출발점이었다. ↩

4. 이 루틴을 흔히 시그마 벡터(sigma vector) 생성기라 부른다. Full CI 코드의 성능은 사실상 이 루틴 하나가 결정하며, 논문 한 편이 통째로 “$\sigma$-벡터를 어떻게 더 빨리 만드는가”에 바쳐지기도 한다. 고유값 알고리즘은 그 위에 얹힌 얇은 껍데기에 가깝다. ↩