データの傾きを求めたい ==> scipy.interpolate.Rbf を利用した微分

1. 概要

得られたデータを補間しておおよその分布や傾向を把握するという事を, 多くの人は行ったことがあると思います.
(一番身近な例はExcelの機能で近似曲線だと思います.)
そして時々, データを結び付けた後の線(1次元補間結果)や面(2次元補間結果)の傾きを求める必要に遭遇します.
データ点の傾きを求めるには色々なアプローチがあると思いますが, 本稿ではRbf補間を利用する方法を紹介します.
ライブラリは scipy.interpolate.Rbf を利用します.

2. Rbf 補間とは ?

Rbf補間は放射基底関数(Radial Basis Function, RBF)という関数を足し合わせることでデータの補間を行う方法です.
放射基底関数とは, ある適当な点からの距離にだけ依存する関数のことで, Rbf補間の場合は得られたデータ点がこの適当な点になります.

たとえば, データ点 $(x_1, y_1), (x_2, y_2), (x_3, y_3), (x_4, y_4)$ という4点を補間する関数$f(x)$は

f(x)=w_1\phi(r_1) + w_2\phi(r_2)+w_3\phi(r_3)+w_4\phi(r_4)\\
r_j = \sqrt{(x-x_j)^2}

で得られます. 式中の $\phi(r)$ が放射基底関数です. また $w_1, w_2, w_3, w_4$ は

y_1=w_1\phi(r_{11}) + w_2\phi(r_{12})+w_3\phi(r_{13})+w_4\phi(r_{14})\\
y_2=w_1\phi(r_{21}) + w_2\phi(r_{22})+w_3\phi(r_{23})+w_4\phi(r_{24})\\ 
y_3=w_1\phi(r_{31}) + w_2\phi(r_{32})+w_3\phi(r_{33})+w_4\phi(r_{34})\\
y_4=w_1\phi(r_{41}) + w_2\phi(r_{42})+w_3\phi(r_{43})+w_4\phi(r_{44})\\
r_{ij}=\sqrt{(x_i-x_j)^2}

を満たす値になります.
($w_1, w_2, w_3, w_4$という4変数に対して方程式が４個なので自動的に定まります. 実際に $w$ を求める計算はscipy.interpolate.Rbfが自動的に行うので割愛します.)

$\phi(r)$ には様々な種類があります. scipy.interpolate.Rbfの場合は引数functionで指定可能です.

function	$\phi(r)$
multiquadric	$\sqrt{(\frac{r}{\varepsilon})^2 + 1} $
inverse	$\frac{1}{\sqrt{(\frac{r}{\varepsilon})^2 + 1}}$
gausssian	$\exp(-(\frac{r}{\varepsilon})^2)$
linear	$r$
cubic	$r^3$
quintic	$r^5$
thin plate	$r^2\log(r)$

補足しておくと, multiquadric, inverse, gaussianには$\varepsilon$という設定値があります.
特に指定しなければ, scipy.interpolate.Rbfは補間したときに使用した$x_1, x_2,...,x_n$の平均距離をデフォルト値にします.

$$
\varepsilon=\frac{\max(x_1,x_2,...,x_n)-\min(x_1,x_2,...,x_n)}{n}
$$

(平均距離であれば, 分母は$n$ではなく$n-1$の方が正しいはずですが, scipy.interpolate.Rbfでは分母を$n$としています.)

3. Rbf 補間による傾きの導出

一例として下記サンプル関数 $g(x)$: sample_func_g に対して, $x_i=10, 20, 30, 40, ...,90$ の点が得られている場合を考えます.
Rbf 補間によって導出される傾きの値は $dg/dx$: sample_func_dg_dx と近い値になるでしょうか.

$$
g(x)=x+20 \sin{\left(\frac{2\pi x}{60}\right)}+100 \exp\left[-\left(\frac{x-50}{8}\right)^2\right]
$$

import numpy as np
from scipy.interpolate import Rbf
from matplotlib import pyplot as plt

def sample_func_g(x):
    return x + 20*np.sin(2*np.pi*x/60) + 100*np.exp(-((x-50)/8)**2)

x_j = np.arange(10, 95,10)
y_j = sample_func_g(x_j) 

