ARPACK

1. 개요[편집]

ARPACK ARnoldi PACKage
개발	Lehoucq · Sorensen · Yang (Rice University)
공개	1996 (사용자 안내서 SIAM, 1998)
언어	Fortran 77
알고리즘	암묵적 재시작 아놀디 / 란초스 (IRAM · IRLM)
인터페이스	역통신 (reverse communication)
현행 포크	ARPACK-NG
래퍼	SciPy eigsh/eigs, Octave eigs, MATLAB eigs(구버전)

eigsh(A, k=6) 한 줄 뒤에는 1996년에 멈춘 포트란 77 코드가 돌고 있다.

ARPACK(ARnoldi PACKage)은 대규모 희소행렬의 고유값 몇 개를 구하기 위한 포트란 77 수치 라이브러리로, 암묵적 재시작 아놀디법(IRAM)과 그 대칭판인 암묵적 재시작 란초스법(IRLM)의 표준 구현체다. 라이스 대학의 리쿡(R. B. Lehoucq)·소렌슨(D. C. Sorensen)·양(C. Yang)이 만들었고, 1998년 SIAM에서 나온 사용자 안내서가 사실상 이 분야의 실무 매뉴얼 역할을 했다.

알고리즘 자체의 설명은 아놀디 알고리즘과 란초스 알고리즘 문서에 있다. 이 문서는 라이브러리로서의 ARPACK — 어떤 API를 갖고, 어떤 모드가 있고, 어디서 사람들이 넘어지는지 — 를 다룬다. 30년 가까이 된 코드가 아직도 SciPy·Octave·수많은 상용 코드의 밑바닥에 깔려 있다는 사실 자체가 이 라이브러리의 위상이다.¹

2. 무엇을 계산하는가[편집]

주는 것은 $n \times n$ 행렬(또는 그것을 벡터에 곱하는 방법)이고, 받는 것은 사용자가 지정한 조건을 만족하는 고유쌍 nev개다. $n$은 수백만이어도 되고, nev는 보통 1~수십이다. 전부를 구하려는 시도는 설계 목적 밖이며, 그건 조밀 행렬용 LAPACK의 일이다.

핵심 동작은 크기 ncv의 아놀디/란초스 분해를 만든 뒤, 원하지 않는 리츠 값들을 시프트로 삼아 암묵적 QR 스텝을 걸어 그 성분을 걸러내고 nev 크기로 접어 넣는 것이다. 이 재시작 조작이 시작 벡터에 다항식 필터를 곱하는 것과 정확히 동등하면서 행렬 곱을 한 번도 추가로 요구하지 않는다는 것이 IRAM의 요지다. 사용자 입장에서 이것이 의미하는 바는 단순하다 — 메모리 사용량이 ncv개의 $n$차원 벡터로 고정된다. 수렴할 때까지 부분공간이 자라는 구조가 아니다.

수렴한 리츠 값은 잠금(locking)되고, 원하지 않는 것은 정화(purging)로 걷어낸다. 이 두 조작이 수치적으로 까다로웠던 점이 후속 크릴로프-슈어 방법이 나온 배경이다.

3. 역통신 인터페이스[편집]

ARPACK을 처음 만난 사람이 가장 당황하는 부분이 역통신(reverse communication)이다. 보통의 라이브러리는 사용자가 행렬이나 콜백 함수 포인터를 넘기면 라이브러리가 그걸 부른다. ARPACK은 반대로 한다. 루틴이 중간에 리턴하면서 “이 벡터에 $A$를 곱해서 다시 불러줘”라고 요청한다.

10 continue
   call dsaupd(ido, bmat, n, which, nev, tol, resid,
  &            ncv, v, ldv, iparam, ipntr, workd, workl,
  &            lworkl, info)
   if (ido .eq. -1 .or. ido .eq. 1) then
      call matvec(n, workd(ipntr(1)), workd(ipntr(2)))
      go to 10
   end if
   call dseupd(...)

ido가 상태 플래그다. $-1$·$1$이면 $\mathbf{y} \leftarrow OP\,\mathbf{x}$를, $2$면 $\mathbf{y} \leftarrow B\mathbf{x}$를 계산해 다시 호출하고, $3$이면 시프트를 직접 제공하며, $99$면 끝났다는 뜻이다. ipntr 배열이 작업 배열 workd 안에서 입력·출력 벡터의 오프셋을 알려준다.

