回転式探査レーダー・拡散式パワーレーダー

回転式探査レーダー・拡散式パワーレーダー

回転式探査レーダー付勾配降下法

 前回の続きです。今回は 2 変数関数 $z=f(x, y)$ のレーダー勾配降下(←なんかもう名称を略してる)をやります。基本的な考え方自体は 1 変数のときと同じですが、3 次元空間を探査する必要があるので、搭載するレーダーもそれなりに工夫を凝らさなくてはなりません。

回転式探査レーダー

 今回紹介する 回転式レーダー は下図のようなイメージです。

 ビームの方角変更

 $1$ ステップごとに $4$ 本のビームが同時照射されます。
 それぞれのビームの射出方向は $90^\circ$ の角度をなしていて、ステップごとに照射軸が反時計回りに回転しつつ、さらに動径方向のビームの長さ $r(t)$ も周期的に変化させます ($t$ はステップ数)。個々のビームの長さと方角 $\boldsymbol{R}(t)$ はベクトルを使って
 
\[\boldsymbol{R}(t)=\begin{bmatrix}
r(t)\cos \theta(t)\\ r(t)\sin \theta(t)\end{bmatrix}\]
のように表されます。$r(t)$ の設定方法は色々考えられますが、今回は三角関数を使って
 
\[r(t)=X\sin a(t)\]
のように設定しておきます。$X$ は $x$ の最大値 $x_{\mathrm{max}}$ と最小値 $x_{\mathrm{min}}$ の差
 
\[X=x_{\mathrm{max}}-x_{\mathrm{min}}\]
によって定義されます。これが動径方向の振幅となります。
 

テストコード

 今回はテスト関数として
\[z=2\sin x+\sin\sqrt{|x|}+\sin\frac{y}{2}\]
を使うことにします。この関数を等高線で図示すると次のようになります。

 Python テストに使用する関数の等高線図

 点線で描かれた部分が $z$ が負の値をとる領域です。浅い窪地に捕まらないように、目的地にできるだけ早く達することを試みます。

# https://python.atelierkobato.com/rotation/
# リストR2-A-1
# 回転式レーダー付勾配降下法

# 必要なモジュールをインポート
import numpy as np
import matplotlib.pyplot as plt

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

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

# ★★★★★★★★★★

# 等高線を描画
n = 129
x = np.linspace(0, 12, n)
y = np.linspace(0, 12, n)
x_max = np.max(x)
x_min = np.min(x)
y_max = np.max(y)
y_min = np.min(y)
X, Y = np.meshgrid(x, y)
Z = func(X, Y)

# FigureとAxesの設定
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([x_min, x_max])
ax.set_ylim([y_min, y_max])

# 等高線を描画
ct = ax.contour(X, Y, Z, 12, colors = "black")
ax.clabel(ct, fmt = "%0.1f", fontsize = 10)

# ★★★★★★★★★★

# 回転式探査レーダー勾配降下法による最小値の探索

# 初期値の設定
x_init, y_init = (2, 10)

# x,yに初期値を代入
x, y = (x_init, y_init)

# 学習率α
alpha = 0.01

# 収束判定条件
eps = 0.01

# ループの最大数と最小数
t_max = 2000
t_min = 500

# 回転レーダーの各種設定
a = 0.2
u = 1.0
du = 1.0
angle_0 = np.arange(0, np.pi, np.pi/4)
d_angle = np.pi / 16
angle = angle_0

xs_list = np.zeros(4)
ys_list = np.zeros(4)

for t in range(1, t_max):
    u += du
    angle += d_angle
    x -= alpha * grad(x, y)[0]
    y -= alpha * grad(x, y)[1]
    
    r = (x_max - x_min) * np.sin(a * u)
    xs_list = x + r * np.cos(angle)
    ys_list = y + r * np.sin(angle)
    zs_list = func(xs_list, ys_list)
    zs_min = np.min(zs_list)
    xs = xs_list[zs_list == zs_min]
    ys = ys_list[zs_list == zs_min]

    # 転移条件
    tf_1 = func(xs, ys) < func(x, y)
    tf_2 = x_min < xs[0] < x_max
    tf_3 = y_min < ys[0] < y_max

    # 条件をすべて満たせば転移して、ビーム長さを初期化
    if tf_1 and tf_2 and tf_3:
        u = 0
        x = xs[0]
        y = ys[0]

    # (x, y)をプロット
    ax.plot(x, y, color = "red",
            marker = "o", markersize = 1,
            zorder = 1, linestyle = "none")

    # 勾配の大きさの計算
    d = np.sqrt(grad(x, y)[0]**2 + grad(x, y)[1]**2)
    
    # 条件を満たせばループ終了
    if t > t_min and d < eps:
        break

print("t = {}".format(t))
print("(x, y, z) = ({0:.3f}, {1:.3f}, {2:.3f})".format(x, y, func(x, y)))
t = 1080
(x, y, z) = (11.000, 9.440, -3.174)

 窪地のレーダー探査(等高線)

 無事に目的地に到達することができました。勾配法ってゲームみたいで面白いですよね。難しい関数を用意して、それを攻略するようなレーダーを開発するみたいなことを繰り返せば、楽しみながらコードを改良できます。
 

ステップ数短縮効果と処理時間増加のトレードオフ

 図の赤い軌跡を見ると、窪地を見つけるまでのステップ数は多くなさそうです(レーダーによって素早く補足されています)が、その窪みの斜面を下ること(つまり普通の勾配降下)にステップ数の大半を費やしています。前回記事の最後で「ビームが近い場所を探査するときにショートカットすることがある」と言いましたが、これを意図的に実行するために、非常に短いビームを照射するショートレンジレーダーを搭載することも考えられます。とはいえ、レーダーの搭載数を増やすと負荷がかかって処理速度が遅くなることがあるので、慎重にステップ短縮効果と処理時間増加の兼ね合いを見極める必要があります(Jupyter Notebook における処理時間の計測方法については、こちらの記事を参照してください)。でもそんな細かい事をやる前に、もっと派手なレーダーを試してみたいのですよね?(←ただ自分がやりたいだけ)。というわけで、次は拡散式パワーレーダーについて解説します。
 

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

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

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

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 つ以上の変数をもつ関数にも応用することができます。次回はそのあたりのことを解説します。