DumPy는 루프와 인덱스 문법을 되살리면서도 GPU에서 벡터화된 배열 연산을 수행하도록 제안하는 NumPy 대안입니다.
최고
주제
팔로우
소개
내가 배열 언어에 바라는 것은 다음과 같다.
NumPy는 이 중 세 가지를 놓친다고 생각한다. 그래서 나는 불필요한 생각의 90%를 없애면서도 기능 손실은 없는 “수정안”을 제안하고 싶다. 이는 모든 머신러닝 라이브러리처럼 NumPy에 기반한 모든 것도 고칠 수 있다.
거창하게 들린다는 것을 안다. 아마도 당신은 늙은 dynomight가 마침내 제정신을 잃었다고 생각할 수 있다. 그러니 지금 경고하겠다. 내 해결책은 완전히 영리하지 않다. 여기서 영리한 점이 있다면, 영리함을 한결같이 거부한다는 점이다.
수정안의 동기를 설명하기 위해 NumPy가 어떻게 잘못되었는지에 대한 내 이야기를 해 보겠다. NumPy는 배열 연산과 선형대수를 위한 작고 좋은 라이브러리로 시작했다. 모든 것이 2차원 이하라면 훌륭하다. 하지만 어느 시점에 누군가 고차원 배열을 들고 나타났다. Python에서 루프가 빨랐다면 NumPy는 이렇게 말했을 것이다. “3차원 이상을 가진 분 안녕하세요, 제가 깔끔하고 단순하게 유지될 수 있도록 루프에서 제 2차원 이하 함수를 호출해 주세요, xox, NumPy.”
하지만 루프가 느렸기 때문에 NumPy는 보통 루프로 처리했을 모든 복잡성을 개별 함수 안으로 밀어 넣었다. 나는 이것이 재앙이었다고 생각한다. np.func(A,B) 같은 함수 호출을 볼 때마다 다음을 생각해야 하기 때문이다.
np.func는 그 형상을 보면 무엇을 하는가?함수마다 규칙이 다르다. 때로는 당혹스럽다. 즉, 계속 생각하고 특정 함수의 변덕을 달래기 위해 차원을 이리저리 옮겨야 한다. 함수들 이 당신의 변덕을 달래야 한다!
A*B나 A[B,C]처럼 단순해 보이는 것조차 시작 형상에 따라 아주 다른 일을 한다. 그리고 그 시작 형상은 흔히 이전 함수의 출력이므로 복잡성이 소용돌이친다.
무엇보다 나쁜 점은, 새로운 2차원 이하 함수를 작성하면 고차원 배열이 당신의 문제가 된다는 것이다. 규칙을 무엇으로 할지 당신이 결정해야 하고, 그러고 나서 함수를 훨씬 더 복잡한 방식으로 다시 작성해야 한다—
뒤쪽에서 들리는 목소리: Python은 형편없어! 진짜 언어를 썼다면 루프가 빨랐을 거야! 이 문제는 멍청해!
10년 전에는 강력한 주장이다. 하지만 이제는 모든 것이 GPU이고, GPU는 루프를 싫어한다. 오늘날 배열 패키지는 Python 같은 모습이지만 실제로는 모든 것을 배열 전체에서 병렬로 실행되는 특수 GPU 명령으로 몰래 컴파일하는 내장 언어인 쾌활한 인터페이스다. 큰 배열에는 GPU가 필요하다. 따라서 이제 호스트 언어의 속도는 그다지 중요하지 않다고 생각한다.
역설적으로 Python의 느림은 GPU가 장악하기도 전에 모든 것을 루프 없이 작동하도록 설계하게 강제했으므로 장점 이 되었을지도 모른다.
그래도 생각하는 것은 나쁘고 NumPy는 나를 생각하게 하므로, 나는 NumPy를 좋아하지 않는다.
여기 내 지극히 영리하지 않은 생각이 있다. 루프가 더 좋았다는 것을 그냥 인정하자. 고차원에서는 아직 누구도 루프와 인덱스를 이기는 표기법을 만들어 내지 못했다. 그러니 이렇게 하자.
그것이 기본적으로 전부다. 이 세 가지 항목을 받아들이면 아래에서 내가 하는 모든 일을 아마 다시 도출할 수 있을 것이다. 영리하지 않다고 했잖는가.
X와 Y가 2차원 배열이고 A가 4차원 배열이라고 하자. 그리고 Zij = (Yj)T (Aij)-1 Xi인 2차원 배열 Z를 찾고 싶다고 하자. 루프를 작성할 수 있다면 쉽다.
import numpy as np
Z = np.empty((X.shape[0], Y.shape[0]))
for i in range(X.shape[0]):
for j in range(Y.shape[0]):
Z[i,j] = Y[j] @ np.linalg.solve(A[i,j], X[i])
예쁘지는 않다. 짧지도 빠르지도 않다. 하지만 쉽다!
한편 NumPy에서 이를 효율적으로 하려면 어떻게 해야 할까? 이렇게 한다.
AiX = np.linalg.solve(A.transpose(1,0,2,3),
X[None,...,None])[...,0]
Z = np.sum(AiX * Y[:,None], axis=-1).T
NumPy 오타쿠가 아니라면 이것이 아웃사이더 아트처럼 보일 수 있다. 안심하라. 방금 내가 썼지만 내게도 그렇게 보인다. 왜 이토록 혼란스러운가? 높은 수준에서 보면 np.linalg.solve, np.sum, 그리고 곱셈(*)에는 복잡한 규칙이 있고 이 특정 문제를 잘 풀도록 함께 설계된 것이 아니기 때문이다. 문제는 무한히 많으므로 그것은 불가능했을 것이다. 그러므로 그 함수들을 만족시키려면 배열을 아주 많이 으깨고 옮겨야 한다.
더 말할 것도 없이, DumPy(겉보기에는 Dynomight NumPy)로 이 문제를 푸는 방법은 다음과 같다.
import dumpy as dp
A = dp.Array(A)
X = dp.Array(X)
Y = dp.Array(Y)
Z = dp.Slot()
Z['i','j'] = Y['j',:] @ dp.linalg.solve(A['i','j',:,:], X['i',:])
그렇다! 원한다면 다음 동등한 문법도 사용할 수 있다.
Z = dp.Slot()
with dp.Range(X.shape[0]) as i:
with dp.Range(Y.shape[0]) as j:
Z[i,j] = Y[j,:] @ dp.linalg.solve(A[i,j,:,:], X[i,:])
둘 다 완전히 벡터화되어 있다. 뒤에서 루프가 실행되지 않는다. GPU가 있다면 GPU에서 실행된다.
마법처럼 보이지만 실제 작동 방식은 꽤 단순하다.
DumPy 배열을 문자열(또는 dp.Range 객체)로 인덱싱하면, 차원이 더 적은 척하는 특별한 “매핑된” 배열을 만든다.
DumPy 함수(예: dp.linalg.solve 또는 @로 호출하는 dp.matmul)가 호출되면, 인수 중 매핑된 차원을 가진 것이 있는지 확인한다. 있다면 레이블을 공유하는 매핑 차원을 맞추어 계산을 자동으로 벡터화한다.
매핑 차원을 가진 배열을 dp.Slot에 할당하면 지정한 위치로 “매핑 해제”한다.
사악한 메타프로그래밍 추상 구문 트리 매크로 바이트코드 가로채기는 필요 없다. 다음 코드를 실행하면,
Z = dp.Slot()
Z['i','j'] = Y['j',:] @ dp.linalg.solve(A['i','j',:,:], X['i',:])
뒤에서는 다음이 일어난다.
a = A.map_axes([0, 1], ['i', 'j'])
x = X.map_axes([0], ['i'])
y = Y.map_axes([0], ['j'])
z = y @ dp.linalg.solve(a, x)
Z = z.unmap('i','j')
(수많은 assert와 주석이 들어간 버전은 여기를 클릭.)
# 먼저 A를 매핑한다
a = A.map_axes([0, 1], ['i', 'j'])
assert A.ndim == 4
assert a.ndim == 2 # 더 적은 차원을 가진 척한다
assert a.data.shape == A.shape # 비밀 매핑 데이터
assert a.axes == ('i', 'j', None, None) # 비밀 매핑 축
assert a.shape == (a.data.shape[2], a.data.shape[3])
# 형상은 매핑되지 않은 (None) 축으로 결정된다
# 이제 X를 매핑한다
x = X.map_axes([0], ['i'])
assert X.ndim == 2
assert x.ndim == 1
assert x.data.shape == X.shape
assert x.axes == ('i', None)
assert x.shape == (x.data.shape[1], )
# 이제 Y를 매핑한다
y = Y.map_axes([0], ['j'])
assert Y.ndim == 2
assert y.ndim == 1
assert y.shape == (Y.shape[1],)
assert y.axes == ('j', None)
assert y.data.shape == Y.shape
assert y.shape == (y.data.shape[1],)
# 실제 계산을 수행한다. 우연히도 'j' 차원이 먼저 저장되는데,
# (y에서 먼저 발견되기 때문이다.)
# 하지만 그것을 생각할 필요는 전혀 없다!
z = y @ dp.linalg.solve(a, x)
assert z.ndim == 0
assert z.shape == ()
assert z.axes == ('j','i')
assert z.data.shape == (Y.shape[0], X.shape[0])
# 매핑된 축을 매핑 해제한다
Z = z.unmap('i','j')
assert Z.ndim == 2
assert Z.shape == (X.shape[0], Y.shape[0])
어려운 부분을 건너뛴 것처럼 보일 수 있다. dp.linalg.solve는 입력 차원의 어떤 조합에 대해서도 어떻게 벡터화할지 어떻게 아는가? DumPy에 포함되는 모든 함수마다 이를 해야 하는 것 아닌가? 어렵지 않은가?
어렵다. 하지만 jax.vmap이 이미 해냈다. 이것은 (JAX의 NumPy 버전으로) 정의한 함수를 입력 차원의 어떤 집합에 대해서도 벡터화한다. DumPy는 실제 벡터화를 모두 수행하기 위해 이에 의존한다. (vmap이 엉성하고 고장 난 편이 좋다면 PyTorch의 torch.vmap을 진심으로 추천한다.)
하지만 잠깐. vmap이 이미 있다면 왜 DumPy가 필요한가? 이유는 다음과 같다.
import jax
from jax import numpy as jnp
Z = jax.vmap(
jax.vmap(
lambda x, y, a: y @ jnp.linalg.solve(a, x),
in_axes=[None, 0, 0]
),
in_axes=[0, None, 0]
)(X, Y, A)
이것이 vmap으로 같은 문제를 푸는 방법이다. (그리고 기본적으로 DumPy가 뒤에서 하는 일이다.)
나는 vmap이 NumPy 생태계의 최고 요소 중 하나라고 생각한다. 위 코드는 기본 NumPy 버전보다 진정으로 더 나아 보인다. 하지만 여전히 생각이 많이 필요하다! 왜 내부 vmap에는 in_axes=[None, 0, 0]을 넣고 외부 것에는 in_axes=[0, None, 0]을 넣는가? A의 두 번째 차원에 대해 벡터화해야 하는데도 왜 모든 축이 0인가? 답은 있지만 생각이 필요하다. 루프와 인덱스가 더 낫다.
좋다, 한 가지는 조금 영리하게 했다. Hij = 1/(1+i+j)인 Hilbert 행렬을 만들고 싶다고 하자. 기본 NumPy에서는 이렇게 해야 한다.
X = 1 / (1 + np.arange(5)[:,None] + np.arange(5)[None,:]) # 으응?
DumPy에서는 그냥 이렇게 쓸 수 있다.
X = dp.Slot()
with dp.Range(5) as i:
with dp.Range(5) as j:
X[i,j] = 1 / (i + j + 1)
그렇다! 작동한다! dp.Range가 문자열처럼도, 그 문자열을 따라 매핑된 배열처럼도 작동하기 때문이다. 따라서 위 코드는 대략 다음과 동등하다.
I = dp.Array([0,1,2,3,4])
J = dp.Array([0,1,2,3,4])
X = dp.Slot()
X['i','j'] = 1 / (1 + I['i'] + J['j'])
보라. 여전히 마법은 없다.
실제로 dp.Range는 무작위 문자열을 선택한다. (클래스는 충돌을 막기 위해 활성 범위 스택을 유지한다.) 따라서 더 자세히 말하면 위 코드는 다음과 비슷해진다.
I = dp.Array([0,1,2,3,4])
J = dp.Array([0,1,2,3,4])
i = I.map_axes([0],'range_EZaW')
j = J.map_axes([0],'range_ailw')
x = 1 / (1 + i + j) # 벡터화됨
X = x.unmap('range_EZaW','range_ailw')
DumPy가 실제로 더 나은지 시험하기 위해, 복잡도가 증가하는 여섯 문제를 골라 루프, NumPy, JAX(vmap 사용), DumPy로 각각 구현했다.
이 예시들에서는 항상 입력 배열이 사용 중인 시스템의 클래스라고 가정한다. 직접 실행하려면 np.array / jnp.array / dp.Array를 사용한 변환을 조금 추가해야 한다.
# 루프
H = np.empty((N, N))
for i in range(N):
for j in range(N):
H[i, j] = 1 / (i + j + 1)
# NumPy
i = np.arange(N)
j = np.arange(N)
H = 1 / (i[:, None] + j[None, :] + 1)
# JAX
indices = jnp.arange(N)
H = jax.vmap(
jax.vmap(
lambda i, j: 1 / (i + j + 1),
[0, None]),
[None, 0]
)(indices, indices)
# DumPy
H = dp.Slot()
with dp.Range(N) as i:
with dp.Range(N) as j:
H[i, j] = 1 / (i + j + 1) # 그렇다! 작동한다!
# 루프
C = np.zeros((X.shape[0],X.shape[1],X.shape[1]))
for n in range(X.shape[0]):
C[n] = np.cov(X[n])
# NumPy
mu = np.mean(X, axis=2)[:, :, None] # 으응?
C = np.sum((X - mu)[:, None, :, :] *
(X - mu)[:, :, None, :],
axis=3) / (X.shape[2] - 1) # 으으응??
# JAX
C = jax.vmap(jnp.cov)(X)
# DumPy
C = dp.Slot()
with dp.Range(X.shape[0]) as n:
C[n, :, :] = dp.cov(X[n, :, :])
# DumPy (대안)
C = dp.Slot()
C['n',:,:] = dp.cov(X['n',:,:])
(numpy.lib.stride_tricks은 존재하지 않는 척하자.)
# 루프
B = np.empty(N - window + 1)
for i in range(N - window + 1):
B[i] = np.mean(A[i:i + window])
# NumPy
i = np.arange(N - window + 1)[:, None]
j = np.arange(window)[None, :]
B = np.mean(A[i+j], axis=-1)
# JAX
idx = jnp.arange(window)
B = jax.vmap(
lambda i: jnp.mean(A[i+idx]),
)(jnp.arange(N - window + 1))
# DumPy
windowed = dp.Slot()
B = dp.Slot()
with dp.Range(N - window + 1) as i:
with dp.Range(window) as j:
windowed[i, j] = A[i + j]
B[i] = dp.mean(windowed[i, :])
# 참고: B[i] = dp.mean(A[i:i+window])는
# dp.Range를 슬라이스에 사용할 수 없으므로 작동하지 않는다
목표는 다음을 만족하는 E를 만드는 것이다.
E[i1, i2, :, i3] = A[B[i1], C[i1, i2], ::2, D[i2, i3]].
# 설정
K = 4
L = 5
M = 6
N = 7
A = np.random.randn(K, L, M, N)
B = np.random.randint(0, K, size=(9,))
C = np.random.randint(0, L, size=(9, 10))
D = np.random.randint(0, N, size=(10, 11))
# 루프
E = np.empty((9, 10, M // 2, 11))
for i1 in range(9):
for i2 in range(10):
for i3 in range(11):
E[i1,i2,:,i3] = A[B[i1],C[i1, i2],::2,D[i2, i3]]
# NumPy
E = A[B[:, None, None],
C[:, :, None],
::2,
D[None, :, :]
].transpose((0, 1, 3, 2))
# JAX
E = jax.vmap(
jax.vmap(
jax.vmap(
lambda b, c, d: A[b, c, ::2, d],
in_axes=[None,None,0]),
in_axes=[None,0,0]),
in_axes=[0,0,None]
)(B,C,D).transpose(0,1,3,2)
# DumPy
E = dp.Slot()
with dp.Range(9) as i1:
with dp.Range(10) as i2:
with dp.Range(11) as i3:
E[i1,i2,:,i3] = A[B[i1],C[i1, i2],::2,D[i2, i3]]
# DumPy (대안)
E = dp.Slot()
E['i1','i2',:,'i3'] = A[B['i1'], C['i1','i2'], ::2, D['i2','i3']]
이 문제의 목표는 벡터 목록, 가우시안 매개변수 목록, 그리고 각 벡터를 매개변수 목록에 매핑하는 배열이 주어졌을 때 각각의 대응하는 벡터/매개변수 조합을 평가하는 것이다. 형식적으로 2차원 X, B, C, means와 3차원 covs가 주어졌을 때 목표는 다음을 만족하는 A를 만드는 것이다.
Aij = log N( Xi | meansBij, covsCij).
# 설정
ndims = 3
ndata = 10
neval = 5
ndist = 7
X = np.random.randn(ndata, ndims)
B = np.random.randint(0, ndist, size=(ndata, neval))
C = np.random.randint(0, ndist, size=(ndata, neval))
means = np.random.randn(ndist, ndims)
scales = np.array(np.random.randn(ndist, ndims, ndims))
covs = np.array([scale @ scale.T for scale in scales])
# 루프
def log_prob(x, mean, cov):
diff = x - mean
y = np.linalg.solve(cov, diff)
quad = diff @ y
logdet = np.linalg.slogdet(2 * np.pi * cov)[1]
return -0.5 * (quad + logdet)
A = np.empty((ndata, neval))
for i in range(ndata):
for j in range(neval):
A[i, j] = log_prob(X[i, :],
means[B[i, j], :],
covs[C[i, j], :, :])
# NumPy
diff = X[:, None, :] - means[B]
y = np.linalg.solve(covs[C], diff[..., None])
quad = np.sum(diff * y[..., 0], axis=-1)
logdet = np.linalg.slogdet(2 * np.pi * covs[C])[1]
A = -0.5 * (quad + logdet)
# JAX
A = jax.vmap(
jax.vmap(
log_prob_gauss,
in_axes=[None, 0, 0]
),
)(X, means[B], covs[C])
# DumPy
def log_prob(x, mean, cov):
diff = x - mean
quad = diff @ dp.linalg.solve(cov, diff)
logdet = dp.linalg.slogdet(2 * jnp.pi * cov)[1]
A = dp.Slot()
with dp.Range(ndata) as i:
with dp.Range(neval) as j:
A[i, j] = log_prob(X[i,:],
means[B[i,j],:],
covs[C[i,j],:,:])
# DumPy (대안)
A = dp.Slot()
with dp.Range(ndata) as i:
with dp.Range(neval) as j:
mean = means[B[i,j],:]
cov = covs[C[i,j],:,:]
diff = X[i,:] - mean
quad = diff @ dp.linalg.solve(cov, diff)
logdet = dp.linalg.slogdet(2 * jnp.pi * cov)[1]
A[i,j] = -0.5 * (quad + logdet)
이전 글의 논의도 보라.
# 설정
input_dim = 4
seq_len = 4
d_k = 5
d_v = input_dim
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 softmax_numpy(x, axis=-1):
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):
Q = X @ W_q
K = X @ W_k
V = X @ W_v
scores = Q @ K.T / np.sqrt(d_k)
attention_weights = softmax_numpy(scores, axis=-1)
output = attention_weights @ V
return output
def multi_head_attention_loops(X, W_q, W_k, W_v, W_o):
projected = []
for n in range(n_head):
my_output = attention(X,
W_q[n, :, :],
W_k[n, :, :],
W_v[n, :, :])
my_proj = my_output @ W_o[n, :, :]
projected.append(my_proj)
projected = np.array(projected)
final = []
for i in range(seq_len):
my_final = np.ravel(projected[:, i, :])
final.append(my_final)
return np.array(final)
# NumPy
def softmax_numpy(x, axis=-1): # 반복
e_x = np.exp(x - np.max(x, axis=axis, keepdims=True))
return e_x / np.sum(e_x, axis=axis, keepdims=True)
def multi_head_attention_numpy(X, W_q, W_k, W_v, W_o):
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_numpy(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)
# JAX
def softmax_jax(x, axis=-1):
e_x = jnp.exp(x - jnp.max(x, axis=axis, keepdims=True))
return e_x / jnp.sum(e_x, axis=axis, keepdims=True)
def attention_jax(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 / jnp.sqrt(d_k)
attention_weights = softmax_jax(scores, axis=-1)
output = attention_weights @ V
return output
def multi_head_attention_jax(X, W_q, W_k, W_v, W_o):
def myfun(X, w_q, w_k, w_v, w_o):
return attention_jax(X, w_q, w_k, w_v) @ w_o
projected = jax.vmap(myfun,
in_axes=[None, 0, 0, 0, 0]
)(X, W_q, W_k, W_v, W_o)
return jax.vmap(jnp.ravel, in_axes=1)(projected)
# DumPy
def softmax_dumpy(x):
assert x.ndim == 1 # 차원을 생각할 필요가 없다!
e_x = dp.exp(x - dp.max(x))
return e_x / dp.sum(e_x)
@dp.wrap # Slots가 있는 함수를 자동 벡터화하려면 필요하다
def attention_dumpy(X, W_q, W_k, W_v):
Q = X @ W_q
K = X @ W_k
V = X @ W_v
scores = Q @ K.T / np.sqrt(d_k)
attention_weights = dp.Slot()
with dp.Range(seq_len) as i:
attention_weights[i, :] = softmax_dumpy(scores[i, :])
output = attention_weights @ V
return output
def multi_head_attention_dumpy(X, W_q, W_k, W_v, W_o):
output = dp.Slot()
projected = dp.Slot()
final = dp.Slot()
with dp.Range(n_head) as n:
output[n, :, :] = attention_dumpy(X,
W_q[n, :, :],
W_k[n, :, :],
W_v[n, :, :])
projected[n, :, :] = output[n, :, :] @ W_o[n, :, :]
with dp.Range(seq_len) as i:
final[i, :] = dp.ravel(projected[:, i, :])
return final
# DumPy (대안)
def multi_head_attention(X, W_q, W_k, W_v, W_o):
attn_weights = dp.Slot()
projected = dp.Slot()
final = dp.Slot()
with dp.Range(n_head) as n:
Q = X @ W_q[n, :, :]
K = X @ W_k[n, :, :]
V = X @ W_v[n, :, :]
scores = Q @ K.T / np.sqrt(d_k)
with dp.Range(seq_len) as i:
attn_weights[n, i, :] = softmax_dumpy(scores[i, :])
projected[n, :, :] = attn_weights[n, :, :] @ V @ W_o[n, :, :]
with dp.Range(seq_len) as i:
final[i, :] = dp.ravel(projected[:, i, :])
return final
각 구현에 1~10 척도의 주관적 “좋음” 점수를 매겼다. 각 문제에서 최고 구현에는 항상 10점을 주었고, 필요한 생각의 양에 따라 다른 것들의 점수를 깎았다.
| 문제 | 루프 | Numpy | JAX (vmap) | DumPy |
|---|---|---|---|---|
| Hilbert 행렬 | 10 | 7 | 7 | 10 |
| 공분산 | 9 | 4 | 10 | 9 |
| 이동 평균 | 10 | 6 | 6 | 8 |
| 인덱싱 | 10 | 5 | 4 | 10 |
| 가우시안 | 10 | 3 | 6 | 10 |
| 어텐션 | 10 | 1 | 8 | 10 |
| 평균 | 9.8 | 4.3 | 6.8 | 9.5 |
이 의심스러운 방법론과 지어낸 숫자들에 따르면 DumPy는 루프의 96.93877%만큼 좋다! 아는 것이 힘이다! 하지만 진지하게 말하면, 주관적이기는 해도 내 점수가 너무 논쟁적이라고는 생각하지 않는다. 가장 논쟁적인 것은 아마 JAX의 어텐션 점수일 것이다.
DumPy가 NumPy에 추가하는 유일한 것은 인덱스를 위한 멋진 표기법이다. 그게 전부다.
DumPy를 좋게 만드는 것은 많은 것을 제거한다는 점이기도 하다고 생각한다. 대략적으로 루프가 없기 때문에 존재하는, 혼란스러운 모든 것을 제거하려 했다. 정확히 올바른 위치에 선을 그었는지는 모르겠지만, 뭔가를 제거하는 방향이 옳다는 데는 자신이 있다.
NumPy에서 A*B는 A와 B가 모두 스칼라이면 작동한다. 또는 A가 5×1×6이고 B가 5×1×6×1이어도 된다. 하지만 A가 1×5×6이고 B가 1×5×6×1이면 안 된다. 응?
사실 곱셈 같은 스칼라 연산의 브로드캐스팅 규칙은 그렇게 복잡하지는 않다. 그래도 나는 좋아하지 않는다. A*B를 볼 때마다 그 형상이 무엇인지와 계산이 무엇을 할지 걱정해야 하기 때문이다.
그래서 제거했다. DumPy에서는 A 또는 B 중 하나가 스칼라이거나 A와 B가 정확히 같은 형상일 때만 A*B를 할 수 있다. 그게 전부고, 나머지는 오류를 낸다. 대신 인덱스를 사용하여 무엇을 하는지 명확히 하라. 다음 대신,
C = A[...,None] * B[None]
이렇게 써라.
C['i','j','k'] = A['i','j'] * B['j','k']
NumPy의 인덱싱은 터무니없이 복잡하다. A[B,C,D]를 쓰면 모든 형상에 따라 많은 다른 일을 할 수 있다. DumPy에서는 스칼라 인덱스만 허용하고 싶기도 했다. 그래도 루프로 고급 작업을 할 수 있으므로 그렇게 나쁘지는 않았을 것이다. 하지만 A와 B가 단순한 1차원 배열일 때 A[B]를 쓰지 못하면 꽤 짜증 난다.
그래서 잠정적으로 좀 더 실용적으로 결정했다. DumPy에서는 정수, 슬라이스, 또는 (매핑될 수도 있는) Array로 인덱싱할 수 있다. 하지만 스칼라가 아닌 Array 인덱스는 하나만 가능하다. 생각이 필요 없는 가장 일반적인 문법이기 때문에 이렇게 정했다.
무슨 뜻인지 보여 주겠다. 다음을 보면,
A[1, 1:6, C, 2:10] # numpy와 dumpy 모두에서 합법
출력 형상이 무엇인지 “분명하다”. (1:6의 형상, 그다음 C의 형상, 그다음 2:10의 형상.) 충분히 단순하다. 하지만 다음처럼 다차원 배열 입력이 두 개 있으면,
A[B, 1:6, C, 2:10] # numpy에서는 합법, dumpy에서는 금지
갑자기 모든 것이 난장판이 된다. B와 C 사이의 브로드캐스팅, 직교 인덱스와 점별 인덱스, 배열과 다르게 동작하는 슬라이스, 출력 차원이 어디로 가는지에 관한 특이점을 생각해야 한다. 그래서 DumPy는 이를 금지한다. 대신 다음 중 하나를 써야 한다.
D['i',:,:] = A[B['i'], 1:6, C['i'], 2:10] # (1)
D[:,:,'i'] = A[B['i'], 1:6, C['i'], 2:10] # (2)
D['i','j',:,:] = A[B['i'], 1:6, C['j'], 2:10] # (3)
D['i','j',:,:] = A[B['i','j'], 1:6, C['i'], 2:10] # (4)
D['i','j',:,:] = A[B['i','j'], 1:6, C['i','j'], 2:10] # (5)
모두 보이는 그대로 정확히 작동한다.
아, 그리고 한 가지 더! DumPy에서는 모든 차원을 인덱싱해야 한다. NumPy에서 A가 3차원이면 A[2]는 A[2,:,:]와 동등하다. 때로는 편리하지만, A[2]를 볼 때마다 A가 몇 차원인지 걱정해야 한다는 뜻이다. DumPy에서는 배열을 인덱싱하거나 dp.Slot에 할당할 때마다 모든 인덱스가 포함되었는지 검사한다. 따라서 위의 옵션 (4)를 보면 다음을 알 수 있다.
A는 4차원이다.B는 2차원이다.C는 1차원이다.D는 4차원이다.언제나, 언제나, 언제나. 경우의 수 없고, 생각도 없다.
다시 말하지만, 많은 NumPy 함수에는 벡터화를 위한 복잡한 관례가 있다. np.linalg.solve는 대략 이렇게 말한다. “입력이 2차원 이하이면 자명한 일을 하라. 그렇지 않으면 극도로 혼란스러운 브로드캐스팅 일을 하라.” DumPy는 혼란스러운 브로드캐스팅 일을 제거한다. dp.linalg.solve(A,B)를 보면 A와 B가 2차원 이하이고, 따라서 까다로운 일이 일어나지 않는다는 것을 안다.
마찬가지로 NumPy에서 A @ B는 np.matmul(A,B)``와 동등하다. 두 입력이 2차원 이하일 때 이는 “자명한 일”을 한다. (내적 또는 어떤 종류의 행렬/벡터 곱셈이다.) 그렇지 않으면 브로드캐스팅하거나 벡터화하거나 뭔가를 한다? 나는 절대 기억할 수 없다. DumPy에서는 A @ B`를 1차원 또는 2차원 배열로만 제한하므로 그런 문제가 없다.
더 많은 차원이 필요하다면 문제없다. 인덱스를 사용하라.
기능을 제거하는 것은 성가롭게 보일 수 있지만, 말하겠다. 그냥 해 보라. 이런 방식으로 프로그래밍하면 가능한 오류의 부류가 하나씩 사라지면서 멋진 평온함이 찾아온다.
다르게 말하면, 왜 멋진 기능을 선택 사항으로 남기지 않는가? 선택 사항은 생각을 뜻하기 때문이다! 나는 단순한 방식으로 프로그래밍하고 싶다. 실수로 혼란스러운 브로드캐스팅 광기를 촉발하지는 않을지 걱정하고 싶지 않다. 그것은 실수일 것이기 때문이다. 컴퓨터가 실수를 잡아주길 원하지, 의도하지 않은 이상한 일을 조용히 하기를 원하지 않는다.
원칙적으로 모든 혼란스러운 배치 처리를 보존하는 evil_solve 메서드가 있어도 괜찮다. 정말로 원한다면 직접 만들 수 있다.
evil_solve = dp.MappedFunction(jnp.linalg.solve) # 권장하지 않음
같은 래퍼를 사용하여 어떤 JAX NumPy 함수든 DumPy와 작동하도록 변환할 수 있다.
수학을 생각해 보라. 2차원 이하에서 좌표 없는 선형대수 표기법은 훌륭하다. 하지만 고차원 텐서에는 경우가 너무 많아서 대부분의 물리학자는 그냥 좌표를 사용한다.
그래서 이 해결책은 내게 아주 자명해 보인다. 솔직히 왜 아직 표준이 아닌지 약간 혼란스럽다. 내가 뭔가 놓치고 있는가?
내가 NumPy에 불평하면 많은 사람이 A, J, K, Q 같은 APL 계열 언어를 살펴보라고 제안한다. (한 글자 언어는 C, D, F, R, T, X 및 수많은 다른 것을 제외하면 모두 APL 비슷하다. 편리하지 않은가?) 이들의 명백한 단점은 다음과 같다.
그중 어느 것도 나를 괴롭히지 않는다. 언어가 더 낫다면, 우리는 그것들을 사용하는 법을 배우고 GPU에서 자동미분을 하게 만들어야 한다. 하지만 그것들이 더 낫다 고 확신하지는 못하겠다. 실제로 이 언어들을 배워 보면, 기호로 된 횡설수설은 기본적으로 앞에서 NumPy에서 본 것과 같은 종류의 차원 으깨기를 하는 것이라는 사실을 알게 된다.
AiX = np.linalg.solve(A.transpose(1,0,2,3),
X[None,...,None])[...,0]
Z = np.sum(AiX * Y[:,None], axis=-1).T
이유는 NumPy와 vmap처럼 이 언어들도 차원을 이름이 아니라 위치 로 정렬하기 때문이다. 차원을 반드시 으깨야 한다면 최고의 도구를 쓰고 싶다. 하지만 차원을 전혀 으깨지 않는 편이 더 좋다.
사람들은 xarray처럼 “이름 있는 차원을 가진 NumPy”도 흔히 제안한다. (PyTorch도 성의 없는 구현을 갖고 있다.) 물론 DumPy도 이름 있는 차원을 사용하지만, 결정적인 차이가 있다. xarray에서는 이름이 배열 자체의 일부인 반면 DumPy에서는 배열 바깥에 존재한다.
어떤 경우에는 영구적인 이름 있는 차원이 아주 좋다. 하지만 선형대수에서는 혼란스럽다. 예를 들어 A가 이름 있는 차원 `