희소 행렬 연산

14.1 배경: 대부분이 0인 행렬을 다루는 법

과학 계산이나 그래프 알고리즘을 조금만 깊게 파다 보면 한 가지 사실에 부딪힌다. 현실 문제의 행렬은 거의 다 희소하다는 점이다. 유한 요소법(FEM)으로 풀어낸 강성 행렬, 회로 시뮬레이션의 노드 행렬, 소셜 네트워크의 인접 행렬, 검색 엔진의 PageRank 행렬, 추천 시스템의 사용자-아이템 행렬 모두 비제로(nonzero) 비율이 1%도 채 안 되는 경우가 많다. 100만 × 100만 행렬에 비제로가 1000만 개라면 밀도(density)는 0.001%. 이걸 빽빽한 2차원 배열로 저장하면 메모리는 4 TB(8 byte 기준), 정작 의미 있는 값은 80 MB. 99.998%가 그냥 0을 들고 있는 셈이다.

희소 행렬을 다룬다는 건 두 가지를 동시에 푼다는 뜻이다. 첫째, 저장 형식(format)을 어떻게 잡을지. 비제로만 챙기되 위치 정보를 어디에 어떻게 끼워 넣을지가 관건이다. 둘째, 그 형식 위에서 희소 행렬-벡터 곱(SpMV)이나 희소 행렬-행렬 곱(SpMM) 같은 연산을 어떻게 효율적으로 수행할지. 특히 SpMV는 거의 모든 반복 솔버(conjugate gradient, GMRES, BiCGStab)의 핵심 워크호스이고, 그래프 알고리즘의 한 단계이기도 해서 가장 깊게 다뤄진다.

14.2 COO 형식과 단순 SpMV 커널

가장 직관적인 형식은 COO(Coordinate)다. 비제로 하나당 (행, 열, 값) 세 쌍을 그대로 저장한다. 세 개의 평행 배열 row[], col[], val[]를 둔다. 길이는 모두 nnz(number of non-zeros).

// 예: 4x4 행렬, nnz = 5
// [ 1 0 0 0 ]
// [ 0 0 2 3 ]
// [ 0 4 0 0 ]
// [ 0 0 0 5 ]

int   row[5] = {0, 1, 1, 2, 3};
int   col[5] = {0, 2, 3, 1, 3};
float val[5] = {1, 2, 3, 4, 5};

COO의 장점은 단순함이다. 행/열을 신경 쓰지 않고 비제로만 쭉 늘어놓으니, 동적으로 비제로를 추가하거나 제거할 때 자유롭다. 그래서 행렬 어셈블리(예: FEM 요소 행렬을 합성할 때)에 자주 등장한다. 또한 다른 형식으로 변환할 때 중간 매개체로도 쓰인다.

SpMV를 짠다고 하면 자연스러운 매핑은 비제로 한 개당 스레드 한 개다. 각 스레드가 자기 비제로의 val * x[col]을 계산해 y[row]에 더한다. 그런데 같은 행을 가진 비제로는 여러 스레드가 동시에 같은 y[row]에 누적해야 한다. 충돌을 피하려면 atomicAdd가 필요하다.

__global__ void spmvCOO(int nnz,
                         const int* row,
                         const int* col,
                         const float* val,
                         const float* x,
                         float* y) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= nnz) return;
  float contrib = val[i] * x[col[i]];
  atomicAdd(&y[row[i]], contrib);
}

코드는 짧고 부하 균형도 거의 완벽하다(스레드 하나당 비제로 하나, 모두 같은 일). 그러나 두 가지 비효율이 있다. 첫째, atomicAdd가 충돌이 잦은 행에서는 직렬화된다. 둘째, x[col[i]] 액세스가 col 패턴에 따라 무작위라 캐시 효율이 떨어진다. 그래도 매우 불규칙한 분포에서는 부하 균형 덕분에 의외로 잘 통한다.

14.3 CSR 형식: 행 단위로 정리하기

CSR(Compressed Sparse Row)은 가장 널리 쓰이는 형식이다. 비제로를 행 순서로 정렬하고, 행마다 시작 위치만 따로 기록한다. 세 배열을 둔다.

