はじめに

　黄金比の定義や性質、黄金比を利用した図形やアルゴリズムについて、数式やプログラム、図を用いて理解を目指すシリーズです。

　この記事では、黄金分割探索について、図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
黄金分割探索の可視化
参考文献
おわりに

黄金分割探索の可視化

　黄金比(golden ratio)を用いた三分探索(ternary search)である黄金分割探索(golden section search)のアルゴリズムを図で確認します。
　黄金比については「黄金比の計算式の導出 - からっぽのしょこ」、黄金分割については「【R】黄金比の可視化 - からっぽのしょこ」、基本形である三分探索については「【Python】三分探索の可視化 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import numpy as np
import matplotlib.pyplot as plt

関数の設定

　まずは、最小値(または最大値)の探索対象とする関数を設定します。この記事では、最小値の探索のみを扱います。

　自作関数を作成します。

# 関数を作成
def eval_fnc(x):

    # 計算式を指定
    y = 0.5 * x**2 - 4.0 * x + 12
    return y

　下に凸の関数(計算式)を指定します。上に凸の関数の場合は、-y を出力すると(簡易的に)扱えます。
　この例では、 $f(x) = \frac{1}{2} x^2 - 4 x + 12$ とします。

　変数の範囲を指定して、曲線の座標を作成します。

# 変数の範囲を指定
x_min = -5.0
x_max = 15.0

# 曲線用の変数を作成
x_vec = np.linspace(start=x_min, stop=x_max, num=1001)

# 関数曲線を計算
f_x_vec = eval_fnc(x_vec)

　変数 $x$ の範囲を指定して、関数 $f(x)$ を計算します。

　関数のグラフを作成します。

# ラベル用の文字列を作成
fnc_label = '$f(x) = 0.5 x^2 - 4 x + 12$'

# 関数曲線を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.plot(x_vec, f_x_vec, color='black', label=fnc_label) # 関数
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
fig.suptitle('convex function', fontsize=20)
ax.grid()
ax.legend()
plt.show()

　この関数(曲線) $f(x)$ の最小値となる変数 $x$ を求めます。

　ちなみにこの例の関数の場合は、平方完成により解析的に求まります。

式変形(クリックで展開)

$\displaystyle \begin{aligned} f(x) &= \frac{1}{2} x^2 - 4 x + 12 \\ &= \frac{1}{2} ( x^2 - 8 x ) + 12 \\ &= \frac{1}{2} ( x^2 - 2 \cdot 4 x + 4^2 - 4^2 ) + 12 \\ &= \frac{1}{2} ( x^2 - 2 \cdot 4 x + 4^2 ) - 8 + 12 \\ &= \frac{1}{2} (x - 4)^2 + 4 \end{aligned}$

　 $x = 4$ のとき最小値 $f(x) = 4$ になります。

黄金分割の図

　次は、黄金分割の操作と性質を図で示します。

　分割比率を設定します。

# 黄金比の逆数を計算
ratio_g = 0.5 * (np.sqrt(5.0) - 1.0)
print(ratio_g)

0.6180339887498949

　黄金分割探索では、黄金比 $\phi = \frac{1 + \sqrt{5}}{2}$ の逆数 $\frac{1}{\phi} = \frac{\sqrt{5} - 1}{2}$ を用います。逆数は $-1$ 乗 $\phi^{-1} = \frac{1}{\phi}$ でも表わせます。

　1回目の試行における分割点を作成します。

# 探索範囲を設定
x_lower = x_min
x_upper = x_max

# 探索範囲を計算
x_size = x_upper - x_lower
print(x_size)

# 分割点を計算
x1 = x_lower + ratio_g * x_size
x2 = x_upper - ratio_g * x_size
print(x1)
print(x2)

20.0
7.360679774997898
2.639320225002102

