LAPACK

1. 개요

LAPACK Linear Algebra PACKage
첫 공개	1992 (사용자 안내서 SIAM, 3판 1999)
언어	Fortran 90 (초기 F77)
전신	EISPACK(1970s) · LINPACK(1979)
대상	조밀·밴드 행렬 (희소는 대상 아님)
기반	BLAS 레벨 3
저장 방식	열 우선(column-major) + leading dimension
라이선스	수정 BSD

numpy.linalg.solve(A, b) 한 줄 뒤에는 1992년에 이름이 정해진 포트란 서브루틴이 있다.

LAPACK(Linear Algebra PACKage)은 조밀 및 밴드 행렬에 대한 연립방정식 풀이, 최소제곱, 고유값 문제, 특이값 분해의 사실상 표준 구현체다. LU 분해·QR 분해·촐레스키 분해·특이값 분해·고유값 문제를 다루는 모든 상용·오픈소스 코드가 결국 이 API로 내려간다. MATLAB의 백슬래시, NumPy/SciPy의 linalg, R의 solve, 상용 유한요소법 코드의 조밀 블록 처리가 전부 여기서 만난다.

이 문서는 알고리즘이 아니라 라이브러리로서의 LAPACK을 다룬다. 각 분해의 수학은 위 개별 문서에 있고, 여기서는 왜 이런 구조로 짜여 있으며 사용자가 어디서 넘어지는지를 본다. 결론부터 말하면 LAPACK을 이해한다는 것은 대부분 성능이 어디서 나오는가(BLAS 레벨 3) 와 호출 규약이 어디서 배신하는가(열 우선·lda·info) 두 가지를 아는 일이다.

2. BLAS 레벨과 블록 알고리즘

LAPACK의 성능은 사실상 BLAS(Basic Linear Algebra Subprograms) 위에 올라탄 결과다. BLAS는 세 레벨로 나뉜다.

레벨	연산	연산량	데이터	대표
1	벡터-벡터	$O(n)$	$O(n)$	daxpy, ddot
2	행렬-벡터	$O(n^2)$	$O(n^2)$	dgemv, dger
3	행렬-행렬	$O(n^3)$	$O(n^2)$	dgemm

핵심은 마지막 줄이다. 레벨 3만이 데이터 한 바이트당 연산 수가 $n$에 비례하므로, 캐시에 올린 블록을 여러 번 재사용해 부동소수점 유닛을 포화시킬 수 있다. 레벨 1·2는 아무리 잘 짜도 메모리 대역폭에서 막힌다.

족보를 짚으면 이렇다. 1970년대에 고유값 문제를 담당한 EISPACK과 1979년의 LINPACK이 먼저 있었고, 두 갈래를 하나의 규약으로 통합해 다시 쓴 것이 LAPACK이다. 전신인 LINPACK은 레벨 1 BLAS 기반이었고, 캐시 계층이 깊어진 1980년대 기계에서 성능이 무너졌다. LAPACK의 존재 이유가 여기 있다 — 같은 알고리즘을 블록 단위로 다시 짜서 작업량 대부분을 dgemm으로 밀어 넣는 것. 예컨대 dgetrf(LU)는 열 하나씩 소거하지 않고 $n_b$열 패널을 처리한 뒤 나머지를 행렬-행렬 곱으로 갱신한다. 블록 크기 $n_b$는 ILAENV가 기계별로 답한다.¹

3. 세 계층과 이름 읽는 법

루틴은 세 계층으로 나뉜다.

• 드라이버 루틴 — 문제 하나를 통째로 푼다. dgesv(연립방정식), dgels(최소제곱), dsyev(대칭 고유값), dgesvd(SVD).
• 계산 루틴 — 한 단계만 한다. dgetrf(LU 분해), dgetrs(삼각계 풀이), dgeqrf(QR 분해).
• 보조 루틴 — 하우스홀더 반사자 생성 dlarfg, 기계 상수 조회 dlamch 등.

