Recommanded Free YOUTUBE Lecture: <% selectedImage[1] %>

Anatomy of an array

소개

이 문서를 읽기 위해서는 NumPy에 대한 기본적인 경험이 있어야 한다. 그래서 기본적 Numpy 배열의 기본구조 특히 메모리 레이아웃, 뷰, 복사(copy), 데이터타입등을 설명하려 한다.

dtype이 np.float32인 배열의 모든 값을 지우려는 간단한 예제를 생각해보자. 속도를 극대화하기 위해서는 어덯게 해야 할 까. 아래의 구문은 (적어도 NumPy에 익숙한 사람들에게는) 명확하지만, 우리가 원하는 건 속도다.
>>> timeit("Z.view(np.float16)[...] = 0", globals())
100 loops, best of 3: 2.72 msec per loop
>>> timeit("Z.view(np.int16)[...] = 0", globals())
100 loops, best of 3: 2.77 msec per loop
>>> timeit("Z.view(np.int32)[...] = 0", globals())
100 loops, best of 3: 1.29 msec per loop
>>> timeit("Z.view(np.float32)[...] = 0", globals())
100 loops, best of 3: 1.33 msec per loop
>>> timeit("Z.view(np.int64)[...] = 0", globals())
100 loops, best of 3: 874 usec per loop
>>> timeit("Z.view(np.float64)[...] = 0", globals())
100 loops, best of 3: 865 usec per loop
>>> timeit("Z.view(np.complex128)[...] = 0", globals())
100 loops, best of 3: 841 usec per loop
>>> timeit("Z.view(np.int8)[...] = 0", globals())
100 loops, best of 3: 630 usec per loop
흥미롭게도 np.float64와 같은 더 큰 데이터 유형으로 배열을 캐스팅함으로써 25%의 속도요소를 얻었다. 그러나 배열을 바이트배열(np.int8)로 했을 때 50%의 속도향상이 있었다. 이러한 속도향상의 이유는 NumPy의 내부 구현과 컴파일러 최적화의 차이일 것이다. 이 예제는 NumPy의 철학을 보여준다.
위 내용은 원문의 내용이고, 내가 테스트한 환경은 약간 달랐다.
>>> timeit("Z.view(np.float16)[...] = 0", globals())
100 loops, best of 3: 1.59 msec per loop
>>> timeit("Z.view(np.int16)[...] = 0", globals())
100 loops, best of 3: 1.59 msec per loop
>>> timeit("Z.view(np.int32)[...] = 0", globals())
100 loops, best of 3: 1.59 msec per loop
>>> timeit("Z.view(np.float32)[...] = 0", globals())
100 loops, best of 3: 1.59 msec per loop
>>> timeit("Z.view(np.int64)[...] = 0", globals())
100 loops, best of 3: 1.58 msec per loop
>>> timeit("Z.view(np.float64)[...] = 0", globals())
100 loops, best of 3: 1.58 msec per loop
>>> timeit("Z.view(np.complex128)[...] = 0", globals())
100 loops, best of 3: 1.89 msec per loop
>>> timeit("Z.view(np.int8)[...] = 0", globals())
100 loops, best of 3: 1.59 msec per loop
차이가 없다. 그동안 NumPy가 나름대로의 방향으로 최적화했기 때문인 것 같다.

메모리 레이아웃

ndarray는 NumPy의 다차원 배열객체로, N 개의 정수를 각 블록의 위치에 매핑하는 색인 체계를 가지고 있는 1차원 세그먼트로 구성이 된다. 다시 말해서 인덱스를 이용해서 접근 할 수 있는 연속적인 메모리블럭이다. N차원 배열을 1차원으로 환원하기 때문에 훨씬 빠르게 작동한다.

