[Linear Algebra] BLAS: Basic Linear Algebra Subprograms

BLAS(Basic Linear Algebra Subprograms) 는 vector와 matrix 연산을 위한 low-level routine의 표준 interface이다. 중요한 점은 BLAS 자체가 하나의 특정 library 이름이 아니라는 것이다.

BLAS는 다음과 같은 약속에 가깝다.

벡터 내적은 dot이라는 형태로 제공한다.
행렬-벡터 곱은 gemv라는 형태로 제공한다.
행렬-행렬 곱은 gemm이라는 형태로 제공한다.

실제 구현체는 따로 있다.

구현체	특징
OpenBLAS	널리 쓰이는 open-source BLAS 구현체
Intel MKL	Intel CPU에 강하게 최적화된 수치 연산 library
BLIS	BLAS-like framework와 구현체
cuBLAS	NVIDIA GPU용 BLAS 구현체
Apple Accelerate	macOS/iOS에서 제공되는 optimized math framework

NumPy, SciPy, MATLAB, R 같은 수치 계산 도구는 내부적으로 이런 BLAS 구현체를 사용한다. 그래서 Python loop로 직접 행렬 곱을 구현하는 것보다 numpy.matmul이 훨씬 빠르다.

BLAS가 하는 일

선형대수 계산은 결국 몇 가지 기본 연산으로 반복된다.

연산	수식	예
vector scaling	$x \leftarrow \alpha x$	벡터 크기 조정
vector addition	$y \leftarrow \alpha x + y$	gradient update
dot product	$x^Ty$	유사도, projection
matrix-vector product	$y \leftarrow Ax$	linear model
matrix-matrix product	$C \leftarrow AB$	neural network, simulation

BLAS는 이런 연산을 hardware에 맞게 최적화된 routine으로 제공한다. 사용자는 수학적 연산 단위로 코드를 작성하고, 실제 loop ordering, cache blocking, vectorization은 BLAS 구현체가 처리한다.

Level 1 BLAS

Level 1 BLAS는 vector-vector operation을 다룬다. 주로 $O(n)$ 연산이다.

AXPY

대표 연산은 axpy이다. 이름은 다음 수식에서 왔다.

y \leftarrow \alpha x + y

여기서 $x,y\in\mathbb{R}^n$ 이고 $\alpha$ 는 scalar이다.

Gradient descent update도 비슷한 형태로 볼 수 있다.

w \leftarrow w - \eta \nabla L(w)

이를 axpy 형태로 쓰면,

w \leftarrow (-\eta)\nabla L(w) + w

이다.

DOT

두 vector의 내적은 다음과 같다.

\operatorname{dot}(x,y) = x^Ty = \sum_{i=1}^{n}x_i y_i

이는 cosine similarity, projection, least squares 등에서 계속 등장한다.

SCAL

Vector scaling은 다음 연산이다.

x \leftarrow \alpha x

NRM2

Euclidean norm은 다음과 같다.

\|x\|_2 = \sqrt{x^Tx}

BLAS에서는 이런 vector primitive가 Level 1에 속한다.

Level 2 BLAS

Level 2 BLAS는 matrix-vector operation을 다룬다. 대표적인 복잡도는 $O(n^2)$ 이다.

GEMV

가장 대표적인 연산은 gemv이다. General matrix-vector multiplication을 의미한다.

y \leftarrow \alpha Ax + \beta y

여기서

A\in\mathbb{R}^{m\times n}, \quad x\in\mathbb{R}^{n}, \quad y\in\mathbb{R}^{m}

이다.

Linear model은 이 형태로 볼 수 있다.

\hat{y}=Xw

여기서 $X$ 는 data matrix, $w$ 는 weight vector이다.

GER

ger는 rank-1 update이다.

A \leftarrow \alpha xy^T + A

여기서 $xy^T$ 는 outer product이다.

xy^T = \begin{bmatrix} x_1y_1 & x_1y_2 & \cdots & x_1y_n \\ x_2y_1 & x_2y_2 & \cdots & x_2y_n \\ \vdots & \vdots & \ddots & \vdots \\ x_my_1 & x_my_2 & \cdots & x_my_n \end{bmatrix}

Level 3 BLAS

Level 3 BLAS는 matrix-matrix operation을 다룬다. 대표적인 복잡도는 $O(n^3)$ 이다.

가장 중요한 연산은 gemm이다.

C \leftarrow \alpha AB + \beta C

여기서

A\in\mathbb{R}^{m\times k}, \quad B\in\mathbb{R}^{k\times n}, \quad C\in\mathbb{R}^{m\times n}

이다.

Matrix multiplication은 naive하게 보면 세 겹 loop이다.

