야코비-데이비드슨 방법

편집 역사 토론
수치해석 계산물리 마지막 수정: 2026-08-06 05:03:18

1. 개요[편집]

야코비-데이비드슨 방법
Jacobi–Davidson method
발표Sleijpen & van der Vorst (1996)
계보Jacobi(1846) 직교보정 + Davidson(1975) 부분공간 가속
핵심사영된 보정 방정식을 부정확하게 풀어 부분공간을 확장
강점내부 고유값, 비대칭·일반화·다항 고유값 문제, 전처리 자유
구현SLEPc EPSJD, PRIMME, JDQR/JDQZ

시프트를 정답에 붙일수록 계는 특이해진다. 그러면 특이해지는 방향만 빼고 풀면 되지 않을까?

야코비-데이비드슨 방법(Jacobi–Davidson method)은 슬레이펀과 판 데르 포르스트가 1996년에 제안한 대형 고유값 문제 해법으로, 현재의 근사 고유쌍에 직교하는 방향으로만 보정을 구해 탐색 부분공간을 넓혀 가는 부분공간 가속 뉴턴법이다. 크리로프 부분공간법 계열이 부분공간의 모양을 스스로 정하는 것과 달리, 이쪽은 확장 벡터를 매 스텝 자유롭게 고를 수 있다는 점이 결정적인 차이다.

기본 골격은 레일리-리츠다. 정규직교 기저 VmV_m 이 주어지면 Vm∗AVmV_m^{*}AV_m 의 고유쌍에서 리츠 쌍 (θ,u)(\theta, u), u=Vmsu = V_m s 를 뽑고 잔차 r=Au−θur = Au - \theta u 를 잰다. θ=u∗Au\theta = u^{*}Au (레일리 몫)로 잡으면 r⊥ur \perp u 다. 그다음이 이 방법의 전부인 보정 방정식이다.

(I−uu∗)(A−θI)(I−uu∗) t=−r,t⊥u(I - uu^{*})(A - \theta I)(I - uu^{*})\, t = -r, \qquad t \perp u

여기서 얻은 tt 를 기저에 직교화해 붙여 Vm+1V_{m+1} 을 만들고 반복한다. 한 사이클을 정리하면 이렇다.

  1. VmV_m 에서 Hm=Vm∗AVmH_m = V_m^{*}AV_m 을 만들고 목표에 가장 가까운 리츠 쌍 (θ,u)(\theta, u) 를 뽑는다.
  2. r=Au−θur = Au - \theta u 를 계산하고, ∥r∥\|r\| 이 허용오차 이하면 수렴 처리 후 수축한다.
  3. 보정 방정식을 (부정확하게) 풀어 t⊥ut \perp u 를 얻는다.
  4. tt 를 VmV_m 에 대해 직교화(수치적으로는 반복 그람-슈미트)해 새 기저 벡터로 추가한다.
  5. 기저 크기가 mmax⁡m_{\max} 에 닿으면 좋은 리츠 벡터 mmin⁡m_{\min} 개만 남기고 재시작한다.

스텝당 비용은 행렬-벡터 곱 몇 번(내부 반복 횟수만큼)과 O(nm)O(nm) 규모의 직교화이고, 메모리는 기저 VmV_m 과 AVmAV_m 을 들고 있는 O(nm)O(nm) 이다.

2. 왜 사영자로 감싸는가[편집]

정답 λ\lambda 에 가까워질수록 A−θIA - \theta I 는 특이행렬에 접근한다. 이것이 역반복법이나 레일리 몫 반복이 매 스텝 마주하는 딜레마다 — 시프트를 정확히 잡을수록 수렴은 빨라지지만 풀어야 할 선형계는 나빠진다.1 직접법으로 LU를 뜨면 이 특이성이 오히려 원하는 고유벡터 방향으로 오차를 증폭시켜 주기 때문에 실질적으로 무해하지만, 대형 희소 문제에서 반복법으로 풀어야 하는 순간 이야기가 달라진다. 전처리기를 붙인 GMRES에게 “거의 특이한 계를 풀어라”라고 시키면 그대로 발산한다.