　変数 $x$ の下限を $x_{\mathrm{lower}}$ 、上限を $x_{\mathrm{upper}}$ としたとき、探索範囲のサイズは $x_{\mathrm{upper}} - x_{\mathrm{lower}}$ です。
　探索範囲全体(上限と下限の差)を $1$ として、下限と上限それぞれからのサイズ比が $\frac{1}{\phi}$ になるような2点で3分割します。2点は次の式で計算できます。

$\displaystyle \begin{aligned} x_1 &= x_{\mathrm{lower}} + \frac{1}{\phi} (x_{\mathrm{upper}} - x_{\mathrm{lower}}) \\ x_2 &= x_{\mathrm{upper}} - \frac{1}{\phi} (x_{\mathrm{upper}} - x_{\mathrm{lower}}) \end{aligned}$

　2点は $x_2 \lt x_1$ の関係(図で言うと $x_1$ の点が右側)です。

　1回目の試行における黄金分割をグラフで確認します。

作図コード(クリックで展開)

# グラフサイズを設定
u = 10
y_max = np.ceil(f_x_vec.max() /u)*u # u単位で切り上げ
y_min = -0.2 * y_max

# 線分の表示位置の調整値を指定
d = 2.0

# ラベル用の文字列を作成
fml_label = '$\\phi^{-1} = \\frac{\\sqrt{5} - 1}{2} = 0.618...$'

# 黄金分割を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.plot(x_vec, f_x_vec, color='black') # 関数
ax.vlines(x=[x_min, x_max], ymin=y_min, ymax=y_max, 
          color='black', linestyle='dashed') # 探索範囲
ax.vlines(x=[x1, x2], ymin=y_min, ymax=y_max, 
          color=['red', 'purple'], linestyle='dashed') # 分割値
ax.hlines(y=0.0, xmin=x_min, xmax=x_max, 
          color='black', label='$1$') # 線分 1
ax.hlines(y=[-d, -d*2], xmin=[x_min, x2], xmax=[x1, x_max], 
          color='red', label='$\\phi^{-1}$') # 線分 1/φ
ax.hlines(y=[-d*3, -d*3], xmin=[x_min, x1], xmax=[x2, x_max], 
          color='purple', label='$\\phi^{-2}$') # 線分 1/φ^2
ax.hlines(y=-d*4, xmin=x2, xmax=x1, 
          color='blue', label='$\\phi^{-3}$') # 線分 1/φ^3
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title(fml_label, loc='left')
fig.suptitle('golden section', fontsize=20)
ax.grid()
ax.legend(loc='upper left')
ax.set_ylim(ymin=y_min, ymax=y_max)
plt.show()

　探索範囲における変数 $x$ の下限 $x_{\mathrm{lower}}$ と上限 $x_{\mathrm{upper}}$ を黒色の破線、大きい方の分割点 $x_1$ を赤色の破線、小さい方の分割点 $x_2$ を紫色の破線で示します。

　探索範囲全体のサイズ( $x_{\mathrm{lower}}$ から $x_{\mathrm{upper}}$ までのサイズ・黒色の実線)を $1$ とすると、下限または上限から大きい(遠い)方の分割サイズ( $x_{\mathrm{lower}}$ から $x_1$ までのサイズと $x_{\mathrm{upper}}$ から $x_2$ までのサイズ・赤色の実線)との比が $1 : \frac{1}{\phi}$ 、小さい(近い)方の分割サイズ( $x_{\mathrm{lower}}$ から $x_2$ までのサイズと $x_{\mathrm{upper}}$ から $x_1$ までのサイズ・紫色の実線)との比が $1 : \frac{1}{\phi^2}$ になります。
　また、分割点間のサイズ( $x_1$ から $x_2$ までのサイズ・青色の実線)との比は $1 : \frac{1}{\phi^3}$ になります。
　つまり、4つのサイズ比は大きい順に $1 : \frac{1}{\phi} : \frac{1}{\phi^2} : \frac{1}{\phi^3}$ であり、それぞれの比が $1 : \frac{1}{\phi}$ です。この関係は黄金比の性質によるものです。

　2回目の試行における分割点を作成します。

# 探索範囲の上限を更新
x_upper = x1
#x_lower = x2