아래 코드는 9개의 int16 타입의 데이터 9개를 저장 할 수 있는, 3x3 2차원 배열을 만들었다.
Z=np.arange(9).reshape(3,3).astype(np.int16)
Z 객체의 정보를 살펴보자.
>>> print(Z.itemsize)
2
>>> print(Z.shape)
(3, 3)
>>> print(Z.ndim)
2
>>> print(Z)
[[0 1 2]
 [3 4 5]
 [6 7 8]]

2차원 이상의 배열을 탐색 할 때, 각 차원으로 이동 할 수 있는 정보를 담은 튜플이 필요하다. 이것을 strides라고 한다.

Z는 2차원배열이다. 이때 배열에 있는 데이터의 타입은 int16/2byte로 모두 동일하기 때문에, 1차원 세그먼트로 재 구성했을 때, 값을 빠르게 찾을 수 있다. 아이템의 크기가 2byte 이며, 하나의 배열이 3개의 아이템을 가지고 있다는 것만 알고 있으면 된다. 즉 이 경우 strides는 (6,2)가 된다.
  • 0번째 배열의 2번째 값의 위치 : (0*6)+(1*2) = 2
  • 1번째 배열의 3번째 값의 위치 : (1*6)+(2*2) = 10
직접 Z 객체에 대한 strides를 구해보자.
>>> strides = Z.shape[1]*Z.itemsize, Z.itemsize
>>> print(strides)
(6, 2)
>>> print(Z.strides)
(6, 2)
예상과 동일한 값이 나왔다. 이 정보를 이용해서 offset를 계산해서 특정 아이템에 접근하는 코드를 만들 수 있다.
offset_start = 0
for i in range(Z.ndim):
    offset_start += Z.strides[i] * index[i]
offset_end = offset_start + Z.itemsize
tobye 변환 방법을 이용해서 올바른 방법인지 테스트 할 수 있다.
>>> Z = np.arange(9).reshape(3, 3).astype(np.int16)
>>> index = 1, 1
>>> print(Z[index].tobytes())
b'\x04\x00'
>>> offset_start = 0
>>> for i in range(Z.ndim):
...     offset_start += Z.strides[i] * index[i]
>>> offset_end = offset_start + Z.itemsize
>>> print(Z.tobytes()[offset_start:offset_end]
b'\x04\x00'

이해하기 쉽게 그림으로 묘사했다.

 memory layout - 1

메모리의 모습은 아래와 같이 묘사 할 수 있다.

 Layout

이제 Z의 슬라이스(slice)를 가져오면 아래와 같은 배열 Z의 뷰(view)를 만들 수 있다.
V = Z[::2,::2]

이러한 뷰는 dtype과 strides를 이용해서 설정할 수 있다. strides와 dtype만으로 전체적인 뷰의 추론이 가능하기 때문이다.

Item layout

 Item layout

Flattened item layout

 Flttened item layout

Memory layout (C order, big endian)

 Memory layout c order

View and copies

뷰와 카피는 수치계산 최적화에 중요한 개념이다. 앞서 최적화와 관련된 내용을 어느 정도 다루기는 했지만, 그 보다 더 복잡한 문제들이 있다.

참고

numpy.arange

지정된 간격 내에서 균등하게 간격을 둔 값을 반환한다. 0부터 100까지 2간격으로 있는 값들을 출력해보자.
>>> np.arange(0,100,2)
array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32,
       34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66,
       68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98])

View and Copies

View와 Copy는 최적화와 관련된 중요한 옵션이다. 앞 장에서 간단히 살펴보긴 했는데, 그 이상의 복잡한 내용들이 있다.

