6.156 병렬 행렬 연산과 GPU 가속

1. 병렬 연산의 필요성

로봇공학에서 실시간 제어, 시뮬레이션, 경로 계획 등의 연산 요구가 증가함에 따라, 단일 프로세서의 성능만으로는 충분하지 않은 경우가 빈번히 발생한다. 특히 다수의 로봇을 동시에 시뮬레이션하거나, 고차원 최적화 문제를 풀거나, 접촉이 포함된 물리 시뮬레이션을 수행할 때 대규모 행렬 연산이 병목이 된다. 병렬 연산(parallel computing)은 여러 처리 유닛이 동시에 작업을 수행하여 전체 계산 시간을 단축하는 방법이다.

2. 병렬 연산의 기본 개념

2.1 병렬화의 유형

행렬 연산의 병렬화는 크게 두 가지 수준에서 이루어진다.

데이터 병렬성(data parallelism): 동일한 연산을 여러 데이터 원소에 동시에 적용한다. 행렬-벡터 곱에서 각 결과 원소의 계산이 독립적이므로, 이를 병렬로 수행할 수 있다. GPU의 SIMT(Single Instruction, Multiple Threads) 아키텍처가 이에 최적화되어 있다.

태스크 병렬성(task parallelism): 서로 다른 연산을 동시에 수행한다. 다수의 독립적인 로봇의 동역학을 동시에 계산하거나, 파이프라인(pipeline) 방식으로 서로 다른 단계의 연산을 중첩하는 것이 이에 해당한다.

2.2 암달의 법칙과 구스타프슨의 법칙

p개의 프로세서를 사용할 때의 이론적 속도 향상(speedup)은 암달의 법칙(Amdahl’s law)으로 상한이 결정된다. 전체 연산 중 병렬화 가능한 비율을 f라 하면 다음과 같다.

S(p) = \frac{1}{(1 - f) + f/p} \leq \frac{1}{1 - f}

이 식에서 직렬 부분 (1 - f)가 0이 아닌 한 속도 향상에 상한이 존재한다. 예를 들어, 전체 연산의 95%가 병렬화 가능하다면 최대 속도 향상은 1/0.05 = 20배이다.

구스타프슨의 법칙(Gustafson’s law)은 문제 크기가 프로세서 수에 비례하여 증가하는 약한 스케일링(weak scaling) 관점을 제시한다.

S(p) = p - (1 - f)(p - 1) = 1 + f(p - 1)

이 법칙에 따르면, 문제 크기를 함께 키울 경우 프로세서 수에 거의 선형적인 속도 향상을 얻을 수 있다.

2.3 병렬 효율

병렬 효율(parallel efficiency)은 속도 향상을 프로세서 수로 나눈 값이다.

E(p) = \frac{S(p)}{p}

이상적인 경우 E(p) = 1이며, 통신 오버헤드, 부하 불균형, 직렬 구간 등으로 인해 실제로는 E(p) < 1이다.

CPU 기반 병렬 행렬 연산

공유 메모리 병렬화

멀티코어 CPU에서의 병렬화는 공유 메모리(shared memory) 모델에 기반하며, OpenMP가 대표적인 프로그래밍 인터페이스이다.

행렬-벡터 곱 \mathbf{y} = \mathbf{A}\mathbf{x}의 병렬화에서, 결과 벡터의 각 성분 y_i = \sum_{j=1}^{n} a_{ij} x_j는 독립적으로 계산할 수 있으므로 행 기반 분할(row-wise partitioning)이 자연스럽다. p개의 스레드가 각각 n/p개의 행을 담당하면, 이상적인 속도 향상은 p배이다.

행렬-행렬 곱 \mathbf{C} = \mathbf{A}\mathbf{B}에서는 결과 행렬의 각 블록을 독립적으로 계산할 수 있으며, 2차원 블록 분할(2D block partitioning)이 캐시 효율성을 높인다. \mathbf{C}를 p_r \times p_c 블록으로 분할하면 각 블록 \mathbf{C}_{ij}는 다음과 같이 계산된다.

\mathbf{C}_{ij} = \sum_{k=1}^{K} \mathbf{A}_{ik}\mathbf{B}_{kj}

2.4 분산 메모리 병렬화

대규모 클러스터에서는 분산 메모리(distributed memory) 모델을 사용하며, MPI(Message Passing Interface)가 표준 통신 라이브러리이다. 행렬을 여러 프로세스에 분배할 때, 1차원 행 분할, 1차원 열 분할, 2차원 블록 순환 분할(2D block-cyclic distribution) 등의 전략이 사용된다.

ScaLAPACK 라이브러리는 2차원 블록 순환 분할을 기반으로 하며, 행렬을 p_r \times p_c 프로세서 격자에 분배한다. 이 분할에서 행렬의 (i, j) 블록은 프로세서 (i \bmod p_r, \; j \bmod p_c)에 할당된다.

3. GPU 아키텍처와 행렬 연산

3.1 GPU의 구조적 특성