이 설계는 포트란 77에 함수 포인터가 없다는 시대적 제약의 산물이지만, 결과적으로 매우 좋은 성질을 낳았다. 행렬이 무엇인지 ARPACK은 전혀 알 필요가 없다. 행렬이 디스크에 있든, 격자 위 스텐실 연산이든, MPI로 분산돼 있든, 심지어 물리 시뮬레이션 한 스텝을 돌리는 함수든 상관없다. 병렬 컴퓨팅 확장판인 PARPACK이 사용자 코드에 통신을 맡기고 자기는 내적 축약만 신경 쓰면 되는 구조도 여기서 나온다. matrix-free 고유값 해석의 원형이라 할 만하다.

대가도 있다. 원본 ARPACK은 호출 사이에 상태를 SAVE 변수에 보관하기 때문에 재진입 가능(reentrant)하지 않고, 따라서 스레드 안전하지 않다. 한 프로세스에서 여러 고유값 문제를 동시에 돌리려다 조용히 값이 섞이는 사고가 여기서 나온다.

4. 루틴 계열과 모드[편집]

이름 규칙은 [정밀도][문제유형]aupd / ...eupd다. aupd는 아놀디 반복을 돌리는 루틴, eupd는 수렴 후 고유값·고유벡터를 뽑아내는 후처리 루틴이다.

루틴	정밀도·유형	대상
dsaupd / dseupd	배정도 실수, 대칭	대칭·에르미트 실행렬 (IRLM)
dnaupd / dneupd	배정도 실수, 비대칭	일반 실행렬 (IRAM)
znaupd / zneupd	배정도 복소	복소 행렬
ssaupd · snaupd · cnaupd	단정도 대응판	위와 동일

무엇을 구할지는 which 문자열로 지정한다. 대칭은 'LA'(최대 대수값)·'SA'(최소 대수값)·'LM'·'SM'(최대·최소 절댓값)·'BE'(양 끝에서 절반씩), 비대칭은 'LM'·'SM'·'LR'·'SR'(실부)·'LI'·'SI'(허부)다.

iparam(7)로 지정하는 모드가 실무적으로 가장 중요하다.

• 모드 1 — 정규: 표준 문제 $A\mathbf{x} = \lambda\mathbf{x}$, $OP = A$. 사용자는 그냥 곱셈만 제공한다.
• 모드 2 — 정규 역: 일반화 문제 $A\mathbf{x} = \lambda M\mathbf{x}$에서 $OP = M^{-1}A$, $B = M$. $M$이 양정치여야 하며, 보통 촐레스키 분해를 한 번 해두고 재사용한다.
• 모드 3 — 이동-역변환: $OP = (A - \sigma M)^{-1}M$. 스펙트럼 내부를 노릴 때의 정답이며, 원리는 역반복법과 같다. 매 반복 선형계를 풀어야 하므로 희소 LU 분해를 미리 한 방 해두는 것이 전제다.
• 모드 4 — 좌굴: $K\mathbf{x} = \lambda K_G\mathbf{x}$ 꼴에 특화된 변환. 좌굴 고유치 해석용이다.
• 모드 5 — 케일리 변환: $(A - \sigma M)^{-1}(A + \sigma M)$. 특정 상황에서 모드 3보다 스펙트럼 사상이 유리하다.

5. 실무 함정[편집]

which='SM'은 함정이다. 크릴로프 부분공간법은 스펙트럼 바깥쪽부터 수렴한다. 절댓값이 가장 작은 고유값은 대개 스펙트럼 안쪽에 뭉쳐 있어서, 'SM'을 그대로 쓰면 수천 번을 돌아도 수렴하지 않거나 극도로 느리다. 올바른 처방은 $\sigma = 0$으로 이동-역변환을 켜고 'LM'을 요구하는 것 — SciPy 기준 eigsh(A, k, sigma=0, which='LM')이다. 이러면 $OP = A^{-1}$의 최대 고유값을 찾게 되어 원래 문제의 최소 고유값이 나온다.²

수렴하지 않으면 ncv를 늘린다. ncv는 유지할 기저 벡터 수이고 nev < ncv ≤ n이어야 한다. 권장은 대칭에서 ncv ≥ 2*nev, 비대칭에서 ncv ≥ 2*nev+1이다. 재시작 필터가 걸러낼 여유분 $p = \texttt{ncv} - \texttt{nev}$가 작으면 필터 다항식의 차수가 낮아져 수렴이 정체된다. ncv를 키우면 반복당 직교화 비용과 메모리가 늘지만, 안 끝나는 것보다는 낫다. 고유값이 뭉쳐 있거나 축퇴가 있으면 nev를 필요한 것보다 몇 개 더 요구하는 것도 표준 요령이다.

tol의 의미. 리츠 값의 상대 정확도이며, tol ≤ 0이면 기계정밀도로 해석된다. 관성항법이나 모드 해석에서 진동수만 필요하다면 tol=1e-8 정도로 느슨하게 잡아 반복을 크게 줄일 수 있다 — 레일리 몫의 제곱 정확도 덕에 고유벡터를 대충 얻어도 고유값은 훨씬 정확하다.