Direct and indirect access
먼저 indexing과 fancy indexing의 차이점을 살펴보겠다. 배열이 리턴을 할 경우 처음에는 항상 view를 리턴하고 그 다음에는 복사(copy)를 반환한다. 첫번째를 수정하면 원래 배열이 수정되지만 두 번째 반환 배열을 수정 할 경우 (원본은 그대로이고)복사본이 수정되기 때문에 숙지를 하고 있어야 한다.
>>> Z = np.zeros(9)
>>> Z_view = Z[:3]
>>> Z_view[...] = 1
>>> print(Z)
[ 1.  1.  1.  0.  0.  0.  0.  0.  0.]
>>> Z = np.zeros(9)
>>> Z_copy = Z[[0,1,2]]
>>> Z_copy[...] = 1
>>> print(Z)
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
따라서 fancy indexing이 필요한 경우(특히 복잡한 연산이 필요한 작업)라면 fancy index의 복사본으로 작업하는게 나을 것이다.
>>> Z = np.zeros(9)
>>> index = [0,1,2]
>>> Z[index] = 1
>>> print(Z)
[ 1.  1.  1.  0.  0.  0.  0.  0.  0.]
만약 작업하려는 indexing이 뷰인지, 복사인지 확신하지 못하겠다면 base를 이용해서 뷰인지 확인 할 수 있다.
>>> Z = np.random.uniform(0,1,(5,5))
>>> Z1 = Z[:3,:]
>>> Z2 = Z[[0,1,2], :]
>>> print(np.allclose(Z1,Z2))
True
>>> print(Z1.base is Z)
True
>>> print(Z2.base is Z)
False
>>> print(Z2.base is None)
True
NumPy 함수들은 가능한 뷰를(ex. ravel) 반환하지만, 어떤 경우(ex. flatten) 복사를 반환하기도 한다. 주의해서 사용해야 한다.
>>> Z = np.zeros((5,5))
>>> Z.ravel().base is Z
True
>>> Z[::2,::2].ravel().base is Z
False
>>> Z.flatten().base is Z
False

Temporary copy
복사는 명시적으로 만들 수 있디만, 이반적으로는 암시적으로 중간 복사본(intermediate copies)을 만들어서 사용하는 경우가 많다. 배열을 이용한 산술연산을 하는 경우다.
>>> X = np.ones(10, dtype=np.int)
>>> Y = np.ones(10, dtype=np.int)
>>> A = 2*X + 2*Y
위 예제는 3개의 intermediate 배열을 만들었다. 2*X의 결과와 2*Y의 결과가 저장되며, 마지막으로 2*X + 2*Y의 결과가 반환됐다. 이 경우에는 배열이 충분히 작아서 큰 차이가 나지 않는다. 하지만 배열이 큰 경우 이러한 식은 신중하게 사용해야한다. intermediate 배열은 최종 결과만 중요하고 나중에 X, Y가 필요하지 않은 경우에 적당하다.

>>> X = np.ones(10, dtype=np.int)
>>> Y = np.ones(10, dtype=np.int)
>>> np.multiply(X, 2, out=X)
>>> np.multiply(Y, 2, out=Y)
>>> np.add(X, Y, out=X)
위 예제는 임시 배열을 만들지 않고 있다. 당연하지만 사본의 사용은 속도에 나쁜 영향을 미친다.
>>> X = np.ones(1000000000, dtype=np.int)
>>> Y = np.ones(1000000000, dtype=np.int)
>>> timeit("X = X + 2.0*Y", globals())
100 loops, best of 3: 3.61 ms per loop
>>> timeit("X = X + 2*Y", globals())
100 loops, best of 3: 3.47 ms per loop
>>> timeit("X += 2*Y", globals())
100 loops, best of 3: 2.79 ms per loop
>>> timeit("np.add(X, Y, out=X); np.add(X, Y, out=X)", globals())
1000 loops, best of 3: 1.57 ms per loop

마무리

연습문제를 푸는 것으로 이번 장을 마무리 하겠다. Z1, Z2 두개의 벡터가 있다. Z2가 Z1의 뷰인지 아닌지 확인해 보자.
>>> Z1 = np.arange(10)
>>> Z2 = Z1[1:-1:2]

 View or Copy

Z2가 Z1의 뷰인지 확인해 보자.
>>> print(Z2.base is Z1)
True

Z2가 Z1의 뷰라는 걸 확인 할 수 있다. 즉 Z2는 Z1의 start:stop:step으로 표현할 수 있다.