행렬과 벡터 연산 (2) - 확산 방정식 순수 파이썬 구현

 

행렬과 벡터 연산 (1) - 확산 방정식 예제

---

지난 포스트의 의사 코드를 기반으로 순수 파이썬으로 (특히 리스트) 구현해 보겠습니다. 먼저 행렬을 받아 변화된 상태를 반환하는 evolve 함수를 구현합니다.

grid_shape = (640,640)

def evolve(grid, dt, D=1.):
    xmax, ymax = grid_shape
    new_grid = [[0.,0.] * ymax for x in range(xmax)]
    for i in range(xmax):
        for j in range(ymax):
            grid_xx = grid[(i+1)%xmax][j] + grid[(i-1)%xmax][j] - 2. * grid[i][j]            
            grid_yy = grid[i][(j+1)%ymax] + grid[i][(j-1)%ymax] - 2. * grid[i][j]
            new_grid[i][j] = grid[i][j] + D * (grid_xx + grid_yy) * dt
    return new_grid

• grid_shape 는 전역 변수입니다. 
• 실제 이 코드르 사용하게 위해서는 격자를 초기화 한 이후에 이 함수를 호출해야합니다.

---

이제 격자를 초기화하면서 evolve 함수를 호출하는 함수를 작성합니다.

def run_experiment(num_iterations):
    xmax, ymax = grid_shape
    grid = [[0.,0.] * ymax for x in range(xmax)]

    ## 초기 조건
    block_low = int(grid_shape[0]*0.4)
    block_high = int(grid_shape[0]*0.5)
    for i in range(block_low, block_high):
        for j in range(block_low, block_high):
            grid[i][j] = 0.005

    for i in range(num_iterations):
        grid = evolve(grid, dt=0.1)

 

프로파일링

Line Profiler

실제로 위 코드를 돌리면 굉장히 오래 걸립니다. (num_iterations 를 500으로 주었습니다) line_profiler 를 이용해서 어느 라인에서 병목이 생기는지 살펴보기 위해 evolve 함수 위에 @profiler 데코레이터를 추가한 이후에 다음 명령어를 실행합니다.

>>> kernprof -l -v diffusion.py
Wrote profile results to diffusion.py.lprof
Timer unit: 1e-06 s

Total time: 742.194 s
File: diffusion.py
Function: evolve at line 3

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     3                                           @profile
     4                                           def evolve(grid, dt, D=1.):
     5       500       1302.8      2.6      0.0      xmax, ymax = grid_shape
     6       500     944569.7   1889.1      0.1      new_grid = [[0.,0.] * ymax for x in range(xmax)]
     7    320500     230697.2      0.7      0.0      for i in range(xmax):
     8 205120000  113614441.9      0.6     15.3          for j in range(ymax):
     9 204800000  230679505.0      1.1     31.1              grid_xx = grid[(i+1)%xmax][j] + grid[(i-1)%xmax][j] - 2. * grid[i][j]
    10 204800000  220867632.3      1.1     29.8              grid_yy = grid[i][(j+1)%ymax] + grid[i][(j-1)%ymax] - 2. * grid[i][j]
    11 204800000  175855490.1      0.9     23.7              new_grid[i][j] = grid[i][j] + D * (grid_xx + grid_yy) * dt
    12       500        347.6      0.7      0.0      return new_grid

• "xmax, ymax = grid_shape" 부분에서의 Per Hit는 다른 라인에 비해 큽니다. 이는 함수 안에서 grid_shape 전역 변수를 가져와야 하기 때문입니다.
• new_grid 변수를 생성하는 부분을 보면 Per Hit 가 굉장히 크지만 % Time 은 작습니다. 이것은 이 줄을 실행하는 것 자체는 (Per Hit) 시간이 꽤 걸리나 루프 안의 다른 줄에 비해 훨씬 덜 호출되기 때문입니다. (따라서 전체 실행 시간에서 차지하는 비율 % Time 이 낮게 나옵니다) 이를 만약에 격자 크기를 줄이고 num_iterations 를 늘려 함수 호출 횟수를 늘리게 되면 % Time 또한 높게 나오리라는 것을 예측할 수 있겠죠. 
• "for i in range(xmax)" 부분은 320,500 번 실행되었습니다. xmax 가 640 이고 500번 실행했기 때문인데, 루프 종료 조건을 만나기 위해 640번에 한 번 더 실행되기 때문입니다. 