사영자 I−uu∗I - uu^{*} 로 양쪽을 감싸면 문제가 사라진다. 연산자를 uu 의 직교 여공간 u⊥u^{\perp} 로 제한하면, λ\lambda 가 단순 고유값인 한 θ→λ\theta \to \lambda 여도 그 제한은 비특이인 채로 남고 조건수는 스펙트럼 갭에 지배된다. 특이해지는 성분은 정확히 uu 방향 하나인데 그 방향은 이미 부분공간에 있으므로 잃을 정보도 없다. 그래서 보정 방정식을 정확히 풀면 이 방법은 레일리 몫 반복과 같은 수렴 차수를 보인다 — 일반 행렬에서 2차, 에르미트 행렬에서는 3차.

관점을 바꾸면 이건 제약조건 u∗u=1u^{*}u = 1 을 단 비선형계 Au−θu=0Au - \theta u = 0 에 대한 뉴턴-랩슨법 스텝이다. 실제로 u+tu + t 를 정규화한 것이 다음 근사가 되도록 tt 를 요구하면 (A−θI)t=−r+εu(A - \theta I)t = -r + \varepsilon u 가 나오고, 미지의 스칼라 ε\varepsilon 을 제거하기 위해 양변에 I−uu∗I - uu^{*} 를 곱하면 정확히 보정 방정식이 된다. 다만 야코비-데이비드슨은 뉴턴 스텝을 그대로 취하지 않고 누적된 부분공간 전체에서 다시 레일리-리츠를 한다 — 이 부분공간 가속이 순수 뉴턴보다 훨씬 넓은 수렴 영역을 준다. 초기값이 나쁘면 2차 수렴은 애초에 발동하지 않으므로, 전역 단계에서는 부분공간이, 국소 단계에서는 뉴턴이 일하는 셈이다. 데이비드슨 알고리즘과의 대비가 여기서 선명해진다. 데이비드슨은 확장 벡터를 t=(D−θI)−1rt = (D - \theta I)^{-1} r 로 잡는데(DD 는 AA 의 대각), 이것은 대각 우세한 양자화학 CI 행렬에서만 통하는 전처리다. 게다가 전처리기를 “개선”해서 (A−θI)−1(A - \theta I)^{-1} 을 그대로 쓰면 t∝ut \propto u 가 나와 부분공간이 자라지 않고 정체해 버리는 유명한 역설이 있다. 야코비-데이비드슨은 사영으로 그 성분을 애초에 제거함으로써 이 역설을 해소하고, 대각 전처리를 임의의 전처리기로 일반화한다.2

3. 내부-외부 반복과 전처리[편집]

실무에서 보정 방정식을 정확히 푸는 일은 없다. 어차피 θ\theta 자체가 근사값이므로 정밀하게 푸는 것은 낭비다. 표준 전략은 GMRES/BiCGStab 같은 반복법으로 몇 스텝만 돌리는 부정확 해법이고, 잔차 감소 목표를 외부 반복이 진행될수록 조이는 방식(예: 2−m2^{-m})을 쓴다. 이 내부-외부(inner-outer) 구조 덕에 계산량을 수렴 속도와 맞바꿀 수 있다.

전처리기 K≈A−θIK \approx A - \theta I 를 쓸 때는 그것도 같은 사영자로 감싼 K~=(I−uu∗)K(I−uu∗)\tilde{K} = (I-uu^{*})K(I-uu^{*}) 를 써야 한다. 다행히 K~−1z\tilde{K}^{-1}z 는

K~−1z=K−1z−u∗K−1zu∗K−1u K−1u\tilde{K}^{-1} z = K^{-1}z - \frac{u^{*}K^{-1}z}{u^{*}K^{-1}u}\,K^{-1}u

로 계산되므로, 외부 반복당 K−1uK^{-1}u 를 한 번 구해 두면 내부 반복마다 추가 비용은 K−1zK^{-1}z 한 번과 내적 두 번뿐이다. 불완전 LU든 다중격자든 영역 분할법 전처리든 그대로 꽂을 수 있다는 것이 이 방법의 실전 경쟁력이다.

4. 내부 고유값과 확장성[편집]

란초스 알고리즘이나 아놀디 알고리즘이 만드는 크릴로프 부분공간은 다항식 필터의 성질상 스펙트럼 양 끝의 고유값을 먼저 잡는다. 스펙트럼 한가운데 있는 내부 고유값(interior eigenvalue)을 원하면 이동-역변환으로 인수분해를 해야 하는데, 그게 감당 안 되는 크기가 진짜 문제다. 야코비-데이비드슨은 부분공간을 크릴로프로 고집하지 않으므로 목표값 τ\tau 근처를 직접 겨냥할 수 있고, 리츠 값 대신 (A−τI)−1(A-\tau I)^{-1} 에 대한 리츠 값을 쓰는 조화 리츠값(harmonic Ritz value) 추출을 곁들이면 내부 고유값 근사의 고질적 불안정성까지 잡힌다.

