행렬 곱셈

수치해석 소프트웨어 마지막 수정: 2026-08-31 04:12:41

1. 개요[편집]

행렬 곱셈 Matrix multiplication
정의	$C_{ij} = \sum_k A_{ik} B_{kj}$
고전 연산량	$2n^3 - n^2$ flops
데이터량	$3n^2$ 워드 — 연산/데이터 비가 $O(n)$
실무 표준	BLAS 레벨 3 `gemm`
점근 지수	$O(n^\omega)$, $2 \le \omega < 2.372$
성능 병목	FLOP 수가 아니라 메모리 계층

삼중 for 루프 세 줄이면 끝나는 연산이다. 그 세 줄을 어떤 순서로 쓰느냐로 성능이 50배 갈린다.

행렬 곱셈(matrix multiplication)은 $A \in \mathbb{R}^{m \times k}$ , $B \in \mathbb{R}^{k \times n}$ 에 대해 $C_{ij} = \sum_{p} A_{ip}B_{pj}$ 로 정의되는 연산이며, 수치선형대수 전체의 성능 기준점이다. 정의만 보면 곱셈 $mnk$ 번과 덧셈 $mn(k-1)$ 번, 정사각 $n$ 이면 $2n^3 - n^2$ flops로 끝나는 시시한 물건이다.

그런데 이 문서가 다루는 진짜 질문은 따로 있다. 왜 같은 $2n^3$ flops를 수행하는 두 코드가 실측에서 수십 배 차이가 나는가. 답은 연산량이 아니라 데이터 이동에 있다. 고전 알고리즘의 점근 지수를 깎으려는 시도(슈트라센 이후)는 반세기째 이어졌지만 실무를 바꾸지 못했고, 정작 산업을 바꾼 것은 지수를 그대로 둔 채 캐시를 공략한 블록 알고리즘이었다. 희소행렬끼리의 곱은 비용 모형 자체가 달라지므로 그쪽 문서로 넘긴다 — 여기서는 밀집 행렬만 다룬다.

2. 연산/데이터 비 — 이 연산이 특별한 이유[편집]

$n \times n$ 곱셈은 $2n^3$ 연산을 하면서 데이터는 $3n^2$ 워드만 만진다. 즉 워드당 연산 수가 $O(n)$ 으로 자란다. LAPACK 문서에 나오는 BLAS 레벨 구분이 정확히 이 비율을 재는 것이다.

레벨	연산	연산량	데이터	연산/데이터
1 (`axpy`)	벡터-벡터	$2n$	$3n$	$O(1)$
2 (`gemv`)	행렬-벡터	$2n^2$	$n^2$	$O(1)$
3 (`gemm`)	행렬-행렬	$2n^3$	$3n^2$	$O(n)$

레벨 1·2 는 비율이 상수라, 아무리 잘 짜도 메모리 대역폭에서 막힌다. 요즘 서버 코어 하나가 사이클당 배정밀도 32 flops를 뽑는데 DRAM에서 워드 하나 끌어오는 데 수백 사이클이 걸리므로, 연산/데이터 비가 상수인 커널은 피크 성능의 5% 근처에서 천장을 친다. 이 “연산 강도 대 달성 가능 성능” 그림을 한 장으로 요약한 것이 루프라인 모델이고, gemm 은 그 그림에서 연산 한계 영역에 들어가는 거의 유일한 수치 커널이다. 레벨 3만이 이 벽을 뚫을 수 있다 — 데이터를 캐시에 한 번 올려놓고 $O(n)$ 번 재사용할 여지가 원리적으로 존재하기 때문이다.

“원리적으로 존재한다”와 “실제로 쓴다”는 다르다. 교과서적 삼중 루프

for i: for j: for p: C[i][j] += A[i][p] * B[p][j]

는 안쪽 루프에서 B 를 열 방향으로 훑는다. 행 우선 저장이면 매 접근이 캐시 라인 하나를 새로 끌어오고, $n$ 이 커져 행 하나가 L2에 안 들어가는 순간 재사용률이 0으로 수렴한다. 루프 순서를 i-p-j 로 바꾸는 것만으로도 몇 배가 개선되지만, 그건 시작일 뿐이다.

3. 블록(타일) 알고리즘[편집]

해법은 행렬을 $b \times b$ 블록으로 자르고, 블록 단위로 같은 곱셈을 재귀적으로 수행하는 것이다.

C_{IJ} \;\mathrel{+}=\; \sum_{P} A_{IP} B_{PJ}

