差分法におけるPythonの高速化テクニック

はじめに

Pythonで数値計算をする場合,書きやすさの反面,実行速度がネックになることが多いです.ここでは差分法での高速化のTipsを紹介し,実用可能なレベルの速度であることを示します.

中央差分

偏微分方程式を数値計算で解を導く場合,実装が簡単な差分法で導くことがあります.微分形式を中央差分で近似すると

\begin{align}
\cfrac{\textrm{d} y}{\textrm{d} x} \approx \cfrac{y_{i+1} – y_{i-1}}{x_{i+1} – x_{i-1}}
\end{align}

と表すことができます.このとき,通常なら

for i in range(1, len(x)-1):
    dydx[i] = (y[i+1] - y[i-1])/(x[i+1] - x[i-1])

とfor文で表しますが,numpyを使うことで次のように書くことができます.

dydx = (y[1:] - y[:-1])/(x[1:] - x[:-1])

これだけで相当高速化されます.

三重対角行列

差分法を利用して微分方程式を解く場合,三重対角行列を扱うことが多いです.三重対角行列とは

\begin{align}
A =
\begin{bmatrix}
a_{1,1} & a_{1,2} & 0 & 0 & \dots & 0 \newline
a_{2,1} & a_{2,2} & a_{2,3} & 0 & \dots & 0 \newline
0 & a_{2,2} & a_{2,3} & a_{2,3} & \dots & 0 \newline
\vdots & & & \ddots & & \vdots \newline
& & & & & a_{n-1,n} \newline
0 & \dots & & 0 & a_{n,n-1} & a_{n,n} \newline
\end{bmatrix}
\end{align}

のように,字のごとく対角成分が三重になっている行列です.三重対角行列を生成するとき,

for i in range(0, n-1):
    A[i, i] = a1[i]
    A[i+1, i] = a2[i]
    A[i, i+1] = a3[i]
A[n-1,n-1] = a1[n-1]

という書き方や,numpyの対角行列を生成する関数を使って

A = np.diag(a1) + np.diag(a2, k=-1) + np.diag(a3, k=1)

という書き方がありますが,

i = np.arange(0, n-1)
A[i, i] = a1[i]
A[i+1, i] = a2[i]
A[i, i+1] = a3[i]

A[n-1,n-1] = a1[n-1]

という書き方がこの中では最も早いです.

実験

どのくらい早くなったか計測してみます.今回は以下のような環境で実行時間を計測しました.

OS Gentoo/Linux
Kernel x86_64 Linux 4.4.6-gentoo
CPU Intel Core i7-2600 CPU
Python 3.4.3
Numpy 1.1.11
Scipy 0.18.0
blas reference(20070226-r4)
cblas reference(20030223-r6)


中央差分では1000万の要素を用意し,三重体格行列では25000の正方行列を用意して計測を行った結果は次のようになりました.

中央差分(for) 12.481740[s]
中央差分(numpy) 0.414090[s]
三重対角行列(for) 0.056523[s]
三重対角行列(diag) 18.103698[s]
三重対角行列(numpy) 0.006361[s]


中央差分については30倍,三重対角行列ではfor文との比較では9倍近く高速化されていることがわかります.diagは内部で行列を生成しているのでとても遅くなっています.

まとめ

Pythonの書き方を変えることで十分な高速化を実現することができました.ちょっと書き方を変えるだけで効果は絶大なのでぜひ使ってみてください.最後にソースコードを載せておきます.

コメントを残す

メールアドレスが公開されることはありません。