6.158 BLAS와 LAPACK 라이브러리의 활용
1. BLAS의 개요와 계층 구조
BLAS(Basic Linear Algebra Subprograms)는 기본 선형대수 연산의 표준화된 인터페이스이다. 1979년에 처음 제안된 이후 수치 선형대수학의 핵심 기반으로 자리 잡았으며, 하드웨어 제조사들이 각 프로세서에 최적화된 구현을 제공한다. BLAS는 연산의 복잡도와 데이터 접근 패턴에 따라 세 개의 레벨로 구분된다.
1.1 BLAS Level 1: 벡터-벡터 연산
BLAS Level 1은 벡터 간의 연산을 정의하며, 계산량은 \mathcal{O}(n), 데이터 접근량도 \mathcal{O}(n)이므로 산술 강도가 \mathcal{O}(1)이다. 주요 함수는 다음과 같다.
| 함수 | 연산 | FLOP 수 |
|---|---|---|
xAXPY | \mathbf{y} \leftarrow \alpha\mathbf{x} + \mathbf{y} | 2n |
xDOT | d = \mathbf{x}^\top\mathbf{y} | 2n |
xNRM2 | \lVert\mathbf{x}\rVert_2 | 2n |
xSCAL | \mathbf{x} \leftarrow \alpha\mathbf{x} | n |
xCOPY | \mathbf{y} \leftarrow \mathbf{x} | 0 (데이터 이동만) |
xSWAP | \mathbf{x} \leftrightarrow \mathbf{y} | 0 (데이터 이동만) |
IxAMAX | \arg\max_i \vert x_i \vert | n |
여기서 접두사 x는 데이터 타입을 나타내며, S(단정밀도), D(배정밀도), C(단정밀도 복소수), Z(배정밀도 복소수)가 있다.
1.2 BLAS Level 2: 행렬-벡터 연산
BLAS Level 2는 행렬과 벡터 간의 연산을 정의하며, 계산량은 \mathcal{O}(n^2), 데이터 접근량도 \mathcal{O}(n^2)이므로 산술 강도가 \mathcal{O}(1)이다. 주요 함수는 다음과 같다.
| 함수 | 연산 | FLOP 수 |
|---|---|---|
xGEMV | \mathbf{y} \leftarrow \alpha\mathbf{A}\mathbf{x} + \beta\mathbf{y} | 2mn |
xSYMV | 대칭 행렬-벡터 곱 | 2n^2 |
xTRMV | 삼각 행렬-벡터 곱 | n^2 |
xTRSV | 삼각 시스템 풀이 \mathbf{T}\mathbf{x} = \mathbf{b} | n^2 |
xGER | 순위-1 갱신 \mathbf{A} \leftarrow \alpha\mathbf{x}\mathbf{y}^\top + \mathbf{A} | 2mn |
xSYR | 대칭 순위-1 갱신 | n^2 |
로봇공학에서 xGEMV는 자코비안과 관절 속도의 곱 \dot{\mathbf{x}} = \mathbf{J}\dot{\mathbf{q}}, 관성 행렬과 가속도의 곱 \mathbf{M}\ddot{\mathbf{q}} 등에 직접 활용된다. xTRSV는 Cholesky 분해 후 전방/후방 대입에 사용된다.
1.3 BLAS Level 3: 행렬-행렬 연산
BLAS Level 3은 행렬 간의 연산을 정의하며, 계산량은 \mathcal{O}(n^3), 데이터 접근량은 \mathcal{O}(n^2)이므로 산술 강도가 \mathcal{O}(n)이다. 이 높은 산술 강도 덕분에 캐시 재사용이 극대화되어 피크 성능에 가장 가깝게 도달할 수 있다.
| 함수 | 연산 | FLOP 수 |
|---|---|---|
xGEMM | \mathbf{C} \leftarrow \alpha\mathbf{A}\mathbf{B} + \beta\mathbf{C} | 2mnk |
xSYMM | 대칭 행렬 곱 | 2mn^2 |
xTRMM | 삼각 행렬 곱 | mn^2 |
xTRSM | 다중 우변 삼각 풀이 \mathbf{T}\mathbf{X} = \mathbf{B} | mn^2 |
xSYRK | 대칭 순위-k 갱신 \mathbf{C} \leftarrow \alpha\mathbf{A}\mathbf{A}^\top + \beta\mathbf{C} | kn^2 |
xGEMM은 수치 선형대수학에서 가장 중요한 함수이며, 많은 고수준 알고리즘(LU 분해, Cholesky 분해, QR 분해 등)이 내부적으로 xGEMM을 호출하도록 설계된다. 이를 “BLAS Level 3 기반 알고리즘“이라 한다.
2. BLAS의 최적화 구현
2.1 하드웨어 최적화 BLAS 구현체
표준 BLAS 인터페이스의 최적화된 구현체는 다음과 같다.
| 구현체 | 제공자 | 특성 |
|---|---|---|
| Intel MKL (oneMKL) | Intel | Intel CPU에 최적화, AVX-512 지원 |
| OpenBLAS | 오픈소스 커뮤니티 | 다양한 CPU 지원, GotoBLAS 기반 |
| BLIS | 오픈소스 (UT Austin) | 모듈화된 프레임워크, 이식성 우수 |
| Apple Accelerate | Apple | Apple Silicon에 최적화 |
| ATLAS | 오픈소스 | 자동 튜닝 기반 |
| Arm Performance Libraries | Arm | Arm 프로세서에 최적화 |
최적화된 BLAS 구현은 나이브 구현에 비해 수십 배의 성능 향상을 제공한다. 예를 들어, 1000 \times 1000 행렬의 DGEMM에서 나이브 3중 루프는 피크 성능의 1~5%를 달성하는 반면, Intel MKL은 90% 이상을 달성한다.
2.2 GEMM 최적화의 핵심 기법
xGEMM의 고성능 구현에서 사용되는 핵심 기법은 다음과 같다.
블록 분할(blocking): 행렬을 L1, L2, L3 캐시 크기에 맞는 블록으로 분할하여 데이터 재사용을 극대화한다. 3단계 블록 분할에서 블록 크기 m_c, k_c, n_c는 각 캐시 수준에 최적화된다.
패킹(packing): 블록 데이터를 연속 메모리 영역에 복사하여 캐시 적중률을 높이고, 메모리 접근 패턴을 규칙적으로 만든다.
마이크로커널(micro-kernel): 가장 내부의 연산 루프를 레지스터 수준에서 최적화한 코드 조각이다. 일반적으로 m_r \times n_r 크기의 출력 블록을 레지스터에 유지하면서 순위-1 갱신을 반복 수행한다.
\mathbf{C}_{m_r \times n_r} \leftarrow \mathbf{C}_{m_r \times n_r} + \mathbf{a}_{m_r \times 1} \mathbf{b}_{1 \times n_r}^\top
AVX-512를 사용하는 경우, 배정밀도에서 m_r = 8, n_r = 6 또는 유사한 크기가 선택되어 48개의 출력 원소를 레지스터에 유지한다.
LAPACK의 구조와 주요 루틴
LAPACK(Linear Algebra PACKage)은 BLAS 위에 구축된 고수준 선형대수 라이브러리로, 선형 시스템 풀이, 고유값 문제, 특이값 분해, 최소자승 문제 등의 루틴을 제공한다.
LAPACK의 명명 규칙
LAPACK 함수 이름은 xYYZZZ 형식을 따르며, 각 부분의 의미는 다음과 같다.
x: 데이터 타입 (S,D,C,Z)YY: 행렬 유형 (GE=일반,SY=대칭,PO=양의 정부호,TR=삼각,HE=에르미트,GB=일반 밴드,PB=양의 정부호 밴드)ZZZ: 연산 유형
선형 시스템 풀이
| 함수 | 행렬 유형 | 방법 | 용도 |
|---|---|---|---|
xGESV | 일반 | LU 분해 + 피벗팅 | 일반 선형 시스템 |
xPOSV | SPD | Cholesky 분해 | 관성 행렬 시스템 |
xSYSV | 대칭 (부정치) | Bunch-Kaufman 분해 | KKT 시스템 |
xGBSV | 일반 밴드 | 밴드 LU 분해 | 밴드 구조 시스템 |
xPBSV | SPD 밴드 | 밴드 Cholesky | 밴드 SPD 시스템 |
xGELS | 일반 (과결정/미결정) | QR 또는 LQ 분해 | 최소자승 문제 |
로봇공학에서 가장 빈번히 사용되는 루틴은 xPOSV(관성 행렬 기반 시스템)와 xGELS(역기구학의 최소자승 풀이)이다.
행렬 분해
| 함수 | 분해 | 입력 행렬 유형 |
|---|---|---|
xGETRF | LU 분해 (부분 피벗팅) | 일반 |
xPOTRF | Cholesky 분해 | SPD |
xGEQRF | QR 분해 (Householder) | 일반 |
xSYTRF | \mathbf{L}\mathbf{D}\mathbf{L}^\top 분해 | 대칭 |
xGESVD | 특이값 분해 | 일반 |
xGEEV | 고유값 분해 | 일반 |
xSYEV | 대칭 고유값 분해 | 대칭 |
고유값 문제
로봇공학에서 고유값 분해는 다음과 같은 상황에서 사용된다.
조작성 해석: 자코비안 \mathbf{J}의 특이값 분해 \mathbf{J} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top에서 최소 특이값 \sigma_{\min}은 특이 자세에의 근접도를 나타낸다. 이는 xGESVD로 계산된다.
진동 해석: 로봇 시스템의 고유 진동수는 일반화 고유값 문제 \mathbf{K}\mathbf{v} = \lambda\mathbf{M}\mathbf{v}의 해로 구해진다. \mathbf{M}이 SPD이면 Cholesky 분해 \mathbf{M} = \mathbf{L}\mathbf{L}^\top을 통해 표준 고유값 문제로 변환할 수 있다.
\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\top}(\mathbf{L}^\top\mathbf{v}) = \lambda(\mathbf{L}^\top\mathbf{v})
이 변환된 문제는 xSYEV로 풀 수 있다.
주성분 분석(PCA): 로봇 센서 데이터의 차원 축소나 동작 분석에서 공분산 행렬의 고유값 분해가 사용된다.
3. LAPACK의 분해-풀이 패턴
LAPACK은 행렬 분해와 풀이를 분리하는 설계 패턴을 따른다. 이는 동일한 분해를 여러 우변 벡터에 재사용할 수 있게 해 준다.
예를 들어, Cholesky 분해 기반 풀이는 다음 세 단계로 수행된다.
- 분해:
xPOTRF로 \mathbf{A} = \mathbf{L}\mathbf{L}^\top을 계산한다. (\frac{1}{3}n^3 FLOP) - 풀이:
xPOTRS로 \mathbf{L}\mathbf{L}^\top\mathbf{x} = \mathbf{b}를 풀이한다. (2n^2 FLOP) - 정밀도 향상(선택):
xPORFS로 반복적 정밀도 향상을 수행한다.
로봇 제어에서 관성 행렬 \mathbf{M}(\mathbf{q})가 제어 주기 동안 거의 변하지 않는다면, 분해를 한 번 수행하고 여러 우변에 대해 풀이만 반복할 수 있다. 분해의 비용이 \mathcal{O}(n^3)이고 풀이의 비용이 \mathcal{O}(n^2)이므로, 우변이 여러 개인 경우 큰 절감 효과가 있다.
4. 로봇공학에서의 실제 사용 패턴
4.1 역동역학의 Cholesky 기반 풀이
순방향 동역학에서 관절 가속도를 구하기 위해 다음 시스템을 풀어야 한다.
\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} = \boldsymbol{\tau} - \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} - \mathbf{g}(\mathbf{q})
이 풀이의 LAPACK 기반 구현은 다음과 같은 호출 순서로 이루어진다.
- CRBA로 \mathbf{M}(\mathbf{q})를 계산한다.
- RNEA로 \mathbf{h} = \mathbf{C}\dot{\mathbf{q}} + \mathbf{g}를 계산한다.
- \mathbf{b} = \boldsymbol{\tau} - \mathbf{h}를 계산한다.
DPOTRF로 \mathbf{M} = \mathbf{L}\mathbf{L}^\top을 분해한다.DPOTRS로 \mathbf{L}\mathbf{L}^\top\ddot{\mathbf{q}} = \mathbf{b}를 풀이한다.
최소자승 역기구학
과잉 구동(redundant) 로봇의 역기구학에서 최소 노름 해를 구하기 위해 다음 최소자승 문제를 풀이한다.
\dot{\mathbf{q}} = \arg\min_{\dot{\mathbf{q}}} \|\dot{\mathbf{q}}\|^2 \quad \text{s.t.} \quad \mathbf{J}\dot{\mathbf{q}} = \dot{\mathbf{x}}
이는 DGELS를 호출하여 QR 분해 기반으로 풀 수 있으며, 또는 DGELSD(SVD 기반)를 사용하면 순위 부족(rank-deficient) 상황에서도 안정적인 해를 얻을 수 있다.
4.2 행렬 조건수 추정
실시간 제어에서 행렬의 조건수를 정확히 계산하는 것은 비용이 크므로(\mathcal{O}(n^3)), LAPACK의 조건수 추정 루틴 xPOCON(Cholesky 분해 후) 또는 xGECON(LU 분해 후)을 사용하여 \mathcal{O}(n^2)에 근사값을 얻는 것이 실용적이다. 조건수가 미리 설정한 임계값을 초과하면 특이 자세 회피 전략을 활성화하는 방식으로 활용할 수 있다.
5. 성능 고려 사항
5.1 행렬 저장 형식
LAPACK은 포트란(Fortran) 기반이므로 열 우선(column-major) 저장 순서를 사용한다. C/C++에서 호출할 때는 이 점에 주의하여야 하며, LAPACKE(C 인터페이스)는 행 우선(row-major) 옵션을 제공한다. 저장 순서의 불일치는 불필요한 전치 연산이나 캐시 비효율을 초래할 수 있다.
5.2 소규모 행렬에서의 오버헤드
BLAS와 LAPACK은 중~대규모 행렬에 최적화되어 있으며, 소규모 행렬(n < 10)에서는 함수 호출 오버헤드, 스레드 생성 비용 등이 실제 연산 비용보다 클 수 있다. 6자유도 로봇의 6 \times 6 관성 행렬에 대해 LAPACK을 호출하는 것은 비효율적일 수 있으며, 이 경우 수동으로 최적화된 인라인 코드가 더 빠르다.
| 행렬 크기 | LAPACK DPOTRF | 인라인 Cholesky | 비율 |
|---|---|---|---|
| 6 \times 6 | \approx 1 \; \mu\text{s} | \approx 0.05 \; \mu\text{s} | 20:1 |
| 30 \times 30 | \approx 3 \; \mu\text{s} | \approx 2 \; \mu\text{s} | 1.5:1 |
| 100 \times 100 | \approx 15 \; \mu\text{s} | \approx 50 \; \mu\text{s} | 0.3:1 |
이 표에서 볼 수 있듯이, n \geq 30 정도부터 LAPACK의 최적화가 효과를 발휘하기 시작하며, n \geq 100에서는 압도적인 성능 차이를 보인다.
참고문헌
- Anderson, E., et al. (1999). LAPACK Users’ Guide (3rd ed.). SIAM.
- Dongarra, J., et al. (1990). “A set of level 3 basic linear algebra subprograms.” ACM Transactions on Mathematical Software, 16(1), 1–17.
- Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
- Goto, K., & Van De Geijn, R. A. (2008). “Anatomy of high-performance matrix multiplication.” ACM Transactions on Mathematical Software, 34(3), 1–25.
- Van Zee, F. G., & Van De Geijn, R. A. (2015). “BLIS: A Framework for Rapidly Instantiating BLAS Functionality.” ACM Transactions on Mathematical Software, 41(3), 1–33.
v 0.1