개선 - 불필요한 변수 생성 없애기

일단 개선할 수 있는 지점은 new_grid 변수를 생성하는 부분일 겁니다. new_grid 리스트는 모양과 크기 저장값이 같기 때문에 evolve 함수를 호출할 때마다 생성하는 것이 아니라 처음부터 한 번만 저장하고 재활용하는 방법으로 최적화 할 수 있겠죠. 즉, 루프 밖으로 빼야 합니다. 다음 예를 보면 더 명확해집니다.

from math import sin

def loop_slow(num_iterations):
    result = 0
    for i in range(num_iterations):
        result += i * sin(num_iterations)
    return result

def loop_fast(num_iterations):
    result = 0
    factor = sin(num_iterations)
    for i in range(num_iterations):
        result += i
    return result * factor
    
>>> %timeit loop_slow(10000)
1.62 ms ± 41.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

>>> %timeit loop_fast(10000)
702 µs ± 86.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

• 루프 안에서 sin 값을 매번 계산할 필요가 없으므로 속도가 훨씬 빨라집니다.

---

위의 예제를 참고하면 evolve, run_experiment 함수를 다음과 같이 개선할 수 있습니다.

def evolve(grid, grid_out, dt, D=1.):
    xmax, ymax = grid_shape
    for i in range(xmax):
        for j in range(ymax):
            grid_xx = grid[(i+1)%xmax][j] + grid[(i-1)%xmax][j] - 2. * grid[i][j]
            grid_yy = grid[i][(j+1)%ymax] + grid[i][(j-1)%ymax] - 2. * grid[i][j]
            grid_out[i][j] = grid[i][j] + D * (grid_xx + grid_yy) * dt
            
            
def run_experiment(num_iterations):
    xmax, ymax = grid_shape
    next_grid = [[0.,0.] * ymax for x in range(xmax)]
    grid = [[0.,0.] * ymax for x in range(xmax)]

    ## 초기 조건
    block_low = int(grid_shape[0]*0.4)
    block_high = int(grid_shape[0]*0.5)
    for i in range(block_low, block_high):
        for j in range(block_low, block_high):
            grid[i][j] = 0.005

    for i in range(num_iterations):
        evolve(grid, next_grid, dt=0.1)
        grid, next_grid = next_grid, grid  # next_grid 갱신하였으므로 grid는 재활용

---

이를 다시 프로파일링하면 다음과 같습니다. 매 루프마다 불필요한 리스트 할당을 없애니 속도가 15% 정도 빨라졌습니다. 변수나 리스트를 저장하려고 메모리를 요청할 때마다 파이썬은 운영체제에 새로운 공간 할당을 요청해야하고 이 과정이 속도적인 측면에서 병목이 되기 때문입니다. 따라서 가능하다면 이미 할당된 공간을 재사용하는 것이 성능 향상에 도움이 됩니다.

>>> kernprof -l -v diffusion.py
Wrote profile results to diffusion.py.lprof
Timer unit: 1e-06 s

Total time: 638.186 s
File: diffusion.py
Function: evolve at line 3

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     3                                           @profile
     4                                           def evolve(grid, grid_out, dt, D=1.):
     5       500        601.4      1.2      0.0      xmax, ymax = grid_shape
     6    320500     168701.9      0.5      0.0      for i in range(xmax):
     7 205120000   88100561.4      0.4     13.8          for j in range(ymax):
     8 204800000  202270296.0      1.0     31.7              grid_xx = grid[(i + 1) % xmax][j] + grid[(i - 1) % xmax][j] - 2. * grid[i][j]
     9 204800000  193962610.2      0.9     30.4              grid_yy = grid[i][(j + 1) % ymax] + grid[i][(j - 1) % ymax] - 2. * grid[i][j]
    10 204800000  153682800.6      0.8     24.1              grid_out[i][j] = grid[i][j] + D * (grid_xx + grid_yy) * dt

 

---

행렬과 벡터 연산 (3) - 파이썬 리스트와 numpy 연산 속도 차이

행렬과 벡터 연산 (4) - numpy 배열을 이용한 확산 방정식

행렬과 벡터 연산 (5) - numpy 배열 메모리 최적화

행렬과 벡터 연산 (6) - numexpr 모듈 이용하기