for i in range(m):
    for j in range(n):
        for p in range(k):
            C[i, j] += A[i, p] * B[p, j]

하지만 이 코드는 보통 최적화된 BLAS gemm보다 훨씬 느리다. 이유는 단순히 loop가 Python이라서만은 아니다. Compiled language로 작성해도 memory hierarchy를 고려하지 않으면 성능이 낮다.

왜 GEMM이 중요한가

현대 CPU/GPU에서 빠른 계산을 하려면 연산량만 보는 것이 부족하다. Memory에서 데이터를 가져오는 비용이 크기 때문이다.

Matrix multiplication은 data reuse가 좋다. 예를 들어 $C=AB$ 에서 $A$ 의 한 block과 $B$ 의 한 block은 여러 output element 계산에 반복 사용된다.

BLAS 구현체는 이 특성을 이용해 다음 최적화를 수행한다.

최적화	의미
Blocking/Tiling	matrix를 cache에 잘 맞는 block으로 나눈다.
SIMD vectorization	CPU vector instruction을 사용한다.
Loop unrolling	loop overhead를 줄이고 instruction scheduling을 개선한다.
Multi-threading	여러 core에 작업을 나눈다.
Packing	memory access pattern이 좋도록 data layout을 재배치한다.

그래서 deep learning, computer graphics, scientific computing에서 matrix multiplication 성능은 전체 성능에 직접적인 영향을 준다.

NumPy와 BLAS

NumPy array 연산은 Python에서 호출하지만, 큰 선형대수 연산은 C/Fortran으로 구현된 optimized routine으로 내려간다.

예를 들어 다음 연산은 Python loop가 아니라 BLAS로 처리될 수 있다.

import numpy as np

A = np.random.randn(1000, 1000)
B = np.random.randn(1000, 1000)

C = A @ B

반면 직접 loop를 작성하면 Python interpreter overhead가 커진다.

C = np.zeros((1000, 1000))
for i in range(1000):
    for j in range(1000):
        for k in range(1000):
            C[i, j] += A[i, k] * B[k, j]

수치 계산 코드에서 중요한 원칙은 다음이다.

가능하면 작은 scalar 연산을 Python loop로 반복하지 말고, 큰 vector/matrix 연산으로 묶어서 BLAS가 처리하게 만든다.

Memory Layout

BLAS는 원래 Fortran ecosystem에서 발전했기 때문에 column-major layout과 깊은 관련이 있다. 반면 C와 NumPy 기본 배열은 row-major layout을 많이 쓴다.

Layout	연속적으로 저장되는 방향
Row-major	같은 row의 원소들이 연속
Column-major	같은 column의 원소들이 연속

Row-major matrix는 다음 순서로 저장된다.

[[1, 2, 3],
 [4, 5, 6]]

memory: 1, 2, 3, 4, 5, 6

Column-major matrix는 다음 순서로 저장된다.

memory: 1, 4, 2, 5, 3, 6

Array가 contiguous하지 않거나 stride가 복잡하면 BLAS 호출 전에 copy가 발생할 수 있다. 따라서 성능이 중요한 코드에서는 shape뿐 아니라 memory layout도 고려해야 한다.

np.ascontiguousarray(x)
np.asfortranarray(x)

BLAS Level 정리

BLAS는 세 level로 나누어 이해하면 쉽다.

Level	연산	대표 routine	복잡도
Level 1	vector-vector	`axpy`, `dot`, `scal`, `nrm2`	$O(n)$
Level 2	matrix-vector	`gemv`, `ger`	$O(n^2)$
Level 3	matrix-matrix	`gemm`, `syrk`, `trmm`	$O(n^3)$

성능 관점에서는 Level 3가 특히 중요하다. 연산량은 크지만 data reuse가 좋아 cache와 vector unit을 잘 활용할 수 있다.

정리

BLAS는 선형대수 연산을 빠르게 하기 위한 표준 interface이다. 우리가 직접 세 겹 loop를 작성하지 않아도 되는 이유는, 오래 최적화된 BLAS 구현체들이 hardware에 맞는 방식으로 계산을 수행하기 때문이다.

핵심은 다음 한 줄이다.

\boxed{ \text{작은 scalar loop를 큰 matrix operation으로 바꾸면 BLAS가 일을 한다.} }

그래서 NumPy를 빠르게 쓰는 방법도 결국 BLAS 친화적인 코드를 쓰는 것이다.

# 느리기 쉬운 방식
for i in range(n):
    y[i] = np.dot(A[i], x)

# BLAS가 잘 처리하는 방식
y = A @ x

수치 계산에서 vectorization은 단순한 문법 취향이 아니라, optimized linear algebra backend를 활용하기 위한 구조적 선택이다.