[고성능 파이썬] 6.1 행렬과 벡터연산

다룰 내용

• 벡터 연산 병목 지점의 위치
• 계산 과정 중에 CPU를 효율적으로 사용하는지 알아볼 수 있는 도구 - perf
• 산술 계산에서 numpy가 순수 파이썬보다 더 빠른 이유 - 메모리 단편화와 메모리 지역성
• 캐시 미스와 페이지 폴트
• 코드 내에서 메모리 할당을 추적하는 방법

여러개의 연산을 한번에 수행할 수 있는 벡터 연산은 프로그램을 빠르게 실행하기 위한 필수 불가결한 요소이다.

이번장에서는 확산 방정식을 연산하는 동안에 CPU의 동작 과정을 이해함으로써 벡터 연산에 대해 알아보고 numpy가 어떻게 돌아가는지, numpy위에 구축된 pandas를 더 빠르게 실행할 수 있는 방법은 무엇인지 알아본다.

 

예제에서 사용한 확산 방정식

이부분은 사실 예제에서 사용하는 확산 방정식을 설명하는 부분으로 넘어가도 되는 부분이다.

(다 작성하고 빼도 될 것 같다는 생각을 했지만 노력이 아까워서 넣는다...)

 

확산 방정식은 아래와 같다.

$$\cfrac{\partial}{\partial t} u(x, t) = D \cdot \cfrac{\partial^2}{\partial x^2} u(x, t)$$

 

$u$는 확산하는 유체의 양을 나타내는 벡터이다. 예를 들어 이 값이 0이면 물만 존재하고, 1이면 물감만 존재하는 것이다.

(일반적으로 u는 3차원 공간에 든 유체를 나타내는 3차원 행렬이지만 이 식은 1차원에 해당하는 식이다. 먼저 1차원에 대한 식을 이해하고 예제에서는 이를 확장해 $x$와 $y$를 사용하는 2차원 식을 사용할 예정이다.)

 

$D$는 확산계수로, $D$값이 클수록 쉽게 확산한다. 여기서는 계산이 너무 복잡해지지 않게 1로 지정한다.

 

이제 미분방정식을 오일러 방법을 이용해 불연속적인 공간과 시간으로 근사값을 구한다.

$$\cfrac{\partial}{\partial t} u(x, t) \approx \cfrac{u(x, t + dt) - u(x, t)}{dt}$$

여기서 $dt$는 이 방정식을 풀 시간의 정밀도를 나타낸다. 더 정확한 값을 얻고 싶으면 $dt$를 더 작은 값으로 사용하면 되는 것이다.

 

이렇게 식을 변형하면 초기 상태 $u(x,0)$에서 미래 특정의 상태 $u(x, dt)$로 변하는 과정을 나타낼 수 있다.

이런 문제를 초기값 문제(initial value problem) 또는 코시 문제(Cauchy problem)라고 한다.

 

이제 미분방정식을 유한 차분 근사법을 이용해서 최종 방정식을 유도하자.

$$u(x, t + dt) = u(x, t) + dt * D * \cfrac {u(x + dx, t) + u(x - dx, t) - 2 \cdot u(x,t)}{dx^2}$$

 

$dt$가 초당 프레임 수를 나타낸다면 $dx$는 동영상의 해상도를 나타낸다고 할 수 있다. ($dx$ 값이 작아질수록 행렬의 셀 크기도 점점 작아지므로 더 정확한 값을 얻을 수 있다.)

여기서는 계산을 간단히 하기 위해 $dx = 1$로 정의한다.

 

그런데 몇가지 주의해야 할 사항이 있다.

먼저 $u$로 표시한 공간 색인($x$ 매개변수)은 행렬의 색인을 나타낸다.

그러면 $x$가 행렬의 시작점일 때 $x - dx$값은 어떻게 구해야 하는가? 이 문제를 경계 조건(boundary condition)이라 한다.

이런 경계 조건은 '행렬의 경계를 벗어나는 모든 값은 0으로 둔다'처럼 설정할 수 있다.

이 예제에서는 단순히 0으로 두지 않고 (i+1)%N 과 같은 주기적 경계 조건으로 설정한다.

 

시간 변화에 따른 $u$값을 어떻게 저장할지도 고려해야 한다.