# 探索範囲を計算
x_size = x_upper - x_lower
print(x_size)

# 分割点を計算
x3 = x_upper - ratio_g * x_size
#x3 = x_lower + ratio_g * x_size
print(x3)

12.360679774997898
-0.27864045000420656

　 $x_1 \gt x_2, f(x_1) \gt f(x_2)$ の場合は、上限を $x_1$ に更新 $x_{\mathrm{upper}} \leftarrow x_1$ します。この例だと、こちらの操作を行います。
　1試行目と同様に、更新した探索範囲全体を $1$ として、下限と上限それぞれからのサイズ比が $\frac{1}{\phi}$ になるような2点で3分割します。2点は先ほどの式で計算できます。

$\displaystyle \begin{aligned} x_2 &= x_{\mathrm{lower}} + \frac{1}{\phi} (x_{\mathrm{upper}} - x_{\mathrm{lower}}) \\ x_3 &= x_{\mathrm{upper}} - \frac{1}{\phi} (x_{\mathrm{upper}} - x_{\mathrm{lower}}) \end{aligned}$

　2点は $x_3 \lt x_2$ の関係(図で言うと $x_3$ の点が左側)です。

　 $x_1 \gt x_2, f(x_1) \lt f(x_2)$ の場合は、下限を $x_2$ に更新 $x_{\mathrm{lower}} \leftarrow x_2$ して、同様に分割します。

$\displaystyle \begin{aligned} x_3 &= x_{\mathrm{lower}} + \frac{1}{\phi} (x_{\mathrm{upper}} - x_{\mathrm{lower}}) \\ x_1 &= x_{\mathrm{upper}} - \frac{1}{\phi} (x_{\mathrm{upper}} - x_{\mathrm{lower}}) \end{aligned}$

　2点は $x_1 \lt x_3$ の関係(図で言うと $x_3$ の点が右側)です。

　どちらの場合でも、2点のうちの1点が更新しなかった方の点と一致します。つまり、更新されなかった点(関数値が小さい方の点)を再利用し、もう一方の点(式)のみを計算して探索を続けられます。

　2回目の試行における黄金分割をグラフで確認します。

作図コード(クリックで展開)

### 資料作成用:(注: 上限を更新した場合のみ対応)

# 線分の表示位置の調整値を指定
d = 2.0

# 黄金分割を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.plot(x_vec, f_x_vec, color='black') # 関数
ax.vlines(x=[x_min, x1], ymin=y_min, ymax=y_max, 
          color='black', linestyle='dashed') # 探索範囲
ax.vlines(x=[x2, x3], ymin=y_min, ymax=y_max, 
          color=['purple', 'blue'], linestyle='dashed') # 分割値
ax.hlines(y=0.0, xmin=x_min, xmax=x_max, 
          color='black', label='$1$') # 線分 1
ax.hlines(y=-d, xmin=x_min, xmax=x1, 
          color='red', label='$\\phi^{-1}$') # 線分 1/φ
ax.hlines(y=[-d*2, -d*3], xmin=[x_min, x3], xmax=[x2, x1], 
          color='purple', label='$\\phi^{-2}$') # 線分 1/φ^2
ax.hlines(y=[-d*4, -d*4], xmin=[x_min, x2], xmax=[x3, x1], 
          color='blue', label='$\\phi^{-3}$') # 線分 1/φ^3
ax.hlines(y=-d*5, xmin=x3, xmax=x2, 
          color='gold', label='$\\phi^{-4}$') # 線分 1/φ^4
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title(fml_label, loc='left')
fig.suptitle('golden section', fontsize=20)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.grid()
ax.legend(loc='upper left')
plt.show()

　この試行の探索範囲における変数 $x$ の下限 $x_{\mathrm{lower}}$ と上限 $x_{\mathrm{upper}} = x_1$ を黒色の破線、大きい方の分割点 $x_2$ を紫色の破線、小さい方の分割点 $x_3$ を青色の破線で示します。