부분공간이 커지면 메모리와 직교화 비용이 선형·제곱으로 늘어나므로 mmax⁡m_{\max} 에 도달하면 좋은 리츠 벡터 몇 개만 남기고 재시작한다(thick restart). 이미 수렴한 고유쌍은 수축(deflation)으로 계에서 떼어내되, 고유벡터가 아니라 슈어 분해의 슈어 벡터를 누적해 부분 슈어 형 AQk=QkRkAQ_k = Q_kR_k 를 키워 가는 것이 표준이다. 결함 행렬이나 중근에서도 무너지지 않기 때문이다. 이 구현이 JDQR이고, 일반화 문제 Ax=λBxAx = \lambda Bx 에 대해 부분 일반화 슈어 형(QZ)을 키우는 판본이 JDQZ다. 이차 고유값 문제 (λ2M+λC+K)x=0(\lambda^2 M + \lambda C + K)x = 0 처럼 선형화하면 크기가 두 배로 뛰는 문제에서도 사영 구조를 그대로 유지할 수 있어, 감쇠 진동계의 모드 해석이나 유체 안정성 해석에서 자주 쓰인다.3

쓰이는 곳도 이 성질을 그대로 따라간다. 전자구조 계산에서 밀도범함수이론 해밀토니안의 최저 수백 개 고유쌍을 구하는 것은 데이비드슨 계열의 전통적 영역이고, 배치상호작용 행렬처럼 대각 우세가 강한 경우는 원조 데이비드슨으로도 충분하다. 반면 자기유체역학 안정성 해석이나 오어-조머펠트 방정식 이산화처럼 비대칭·비정규 스펙트럼 한가운데를 뒤져야 하는 문제, 감쇠 구조물의 복소 모드 추출이 야코비-데이비드슨의 무대다.

실무 구현으로는 SLEPc의 EPSJD/EPSGD, PRIMME, Trilinos의 Anasazi가 대표적이다. 대칭 양정부호 문제에 한정하면 LOBPCG가 더 간결하고 빠른 경우가 많으므로, 야코비-데이비드슨이 힘을 발휘하는 곳은 비대칭·내부·일반화·다항처럼 크릴로프 계열이 불편해지는 영역이다.4

5. 관련 문서[편집]

6. Footnotes[편집]

  1. 이 딜레마를 처음 보면 “그럼 역반복법은 왜 되는 거냐”는 의문이 든다. 답은 계가 나빠지면서 생기는 오차가 하필 최소 특이벡터, 즉 우리가 원하는 고유벡터 방향이라는 것. 세상에서 가장 유용한 수치적 재앙이다. 다만 이 행운은 직접 인수분해에만 해당하고 반복 해법에는 적용되지 않는다. ↩

  2. 이름의 “야코비” 쪽이 여기서 나온다. 야코비가 1846년 논문에서 쓴 것은 회전 대각화만이 아니라, 현재 근사 고유벡터에 직교하는 성분에 대해서만 보정을 푸는 절차였다. 150년 뒤에 그 아이디어가 재발견되어 데이비드슨과 결합한 셈이다. ↩

  3. 다항 고유값 문제를 선형화하면 2n2n 또는 ℓn\ell n 차 문제가 되는데, 크릴로프 계열은 그 확대된 공간에서 기저를 통째로 저장해야 한다. 사영 기반 방법은 확장 벡터를 원래 크기로 되돌려 다룰 수 있어 메모리에서 이긴다. TOAR 같은 구조 보존 크릴로프가 이 격차를 많이 좁히긴 했다. ↩

  4. “그래서 뭘 써야 하냐”에 대한 실무적 답: 대칭 + 양 끝 고유값이면 란초스, 대칭 양정부호 + 전처리기가 좋으면 LOBPCG, 비대칭이면 아놀디/크릴로프-슈어, 목표값 근처 내부 고유값인데 인수분해가 불가능하면 야코비-데이비드슨. 인수분해가 가능하면 그냥 이동-역변환 아놀디가 대체로 이긴다. ↩