현재 상태, 다음 상태를 나타내는 행렬이 최소 2개 필요한데 그러면 두 행렬을 계속 갱신해야 하므로 성능에 미치는 영향은 어마어마하다. (뒤에서 다룸)

 

# Create the initial conditions
u = vector of length N
for i in range(N):
    u = 0 if there is water, 1 if there is dye

# Evolve the initial conditions
D = 1
t = 0
dt = 0.0001
while True:
    print(f"Current time is: {t}")
    unew = vector of size N

    # Update step for every cell
    for i in range(N):
        unew[i] = u[i] + D * dt * (u[(i+1)%N] + u[(i-1)%N] - 2 * u[i])
    # Move the updated solution into u
    u = unew

    visualize(u)

위의 코드는 1차원 확산 의사코드로 while문 내의 식을 보면 우리가 구한 최종방정식을 사용하고 있음을 알 수 있다.

 

2차원 확산 방정식은 다음과 같이 나타낼 수 있다.

$$ \cfrac{\partial}{\partial t} u(x, y, t) = D \cdot \left( \cfrac {\partial^2}{\partial x^2} u(x, y, t) + \cfrac{\partial^2}{\partial y^2} u(x, y, t) \right) $$ 

 

이 식도 위의 1차원식을 유도한 것처럼 유도하면 다음과 같다.

 

$$u(x, y, t + dt) = u(x, y, t) + dt * D * \left( \cfrac {u(x + dx, y, t) + u(x - dx, y, t) - 2 \cdot u(x, y, t)}{dx^2} + \cfrac {u(x, y + dy, t) + u(x, y- dy, t) - 2 \cdot u(x, y, t)}{dy^2} \right)$$

 

이를 코드로 나타내면 다음과 같다.

($dx$, $dy$, $D$는 1이므로 생략했다.)

 

파이썬 리스트로 만든 확산 방정식 코드

for i in range(N):
    for j in range(M):
        unew[i][j] = u[i][j] + dt * (
            (u[(i + 1) % N][j] + u[(i - 1) % N][j] - 2 * u[i][j]) + # d^2 u / dx^2
            (u[i][(j + 1) % M] + u[i][(j - 1) % M] - 2 * u[i][j])   # d^2 u / dy^2
        )

 

메모리 단편화

파이썬은 기본적으로 CPU의 벡터 연산을 지원하지 않는다.

그 이유는 첫째로 파이썬 리스트의 원소들이 실제로는 메모리상에 순서대로 나열되어 있지 않고 여기저기 떨어져 있는 데이터들을 가리키는 포인터를 저장하기 때문이다.

둘째로는 파이썬 바이트 코드는 벡터 연산에 최적화되지 않았기 때문이다.

 

첫번째 이유에서의 메모리가 여기저기 떨어져있다는 것을 메모리 단편화라고 한다.

데이터가 단편화되었다면 전체 블록을 한 번에 읽는 대신 데이터 조각을 하나씩 읽어야 하기에 메모리 복사 오버헤드가 더 많이 발생하고 필요한 데이터를 모두 복사할 때까지 CPU가 대기해야 한다.

 

또한, CPU에 데이터가 필요할 때 이를 공급하는 과정에서 데이터가 처리되는 속도보다 공급되는 속도가 느린 폰 노이만 병목 현상이 일어난다. (이는 메모리와 CPU 사이의 대역폭이 제한적이어서 발생하는 현상이다.)

 

그렇기에 CPU는 분기 예측과 파이프라이닝을 통해 다음 연산에 필요한 데이터를 예측해서 캐시에 미리 올려놓는데, 이 예측이 틀렸을 경우를 캐시 미스라고 한다. 이런 캐시 미스는 perf라는 도구로 측정할 수 있다.

파이썬 리스트를 사용하는 확산 방정식 코드를 perf로 측정해보면 약 50% 이상의 캐시가 캐시 미스가 발생하는 것을 확인할 수 있다.

메모리가 단편화되어 있기 때문에 이만큼의 캐시 미스가 발생하는 것이다.

 

이렇게 파이썬의 리스트는 필요한 데이터를 한번에 CPU 캐시에 담아둘 수 없기 때문에 연산을 벡터화할 수 없다.