블록 세 개 $A_{IP}, B_{PJ}, C_{IJ}$ 가 빠른 메모리에 동시에 들어가면( $3b^2 \le M$ ), 그 안에서 $2b^3$ flops가 데이터 재적재 없이 돌아간다. 전체 데이터 이동량은 블록 로드 횟수 $\times$ 블록 크기 $= \Theta(n^3/b) = \Theta(n^3/\sqrt{M})$ 이다. 캐시를 키우면 이동량이 그 제곱근에 반비례해 줄어든다는 뜻.

이 $\Theta(n^3/\sqrt{M})$ 가 개선 여지가 없는 하한이라는 것이 홍과 쿵의 1981년 결과이고, 그래서 블록 알고리즘은 “잘 하면 되는 요령”이 아니라 최적이다. 하한 이론과 분산 메모리 확장은 통신 회피 알고리즘 문서에 있다.

현대 구현(GotoBLAS 계보 — OpenBLAS, BLIS)은 이 아이디어를 캐시 계층마다 한 겹씩, 다섯 겹의 루프로 편다. 대략 이렇다.

$B$ 의 $k_c \times n_c$ 패널을 연속 메모리로 재포장해 L3에 상주시킨다.
$A$ 의 $m_c \times k_c$ 블록을 재포장해 L2에 상주시킨다.
가장 안쪽 마이크로커널은 $m_r \times n_r$ (AVX2 배정밀도면 흔히 $6\times8$ 이나 $8\times6$ ) 크기의 $C$ 조각을 레지스터에 통째로 올려놓고, $k_c$ 번의 랭크-1 갱신을 SIMD FMA 명령으로 때린다.

포장(packing)이 핵심이다. 복사 비용 $O(n^2)$ 를 지불하는 대신 안쪽 루프의 모든 접근이 순차 스트림이 되어 프리페처가 완벽하게 동작하고, TLB 미스도 사라진다. 이 구조로 짠 dgemm 은 실측에서 코어 피크의 85~95%를 낸다 — 수치 커널 중 이 수준에 도달하는 것은 사실상 이것뿐이다.¹

같은 C=AB 를 네 가지 순회 일정으로 브라우저에서 직접 재서 MFLOPS 막대로 쌓는다 — 전부 정확히 2N³ FLOP 인데 4×4 레지스터 블록이 ijk 삼중 루프의 3.7배(N=192 에서 1904 대 512 MF/s)다. 타일 크기는 8부터 64까지 밀어도 0.80~0.95배에 머문다: 이 규모에서 처리량을 정하는 것은 타일이 아니라 MAC 당 배열 로드 수(2.00 → 0.50)다.

위 실측이 보여 주는 것은 계층마다 병목이 다르다는 사실이다. $N \le 256$ 이면 $B$ 전체가 이미 L2에 들어앉아 캐시 타일링이 없앨 용량 미스가 아예 없고, 그래서 타일 크기를 8부터 64까지 밀어도 처리량이 거의 안 움직인다. 반면 가장 안쪽 레지스터 블록을 $1\times1$ 에서 $4\times4$ 로 키우면 MAC 당 배열 로드가 2.00에서 0.50으로 떨어지면서 3.7배가 나온다. 캐시 블로킹이 의미를 갖는 것은 작업 집합이 캐시를 실제로 넘칠 때이며, 그 전까지 이득을 내는 층은 레지스터다.

블록 크기는 공짜로 커지지 않는다. $b$ 를 키우면 재사용은 늘지만 작업 집합이 캐시를 넘어가는 순간 성능이 절벽처럼 떨어지고, 연관도(associativity) 충돌 때문에 특정 $b$ 에서만 국소적으로 무너지기도 한다. LAPACK이 블록 크기를 ILAENV 라는 기계 의존 함수로 물어보는 이유, BLIS가 캐시 크기·연관도·레지스터 수로부터 $m_c, k_c, n_c$ 를 해석적으로 유도하는 논문을 따로 쓴 이유가 여기 있다.

4. GEMM 위에 쌓인 세계[편집]

gemm 이 기준점이 되면 나머지 알고리즘의 설계 목표가 바뀐다 — 자기 연산량의 최대한을 gemm 호출로 밀어 넣는 것이 된다.