GPU(Graphics Processing Unit)는 수천 개의 경량 코어(core)로 구성되어 대규모 데이터 병렬 연산에 최적화된 프로세서이다. NVIDIA GPU의 경우 다음과 같은 계층적 구조를 가진다.

GPU와 CPU의 주요 차이점은 다음과 같다.

3.2 GPU 메모리 계층

GPU의 메모리 계층은 성능에 직접적인 영향을 미친다.

• 전역 메모리(global memory): 용량이 크지만(수~수십 GB) 접근 지연이 크다(수백 클럭 사이클). 행렬 데이터의 주 저장소이다.
• 공유 메모리(shared memory): 블록 내 스레드 간 공유되며, 접근 속도가 빠르다(수 클럭 사이클). 타일 기반 행렬 곱셈에서 핵심적으로 활용된다.
• 레지스터(register): 각 스레드에 할당된 가장 빠른 저장소이다.

4. GPU 기반 행렬-행렬 곱의 최적화

4.1 타일 기반 행렬 곱셈

GPU에서 행렬-행렬 곱 \mathbf{C} = \mathbf{A}\mathbf{B}를 효율적으로 수행하기 위해 타일(tile) 기반 알고리즘을 사용한다. 행렬을 T \times T 크기의 타일로 분할하고, 각 스레드 블록이 하나의 출력 타일 \mathbf{C}_{ij}를 계산한다.

\mathbf{C}_{ij} = \sum_{k=0}^{n/T - 1} \mathbf{A}_{ik}\mathbf{B}_{kj}

각 단계에서 \mathbf{A}_{ik}와 \mathbf{B}_{kj}를 공유 메모리에 적재한 후, 공유 메모리에서 부분 곱을 계산한다. 이 방식으로 전역 메모리 접근 횟수를 T배 줄일 수 있다.

구체적으로, 나이브 구현에서의 전역 메모리 접근 횟수가 2n^3이라면, 타일 기반 구현에서는 2n^3/T로 감소한다. 산술 강도(arithmetic intensity)는 다음과 같이 개선된다.

\text{AI}_{\text{naive}} = \frac{2n^3}{2 \cdot 8n^3} = \frac{1}{8} \; \text{FLOP/바이트}

\text{AI}_{\text{tiled}} = \frac{2n^3}{2 \cdot 8n^3/T} = \frac{T}{8} \; \text{FLOP/바이트}

T = 32이면 산술 강도가 32배 증가하여 메모리 대역폭 제한에서 벗어날 수 있다.

텐서 코어의 활용

최신 NVIDIA GPU(Volta 이후)에는 텐서 코어(Tensor Core)가 탑재되어 있으며, 4 \times 4 행렬의 FMA 연산을 단일 사이클에 수행한다.

\mathbf{D} = \mathbf{A}\mathbf{B} + \mathbf{C}, \quad \mathbf{A}, \mathbf{B} \in \mathbb{R}^{4 \times 4}

텐서 코어는 혼합 정밀도(mixed precision)를 지원하여, FP16 입력과 FP32 누적으로 높은 처리량을 달성한다. 로봇공학에서 시뮬레이션의 일부 단계(예: 신경망 추론, 대략적 동역학 계산)에 혼합 정밀도를 적용하면 성능을 크게 향상시킬 수 있다.

5. cuBLAS와 cuSPARSE

5.1 cuBLAS

NVIDIA의 cuBLAS 라이브러리는 GPU에서의 BLAS(Basic Linear Algebra Subprograms) 연산을 제공한다. 주요 함수는 다음과 같다.

cuBLAS의 Dgemm은 최신 GPU에서 피크 성능의 90% 이상을 달성하며, 배정밀도 기준으로 수 TFLOP/s의 처리량을 보인다.

5.2 cuSPARSE

희소 행렬 연산을 위한 cuSPARSE 라이브러리는 CSR(Compressed Sparse Row), CSC(Compressed Sparse Column), COO(Coordinate) 등의 희소 행렬 형식을 지원한다.

희소 행렬-벡터 곱(SpMV) \mathbf{y} = \mathbf{A}\mathbf{x}는 CSR 형식에서 다음과 같이 계산된다. 행 i의 비영 원소가 열 인덱스 j_1, j_2, \ldots, j_{m_i}에 위치하면 다음과 같다.

y_i = \sum_{l=1}^{m_i} a_{ij_l} x_{j_l}

GPU에서 SpMV의 병렬화는 각 행을 하나의 워프에 할당하는 방식이 일반적이다. 그러나 행별 비영 원소 수의 편차가 크면 부하 불균형이 발생하므로, 적응적(adaptive) 분할 전략이 필요하다.

병렬 행렬 분해

병렬 Cholesky 분해

대칭 양의 정부호 행렬의 Cholesky 분해 \mathbf{A} = \mathbf{L}\mathbf{L}^\top을 블록 단위로 병렬화할 수 있다. 행렬을 N_b \times N_b 블록으로 분할하면, 각 단계는 다음 네 가지 연산으로 구성된다.

