[研究レポート] 探査レーダー付勾配降下法

[研究レポート] 探査レーダー付勾配降下法

研究レポートシリーズをはじめました

 研究レポートシリーズは、BlogCat の個人的な研究ノートを公開してしまおうという大胆な企画です。質の良し悪しはともかくとして、一般の書籍などには載っていないようなアルゴリズムを見れちゃったりします。
「そんなの公開しちゃっていいの?」
と思われるかもしれませんが、今はオープンソースな時代(?)だし、ノートを見た皆さんから
「ここは、こうしたほうがいいんじゃない?」
みたいな意見をお寄せいただければ、私自身のコードを改良することもできるし、何よりこのサイトの評判が上がれば収入も増え ...... じゃなくて、公開によって多少なりとも日本の技術発展に貢献できるのではないかと考えております(←そんな大したものではない)。

 思いつきで作ったようなコードもぱっと載せたりするので、それを叩き台にして皆さんと共同作業で開発できたらいいなと考えております。当サイトの他の記事にあるような堅苦しい口調はやめて、ざっくばらんな口調でお話します。

探査レーダー付勾配降下法

 記念すべき(?)第 1 回目は、勾配降下法の新しい手法である 探査レーダー付勾配降下法 (GDSR : Gradient Descent with Search Radar) を公開します(≫ 勾配降下法の基本アルゴリズムについては、こちらの記事をご確認ください)。

 勾配降下法というのは、(斜面の上に置かれた点にとって)自分の足元の情報しかないのですね。曲面全体の様子なんて全くわからない。暗闇の中で斜面をそろりそろりと降りていく感じです。だから窪みの底に達してしまうと、そこが本当に世界で一番深いところなのかどうかわからない。他にもあるかもしれないけど、まさか全部歩いて調べるわけにもいかない(処理時間が膨大になる )。これが勾配降下法の弱点なのです。

 というわけで、もっと周囲の状況がわかるように、レーダー機能を備えてあげようというのが 探査レーダー付勾配降下法 の考え方です。1変数関数の場合は下図のようなイメージになります。

Python 窪地の探査イラスト01

 ビームはステップごとに斜面に置かれた点から水平に発射されて深さ (曲線の y 座標) を測ります。そこが現在居る場所よりも深ければ、その点まで転移します。

 どうでもいいような窪みをさっさと放棄できるので、探査時間を節約できます。搭載するレーダーの細かな仕様は色々考えられますが、基本的にビームの長さは(三角関数などを使って)ステップごとに少しずつ変動させるようにしておきます。複数のレーダーを上手く組合わせると、さらに効率を上げることができます(組合せ方は結構難しいので、あれこれ試行錯誤も必要です)。ただし、あまり色々と詰め込むと処理に負担がかかるので、2 個か 3 個ぐらいにしておいたほうがいいでしょう。今回は 2 種類のレーダーを搭載します。
 

概周期関数

 テストに用いる関数は $y=\sin x+\sin\sqrt{x}$ です。
 これは概周期関数とよばれる関数の一種で、おおよそ周期的だけど、厳密には周期関数ではありません。様々な深さの窪み(極小値)があるので、テストにぴったりです。

Python 概周期関数のグラフ

 ちなみに、このような奇妙な曲線は Exel VBA 数学実験室 というサイトにたくさん載っているので、興味のある方はぜひお立ち寄りください(←宣伝)。
 

テストコード

 それでは下のコードを動かして、この関数の最小値を求めてみましょう。

# https://python.atelierkobato.com/radar/ ‎
# リストR1-A-1

# 探査レーダー付勾配降下法のテストコード
# Gradient Descent with Search Radar(GDSR)

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

# 関数を定義
def func(x):
    return np.sin(x) + np.sin(np.sqrt(x))

# 導関数
def dfunc(x):
    return np.cos(x) + 0.5 * np.cos(np.sqrt(x)) / np.sqrt(x)

# ★★★★★★★★★★

# データを用意
x = np.linspace(0, 16*np.pi, 257)
x_max = np.max(x)
x_min = np.min(x)

# ★★★★★★★★★★

# データを可視化
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.set_title("y = sinx + sin√x", fontsize = 16)
ax.set_xlabel("x", fontsize = 16)
ax.set_ylabel("y", fontsize = 16)
ax.set_xlim([0, 16*np.pi])
ax.plot(x, func(x), color = "gray", zorder = 1)

# ★★★★★★★★★★

# Search Radar付勾配降下法

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

# ループ最大数
k_max = 1000

# xの初期値
x = 45

# レーダー1 ビーム長変動周波数
r = 0.2

# レーダー2 ビーム長初期設定
s = 1.0

# レーダー2 ビーム長変動定数
ds = 0.2

for k in range(1, k_max):
    s += ds
    x1 = x + (x_max - x_min) * np.sin(r * k)
    x2 = x + (-1)**k * s

    if func(x1) < func(x2):
        xt = x1
    else:
        xt = x2

    if func(xt) < func(x) and x_min < xt < x_max:
        s = 0
        x = xt

    x -= alpha * dfunc(x)
    ax.scatter(x, func(x), s = 5, color = "red", zorder = 2)
    
    # 勾配の大きさを計算
    d = abs(dfunc(x))

    # 勾配の大きさが収束判定条件以下ならループを抜ける
    if d < eps:
        break

# x,yの数値を丸める
x = np.round(x, 3)
y = np.round(func(x), 3)

print("k = {}".format(k))
print("(x, y) = ({}, {})".format(x, y))
k = 303
(x, y) = (23.557, -1.99)

Python 探査レーダーテスト修正版2

 303 ステップで最小値を発見することができました。
 レーダー 1 から射出されるビームの長さは
 
\[L_1=X\sin rk\]
で表されます。$X$ は現在考えている $x$ 軸の範囲の幅 ($16\pi$)、$k$ はステップ数、$r$ は「ステップごとのビームの長さの変化」を調整する周波数です。

 下の図にあるように、レーダー 2 のビームはステップごとに射出方向と長さを変えます。

 Python レーダーによる窪地の探査イラスト02

 ただし、転移が起こるとビームの長さは初期化されるようになっています。
 コードの中にある r, s, ds などは状況に応じて設定し直してください。
 ちなみに転移は現在いる窪みの中でも起こり得るので、勾配降下を加速させる効果も併せ持っています。

 Python 転移(テレポーテーション)

 後で紹介する 拡散式パワーレーダー を使ったほうが速いと思いますが、1 変数関数の場合はベーシックな探査レーダーで十分だと思います。ていうか、1 変数だと処理速度に大した差はないかもしれません。試したことないからわからないけど(←いい加減)。
 

今後の研究課題

 今後の研究課題としては、幅の狭い窪地をもつような難易度の高い関数でも成功するかどうかです。今のところ失敗例はありませんが、コメントで難しいテスト関数をお寄せくだされば、それで試してみて、成否に関わらず、また記事で発表したいと思います。レーダーの機能についても、まだまだ色々なものがありそうなので、皆さんから提案をお待ちしております。次回はレーダーを使って 2 変数関数の最小値を探します。