LU 분해의 블록판(dgetrf)이 표본이다. 열 하나씩 소거하는 대신 $n_b$ 열 패널을 먼저 분해하고( $O(n^2 n_b)$ , 레벨 2), 남은 부분행렬을 $C \mathrel{-}= A B$ 한 번으로 갱신한다(레벨 3). 전체 $\frac{2}{3}n^3$ flops 중 패널에 남는 몫은 $O(n^2 n_b)$ 뿐이라, $n \gg n_b$ 면 연산의 거의 전부가 gemm 안에서 벌어진다. QR 분해의 WY 표현, 촐레스키 분해의 블록판, 슈어 분해의 다중 시프트 QR, 삼중대각화가 전부 같은 전략이다.

전신인 LINPACK이 레벨 1 BLAS 기반이었다가 캐시 계층이 깊어진 기계에서 성능이 무너졌고, 그것을 블록 알고리즘으로 다시 쓴 것이 LAPACK이라는 족보도 결국 이 이야기다. 그래서 “라이브러리를 갈아 끼웠더니 10배 빨라졌다” 는 흔한 미스터리의 정답은 거의 항상 링크된 BLAS 구현이지 상위 코드가 아니다.

응용 쪽에서도 마찬가지다. 합성곱 신경망의 합성곱을 im2col로 펴서 gemm 으로 돌리는 것, 어텐션 메커니즘이 가산형 대신 내적형으로 수렴한 것, 혼합 정밀도 학습이 텐서 코어의 반정밀도 gemm 처리량을 노리는 것 — 전부 “이 하드웨어가 제일 잘하는 연산 모양으로 문제를 번역한다”는 같은 동작이다. GPU 컴퓨팅이 과학계산의 왕좌에 앉은 것도 결국 이 커널 하나 덕이다.

5. 곱셈 순서 — 결합법칙은 공짜가 아니다[편집]

행렬 곱은 교환법칙은 없지만 결합법칙은 성립한다. 그런데 결합 순서에 따라 연산량이 자릿수로 달라진다. $A \in \mathbb{R}^{n\times 1}$ , $B \in \mathbb{R}^{1 \times n}$ , $x \in \mathbb{R}^n$ 일 때 $(AB)x$ 는 $O(n^2)$ 이지만 $A(Bx)$ 는 $O(n)$ 이다. 랭크-1 행렬을 절대 조립하지 말라는 저랭크 근사의 실무 수칙이 이 한 줄에서 나온다.

여러 개를 연쇄로 곱할 때 최적 괄호 배치를 찾는 것은 동적 계획법의 교과서 예제(행렬 연쇄 곱셈)다. 그리고 이 고전적 교훈이 딥러닝에서 재발견됐다 — 선형 어텐션이 $\mathrm{softmax}$ 를 버린 대가로 얻은 것은 정확히 $(QK^\top)V$ 를 $Q(K^\top V)$ 로 묶을 자유이고, 그 결과 비용이 $O(n^2 d)$ 에서 $O(n d^2)$ 로 떨어진다.

6. 점근 지수 $\omega$ — 은하계 알고리즘의 세계[편집]

$2\times2$ 블록 곱셈은 정의대로면 곱셈 8번이다. 1969년 슈트라센은 덧셈을 18번 더 쓰는 대신 곱셈 7번으로 끝내는 항등식을 찾아냈다. 재귀로 펼치면

T(n) = 7\,T(n/2) + \Theta(n^2) \;\Longrightarrow\; T(n) = \Theta\!\left(n^{\log_2 7}\right) = \Theta(n^{2.807})

$O(n^3)$ 이 최적이 아님을 처음 보인 결과이고, 이후 반세기 동안 지수 $\omega$ 를 깎는 경주가 이어졌다. 쿠퍼스미스-위노그라드(1990)가 $\omega < 2.376$ , 이후 레이저 방법을 정교화한 후속 연구들이 $2.372$ 아래까지 밀어 내렸다. 하한은 자명하게 $\omega \ge 2$ (출력만 $n^2$ 개)이고, $\omega = 2$ 인지는 열린 문제다.

그런데 실무 코드는 거의 전부 $O(n^3)$ 을 돌린다. 이유가 셋이다.