이름은 암호가 아니라 규약이다. XYYZZZ 꼴로, X는 정밀도(s 단정도, d 배정도, c/z 복소), YY는 행렬 종류(ge 일반, sy 대칭, po 대칭 양정부호, gb 일반 밴드, tr 삼각, he 에르미트), ZZZ는 계산 내용(sv 풀이, trf 분해, trs 대입, ev 고유값, qrf QR)이다.

• dgesv = 배정도 + 일반 + 풀이
• dsyev = 배정도 + 대칭 + 고유값
• dpotrf = 배정도 + 양정부호 + 분해(= 촐레스키 분해)

드라이버에는 간이판과 확장판이 따로 있다. dgesv가 그냥 푸는 데 반해 dgesvx는 행·열 스케일링, 조건수 추정, 반복 개선, 전진·후진 오차 한계까지 돌려준다. 답만 필요한지 답의 신뢰도까지 필요한지에 따라 고르면 된다. 대칭 고유값 하나만 해도 dsyev(QR), dsyevd(분할 정복), dsyevr(MRRR), dsyevx(선택적 계산)가 다 있는데, 전부 필요하면 dsyevd, 일부 구간만 필요하면 dsyevr가 대체로 정답이다.

4. C·파이썬 사용자를 물어뜯는 지점

열 우선 저장. 포트란 관례를 그대로 쓴다. C의 2차원 배열은 행 우선이므로 그대로 넘기면 조용히 전치된 행렬을 푼다. 대칭 행렬로 시험하면 멀쩡히 통과했다가 비대칭에서 터지는 최악의 버그 패턴이 여기서 나온다.²

leading dimension(lda). 배열의 열 간 보폭이다. lda ≥ max(1, m)이면 되고, 큰 배열의 부분행렬을 복사 없이 넘길 수 있게 해 주는 장치다. 부분행렬을 넘기면서 lda를 부분행렬 크기로 준 순간 엉뚱한 메모리를 읽는다.

작업배열 조회. 다수의 루틴이 work 배열을 요구한다. lwork = -1로 한 번 부르면 실제 계산은 하지 않고 최적 크기를 work(1)에 실수로 돌려준다. 이걸 정수로 캐스팅해 할당한 뒤 진짜로 부르는 2단 호출이 관례다. 최소 크기만 주면 돌긴 하지만 언블록 경로로 떨어져 성능이 몇 배 손해다.

info 해석. info = 0이면 성공, info = -i면 $i$번째 인자가 잘못됐다는 뜻이며, info > 0은 루틴마다 의미가 다르다. dgesv에서 info = i는 $U_{ii}$가 정확히 0이라 특이하다는 뜻이고, dsyev에서는 수렴하지 못한 비대각 원소 개수다. 양수 info를 무시하고 결과 배열을 그대로 쓰는 것이 실무에서 가장 흔한 사고다.

5. 연산량과 오차 한계 — 직접 검산

LAPACK을 고를 때 실제로 쓰이는 숫자는 두 종류다. 연산량(얼마나 걸리나)과 오차 한계(답을 믿어도 되나). 둘 다 교과서 공식이 있지만, 그 공식이 언제부터 맞는지를 직접 세어 보면 쓸모가 달라진다.

5.1. 연산량 공식은 «점근적»이다

비블록 LU 분해와 촐레스키 분해의 곱셈·덧셈·나눗셈을 루프 안에서 하나씩 세면 이렇게 나온다.

n	LU 실측 연산	2n³/3	비	촐레스키 실측	n³/3	비
10	615	667	0.92	385	333	1.16
100	661,650	666,667	0.99	338,350	333,333	1.02
1000	666,166,500	666,666,667	0.999	333,833,500	333,333,333	1.002