x = np.arange(0, 101,1)
y = sample_func_g(x)

fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111)
ax.plot(x, y, color='magenta', lw=3, label=r'$g(x)$')
ax.scatter(x_j, y_j, color='blue', s=100, zorder=2, label=r'$x_j$')
ax.grid()
ax.set_xlim(0,100)
ax.set_ylim(0,140)
ax.legend(fontsize=16)
plt.show()

参考として, 補間した結果をプロットすると次のようになります.
放射基底関数の種類によって微分結果に違いがあるのが分かります.

function_list=['multiquadric', 'inverse', 'gaussian', 'linear', 'cubic', 'quintic', 'thin_plate']
fig, axes = plt.subplots(figsize=(32,10) ,nrows=2, ncols=4)
for i, function in enumerate(function_list):
    row_num = i//axes.shape[1]
    col_num = i%axes.shape[1]
    ax = axes[row_num, col_num]

    interp_model = Rbf(x_j, y_j, function=function)

    ax.scatter(x_j, y_j, color='blue', s=50, zorder=2, label='points for interpolation')
    ax.plot(x, interp_model(x), color='green', linestyle='dashed', lw=3, label=r'$f(x)$')
    ax.plot(x,y, color='magenta', lw=3, label=r'$g(x)$')
    ax.set_title(function, fontsize=16)
    ax.grid()
    if i == 0:
        ax.legend(fontsize=12)
    ax.set_xlim(0,100)
    ax.set_ylim(0,140)

補間結果を見ると何となく分かりますが, $f(x)$を微分すると$g(x)$の微分と大体同じになりそうです.
scipy.interpolate.Rbfに微分する機能は現段階(2023年3月現在)ではありません.
本稿では, scipy.interpolate.Rbfの属性値を利用して微分結果を取得する方法を説明します.

まず, $f(x)$は
$$
f(x)=\sum^N_j{w_j\phi(r_j)}\
r_j=\sqrt{(x-x_j)^2}
$$
なので, その微分$\frac{df}{dx}$は,

\begin{align}
\frac{df}{dx} &= \sum^N_j{w_j \frac{\partial \phi(r_j)}{\partial r} \frac{dr_j}{dx}} \\
&= \sum^N_j{w_j \frac{\partial \phi(r_j)}{\partial r} \frac{x-x_j}{r_j}}
\end{align}

となります.

各$\frac{\partial \phi(r)}{\partial r}$: phi_partial_diff_functionは以下のようになります.

function	$\phi(r)$	$\frac{\partial \phi(r)}{\partial r}$
multiquadric	$\sqrt{(\frac{r}{\varepsilon})^2 + 1}$	$\frac{r}{\varepsilon^2 \sqrt{(\frac{r}{\varepsilon})^2 + 1}}$
inverse	$\frac{1}{\sqrt{(\frac{r}{\varepsilon})^2 + 1}} $	$ -\frac{r}{\varepsilon^2 \left[(\frac{r}{\varepsilon})^2 + 1\right] \sqrt{(\frac{r}{\varepsilon})^2 + 1}}$
gausssian	$ \exp\left[-(\frac{r}{\varepsilon})^2\right]$	$ -\frac{2r}{\varepsilon^2} \exp\left[-(\frac{r}{\varepsilon})^2\right]$
linear	$ r $	$ 1 $
cubic	$ r^3 $	$ 3r^2$
quintic	$r^5$	$5r^4$
thin plate	$r^2\log(r)$	$2r\log(r)+r$

となります.

from scipy.special import xlogy

def phi_partial_diff_function(r, function, epsilon=None):
    if function == 'multiquadric':
        return r/(epsilon**2*np.sqrt((r/epsilon)**2+1))
    elif function == 'inverse':
        return -r/(epsilon**2*((r/epsilon)**2+1)**(3/2))
    elif function == 'gaussian':
        return -2*r/epsilon**2*np.exp(-(r/epsilon)**2)
    elif function=='linear':
        return 1
    elif function == 'cubic':
        return 3*r**2
    elif function == 'quintic':
        return 5*r**4
    elif function == 'thin_plate':
        return 2*xlogy(r,r) + r