　元の探索範囲のサイズ(黒色の実線)を $1$ とすると、探索範囲全体のサイズ( $x_{\mathrm{lower}}$ から $x_{\mathrm{upper}}$ までのサイズ・赤色の実線)との比は $1 : \frac{1}{\phi}$ です。
　下限または上限から大きい(遠い)方の分割サイズ( $x_{\mathrm{lower}}$ から $x_2$ までのサイズと $x_{\mathrm{upper}}$ から $x_3$ までのサイズ・紫色の実線)との比が $1 : \frac{1}{\phi^2}$ 、小さい(近い)方の分割サイズ( $x_{\mathrm{lower}}$ から $x_3$ までのサイズと $x_{\mathrm{upper}}$ から $x_2$ までのサイズ・青色の実線)との比が $1 : \frac{1}{\phi^3}$ になります。
　また、分割点間のサイズ( $x_2$ から $x_3$ までのサイズ・黄色の実線)との比は $1 : \frac{1}{\phi^4}$ になります。
　つまり、5つのサイズ比は大きい順に $1 : \frac{1}{\phi} : \frac{1}{\phi^3} : \frac{1}{\phi^3} : \frac{1}{\phi^4}$ であり、それぞれの比が $1 : \frac{1}{\phi}$ であることが変わりません。

　黄金比の性質を利用することで、三等分するよりも効率よく分割を行いながら探索を行えます。

黄金分割探索の図

　続いて、黄金分割探索の操作を図で示します。

　1回目の試行における曲線上の点の座標を作成します。

# 関数値を計算
f_x1 = eval_fnc(x1)
f_x2 = eval_fnc(x2)
print(f_x1)
print(f_x2)

9.647084275039962
4.925724725044166

　2つの分割点 $x_1, x_2$ の関数の値 $f(x_1), f(x_2)$ を計算します。

　1回目の試行における黄金分割探索をグラフで確認します。

作図コード(クリックで展開)

# グラフサイズを設定
u = 10
y_min = 0.0
y_max = np.ceil(f_x_vec.max() /u)*u # u単位で切り上げ

# ラベル用の文字列を作成
fnc_label = f'$f(x_1) = {f_x1:.5f}, f(x_2) = {f_x2:.5f}$'

# 黄金分割探索を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.plot(x_vec, f_x_vec, color='black') # 関数
ax.vlines(x=[x_min, x_max], ymin=y_min, ymax=y_max, 
          color='black', linestyle='dashed') # 探索範囲
ax.vlines(x=[x1, x2], ymin=y_min, ymax=y_max, 
          color=['red', 'purple'], linestyle='dashed') # 分割値
ax.scatter(x=x1, y=f_x1, 
           s=100, color='red', label='$(x_1, f(x_1))$') # 分割点 1
ax.scatter(x=x2, y=f_x2, 
           s=100, color='purple', label='$(x_2, f(x_2))$') # 分割点 2
ax.hlines(y=[f_x1, f_x2], xmin=x_min, xmax=[x1, x2], 
          color=['red', 'purple'], linestyle='dotted') # 分割点の値
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title(fnc_label, loc='left')
fig.suptitle('golden section search', fontsize=20)
ax.set_xlim(xmin=x_min, xmax=x_max)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.grid()
ax.legend(loc='upper left')
plt.show()

　2つの分割点 $x_1, x_2$ に対応する曲線上の点 $(x_1, f(x_1)), (x_2, f(x_2))$ を比較して、小さい方の点を残し、大きい方の点を次の試行の下限または上限(この例だと $x_1$ を上限)とします。

　2回目の試行における曲線上の点の座標を作成します。

# 関数値を計算
f_x3 = eval_fnc(x3)
print(f_x3)

13.1533820502061

　新たな分割点 $x_3$ の関数の値 $f(x_3)$ を計算します。

　2回目の試行における黄金分割探索をグラフで確認します。

作図コード(クリックで展開)

### 資料作成用:(注: 上限を更新した場合のみ対応)