이런 문제를 해결하려면 리스트가 아닌 데이터를 연속된 메모리 공간에 저장하는 array 모듈을 사용해야 한다.

 

하지만 array를 사용해도 파이썬은 여전히 루프를 벡터화해 계산하는 법을 모르므로 이것만으로는 완벽히 해결할 수 없다.

또한 array는 아주 저수준의 형태로 수를 저장하므로 사용할 때 파이썬에 호환되는 상태로 변경한 뒤 사용해야 하는데, 이로 인해 array를 사용할 때 오버헤드가 발생해 list보다 더 오래걸린다.

 

고로 array는 벡터 연산의 용도보다는 데이터를 저장하고 전송하는데에 더 적합하다고 할 수 있다.

 

Numpy

범용적인 목적으로 만들어진 리스트와 달리 수 배열을 다룬다는 특수한 목적으로 정교하게 만들어진 numpy 배열은 데이터를 연속된 메모리 공간에 저장하고 이에 대한 벡터 연산도 지원한다.

 

numpy는 내부적으로 잘 최적화된 C코드를 통해 CPU에서 지원하는 벡터화 기능의 장점을 활용한다.

또한, numpy 배열은 array 객체와 마찬가지로 저수준 형태로 연속된 메모리 공간에 수를 저장한다.

 

다음은 numpy로 구현한 확산 방정식 코드이다.

from numpy import zeros, roll

def evolve(grid, dt, D=1):
    laplacian = roll(grid, +1, 0) +
                roll(grid, -1, 0) +
                roll(grid, +1, 1) +
                roll(grid, -1, 1) -
                4 * grid
                
    return grid + dt * D * laplacian

파이썬 리스트로 구현한 코드는 2번의 for루프를 통해 구현되었지만 numpy로 구현한 코드는 1번의 루프도 없이 내부적으로 동작하고 있음을 볼 수 있다.

 

그 이유는 리스트는 각 원소별로 연산을 하지만, numpy의 배열은 행렬과 같이 한번에 연산을 하기 때문이다.

 

perf를 이용해 두 코드가 CPU에서 어떻게 작동되는지 살펴보면 약 70배 정도 차이남을 볼 수 있다.

파이썬 리스트

$ perf stat -e cycles,instructions,\
    cache-references,cache-misses,branches,branch-misses,task-clock,faults,\
    minor-faults,cs,migrations python diffusion_python_memory.py


 Performance counter stats for 'python diffusion_python_memory.py':

   415,864,974,126      cycles                    #    2.889 GHz                      
 1,210,522,769,388      instructions              #    2.91  insn per cycle           
       656,345,027      cache-references          #    4.560 M/sec                    
       349,562,390      cache-misses              #   53.259 % of all cache refs      
   251,537,944,600      branches                  # 1747.583 M/sec                    
     1,970,031,461      branch-misses             #    0.78% of all branches          
     143934.730837      task-clock (msec)         #    1.000 CPUs utilized          
            12,791      faults                    #    0.089 K/sec                  
            12,791      minor-faults              #    0.089 K/sec                  
               117      cs                        #    0.001 K/sec                  
                 6      migrations                #    0.000 K/sec                  

     143.935522122 seconds time elapsed

 

numpy 배열

$ perf stat -e cycles,instructions,\
    cache-references,cache-misses,branches,branch-misses,task-clock,faults,\
    minor-faults,cs,migrations python diffusion_numpy.py


 Performance counter stats for 'python diffusion_numpy.py':

     8,432,416,866      cycles                    #    2.886 GHz                      
     7,114,758,602      instructions              #    0.84  insn per cycle           
     1,040,831,469      cache-references          #  356.176 M/sec                    
       216,490,683      cache-misses              #   20.800 % of all cache refs      
     1,252,928,847      branches                  #  428.756 M/sec                    
         8,174,531      branch-misses             #    0.65% of all branches          
       2922.239426      task-clock (msec)         #    1.285 CPUs utilized          
           403,282      faults                    #    0.138 M/sec                  
           403,282      minor-faults              #    0.138 M/sec                  
                96      cs                        #    0.033 K/sec                  
                 5      migrations                #    0.002 K/sec                  

       2.274377105 seconds time elapsed

 

 