データ点群 $x_j$: x_j, $y_j$: y_jを補間した関数の $x=x_i$:x_i における傾き $df/dx$: df_dx_function は以下のようになる.

def div_function(a, b):
    boolean = (a==b) & (b==0)
    b = b + 1*boolean
    return (a/b)*(boolean==False) + 1*boolean
    

def df_dx_function(x_j, y_j, x_i, function):
    interp_model = Rbf(x_j, y_j, function=function)
    def slope_one_point(X):
        r_j = np.sqrt((X - interp_model.xi[0])**2)
        phi_partial_diff=phi_partial_diff_function(r_j, function, interp_model.epsilon)
        return_slope = np.sum(interp_model.nodes*phi_partial_diff*div_function(X-interp_model.xi[0], r_j))
        return return_slope
    vslope_one_point = np.vectorize(slope_one_point)

    return vslope_one_point(x_i)

補足しておくと, interpolate.Rbf の属性値はそれぞれ以下のように対応しています.
$x_j$: Rbf.xi[0]
$y_j$: Rbf.di
$w_j$: Rbf.nodes

df_dx_functionによって求まる傾きが実際の傾きとどの程度近いかを見てみます.
先ほど使用したサンプル関数$g(x)$の微分$dg/dx$:sample_func_dg_dxは下記のように求められます.

$$
g(x)=x+20 \sin{\left(\frac{2\pi x}{60}\right)}+100 \exp\left[-\left(\frac{x-50}{8}\right)^2\right]
$$

$$
\frac{dg}{dx}=1+20\left(\frac{2\pi}{60}\right)\cos\left(\frac{2\pi x}{60}\right) + 100 \exp\left[-\left(\frac{x-50}{8}\right)^2\right]\left[ -2\left(\frac{x-50}{8}\right)\left(\frac{1}{8}\right) \right]
$$

def sample_func_dg_dx(x):
    return 1 + 20*(2*np.pi/60)*np.cos(2*np.pi*x/60) + 100*np.exp(-((x-50)/8)**2)*(-2*((x-50)/8)*(1/8))

それぞれの方法で補間して傾きを求めた結果を描画してみました.
おおよその傾きを再現できていることが分かります.

function_list=['multiquadric', 'inverse', 'gaussian', 'linear', 'cubic', 'quintic', 'thin_plate']
fig, axes = plt.subplots(figsize=(32,10) ,nrows=2, ncols=4)
for i, function in enumerate(function_list):
    row_num = i//axes.shape[1]
    col_num = i%axes.shape[1]
    ax = axes[row_num, col_num]

    interp_model = Rbf(x_j, y_j, function=function)

    ax.plot(x, df_dx_function(x_j, y_j, x, function=function), color='green', lw=3, linestyle='dashed', label=r'$df/dx$')
    ax.plot(x,sample_func_dg_dx(x), color='magenta', lw=3, label=r'$dg/dx$')
    #ax.scatter(x_j, y_j, color='blue', s=50, zorder=2)
    ax.set_title(function, fontsize=16)
    #ax.plot(x, interp_model(x), color='green', linestyle='dashed', lw=3)
    ax.grid()
    if i==0:
        ax.legend(fontsize=12)
    ax.set_xlim(0,100)
    ax.set_ylim(-10,15)
plt.show()

4. Rbf 補間による傾きの導出, 2次元以上の場合

2次元以上の場合も実は補間関数$f(x, y, z,...)$はあまり変わらず,

f(x, y, z,...)=\sum^N_j{w_j\phi(r_j)}\\
r_j=\sqrt{(x-x_j)^2+(y-y_j)^2+(z-z_j)^2+...}

なので, その微分 $df/dx, df/dy, df/dz...$ は

\frac{df}{dx} = \sum^N_j{w_j \frac{\partial \phi(r_j)}{\partial r} \frac{x-x_j}{r_j}}\\
\frac{df}{dy} = \sum^N_j{w_j \frac{\partial \phi(r_j)}{\partial r} \frac{y-y_j}{r_j}}\\
\frac{df}{dz} = \sum^N_j{w_j \frac{\partial \phi(r_j)}{\partial r} \frac{z-z_j}{r_j}}\\
...

となります.

