[研究レポート] 三次元空間ランダム座標立方格子

三次元空間におけるランダム座標の抽出方法

　早いもので（別にそんなに早くないかもしれないけど）、Python 研究レポートも 4 回目を迎えました。前回記事で扱った拡散式パワーレーダーは 3 変数以上の関数にもほとんどそのままの形で応用することができます。たとえば、3 変数関数 $w=f(x,y,z)$ の最小値を探索する場合を考えてみましょう。

立方格子生成関数

　最初に指定範囲内の空間から 無作為（ランダム）に選んだ座標を格納するための立方格子 を生成する r_mesh_3d() 関数を次のように定義します。

# リストR4-A-1

import numpy as np

# ランダム座標を格納する立方格子を生成する関数
def r_mesh_3d(xr, yr, zr, n):
    x = np.random.uniform(xr[0], xr[1], (n, n, n))
    y = np.random.uniform(yr[0], yr[1], (n, n, n))
    z = np.random.uniform(zr[0], zr[1], (n, n, n))
    X, Y, Z = np.array([x, y, z])
    return X, Y, Z

　r_mesh_3d() は 3 次元配列を 3 個返します。
　それを X, Y, Z という名前の変数に入れて X だけ表示してみます。

# リストR4-A-2

# x,y,zの範囲を指定
xr = [0, 10]
yr = [0, 10]
zr = [0, 10]

# 各軸に沿った格子の数
n = 3

# 立方格子を生成
X, Y, Z = r_mesh_3d(xr, yr, zr, n)

print(X)

[[[5.48813504 7.15189366 6.02763376]
  [5.44883183 4.23654799 6.45894113]
  [4.37587211 8.91773001 9.63662761]]

 [[3.83441519 7.91725038 5.2889492 ]
  [5.68044561 9.25596638 0.71036058]
  [0.871293   0.20218397 8.32619846]]

 [[7.78156751 8.70012148 9.78618342]
  [7.99158564 4.61479362 7.80529176]
  [1.18274426 6.39921021 1.43353287]]]

　Y, Z にも、それぞれ同じように 27 個のデータが入っているので、空間から 27 個の座標をランダムに選んでいることになります。
　

3変数関数の最小値

　3 変数関数
　
\[w=\cos(x+y-z)+\sin(2z)\]
の最小値を求めてみましょう。関数の形から最小値が $-2$ であることは自明ですが、3 変数以上ではもはや等高線を描けないので、予め答えの分かっている問題を設定してプログラムが正しく機能するかどうかを確認することも大切です。あ、それから、拡散式パワーレーダーではランダムな座標を選んで移動していくので、初期地点の選び方はほとんど影響しません。

# リストR4-A-3

def func_3v(x, y, z):
    return np.cos(x + y - z) + np.sin(2 * z)

def gfunc_3v(x, y, z):
    gx = - np.sin(x + y -z)
    gy = - np.sin(x + y -z)
    gz = np.sin(x * y * z) + 2 * np.cos(2 * z)
    return np.array([gx, gy, gz])

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

# x,y,zの範囲を指定
xr = [0, 10]
yr = [0, 10]
zr = [0, 10]

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

# 初期地点の高度
wt = func_3v(pt[0], pt[1], pt[2])

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

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

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

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

    if t % fq == 0:
        Xs, Ys, Zs = r_mesh_3d(xr, yr, zr, 32)
        Ws = func_3v(Xs, Ys, Zs)
        Ws_min = np.min(Ws)
        if wt > Ws_min:
            Xs_min = Xs[Ws == Ws_min][0]
            Ys_min = Ys[Ws == Ws_min][0]
            Zs_min = Zs[Ws == Ws_min][0]
            pt = [Xs_min, Ys_min, Zs_min]

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

wt = func_3v(pt[0], pt[1], pt[2])
pt = np.round(pt, 3)
wt = np.round(wt, 3)

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

t = 1188
(x, y, z, w) = (0.279, 5.214, 8.636, -2.000)

　きちんと最小値が求まりましたね。ちなみに次元が増えると探索範囲はぐぐーっと広くなります。二次元世界（平面）から三次元世界（空間）へぐぐーっと広くなるように、三次元から四次元に変わるときも、ぐぐぐーっと広くなります（←頭が悪いから、こんな表現しかできない）。そんなわけで、探索時間もぐぐぐぐーっと増加します。私の環境では 3.26 秒かかりました。それにしても 3 変数以上になると、Matplotlib で探索過程を可視化できないので少し寂しい気もしますね。