$2n^3/3$과 $n^3/3$은 $n$이 커질 때만 맞는 값이다($n=10$에서는 8~16% 어긋난다). 그래서 작은 행렬 수백만 개를 돌리는 코드에서는 이 공식으로 시간을 추정하면 안 되고, 실제로 그런 영역에서는 호출 오버헤드와 언블록 경로가 지배한다. 반대로 비는 촐레스키가 LU의 정확히 절반이라는 점은 $n=100$부터 소수점 둘째 자리까지 성립한다 — 대칭 양정부호임을 알면서 dgesv를 부르는 것은 계산량을 두 배 쓰는 짓이라는 실무 결론이 여기서 나온다.

5.2. 부분 피벗의 성장 인자: 2^(n−1)은 실존한다

부분 피벗 가우스 소거법의 성장 인자(growth factor) $\rho_n = \max|u_{ij}| / \max|a_{ij}|$는 $2^{n-1}$로 상계가 잡혀 있지만 그 상계가 달성 가능하다. 윌킨슨이 제시한 행렬 — 대각 1, 대각 아래 모두 −1, 마지막 열 모두 1 — 로 직접 재면 정확히 맞는다.

n	max|U| / max|A| 실측	2^(n−1)
5	16	16
10	512	512
20	524,288	524,288
24	8,388,608	8,388,608

$n=24$에서 이미 성장이 $8.4\times10^6$이므로, 배정도의 유효 자릿수 16자리 중 7자리를 이 행렬 하나가 먹는다. 그런데도 실무에서 dgesv를 믿는 이유는, 무작위 행렬에서 측정되는 성장이 $n^{2/3}$ 수준에 머물러 이 상계와는 비교가 안 되게 작기 때문이다. 즉 부분 피벗은 “최악의 경우 불안정하지만 평균적으로는 매우 안정”한 전략이고, 이 간극이 수치선형대수에서 가장 유명한 미해결 간극 중 하나다.³ 진짜 걱정되면 완전 피벗(dgetc2 계열)이나 성장 인자 감시가 답이지, 손으로 짠 소거법은 답이 아니다.

5.3. 오차 한계: 조건수가 전부다

dgesvx가 돌려주는 전진 오차 한계의 뼈대는 $\|\hat x - x\| / \|x\| \lesssim \varepsilon \cdot \kappa(A)$다. 악명 높은 힐베르트 행렬 $H_{ij} = 1/(i+j-1)$로 $Hx = b$($x$는 전부 1)를 배정도로 풀어 보면 이 한계가 어떻게 작동하는지 한눈에 보인다.

n	κ∞(H)	실제 상대 전진 오차	εκ 한계	한계 대비
6	2.9×10⁷	5.3×10⁻¹⁰	3.2×10⁻⁹	0.16
10	3.5×10¹³	3.2×10⁻⁴	3.9×10⁻³	0.08
12	3.8×10¹⁶	0.36	4.3	0.08

읽을 것이 세 가지다. 첫째, 한계는 실제 오차를 10배 안쪽으로 과대평가한다 — 비관적이지만 자릿수는 맞게 잡는다. 둘째, $n=12$에서 $\kappa > 1/\varepsilon$이 되면 상대 오차 한계가 1을 넘고, 실제로 답의 36%가 틀린다. 이 지점을 넘으면 알고리즘을 바꿔도 안 된다 — 문제 자체가 배정도로 표현할 수 있는 정보를 다 썼기 때문이고, 처방은 정규화나 고정밀 연산이다. 셋째, 그래서 dgecon·dgesvx의 조건수 추정은 장식이 아니다. 피벗이 0이 아니라고 info = 0이 떴어도, $\kappa$가 $10^{16}$이면 그 답은 숫자 모양의 잡음이다.⁴

6. 생태계와 경계선