info 코드를 반드시 본다. info = 1은 최대 반복에 도달했고 일부만 수렴했다는 뜻이지 실패가 아니다. SciPy는 이걸 ArpackNoConvergence 예외로 올리면서 부분적으로 수렴한 고유쌍을 예외 객체에 담아 준다. 조용히 틀린 답을 반환하는 것보다 훨씬 정직한 설계지만, try/except로 뭉개고 지나가면 그 정보가 버려진다.

축퇴 모드. 단일 벡터 크릴로프 공간은 원리적으로 축퇴 고유값을 한 번에 하나만 잡는다. 대칭 구조물의 겹치는 모드는 반올림 오차 덕에 늦게라도 나오긴 하지만, 그걸 신뢰하고 설계하면 안 된다.

6. 생태계 — ARPACK-NG와 그 이후[편집]

원본 ARPACK은 2001년경 사실상 유지보수가 멈췄다. 데비안·Octave·SciPy 등 각자 패치를 들고 있던 배포자들이 모여 만든 커뮤니티 포크가 ARPACK-NG이며, 오늘날 리눅스 배포판이 libarpack이라 부르는 것은 대부분 이쪽이다. 주요 변화는 autotools·CMake 빌드 체계, 누적 버그 픽스, MPI 기반 PARPACK 통합, 그리고 ISO_C_BINDING 기반 C 인터페이스(ICB) 추가다. ICB 덕에 포트란 심볼 이름 장식 문제 없이 C/C++에서 직접 부를 수 있게 됐다.

• SciPy: scipy.sparse.linalg.eigsh가 dsaupd, eigs가 dnaupd/znaupd를 감싼다. 기본값은 k=6, ncv = min(n, max(2k+1, 20)), maxiter = 10n. 희소 특이값 분해 svds도 기본 솔버로 ARPACK을 쓴다.
• MATLAB: eigs는 오랫동안 ARPACK 래퍼였다. 이후 버전에서 내부 구현이 교체됐지만, 인자 이름(k, sigma, 'smallestabs' 등)과 이동-역변환 관행은 ARPACK 시절의 문법을 그대로 물려받았다.
• Octave / Julia / R: Octave eigs는 ARPACK-NG, Julia는 Arpack.jl 패키지로 분리돼 있다.
• SLEPc: PETSc 위에 세워진 현대적 고유값 문제 프레임워크로, ARPACK에 의존하지 않고 자체 크릴로프-슈어 방법 구현을 기본값으로 쓴다. 다만 외부 솔버로 ARPACK을 붙일 수 있고, 인터페이스 철학(무엇을 몇 개, 어느 변환으로)은 ARPACK이 정립한 것을 계승한다.

포트란 77 코드가 30년 뒤에도 파이썬 한 줄 뒤에서 돌고 있다는 사실은, 수치 라이브러리 세계에서 “제대로 검증된 구현”의 수명이 얼마나 긴지를 보여주는 표본이다.³

7. 관련 문서[편집]

• 아놀디 알고리즘 · 란초스 알고리즘 · 크리로프 부분공간법
• 고유값 문제 · 역반복법 · 레일리 몫
• 희소행렬 · LU 분해 · 촐레스키 분해
• 모드 해석 · 좌굴 · NASTRAN
• 특이값 분해 · 병렬 컴퓨팅
• LAPACK · SLEPc · 크릴로프-슈어 방법

8. Footnotes[편집]

1. 소프트웨어 고고학의 좋은 표본이다. numpy를 설치하면 딸려 오는 바이너리 안에 GOTO 문과 6자 변수명으로 짜인 1990년대 포트란이 들어 있고, 그게 당신의 그래프 라플라시안 고유벡터를 계산해 준다. 신뢰성은 세월로 증명됐지만, 스택 트레이스가 포트란으로 떨어지는 날의 기분은 별개 문제다. ↩

2. 이 함정에 걸린 사람이 하도 많아서 SciPy 문서가 아예 경고 문구를 박아 뒀다. “작은 고유값을 원하면 which='SM'이 아니라 sigma=0”은 사실상 이 바닥의 국룰. 물론 $A$가 특이하면 $\sigma$를 0에서 살짝 옮겨야 한다. ↩

3. 반대로 말하면 새 고유값 라이브러리가 시장을 뚫기가 그만큼 어렵다는 뜻이기도 하다. Spectra·PRIMME·SLEPc 같은 현대적 대안이 성능·기능에서 앞서는 경우가 많은데도, “일단 eigsh 써보고 안 되면 생각하자”에서 진도가 안 나가는 것이 현실이다. ↩