NumPy의 브로드캐스팅, 인덱싱, 벡터화된 함수 설계가 배열 프로그래밍을 불필요하게 어렵고 읽기 힘들게 만드는 이유.
최고
주제
구독
소개
누군가를 진정으로 미워하려면 먼저 사랑했어야 한다고들 말한다. 이것이 일반 원칙으로 참인지는 모르겠지만, 적어도 NumPy와 나의 관계는 정확히 그렇다.
그런데 NumPy는 Python에서 배열 연산을 수행하는 소프트웨어다. 엄청나게 널리 쓰이며 PyTorch 같은 인기 기계 학습 라이브러리 모두에 막대한 영향을 끼쳤다. 이 라이브러리들도 아래에서 논의할 문제 대부분을 공유하지만, 구체성을 위해 NumPy에만 집중하겠다.
NumPy는 쉬운 일을 쉽게 만든다. A가 5×5 행렬이고 x가 길이 5인 벡터일 때, Ay=x를 만족하는 벡터 _y_를 찾고 싶다고 하자. NumPy에서는 이렇게 하면 된다.
y = np.linalg.solve(A, x)
정말 우아하다! 정말 명확하다!
하지만 상황이 아주 조금만 더 복잡해졌다고 하자. A가 100×5×5 배열로 주어진 5×5 행렬 100개의 스택이라고 하자. 그리고 x가 100×5 배열로 주어진 길이 5 벡터 100개의 스택이라고 하자. 또 모든 1≤i≤100에 대해 Aᵢyᵢ=xᵢ를 풀고 싶다고 하자.
루프를 쓸 수 있다면 쉬울 것이다.
y = np.empty_like(x)
for i in range(100):
y[i,:] = np.linalg.solve(A[i,:,:], x[i,:])
하지만 루프를 쓸 수 없다. 어느 정도는 Python에서 루프가 느리다는 한계 때문이다. 그러나 요즘은 모든 것이 GPU이고 큰 배열이 있다면 어느 언어에서든 루프를 쓰고 싶지는 않을 것이다. 그 모든 트랜지스터를 작동시키려면, 배열을 수많은 작은 조각으로 나누고 병렬 처리하는 특수 GPU 함수를 호출해야 한다.
좋은 소식은 NumPy가 그런 특수 루틴을 안다는 것이다(JAX나 CuPy를 쓴다면 적어도 그렇다). 그리고 np.linalg.solve를 올바르게 호출하면 이를 사용할 것이다.
나쁜 소식은 아무도 그 방법을 모른다는 것이다.
못 믿겠는가? 좋다, 다음 중 어느 것이 맞는가?
y = linalg.solve(A,x)
y = linalg.solve(A,x,axis=0)
y = linalg.solve(A,x,axes=[[1,2],1])
y = linalg.solve(A.T, x.T)
y = linalg.solve(A.T, x).T
y = linalg.solve(A, x[None,:,:])
y = linalg.solve(A,x[:,:,None])
y = linalg.solve(A,x[:,:,None])[:,:,0]
y = linalg.solve(A[:,:,:,None],x[:,None,None,:])
y = linalg.solve(A.transpose([1,2,0]),x[:,:,None]).T
아무도 모른다. 그리고 다른 것도 보여주겠다. 여기 문서가 있다.
읽어 보라. 곰곰이 생각해 보라. 이제 주목하라. 여전히 모든 i에 대해 한 번에 Aᵢyᵢ=xᵢ를 푸는 방법을 모른다. 그게 가능한 일이기는 한가? 내가 가능하다고 말한 것은 거짓말이었나?
내가 알기로 사람들이 실제로 하는 일은 되는 듯 보이는 무작위 변형을 시도하는 것이다.
NumPy는 배열에 연산을 적용하는 것이 전부다. 배열의 차원이 2 이하라면 모든 것이 괜찮다. 그러나 조금이라도 복잡한 일을 하면, 배열 A의 몇몇 차원, 배열 B의 다른 차원, 배열 C의 또 다른 차원에 적용하고 싶은 연산을 피할 수 없이 만나게 된다. 그리고 NumPy에는 이를 표현하는 이론이 없다.
내 뜻을 보여주겠다. 다음을 가정하자.
A는 K×L×M 배열이다.B는 L×N 배열이다.C는 K×M 배열이다.그리고 각 k와 n에 대해 L 및 M 차원에 걸친 평균을 계산하고 싶다고 하자. 즉, 다음을 원한다.
Dkn = 1/(LM) × ∑lm Aklm Bln Ckm.
이를 하는 방법은 두 가지다. 첫 번째는 흉측한 차원 정렬 기법을 쓰는 것이다.
D = np.mean(
np.mean(
A[:,:,:,None] *
B[None,:,None,:] *
C[:,None,:,None],
axis=1),
axis=1)
이게 대체 뭐냐고 묻게 될 것이다. 왜 None이 사방에 있는가? NumPy에서 배열을 인덱싱할 때는 None을 써서 새 차원을 삽입할 수 있다. A는 K×L×M이지만 A[:,:,:,None]은 K×L×M×1이다. 마찬가지로 B[None,:,None,:]는 1×L×1×N이고, C[:,None,:,None]은 K×1×M×**1**이다. 이것들을 곱하면 NumPy는 크기 1인 모든 차원을 “브로드캐스트”하여 K×L×M×N 배열을 만든다. 그런 다음 np.mean 호출이 L과 M 차원에 대해 평균을 낸다.
나는 이것이 나쁘다고 생각한다. NumPy를 수년간 써 왔지만, 이런 코드를 항상 실수 없이 작성하는 것은 여전히 불가능하다고 느낀다.
읽기도 거의 불가능하다. 이를 증명하기 위해 방금 동전을 던져서, 앞면이 아닌 경우에만 위에 버그를 넣었다. 버그가 있는가? 확실한가? 아무도 모른다.
두 번째 선택지는 필사적으로 영리해지려 애쓰는 것이다. 삶은 짧고 소중하지만, NumPy 문서를 읽으며 많은 시간을 보내면 언젠가는 np.tensordot이라는 함수가 있고 그것이 작업 대부분을 하게 만들 수 있다는 사실을 깨달을지도 모른다.
D = (1/L) * np.mean(
np.tensordot(A, B, axes=[1,0]) *
C[:,:,None],
axis=1)
이것은 맞다. (장담한다.) 하지만 왜 작동하는가? 정확히 np.tensordot은 무엇을 하는가? 다른 맥락에서 저 코드를 봤다면, 무슨 일이 벌어지는지 조금이라도 알 수 있겠는가?
루프를 쓸 수만 있다면 나는 이렇게 했을 것이다.
D = np.zeros((K,N))
for k in range(K):
for n in range(N):
a = A[k,:,:]
b = B[:,n]
c = C[k,:]
assert a.shape == (L,M)
assert b.shape == (L,)
assert c.shape == (M,)
D[k,n] = np.mean(a * b[:,None] * c[None,:])
NumPy를 너무 많이 작성한 사람들은 이것이 투박하다고 느낄지도 모른다. 나는 그것이 약간의 스톡홀름 증후군이라고 의심한다. 하지만 적어도 _명확하다_는 데는 동의할 수 있을 것이다.
실전에서는 상황이 흔히 더 나쁘다. A의 형상이 K×L×M이 아니라 M×K×L이었다고 하자. 루프를 쓴다면 별일 아니다. 그러나 NumPy에서는 A.transpose([1,2,0]) 같은 괴물을 작성해야 한다. 아니면 A.transpose([2,0,1])이어야 하나? 그것들은 어떤 형상을 만드는가? 아무도 모른다.
루프가 더 나았다.
세 번째 선택지가 있다.
D = 1/(L*M) * np.einsum('klm,ln,km->kn', A, B, C)
아인슈타인 합을 한 번도 본 적 없다면 무서워 보일 수도 있다. 하지만 우리의 목표를 기억하자.
D kn = 1/(LM) × ∑lm A klm B ln C km.
위 코드의 문자열은 기본적으로 세 입력 각각의 인덱스(klm,ln,km)와 출력의 목표 인덱스(->kn)에 라벨을 붙인다. 그러면 np.einsum이 입력의 대응 원소들을 곱하고 출력에 없는 모든 인덱스에 대해 합한다.
개인적으로 np.einsum은 실제로 좋은 NumPy의 드문 부분 중 하나라고 생각한다. 문자열이 조금 번거롭지만, 전체 함수는 이해하기 쉽고, 완전히 명시적이며, 상당히 일반적이고 강력하기 때문에 그만한 가치가 있다.
그런데 np.einsum은 어떻게 이 모든 것을 해내는가? 인덱스를 사용한다. 더 정확히 말하면 인덱스에 기반한 작은 _도메인 특화 언어_를 도입한다. NumPy의 일반 규칙을 따르기를 거부하기 때문에 NumPy의 설계 결함을 겪지 않는다.
하지만 np.einsum은 몇 가지 일만 한다. (Einops는 조금 더 한다.) 어떤 배열의 여러 차원에 걸쳐 다른 함수를 적용하고 싶다면 어떻게 하는가? np.linalg.einsolve는 없다. 그리고 자신만의 함수를 만든다면, _그 함수_의 “아인슈타인” 버전은 분명히 없다.
np.einsum이 좋다는 사실은 NumPy가 어딘가에서 잘못된 길로 갔음을 보여준다고 생각한다.
우리 주제와 비슷한 느낌이 드는 그림이 있다.
내가 배열 언어에 바라는 것은 이렇다. 문법에 특별히 까다롭지는 않지만, 다음과 같으면 좋겠다.
좋지 않은가? NumPy는 원죄 때문에 이를 달성하지 못한다고 생각한다. 인덱스를 없애고 브로드캐스팅으로 대체했다는 점이다. 그리고 브로드캐스팅은 인덱스를 대신할 수 없다.
NumPy의 핵심 비법은 브로드캐스팅이다. 다음 코드를 보자.
A = np.array([[1,2],[3,4],[5,6]])
B = np.array([10,20])
C = A * B
print(C)
출력은 다음과 같다.
[[ 10 40]
[ 30 80]
[ 50 120]]
여기서 A는 3×2 배열이고 B는 길이 2 배열이다. 둘을 곱하면 B는 A의 형상으로 “브로드캐스트”된다. 즉 A의 첫 번째 열은 B[0]=10과 곱해지고 두 번째 열은 B[1]=20과 곱해진다.
단순한 경우에는 좋아 보인다. 하지만 나는 좋아하지 않는다. 한 가지 이유는 위에서 보았듯이, 차원을 맞추기 위해 종종 지저분한 일을 해야 하기 때문이다.
또 다른 이유는 그것이 명시적이지도 읽기 쉽지도 않기 때문이다. A*B는 때로 원소별 곱셈을 하고 때로 더 복잡한 일을 한다. 따라서 A*B를 볼 때마다 브로드캐스팅 규칙의 어떤 경우가 발동하는지 알아내야 한다.
하지만 브로드캐스팅의 진짜 문제는 그것이 다른 모든 것에 퍼져 나가는 방식이다. 아래에서 설명하겠다.
수수께끼가 있다. 다음 코드를 보라.
A = np.ones((10,20,30,40))
i = np.array([1,2,3])
j = np.array([[0],[1]])
B = A[:,i,j,:]
B의 형상은 무엇인가?
정답은 10×2×3×40이다. i와 j 인덱스가 2×3 형상으로 브로드캐스트된 뒤에 어쩌고저쩌고 중얼중얼하기 때문이다. 이것이 말이 된다고 스스로 설득해 보라.
끝났는가? 좋다, 이제 다음을 해 보라.
C = A[:,:,i,j]
D = A[:,i,:,j]
E = A[:,1:4,j,:]
이것들의 형상은 무엇인가?
C는 10×20×2×3이다. 위에서 B에 일어난 일을 생각하면 논리적으로 보인다.
D는 어떤가? 2×3×10×30이다. 이제 어떤 이유에서인지 2와 3이 맨 앞으로 간다?
E는 어떤가? Python의 “슬라이스”는 끝점을 제외하므로 1:4는 [1,2,3]과 같고, 이는 i와 같으니 E는 B와 같다. 하하하, 농담이다! E는 10×3×2×1×40이다.
그렇다, 실제로 이렇게 된다. 못 믿겠다면 직접 해 보라! NumPy가 왜 이렇게 하는지는 안다. NumPy 인덱싱이 작동하는 방식을 설명하는 5000단어짜리 문서를 내가 흡수했기 때문이다. 하지만 그 시간을 돌려받고 싶다.
재미로 여러 AI 모델에게 이 배열들의 형상이 무엇인지 알아내 보라고 했다. 결과는 다음과 같았다. 다음 질의를 사용했다.
다음 Python 코드를 보라.
A = np.ones((10,20,30,40))
i = np.array([1,2,3])
j = np.array([[0],[1]])
B = A[:,i,j,:]
C = A[:,:,i,j]
D = A[:,i,:,j]
E = A[:,1:4,j,:]B, C, D, E의 형상은 무엇인가?
Claude 3.7은 “확장 사고”를 사용했다. 다음은 틀린 출력 전부다.
| AI | B | C | D | E |
|---|---|---|---|---|
| GPT 4.1 | 10×2×3×30 | |||
| Grok 3 | 10×3×30×2 | 10×3×2×40 | ||
| Claude 3 Opus | 10×3×2×30 | 10×20×3×2 | 10×3×30×2 | 10×3×2×40 |
| Llama 4 Maverick | 10×3×30×2 | 10×3×2×40 | ||
| o3 | 10×2×3×30 | |||
| Claude 3.7 | 10×3×30×2 | 10×3×2×40 |
| AI | B | C | D | E |
|---|---|---|---|---|
| GPT 4.1 | ✔️ | ✔️ | X | ✔️ |
| Grok 3 | ✔️ | ✔️ | X | X |
| Claude 3 Opus | X | X | X | X |
| Llama 4 Maverick | ✔️ | ✔️ | X | X |
| o3 | ✔️ | ✔️ | X | ✔️ |
| Claude 3.7 | ✔️ | ✔️ | X | X |
| Gemini 2.5 Pro | ✔️ | ✔️ | ✔️ | ✔️ |
| DeepSeek R1 | ✔️ | ✔️ | ✔️ | ✔️ |
(DeepSeek의 사고 연쇄에는 “잠깐”이 76번 등장했다. 처음에는 모든 것을 맞혔지만, 다시 시도했을 때는 어쩐지 B, C, D를 전부 틀리고 E만 맞혔다.)
이건 미친 일이다. 기본 기능을 쓰는 데 미친 논리 퍼즐을 풀 필요가 있어서는 안 된다.
“좋아, 단순한 방식의 인덱싱만 쓰면 되겠네”라고 생각할 수도 있다. 좋다. 하지만 고급 인덱싱이 필요한 때도 있다. 그리고 단순한 일을 하더라도, 그 미친 경우들을 피하려면 여전히 조심해야 한다.
이는 다시 모든 것을 읽기 어렵게 만든다. 단순한 방식으로 인덱싱을 사용하는 코드를 읽는 것뿐이라 해도, 그것이 단순하다는 것을 어떻게 아는가? A[B,C]를 보면 거의 무엇이든 할 수 있다. 이해하려면 A, B, C의 형상을 기억하고 모든 경우를 따져 봐야 한다. 그리고 물론 A, B, C는 흔히 다른 코드가 만들어 낸 것이며, 그 코드도 또한 생각해야 한다…
왜 NumPy는 그렇게 혼란스러운 np.linalg.solve(A,B) 함수를 갖게 되었을까? 아마 처음에는 A가 2차원 배열이고 b가 1차원 또는 2차원 배열일 때 작동하도록 만들었을 것이다. A⁻¹b 또는 A⁻¹B라는 수학 표기와 똑같이 말이다.
여기까지는 좋다. 하지만 그다음 누군가 3차원 배열을 들고 왔을 것이다. 루프를 쓸 수 있다면 해법은 “루프와 함께 기존 함수를 사용하라”이다. 하지만 루프를 쓸 수 없다. 따라서 기본적으로 세 가지 선택지가 있었다.
axes 인수를 넣을 수 있었다. 어쩌면 solve(A,B,axes=[[1,2],1])라고 작성할 수 있었을 것이다.solve_matrix_vector는 한 가지를 하고 solve_tensor_matrix는 다른 것을 할 수 있었다.solve가 내부적으로 차원 정렬을 시도할 방식에 관한 임의의 선택 말이다. 그러면 그 규칙을 알아내고 따르는 일은 사용자의 문제가 된다.이 모든 선택지는 나쁘다. 서로 다른 경우의 수가 조합적으로 많다는 사실을 어느 것도 제대로 다룰 수 없기 때문이다. NumPy의 선택은 이랬다. 전부 다. 어떤 함수에는 axes 인수가 있다. 어떤 함수에는 이름이 다른 버전이 있다. 어떤 함수에는 규칙이 있다. 어떤 함수에는 규칙과 axes 인수가 모두 있다. 그리고 어떤 함수는 벡터화 버전을 아예 제공하지 않는다.
하지만 NumPy의 가장 큰 결함은 이것이다. 어떤 주어진 형상의 배열로 문제를 푸는 함수를 만들었다고 하자. 이제 그 함수를 더 큰 배열의 특정 차원에 어떻게 적용하는가? 답은 이렇다. 훨씬 더 복잡한 방식으로 함수를 처음부터 다시 작성한다. 프로그래밍의 기본 원칙은 추상화, 즉 단순한 문제를 해결하고 그 해결책을 더 복잡한 문제의 구성 요소로 사용하는 것이다. NumPy는 그렇게 하지 못하게 한다.
내가 말하는 바를 보여 줄 마지막 예시 하나다. NumPy에 대해 불평할 때마다 사람들은 현대 언어 모델의 핵심 기법인 셀프 어텐션 예시를 보고 싶어 한다. 좋다. 내가 “self-attention numpy”를 검색해서 찾은 227개 버전 전부보다 낫다고 겸손하게 제안하는 구현이 여기 있다.
# 당신의 친구 dynomight가 만든 셀프 어텐션
input_dim = 4
seq_len = 4
d_k = 5
d_v = input_dim
X = np.random.randn(seq_len, input_dim)
W_q = np.random.randn(input_dim, d_k)
W_k = np.random.randn(input_dim, d_k)
W_v = np.random.randn(input_dim, d_v)
def softmax(x, axis):
e_x = np.exp(x - np.max(x, axis=axis, keepdims=True))
return e_x / np.sum(e_x, axis=axis, keepdims=True)
def attention(X, W_q, W_k, W_v):
d_k = W_k.shape[1]
Q = X @ W_q
K = X @ W_k
V = X @ W_v
scores = Q @ K.T / np.sqrt(d_k)
attention_weights = softmax(scores, axis=-1)
return attention_weights @ V
이것은 괜찮다. axis 관련 부분이 약간 모호하지만, 뭐 어떤가.
하지만 언어 모델에 정말 필요한 것은 멀티헤드 어텐션이다. 어텐션을 병렬로 여러 번 수행한 뒤 결과를 합치는 것이라고 보면 된다. 이것을 어떻게 할까?
먼저 추상화를 쓸 수 있는 제정신인 세상에 살고 있다고 상상해 보자. 그러면 앞선 함수를 루프에서 호출하기만 하면 된다.
# 당신의 친구 dynomight가 만든 멀티헤드 셀프 어텐션
# 루프를 쓸 수만 있다면
n_head = 2
X = np.random.randn(seq_len, input_dim)
W_q = np.random.randn(n_head, input_dim, d_k)
W_k = np.random.randn(n_head, input_dim, d_k)
W_v = np.random.randn(n_head, input_dim, d_v)
W_o = np.random.randn(n_head, d_v, input_dim // n_head)
def multi_head_attention(X, W_q, W_k, W_v, W_o):
projected = []
for n in range(n_head):
output = attention(X,
W_q[n,:,:],
W_k[n,:,:],
W_v[n,:,:])
my_proj = output @ W_o[n,:,:]
projected.append(my_proj)
projected = np.array(projected)
output = []
for i in range(seq_len):
my_output = np.ravel(projected[:,i,:])
output.append(my_output)
return np.array(output)
멍청해 보이는가? 그렇다—고맙다! 영리함은 나쁘다.
하지만 우리는 제정신인 세상에 살지 않는다. 그래서 대신 이렇게 해야 한다.
# 당신의 친구 dynomight가 만든 멀티헤드 셀프 어텐션
# 전부 벡터화했고 어리둥절하다
def multi_head_attention(X, W_q, W_k, W_v, W_o):
d_k = W_k.shape[-1]
Q = np.einsum('si,hij->hsj', X, W_q)
K = np.einsum('si,hik->hsk', X, W_k)
V = np.einsum('si,hiv->hsv', X, W_v)
scores = Q @ K.transpose(0, 2, 1) / np.sqrt(d_k)
weights = softmax(scores, axis=-1)
output = weights @ V
projected = np.einsum('hsv,hvd->hsd', output, W_o)
return projected.transpose(1, 0, 2).reshape(
seq_len, input_dim)
하! 하하하하하!
분명히 해 두자면, 나는 NumPy가 “다른 모든 배열 언어를 제외하면 최악의 배열 언어”라고만 제안하는 것이다. 더 나은 것을 제안할 수 없다면 불평하는 데 무슨 의미가 있겠는가?
그런데 사실 더 나은 제안이 있다. 나는 모든 성능은 유지하면서 모든 날카로운 모서리를 없앤다고 생각하는 “더 나은” NumPy의 프로토타입을 만들었다. 이것은 짧은 동기 부여용 서론일 뿐이라고 생각했지만, 쓰기 시작한 뒤 악마가 나를 사로잡았고 이렇게 3000단어 뒤에 와 있다.
또한 자신의 광적인 논쟁과 건설적인 배열 언어 API 제안 사이에는 어느 정도 거리를 두는 편이 아마 현명할 것이다. 그래서 내 새 물건은 다음 글에서 다루겠다.
이 양식이나 substack 또는 RSS로 구독하라. 이 양식으로 오류를 수정해 달라. substack이나 hacker news, julialang, fortranlang, lemmy에서 댓글을 달라. DumPy: 멍청해도 괜찮은 NumPy · 컴퓨터수학
사전분포는 이전이 아니라 데이터 이후에 만들면 좋다 · 수학예측