# ラベル用の文字列を作成
fnc_label = f'$f(x_2) = {f_x2:.5f}, f(x_3) = {f_x3:.5f}$'

# 黄金分割を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.plot(x_vec, f_x_vec, color='black') # 関数
ax.vlines(x=[x_min, x1], ymin=y_min, ymax=y_max, 
          color='black', linestyle='dashed') # 探索範囲
ax.vlines(x=[x2, x3], ymin=y_min, ymax=y_max, 
          color=['purple', 'blue'], linestyle='dashed') # 分割値
ax.scatter(x=x2, y=f_x2, 
           s=100, color='purple', label='$(x_2, f(x_2))$') # 分割点 2
ax.scatter(x=x3, y=f_x3, 
           s=100, color='blue', label='$(x_3, f(x_3))$') # 分割点 3
ax.hlines(y=[f_x2, f_x3], xmin=x_min, xmax=[x2, x3], 
          color=['purple', 'blue'], linestyle='dotted') # 分割点の値
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_title(fnc_label, loc='left')
fig.suptitle('golden section search', fontsize=20)
ax.set_xlim(xmin=x_min, xmax=x_max)
ax.set_ylim(ymin=y_min, ymax=y_max)
ax.grid()
ax.legend(loc='upper left')
plt.show()

　1試行目と同様に、2つの分割点 $x_2, x_3$ に対応する曲線上の点 $(x_2, f(x_2)), (x_3, f(x_3))$ を比較して、小さい方の点を残し、大きい方の点を次の試行の下限または上限(この例だと $x_3$ を上限)とします。

　黄金比による分割と更新の操作を繰り返すことで、下に凸な関数の最小値を探索します。

探索の推移

　最後は、黄金分割探索の様子をアニメーションで示します。
　アニメーションの作図コードは「Mathematics/golden_ratio/code/golden_section_search.py at main · anemptyarchive/Mathematics · GitHub」を参照してください。

　黄金分割探索の推移をアニメーションで確認します。

　各試行における変数 $x$ の探索範囲の下限 $x_{\mathrm{lower}}$ と上限 $x_{\mathrm{upper}}$ を黒色の破線、分割点 $x_i, x_j$ をそれぞれ色付きの点と破線で示します。

　試行ごとに、 $f(x_i), f(x_j)$ の大きい方の点(変数の値)で下限または上限を更新して、小さい方の点(更新に使わなかった点)を再利用します。更新後の下限から上限までのサイズを $1$ として、更新した下限または上限からのサイズ比が $1 : \phi$ となる点で新たに分割します。
　この操作を繰り返すことで、変数の範囲を一定の比率で3分割しながら最小値の探索を行えます。

　この記事では、黄金分割探索を確認しました。

参考文献

Pythonで学ぶ空間データサイエンス入門: 地域の特徴を発見する方法

作者:植杉　大
Independently published

Amazon

おわりに

　何のきっかけだったかで黄金比についていくつかの内容を書いたのですが、それがここで役立つことになろうとは。今回は空間統計の中で登場した内容です。話の本筋からは逸れるので別シリーズの記事として書きましたが。別の文脈で勉強した内容が不意に繋がると何とも言えない快感があります。

　上限下限のどちらを更新しても再利用する点と新規の点の計算式や位置関係が試行ごとに入れ替わることに気付くまで混乱しました。

　2024年9月14日は、モーニング娘。の結成27周年の記念日です！

　卒業と加入を繰り返しながら進化していくグループことモーニング娘。！と書いて「最強」な曲を貼るなら最大化問題にしておけばよかったですね。
　今日は他にも記念日だったのですが記事を用意できませんでした。それどころか1か月ぶりくらいに更新できた記事です。色々とネタの準備はしているんですが手が回らないこの頃でままなりません。
　最器用になりたい、切実に。

【次の内容】

つづく

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【Python】黄金分割探索の可視化

はじめに

黄金分割探索の可視化

関数の設定

黄金分割の図

黄金分割探索の図

探索の推移

参考文献

おわりに