[研究レポート] 拡散式パワーレーダー

[研究レポート] 拡散式パワーレーダー

拡散式パワーレーダー付降下勾配法

拡散式パワーレーダーの仕組み

 拡散式パワーレーダー は、一定ステップごとに全領域から無作為に複数の点を選んで高度を測り、その中から最も深い場所を選んで転移(テレポート)します。

Python研究レポート 拡散式パワーレーダー

 一度に数多くの地点を探査できる反面、かなりの重装備になるのではと心配したのですが、乱数配列やマスキング操作など、NumPy の高速処理機能をフル活用してみると、意外と軽量にサクサク動きました。NumPy、恐るべしですね。最初に平面上の複数地点をランダムに選ぶために、r_mesh() という関数を定義しておきます。

# リストR3-A-1

import numpy as np

np.random.seed(0)

def r_mesh(xr, yr, n):
    x = np.random.uniform(xr[0], xr[1], (n, n))
    y = np.random.uniform(yr[0], yr[1], (n, n))
    X, Y = np.array([x, y])
    return X, Y

 r_mesh() は無作為地点を抽出して格子点を生成します。格子点のデータ形式は numpy.meshgrid() によって生成されるデータと同じように、2 つの 2 次元配列の組合わせで表されます。たとえば、x と y の範囲をそれぞれ [0, 20], [0, 10] に設定し、9 地点を同時探索する場合は次のようなコードを記述します。

# リストR3-A-2

xr = [0, 20]
yr = [0, 10]
X, Y = r_mesh(xr, yr, 3)

print(X)
print(Y)
[[10.97627008 14.30378733 12.05526752]
 [10.89766366  8.47309599 12.91788226]
 [ 8.75174423 17.83546002 19.27325521]]
[[3.83441519 7.91725038 5.2889492 ]
 [5.68044561 9.25596638 0.71036058]
 [0.871293   0.20218397 8.32619846]]

 2 変数関数を定義して格子点データを入れてみます。

# リストR3-A-3

def test_func(x, y):
    return np.sin(x * y)

Z = test_func(X, Y)
print(Z)
[[-0.94799371  0.14877576  0.8003294 ]
 [-0.80058291  0.1128836   0.24587626]
 [ 0.97397369 -0.44793251 -0.24907247]]

 np.min(Z) で配列 Z の最小要素(探査地点のうち最も深い場所)を取得することができます。

# リストR3-A-4

Z_min = np.min(Z)
print(Z_min)
-0.9479937145939273

 マスキング機能を使うと、Z の最小要素に対応する X, Y の要素を取得できます。

# リストR3-A-5

X_min = X[Z == Z_min][0]
Y_min = Y[Z == Z_min][0]
print("(X_min, Y_min) = ({0:.2f}, {1:.2f})".format(X_min, Y_min))
(X_min, Y_min) = (10.98, 3.83)

 このように絞り込んだ探査地点と現在地点の高度を比較して、探査地点のほうが深ければ、そこへ転移するようにします。
 

拡散式パワーレーダーによる最小値の探索

 とりあえずレーダーの性能を調べるために勾配降下機能は切り離して、拡散式パワーレーダー のみを用いて最小値へ近づいてみます。テスト関数は 前回記事 で使用したものと同じです。

# リストR3-B-1

import numpy as np
import matplotlib.pyplot as plt

# 乱数を初期化
np.random.seed(0)

def mesh_2d(xr, yr, n):
    x = np.linspace(xr[0], xr[1], n)
    y = np.linspace(yr[0], yr[1], n)
    X, Y = np.meshgrid(x, y)
    return X, Y

def r_mesh(xr, yr, n):
    x = np.random.uniform(xr[0], xr[1], (n, n))
    y = np.random.uniform(yr[0], yr[1], (n, n))
    X, Y = np.array([x, y])
    return X, Y

# テスト関数を定義
def func(x ,y):
    z = 2 * np.sin(x) + np.sin(np.sqrt(np.abs(x))) + np.sin(y/2)
    return z

# 等高線描画データ
xr = [0, 12]
yr = [0, 12]
n = 257
X, Y = mesh_2d(xr, yr, n)
Z = func(X, Y)

# Figureと3DAxesの設定
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(111)
ax.set_xlabel("x", size = 16)
ax.set_ylabel("y", size = 16)
ax.set_xlim([xr[0], xr[1]])
ax.set_ylim([yr[0], yr[1]])

# 等高線をプロット
ct = ax.contour(X, Y, Z, 12, colors = "black")
ax.clabel(ct, fmt = "%0.1f", fontsize = 10, zorder = 1)

# ★★★★★★★★★★

# 初期地点を設定
pt = [2, 10]

# 初期地点の高度
zt = func(pt[0], pt[1])

# ループ最大数
t_max = 100

for t in range(t_max):
    Xs, Ys = r_mesh(xr, yr, 32)
    Zs = func(Xs, Ys)
    Zs_min = np.min(Zs)
    if zt > Zs_min:
        Xs_min = Xs[Zs == Zs_min][0]
        Ys_min = Ys[Zs == Zs_min][0]
        pt = [Xs_min, Ys_min]
        ax.plot(pt[0], pt[1], color = "red",
                marker = "o", markersize = 1,
                zorder = 1, linestyle = "none")

zt = func(pt[0], pt[1])
pt = np.round(pt, 3)
zt = np.round(zt, 3)
print("(x, y, z) = ({0:.3f}, {1:.3f}, {2:.3f})".format(pt[0], pt[1], zt))
(x, y, z) = (11.185, 9.591, -3.162)