• row_ptr[]: 길이 (행수 + 1). i번째 행의 비제로는 row_ptr[i]부터 row_ptr[i+1]까지.
• col[]: 길이 nnz. 각 비제로의 열 인덱스.
• val[]: 길이 nnz. 각 비제로의 값.

// 위와 같은 4x4 예시
int   row_ptr[5] = {0, 1, 3, 4, 5};  // 행0: [0,1), 행1: [1,3), ...
int   col[5]     = {0, 2, 3, 1, 3};
float val[5]     = {1, 2, 3, 4, 5};

CSR로 SpMV를 짜는 가장 단순한 방식은 행 한 개당 스레드 한 개다. 한 스레드가 한 행의 모든 비제로를 순차적으로 읽고, 부분합을 모은 뒤 y[i]에 한 번만 쓴다. atomic이 사라진다.

__global__ void spmvCSR(int M,
                         const int* row_ptr,
                         const int* col,
                         const float* val,
                         const float* x,
                         float* y) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= M) return;
  int start = row_ptr[i];
  int end   = row_ptr[i + 1];
  float sum = 0.0f;
  for (int j = start; j < end; ++j) {
    sum += val[j] * x[col[j]];
  }
  y[i] = sum;
}

대신 두 가지 새로운 문제가 나타난다.

1. 제어 흐름 다이버전스(control divergence): 같은 워프의 스레드들이 서로 다른 길이의 행을 처리하면, 가장 긴 행이 끝날 때까지 워프 전체가 기다려야 한다. 한 행이 비제로 1000개, 옆 행이 5개면, 995회는 한 스레드만 일하는 셈이다.
2. 비코얼레싱(uncoalesced) 액세스: 같은 워프의 스레드들이 각자 자기 행의 첫 비제로를 읽는데, 그 비제로들의 글로벌 메모리 위치는 행마다 다르다. 워프가 한 사이클에 32개의 흩어진 주소를 읽으니 메모리 트랜잭션이 여러 번으로 쪼개진다.

이 두 문제는 행의 길이 분포에 강하게 의존한다. 길이 편차가 작으면 CSR이 쾌적하지만, 편차가 큰 그래프 인접 행렬에서는 처참할 수 있다. 해결책으로 한 행에 한 워프를 배정하는 변형(CSR-vector)이 있다. 워프 내 32 스레드가 협력해 한 행을 처리하면, 인접한 비제로를 읽으니 코얼레싱이 회복된다. 다만 행이 짧으면 워프 활용률이 떨어진다.

14.4 ELL 형식: 패딩으로 격자 맞추기

ELL(ELLPACK) 형식은 CSR의 다이버전스 문제를 정공법으로 해결한다. 모든 행을 가장 긴 행 길이 K에 맞춰 0으로 채우는 것이다(0-padding). 그러면 행마다 똑같이 K개의 슬롯이 생기고, 두 개의 2차원 배열 col[M][K], val[M][K]로 나타낼 수 있다.

그림 14.1 — ELL 형식. 모든 행을 K=2에 맞춘다. 짧은 행은 패딩으로 자리만 채운다.

여기서 한 가지 트릭이 더 있다. 메모리 레이아웃을 열-우선(column-major)으로 잡는 것이다. 즉 같은 슬롯 인덱스 k의 값이 메모리에 연속하도록 val[k * M + i] 같은 식으로 인덱싱한다. 그러면 행 i를 처리하는 스레드 i가 슬롯 k=0, 1, 2, ...를 차례로 읽을 때, 같은 워프의 스레드들이 인접한 i에 대해 같은 슬롯 k를 동시에 읽으니 메모리 액세스가 완벽하게 코얼레싱된다.

__global__ void spmvELL(int M, int K,
                         const int* col,    // [K][M] column-major
                         const float* val,
                         const float* x,
                         float* y) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= M) return;
  float sum = 0.0f;
  for (int k = 0; k < K; ++k) {
    int c = col[k * M + i];   // 워프 내 인접 i → coalesced
    sum += val[k * M + i] * x[c];
  }
  y[i] = sum;
}

ELL은 코얼레싱이 깔끔하고 다이버전스도 거의 없다. 다만 치명적인 약점이 하나 있다. K가 가장 긴 행에 의해 결정된다는 점이다. 대부분의 행이 짧고 한 두 행만 길면 메모리가 폭발한다. 평균 비제로 수가 5인데 최장 행이 1000이면 메모리는 200배가 된다. 패딩의 0을 곱셈에 끼우니 연산도 낭비된다.

