はじめに

先日数値計算をしていて、勾配を計算したい場面がありました。

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()を使った計算例