1. 대각 블록의 Cholesky 분해: \mathbf{L}_{kk} = \text{chol}(\mathbf{A}_{kk}) — 직렬
2. 열 블록의 삼각 풀이: \mathbf{L}_{ik} = \mathbf{A}_{ik}\mathbf{L}_{kk}^{-\top}, i > k — 병렬화 가능
3. 대각 블록 갱신: \mathbf{A}_{ii} \leftarrow \mathbf{A}_{ii} - \mathbf{L}_{ik}\mathbf{L}_{ik}^\top, i > k — 병렬화 가능
4. 비대각 블록 갱신: \mathbf{A}_{ij} \leftarrow \mathbf{A}_{ij} - \mathbf{L}_{ik}\mathbf{L}_{jk}^\top, i > j > k — 대규모 병렬화 가능

단계 4의 DGEMM 호출이 전체 연산의 대부분을 차지하며, 이 부분이 가장 효율적으로 병렬화된다.

배치 행렬 연산

로봇공학에서는 동일한 크기의 소규모 행렬 연산을 대량으로 수행해야 하는 경우가 많다. 예를 들어 다음과 같다.

• 수백 개 접촉점 각각에서 3 \times 3 또는 6 \times 6 행렬 연산
• 몬테카를로 시뮬레이션에서 수천 개 샘플의 독립적 동역학 계산
• 경로 최적화에서 각 시간 단계의 독립적 행렬 분해

cuBLAS의 배치(batched) 함수는 이러한 연산을 단일 GPU 커널로 효율적으로 처리한다.

\mathbf{C}_i = \alpha \mathbf{A}_i \mathbf{B}_i + \beta \mathbf{C}_i, \quad i = 1, 2, \ldots, N_{\text{batch}}

cublasDgemmBatched는 N_{\text{batch}}개의 행렬 곱을 동시에 실행하며, 소규모 행렬(n \leq 32)에서 특히 높은 효율을 보인다.

6. 로봇공학 시뮬레이션에서의 GPU 활용

6.1 병렬 강체 동역학

다수의 로봇 또는 다수의 초기 조건에 대한 동역학 시뮬레이션을 GPU에서 병렬로 수행할 수 있다. 각 시뮬레이션 인스턴스의 관성 행렬 계산, Cholesky 분해, 전방/후방 대입을 배치 연산으로 처리하면, 수천 개의 시뮬레이션을 동시에 실행할 수 있다.

NVIDIA의 Isaac Gym과 같은 GPU 기반 물리 시뮬레이터는 이 원리를 활용하여 수천 개의 로봇 환경을 단일 GPU에서 동시에 시뮬레이션하며, 강화학습 기반 로봇 제어 연구에서 혁신적인 성능 향상을 가져왔다.

6.2 GPU에서의 접촉 역학

접촉이 포함된 시뮬레이션에서 접촉력 계산은 대규모 선형 상보성 문제(LCP) 또는 이차 프로그래밍(QP) 문제의 풀이를 필요로 한다. GPU에서의 병렬 Gauss-Seidel 반복법이나 병렬 ADMM(Alternating Direction Method of Multipliers)이 이러한 문제에 적용되고 있다.

6.3 CPU-GPU 전송 오버헤드

GPU 활용에서 주의할 점은 CPU와 GPU 간의 데이터 전송 오버헤드이다. PCIe 버스의 대역폭은 약 16~32 GB/s로, GPU의 메모리 대역폭(수백~수천 GB/s)에 비해 훨씬 낮다. 따라서 가능한 한 데이터를 GPU 메모리에 유지하고, CPU-GPU 간 전송을 최소화하는 설계가 중요하다. 전체 시뮬레이션 파이프라인을 GPU에서 수행하면 이 오버헤드를 크게 줄일 수 있다.

전송 시간 대비 연산 시간의 비율은 다음과 같이 추정할 수 있다.

\frac{T_{\text{compute}}}{T_{\text{transfer}}} = \frac{2n^3 / P_{\text{GPU}}}{2 \cdot 8n^2 / B_{\text{PCIe}}} = \frac{n \cdot B_{\text{PCIe}}}{8 P_{\text{GPU}}}

여기서 P_{\text{GPU}}는 GPU의 연산 성능(FLOP/s), B_{\text{PCIe}}는 PCIe 대역폭(바이트/s)이다. 이 비율이 1보다 커야 GPU 활용이 유리하며, 이는 행렬 크기 n이 충분히 커야 함을 의미한다.

참고문헌

• Kirk, D. B., & Hwu, W. W. (2016). Programming Massively Parallel Processors: A Hands-on Approach (3rd ed.). Morgan Kaufmann.
• Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
• NVIDIA (2024). cuBLAS Library User Guide. NVIDIA Corporation.
• Makoviychuk, V., et al. (2021). “Isaac Gym: High Performance GPU-Based Physics Simulation For Robot Learning.” NeurIPS 2021.
• Dongarra, J., et al. (2017). “MAGMA: a new generation of linear algebra libraries for GPU and multicore architectures.” SC17 Proceedings.

v 0.1