상수. CW 계열의 숨은 상수는 천문학적이라, 이득이 나타나는 $n$ 이 우주에 존재하는 원자 수를 넘는다. 이런 부류를 은하계 알고리즘(galactic algorithm)이라 부른다. 존재 정리로서의 가치는 크지만 컴파일 대상이 아니다.
수치 안정성. 고전 알고리즘은 성분별 후진 오차 한계 $|\Delta C| \le n u\,|A||B|$ 를 만족한다 — $C$ 의 각 성분이 그 성분을 만드는 데 쓰인 항들의 크기로만 오염된다. 슈트라센은 항들을 더하고 빼면서 이 성질을 잃고, 노름 단위 한계( $\|\Delta C\| \le c(n)u\|A\|\|B\|$ )만 남는다. 크기 차가 큰 성분이 섞인 행렬에서 작은 성분의 상대오차가 통째로 날아갈 수 있다는 뜻이다.²
크로스오버와 형상. 슈트라센의 실측 크로스오버는 대략 $n \sim 10^3$ 이고, 그 지점에서도 이득은 한 자릿수 %다. 재귀 임시 버퍼로 메모리를 더 먹고, 홀수 차원 패딩이 필요하고, 무엇보다 얇고 긴 행렬(딥러닝과 블록 알고리즘에 흔한 형상)에서는 이득이 사라진다.

그래서 슈트라센이 실제로 쓰이는 경우는 “재귀를 한두 단만 펼치고 바닥에서 최적화된 gemm 을 부르는” 하이브리드뿐이며, 그것도 밀집 정사각 대형 행렬에 한정된다. 2022년 이후 강화학습으로 새 곱셈 항등식을 탐색한 결과들(작은 크기에서 기존 기록을 몇 개 깬)이 화제가 됐지만, 성격은 여전히 같다 — 조합론적 발견이지 커널 교체가 아니다.³

$(\min,+)$ 반환처럼 뺄셈이 없는 대수에서는 슈트라센류가 아예 성립하지 않는다는 점도 기억할 만하다. 플로이드-워셜 알고리즘이 블록 타일링은 그대로 물려받으면서 빠른 곱셈은 못 쓰는 이유다. 반대로 불 행렬곱은 정수 곱으로 환원되므로 $\omega$ 의 혜택을 받고, 그래서 이행 폐포의 복잡도가 $\omega$ 에 묶여 있다.

7. 실무 체크리스트[편집]

직접 짜지 마라. 세 줄 루프를 손으로 최적화해 벤더 BLAS를 이기는 일은 없다. 진짜로 필요하면 BLIS의 마이크로커널만 자기 ISA로 갈아 끼우는 것이 정공법이다.
gemm 호출 하나로 뭉칠 수 있으면 뭉쳐라. 작은 gemm 100번은 큰 gemm 1번보다 훨씬 느리다. 배치 gemm API가 존재하는 이유.
레이아웃을 먼저 정하라. BLAS는 열 우선이고, 전치 플래그(transa/transb)는 대개 공짜다 — 포장 단계에서 어차피 복사하기 때문. 반대로 코드 바깥에서 전치 복사를 하고 있다면 그건 낭비다.
실측하고 나서 말하라. “이론상 $n^3$ 이니까 $n$ 을 두 배로 하면 8배” 는 캐시에 들어가는 크기에서만 맞는다. 작업 집합이 L2를, 다시 L3를 넘어가는 지점마다 곡선이 꺾인다. 성능은 flops 계산기가 아니라 스톱워치가 정한다.

8. 관련 문서[편집]

그래서 HPC 강의의 첫 과제가 늘 이것이다. 순진한 삼중 루프로 시작해 루프 교환, 타일링, 포장, 레지스터 블로킹, 벡터화를 차례로 얹으며 GFLOPS 숫자가 계단식으로 올라가는 것을 직접 보는 것. 마지막에 벤더 BLAS와 비교하면 대개 절반쯤에서 멈추는데, 그 나머지 절반이 사람들이 20년간 갈아 넣은 것이다. ↩
성분별 한계와 노름 한계의 차이가 실무에서 문제가 되는 순간은 생각보다 드물다. 그래서 “슈트라센은 불안정하다”는 요약은 과장이고, 정확히는 “고전 알고리즘이 주는 강한 보장 하나를 포기한다”가 맞다. 다만 라이브러리는 사용자가 어떤 행렬을 넣을지 모르므로 강한 쪽을 기본값으로 둔다. 안전한 기본값은 원래 재미없다. ↩
화제가 됐던 결과들 중 상당수는 $\mathbb{Z}_2$ 같은 특정 계수체에서만 성립하거나 $4\times4$ · $5\times5$ 같은 작은 기저 크기에 한정된다. 기저 크기 $n_0$ 에서 곱셈 $r$ 번이면 $\omega \le \log_{n_0} r$ 이므로 지수는 확실히 내려가지만, 재귀 상수와 덧셈 개수까지 함께 봐야 실전 이득이 나온다. 반세기 동안 이 계산을 통과한 것은 여전히 슈트라센 하나뿐이다. ↩