14.5 하이브리드 ELL-COO: 양쪽의 좋은 점만

ELL의 약점을 무마하는 자연스러운 아이디어가 하이브리드(Hybrid, HYB) 형식이다. 행 길이의 어떤 임계값 K*를 정해 두고, 각 행에서 비제로 K* 개까지는 ELL에, 그 이상은 COO에 따로 보관한다. K*는 보통 평균 비제로 수 부근으로 잡는다(예: 행 길이의 중앙값).

그림 14.2 — 하이브리드의 분할 전략. 짧은 막대는 ELL로, 긴 꼬리는 COO로.

SpMV는 두 커널을 차례로 실행한다. ELL 부분으로 대부분의 트래픽을 처리하고, 짧은 COO 부분에 atomicAdd를 적용해 끝낸다. 메모리 낭비는 K*가 평균 부근일 때 거의 없고, 코얼레싱은 ELL 쪽에서 회복된다. 실험적으로 그래프 행렬 같은 편향 분포에서 단일 ELL이나 단일 CSR보다 1.5~3배 빠른 경우가 많다.

// 호스트 측 호출
spmvELL<<<..., ...>>>(M, Kstar, ell_col, ell_val, x, y);
spmvCOO<<<..., ...>>>(coo_nnz, coo_row, coo_col, coo_val, x, y);

14.6 JDS: 행을 정렬해서 패딩을 줄이기

또 다른 방향은 JDS(Jagged Diagonal Storage)다. 아이디어는 단순하다. 행을 비제로 수 기준 내림차순으로 정렬한 뒤 ELL과 비슷하게 슬롯별로 저장한다. 정렬했기 때문에 슬롯 k에서 일하는 행의 수가 점점 줄어들고, 이를 그대로 활용해서 패딩 0을 거의 없앤다.

그림 14.3 — JDS는 행을 길이 순으로 재배열해 슬롯별로 자연스럽게 줄어드는 직사각형(jagged) 모양을 만든다.

JDS는 슬롯 k별로 길이가 다른 별도 배열을 둔다. 슬롯 0에는 모든 행의 첫 비제로가, 슬롯 1에는 길이 ≥ 2인 행의 두 번째 비제로만 들어간다. 메모리 낭비는 거의 사라진다. 다만 정렬을 위해 perm[]이라는 행 순서 배열을 함께 두어야 하고, SpMV 결과를 원래 행 순서로 되돌리려면 추가 단계가 필요하다.

// JDS SpMV (대략적인 구조)
__global__ void spmvJDS(int M, int maxK,
                         const int* slot_start,  // 각 슬롯의 시작 오프셋
                         const int* col,
                         const float* val,
                         const int* perm,        // 정렬된 → 원래 행
                         const float* x,
                         float* y) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i >= M) return;
  float sum = 0.0f;
  for (int k = 0; k < maxK; ++k) {
    int base = slot_start[k];
    int rowsInSlot = slot_start[k + 1] - base;
    if (i >= rowsInSlot) break;        // 짧은 행은 일찍 종료
    sum += val[base + i] * x[col[base + i]];
  }
  y[perm[i]] = sum;
}

JDS는 메모리 효율은 최고지만, 한 워프 안에 짧은 행과 긴 행이 섞여 있으면 다이버전스가 다시 살아난다. 워프 단위로 행 길이가 비슷하도록 인접 정렬된 상태를 유지하기에, 슬롯 끝부분에서는 일부 워프만 일하는 것은 어쩔 수 없다.

14.7 형식 비교와 정리

네 형식의 trade-off를 한눈에 보자.

마지막으로, 실무에서는 cuSPARSE나 MAGMA가 위 형식들을 모두 지원한다. cuSPARSE의 경우 CSR이 기본이고, hybrid나 ELL은 deprecated 흐름이다. 대신 새로운 BSR(Block Sparse Row)이나 blocked-ELL이 도입되어 비제로가 작은 블록 단위로 모여 있는 패턴(딥러닝의 sparse attention 등)에 활용된다. 형식 선택은 결국 행렬의 통계(평균 비제로, 분산, 블록 패턴, 대칭성)에 달려 있다.