마이너 페이지 폴트와 제자리 연산

numpy에 배열을 사용한 evolve 함수를 호출하면 할수록 새로운 numpy 배열이 계속 할당된다. 그 이유는 함수 내부에서 laplacian이라는 새로운 numpy 배열을 만들기 때문이다.

 

마이너 페이지 폴트

이때 발생하는 현상을 마이너 페이지 폴트라고 한다.

마이너 폴트는 프로그램이 메모리에 새로 할당된 공간에 처음 접근할 때 발생한다.

 

이 말이 무슨 말이냐면, 새로 만들어진 객체는 처음부터 메모리가 할당되는 것이 아니다.

커널은 메모리 주소만 먼저 주고 처음 사용될 때 그제서야 메모리를 할당하는 지연 할당을 사용한다.

 

따라서 새로 할당된 데이터에 처음 접근하면 프로그램 실행을 잠깐 중지한 후 실제로 요청받은 메모리 공간이 존재하는지 보고 프로그램에서 사용할 수 있는 참조를 생성한다.

 

이렇게 프로그램을 일시 중지하면 캐시에 있던 상태를 잃어버리고 명령어 파이프라이닝을 수행할 수 없게 된다.

결국 다시 메모리를 캐시로 가져와야 하므로, 엄청나게 많은 비용이 드는 것이다.

 

제자리 연산

그래서 새로운 배열을 만들지 않고 재활용할 배열을 만들어 놓으면 메모리도 아낄 수 있고, 메모리를 할당하는 데 필요한 시간도 절약할 수 있다.

 

이를 위해 사용하는 방법이 제자리 연산인데 이는 배열 A가 있을 때 이를 통해 연산을 한 후 그 결과를 다시 배열 A에 저장하는 것을 뜻한다. (제자리 연산 중 가장 간단한 것 : a += 1)

이렇게 하면 새롭게 메모리를 할당하지 않아도 되기 때문에 성능이 향상된다. (하지만 입력배열과 출력배열이 둘다 캐시에 들어갈만큼 작다면 벡터화의 이점을 누릴 수 있기 때문에 이때는 따로 할당하는 것이 더 속도가 빠르다.)

 

아래 코드가 제자리 연산을 위한 out 라는 배열을 사용한 코드이다.

def laplacian(grid, out):
    np.copyto(out, grid)
    out *= -4
    out += np.roll(grid, +1, 0)
    out += np.roll(grid, -1, 0)
    out += np.roll(grid, +1, 1)
    out += np.roll(grid, -1, 1)


def evolve(grid, dt, out, D=1):
    laplacian(grid, out)
    out *= D * dt
    out += grid

 

하지만 이 역시도 np.roll을 수행하면서 임시 배열을 할당하기 때문에 완전히 해결된 코드는 아니다.

 

이 문제는 아래의 코드처럼 가독성이 떨어지지만 np.roll의 로직을 이용하고 임시 배열을 만들지 않는 roll_add 함수를 만들어 해결할 수 있다.

 

import numpy as np


def roll_add(rollee, shift, axis, out):
    
    if shift == 1 and axis == 0:
        out[1:, :] += rollee[:-1, :]
        out[0, :] += rollee[-1, :]
    elif shift == -1 and axis == 0:
        out[:-1, :] += rollee[1:, :]
        out[-1, :] += rollee[0, :]
    elif shift == 1 and axis == 1:
        out[:, 1:] += rollee[:, :-1]
        out[:, 0] += rollee[:, -1]
    elif shift == -1 and axis == 1:
        out[:, :-1] += rollee[:, 1:]
        out[:, -1] += rollee[:, 0]


def test_roll_add():
    rollee = np.asarray([[1, 2], [3, 4]])
    for shift in (-1, +1):
        for axis in (0, 1):
            out = np.asarray([[6, 3], [9, 2]])
            expected_result = np.roll(rollee, shift, axis=axis) + out
            roll_add(rollee, shift, axis, out)
            assert np.all(expected_result == out)


def laplacian(grid, out):
    np.copyto(out, grid)
    out *= -4
    roll_add(grid, +1, 0, out)
    roll_add(grid, -1, 0, out)
    roll_add(grid, +1, 1, out)
    roll_add(grid, -1, 1, out)