はじめに
先日数値計算をしていて、勾配を計算したい場面がありました。
NumPyアレイの勾配を自力で求めるのは、差分取るだけでしょと思ったのですが、いざ考えてみると細々と考えることがあり、大変でした。
その時にnp.gradient()というものがあることを知り、とても便利でした。
自力で勾配(微分)を求めると何が大変なのか、更にnp.gradient()の使い方について解説します。
自分で勾配(微分)を求めた時の問題
最初は簡単に1次元のy,xがあるとします。
import numpy as np
y = np.array([1,2,4,5,7,11,16])
x = np.arange(y.size)
簡単にx=0,1,2,,,と1ずつ増えていくとします。
この時に微分$\frac{dy}{dx}$を求めるにはどうすればいいでしょうか？
素朴に考えると、yとxの各項の差分を取って比をとればよいです。
dy = y[1:] - y[:-1]
dx = x[1:] - x[:-1]
dydx = dy / dx
dydxの要素数を確認します。
print(dydx.size)
6
最初は要素数7だったので、要素数が1個減っています。
また、各要素数の差分を取ったため、対応するxの値は0,1,2,,,ではなく
x = 0 1 2 3 4
\ / \ / \ /
0.5 1.5 2.5 ...
のように各区間の中点に対応させるのが自然です。
つまり、自力で計算すると下記のような問題があります。
- 元データから要素が1個減る
- 微分結果をどのxに対応させるか考える必要がある
こういう時にnp.gradient()が便利です。
np.gradient()の使い方(1次元)
先ほどの問題であれば下記で微分が求められます。
dx = x[1] - x[0]
dydx = np.gradient(y,dx)
print(dydx.size)
7
要素数も最初の値を維持しています。
また対応するxの値も0,1,2,,でありずれることはありません。
np.gradient()の使い方(2次元)
2Dの場合もほぼ同じように計算できます。
f = np.array([
[1, 2, 6],
[3, 4, 5]
])
dx = 2.
dy = 1.
gs_y,gs_x = np.gradient(f, dy, dx)
print(gs_y.shape)
print(gs_x.shape)
(2,3)
(2,3)
fが2Dの微分したいNumPyアレイ、x方向とy方向の刻み幅がそれぞれ2と1としました。
np.gradient()はx方向とy方向のそれぞれでの勾配が計算されます。
そのため返り値は2個なので、それぞれ変数で受けます。
また、返り値のサイズはfと同じです。
ここで少し注意が必要なのが、dx と dy の順番です。
NumPyの2次元配列では、
axis=0 → 行方向 → y方向
axis=1 → 列方向 → x方向
となっています。
そのため、np.gradient() の返り値もgs_y, gs_x の順番になります。
また、np.gradient() に渡す刻み幅もこの順番に対応するため、
np.gradient(f, dy, dx)
となります。
実際にnp.gradinet()を用いてスペックル計算を行った例については、以下の記事でPythonコード付きで解説しています。
フレネル反射を考慮したスペックル計算
不等間隔データにも使える
x = np.array([0., 1., 1.5, 3., 5.])
y = np.array([1., 2., 4., 5., 7.])
dydx = np.gradient(y, x)
のように、x に各データ点の座標を入れたNumPy配列を渡すことができます。
x の間隔が一定でなくても、np.gradient() がそれぞれの間隔を考慮して勾配を計算してくれます。
np.gradient()はどうやって微分している？
np.gradient() は、場所によって異なる差分方法を使って微分を計算しています。
データの内部では、前後のデータを使う中央差分が使われます。
$$
\frac{df}{dx} \approx \frac{f_{i+1} - f_{i-1}}{2h}
$$
しかし、データの端では片側にしかデータがありません。
そのため、左端では前進差分を行います。
$$
\frac{df}{dx} \approx \frac{f_{1} - f_{0}}{h}
$$
右端では後退差分を使います。
$$
\frac{df}{dx} \approx \frac{f_{n} - f_{n-1}}{h}
$$
イメージとしてはこんな感じです。
左端 内部 右端
f0 -- f1 f(i-1) -- fi -- f(i+1) fn-1 -- fn
↑ ↑ ↑
前進差分 中央差分 後退差分
np.gradient() はデータの位置に応じて差分方法を変えることで、元の配列と同じ要素数の勾配を計算しています。
また端以外では中央差分を使っているので、対応する座標のずれもあまり気にしなくていいのです。
まとめ
今回は、NumPyの np.gradient() を使った勾配（微分）の計算方法について解説しました。
np.gradient() を使うことで、
- 差分を使った微分でも、元の配列と同じ要素数で結果を得られる
- 1次元だけでなく、2次元のデータも各方向の勾配を計算できる
- データ点の間隔が一定でなくても、
xの配列を渡すことで計算できる
といったことができます。
自力で計算すると「結果をどの位置に対応させるか」や「データ点の間隔が一定か」といったことを考える必要があります。
こうした処理をまとめて扱ってくれる np.gradient() は、NumPyで数値微分を行うときに便利な関数です。
関連記事
-
np.gradient()を使った計算例