1次元の場合と同様に下記サンプル関数 $h(x, y)$: sample_func_h を $x$ で微分した $\partial h/\partial x$: sample_func_dh_dx をどの程度再現できるかを確認してみます.

$$
h(x) = 20\sin\left(\frac{2\pi x}{40}\right)\cdot 2y + 20\sin\left(\frac{2\pi y}{60}\right)\cdot 3x + 6000\exp\left[-\left(\frac{x^2+y^2}{60^2}\right)\right]
$$

$$
\frac{dh}{dx}=20\left(\frac{2\pi}{40}\right)\cos\left(\frac{2\pi x}{40}\right)\cdot 2y + 20\sin\left(\frac{2\pi y}{60}\right)\cdot 3+6000\exp\left[-\left(\frac{x^2+y^2}{60^2}\right)\right]\left(-\frac{2x}{60^2}\right)
$$

def sample_func_h(x,y):
    return 20*np.sin(2*np.pi*x/40)*2*y + 20*np.sin(2*np.pi*y/60)*3*x + 6000*np.exp(-(x**2+y**2)/(60**2))

def sample_func_dh_dx(x,y):
    return 20*(2*np.pi/40)*np.cos(2*np.pi*x/40)*2*y + 20*np.sin(2*np.pi*y/60)*3 + 6000*np.exp(-(x**2+y**2)/(60**2))*(-2*x/(60**2))

x = np.arange(-50, 51, 1)
y = np.arange(-50, 51, 1)
x_mesh, y_mesh = np.meshgrid(x, y)

x_j = np.arange(-50, 60, 10)
y_j = np.arange(-50, 60, 10)
x_j_mesh, y_j_mesh = np.meshgrid(x_j, y_j)

def df_dx(x_j, y_j, z_j, x_i, y_i, function):
    interp_model = Rbf(x_j, y_j, z_j, function=function)
    def slope_one_point(X, Y):
        r_j = np.sqrt((X - interp_model.xi[0])**2 + (Y - interp_model.xi[1])**2)
        phi_partial_diff=phi_partial_diff_function(r_j, function, interp_model.epsilon)
        return_slope = np.sum(interp_model.nodes*phi_partial_diff*div_function(X-interp_model.xi[0], r_j))
        return return_slope
    vslope_one_point = np.vectorize(slope_one_point)

    return vslope_one_point(x_i, y_i)

サンプル関数 $h(x)$ をプロットすると以下のようになり,

fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)
ax.contourf(x_mesh, y_mesh, sample_func_h(x_mesh, y_mesh), cmap='jet', vmin=-2000, vmax=7000)
ax.set_xlim(-50, 50)
ax.set_ylim(-50, 50)
ax.grid()
plt.show()

$dh/dx$ をプロットすると以下のようになります.

fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)
ax.contourf(x_mesh, y_mesh, sample_func_dh_dx(x_mesh, y_mesh), cmap='jet', vmin=-400, vmax=400)
ax.set_xlim(-50, 50)
ax.set_ylim(-50, 50)
ax.grid()
plt.show()

先ほどの1次元関数の場合と同じように, 各補間関数で補間して $x$ に対する傾きを求めた結果を図示します.

function_list=['multiquadric', 'inverse', 'gaussian', 'linear', 'cubic', 'quintic', 'thin_plate']
fig, axes = plt.subplots(figsize=(32,10) ,nrows=2, ncols=4)
for i, function in enumerate(function_list):
    row_num = i//axes.shape[1]
    col_num = i%axes.shape[1]
    ax = axes[row_num, col_num]

    ax.contourf(x_mesh, y_mesh, df_dx(x_j_mesh, y_j_mesh, sample_func_h(x_j_mesh, y_j_mesh), x_mesh, y_mesh, function), cmap='jet', vmin=-400, vmax=400)
    ax.set_title(function, fontsize=16)
    #ax.plot(x, interp_model(x), color='green', linestyle='dashed', lw=3)
    ax.grid()

    ax.set_xlim(-50, 50)
    ax.set_ylim(-50, 50)
plt.show()

linearの場合は少し怪しいですが, 2次元の場合でも同様に傾きをおおよそ再現できることが確認できました.