API는 하나지만 구현은 여럿이다. 넷립의 레퍼런스 LAPACK이 정의를 쥐고 있고, 실제 속도는 커널을 손으로 최적화한 구현이 낸다 — OpenBLAS, BLIS, 인텔 oneMKL, AMD AOCL, 애플 Accelerate. 분산 메모리는 2차원 블록 순환 분배를 쓰는 ScaLAPACK, 멀티코어 타일 알고리즘은 PLASMA, GPU 컴퓨팅 쪽은 MAGMA·cuSOLVER·rocSOLVER가 같은 이름 규약을 대체로 승계했다. 즉 인터페이스만 붙잡고 있으면 하드웨어를 갈아타도 코드가 산다.

경계선도 분명하다. LAPACK은 희소행렬을 다루지 않는다. 0을 전부 저장하고 채움(fill-in)을 그대로 받아들이는 설계라, 유한요소 강성행렬을 통째로 넘기면 메모리부터 터진다. 그쪽은 SuiteSparse(UMFPACK·CHOLMOD), MUMPS, SuperLU, PETSc의 영역이고, 희소 고유값은 ARPACK이다. 다만 이 도구들도 내부의 조밀 블록·초노드 연산에서는 결국 BLAS 3와 LAPACK을 부른다.

마지막으로 실무 결론 하나. 직접 짜지 마라. 가우스 소거법을 손으로 구현하는 것은 교육적으로 훌륭하지만, 30년 넘게 검증된 피벗 전략·스케일링·후진 안정성·예외값 처리를 재현할 수는 없다.⁵

7. 관련 문서

• LU 분해 · QR 분해 · 촐레스키 분해 · 특이값 분해
• 고유값 문제 · 최소자승법 · 조건수
• 희소행렬 · ARPACK · 반복법
• 병렬 컴퓨팅 · GPU 컴퓨팅 · 부동소수점 연산
• 가우스 소거법 · 기븐스 회전 · 강성행렬

8. Footnotes

1. 그래서 “레퍼런스 LAPACK이 느리다”는 말은 절반만 맞다. 느린 건 대개 함께 링크된 레퍼런스 BLAS 쪽이고, 같은 LAPACK 소스에 OpenBLAS를 물리면 성능이 자릿수 단위로 달라진다. 파이썬 환경 두 개가 같은 코드에서 10배 차이 나는 미스터리는 거의 항상 여기서 나온다. ↩

2. C 표준 인터페이스 LAPACKE는 LAPACK_ROW_MAJOR를 받아 주지만, 내부적으로 전치 복사를 하거나 등가 문제로 바꿔 푸는 것이라 공짜가 아니다. 정공법은 애초에 열 우선으로 데이터를 잡는 것이고, NumPy에서 order='F'를 쓰면 복사 한 번을 아낀다. 물론 그 한 번을 아끼려다 배열 레이아웃을 헷갈려 하루를 쓰는 게 국룰이다. ↩

3. 무작위 행렬의 성장 인자가 왜 그렇게 작은지는 아직 완전히 설명되지 않았다. 러프하게는 “소거가 진행되면서 부호가 섞여 상쇄된다”는 설명이 쓰이는데, 이 상쇄를 엄밀하게 정량화하는 것이 어렵다. 윌킨슨 행렬이 상계를 달성하는 구조를 전부 분류한 것은 히검 부부(Higham & Higham, 1989)다. ↩

4. 그래서 numpy.linalg.solve가 아무 경고 없이 답을 돌려줬다는 사실은 그 답이 맞다는 증거가 전혀 아니다. numpy.linalg.cond를 한 번 찍는 데 드는 비용은 분해 한 번 값이고, 그 한 줄이 “왜 결과가 이상하냐”로 태우는 반나절보다 싸다. ↩

5. 유명한 반례로 LINPACK 벤치마크가 있다. TOP500 슈퍼컴퓨터 순위가 아직도 조밀 연립방정식 풀이 속도로 매겨지는 이유는, 이 문제가 실제 워크로드를 대표해서가 아니라 40년치 비교 가능한 숫자가 쌓여 있어서다. 그 관성이 싫어서 나온 것이 HPCG 벤치마크다. ↩