Python 拡散式パワーレーダーのテスト

 32×32 = 1024 地点の同時探索を 100 回行いました。
 前回記事の「回転式レーダー付勾配降下法」で求めた値が

  (x, y, z) = (11.000, 9.440, -3.174)

だったので、最小値にかなり近づいていることがわかります。

拡散式パワーレーダー付勾配降下法

 とはいえ、レーダーに頼った転移だけでは小回りがきかないので、窪みの底に達するためには、やはり勾配降下法と併用する必要があります。処理にかかる負荷を軽減するために、拡散式レーダーは 4 ループに 1 回の頻度で探査を行なうように設定します。

# リストR3-C-1

import numpy as np
import matplotlib.pyplot as plt

def mesh_2d(xr, yr, n):
    x = np.linspace(xr[0], xr[1], n)
    y = np.linspace(yr[0], yr[1], n)
    X, Y = np.meshgrid(x, y)
    return X, Y

def r_mesh(xr, yr, n):
    x = np.random.uniform(xr[0], xr[1], (n, n))
    y = np.random.uniform(yr[0], yr[1], (n, n))
    X, Y = np.array([x, y])
    return X, Y

# 乱数を初期化
np.random.seed(0)

# テスト関数を定義
def func(x ,y):
    z = 2 * np.sin(x) + np.sin(np.sqrt(np.abs(x))) + np.sin(y/2)
    return z

# テスト関数の勾配(gradient)
def grad(x, y):
    gx = 2 * np.cos(x) + 0.5 * np.cos(x) / np.sqrt(x)
    gy = 0.5 * np.cos(0.5 * y)
    return np.array([gx, gy])

# 等高線描画データ
xr = [0, 12]
yr = [0, 12]
sr = [-0.05, 0.05]
n = 257
X, Y = mesh_2d(xr, yr, n)
Z = func(X, Y)

# Figureと3DAxesの設定
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(111)
ax.set_xlabel("x", size = 16)
ax.set_ylabel("y", size = 16)
ax.set_xlim([xr[0], xr[1]])
ax.set_ylim([yr[0], yr[1]])

# 等高線をプロット
ct = ax.contour(X, Y, Z, 12, colors = "black")
ax.clabel(ct, fmt = "%0.1f", fontsize = 10, zorder = 1)

# ★★★★★★★★★★

# 初期地点を設定
pt = [2, 10]

# 初期地点の高度
zt = func(pt[0], pt[1])

# ループの最大数と最小数
t_max = 1000
t_min = 100

# レーダーの探査頻度
fq = 4

# 学習率と収束判定条件
alpha = 0.01
eps = 0.01

for t in range(t_max):
    pt[0] -= alpha * grad(pt[0], pt[1])[0]
    pt[1] -= alpha * grad(pt[0], pt[1])[1]

    if t % fq == 0:
        Xs, Ys = r_mesh(xr, yr, 32)
        Zs = func(Xs, Ys)
        Zs_min = np.min(Zs)
        if zt > Zs_min:
            Xs_min = Xs[Zs == Zs_min][0]
            Ys_min = Ys[Zs == Zs_min][0]
            pt = [Xs_min, Ys_min]
            ax.plot(pt[0], pt[1], color = "red",
                    marker = "o", markersize = 1,
                    zorder = 1, linestyle = "none")

    # 勾配の大きさを計算
    g = grad(pt[0], pt[1])
    d = np.sqrt(g[0]**2 + g[1]**2)
    
    # 最小ステップ数に達し、かつ勾配の大きさが基準値以下になったらループ終了
    if t > t_min and d < eps:
        break

zt = func(pt[0], pt[1])
pt = np.round(pt, 3)
zt = np.round(zt, 3)

print("t = {}".format(t))
print("(x, y, z) = ({0:.3f}, {1:.3f}, {2:.3f})".format(pt[0], pt[1], zt))

Python 機械学習 拡散式パワーレーダー付勾配降下法テスト01

t = 908
(x, y, z) = (10.997, 9.397, -3.174)

 トータルのループ回数は t = 908 となっています。処理時間は 640 ms でした。これは回転式レーダーを使用した場合と比較すると圧倒的に高速です。次はもう少し難しい関数
 
\[f(x, y)=\sin ^2 x +\cos\left(1+\frac{xy}{4}\right)\sin \sqrt{x}\]
を使ってテストしてみます。リストR3-C-1 の関数定義の部分を以下のように変えて実行してみてください。

# テスト関数を定義
def func(x ,y):
    pi = np.pi
    sqx = np.sqrt(x)
    z = np.sin(x)**2 + np.cos(1 + x*y/4) * np.sin(sqx)
    return z

# テスト関数の勾配(gradient)
def grad(x, y):
    sin = np.sin
    cos = np.cos
    sqx = np.sqrt(x)
    gx = 2*sin(x)*cos(x) - y * sin(1+x*y/4) * sin(sqx) / 4
    + 0.5 * cos(1 + x*y/4) * cos(sqx)
    gy = -x * sin(1 + x*y/4)
    return np.array([gx, gy])

 実行結果は次のようになります。

Python 機械学習 拡散式パワーレーダー付勾配降下法テスト02

t = 488
(x, y, z) = (3.141, 2.732, -0.980)

 処理時間は 530 ms でした。関数によっては、拡散式パワーレーダーを使っても処理に時間がかかる場合もありますが、それはレーダーの性能というよりは、窪みの底まで移動する速度、つまり学習率の設定に問題があります。学習率を調整する方法(アルファ関数の定義)についてはこちらの記事を参照してください

 拡散式パワーレーダーは、格子の次元数を増やすことによって、3 つ以上の変数をもつ関数にも応用することができます。次回はそのあたりのことを解説します。