ブログ

直線探索(line search)の概要とPythonを用いた計算例

最適化問題を解く際に勾配法などがよく用いられますが、多次元空間上における勾配法などを用いる際には取り得る範囲内でどのステップ幅が良いかを検討する場合があります。この際によく用いられるのが直線探索(line search)です。当記事では直線探索の流れとPythonを用いた計算例をまとめました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

直線探索の概要

目的関数と勾配降下法

目的関数$f(\mathbf{x}_{k})$の$p$次元ベクトル$\mathbf{x}_{k}$に関する勾配降下法(Gradient Descent Method)の漸化式は下記のように表される。
$$
\large
\begin{align}
\mathbf{x}_{k+1} &= \mathbf{x}_{k} – \alpha_{k} \nabla f(\mathbf{x}_{k}) \quad (1) \\
\mathbf{x}_{k} & \in \mathbb{R}^{p} \\
\nabla &= \left(\begin{array}{c} \displaystyle \frac{\partial}{\partial x_1} \\ \vdots \\ \displaystyle \frac{\partial}{\partial x_p} \end{array} \right)
\end{align}
$$

直線探索

前項の「目的関数と勾配降下法」の$(1)$式における$\alpha_{k}$には定数や減衰などを表す関数などを用いることが多いが、多変数関数などを取り扱う際の伝統的な手法の直線探索(line search)も問題によっては役立つ。

直線探索は勾配ベクトルが得られた際に、勾配の方向の目的関数$f(\mathbf{x}_{k})$の値を計算し、最適な値を持つ点のステップ幅を選択するという手法である。原理自体はシンプルである一方でプログラムはやや複雑になるので、次節ではPythonを用いた直線探索の計算例を取り扱った。

直線探索の実装方針

直線探索を実装するにあたっては、探索範囲を元に$\alpha_{k}$の取りうる値を配列で作成し、どの値を用いた際に目的関数が最適値を取るかを調べればよい。たとえば目的関数$f(x)=x^2$の最小化にあたって、$x_k=-2$の際を仮定する。

このとき、$f'(x)=2x$であるので、勾配降下法の式は下記のように表される。
$$
\large
\begin{align}
x_{k+1} &= x_{k} – \alpha_{k} f'(x_k) \\
&= -2 + 4 \alpha_{k}
\end{align}
$$

上記に対し下記を実行することで、$0$〜$1$の$0.01$刻みで$\alpha_{k}$を生成し、どの値が目的関数$f(x)=x^2$を最大化するか調べることができる。

import numpy as np

alpha = np.arange(0,1.01,0.01)
x = -2 + 4*alpha
f_x = x_new**2

x_new = x[np.argmin(f_x)]
print("x_new: {:.1f}".format(x_new))

・実行結果

x_new: 0.0

上記では$0.01$刻みで生成した$\alpha_{k}$に対し、取りうる$x_{k+1}$を計算し、$f(x)$を最小にするインデックスをnp.argminを用いて取得した。最小値における詳細は下記を出力することで確認できる。

print(np.argmin(f_x))
print(alpha[np.argmin(f_x)])
print(f_x[np.argmin(f_x)])

・実行結果

50
0.5
0.0

Pythonを用いた直線探索の計算例

問題設定

$$
\large
\begin{align}
f(\mathbf{x}) &= \mathbf{x}^{\mathrm{T}} A \mathbf{x} \\
A &= \left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \\
\mathbf{x} & \in \mathbb{R}^{2}
\end{align}
$$

上記のように目的関数$f(\mathbf{x})$を定義する。このとき、等高線上の点は下記のような計算で得ることができる。
$$
\large
\begin{align}
U \mathbf{x} &= \frac{1}{\sqrt{2}} \left(\begin{array}{cc} 1 & -1 \\ 1 & 1 \end{array} \right) \left(\begin{array}{c} \displaystyle \frac{5}{\sqrt{3}} \cos{\frac{\pi}{3}} \\ \displaystyle 5 \sin{\frac{\pi}{3}} \end{array} \right) \\
&= \frac{5}{2 \sqrt{6}} \left(\begin{array}{cc} 1 & -1 \\ 1 & 1 \end{array} \right) \left(\begin{array}{c} 1 \\ 3 \end{array} \right) \\
&= \frac{5}{\sqrt{6}} \left(\begin{array}{c} -1 \\ 2 \end{array} \right) = \left(\begin{array}{c} -2.041 \cdots \\ 4.082 \cdots \end{array} \right)
\end{align}
$$

上記の点から勾配法と共役勾配法によって最小値の探索を行う。また、下記を実行することで点の描画を行うことができる。

import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0, 2*np.pi, 100)
c = np.arange(1,6)

x = np.zeros([2,len(theta)])
U = np.array([[1,-1],[1,1]])/np.sqrt(2)
ratio = 1/np.sqrt(3)

plt.figure(figsize=(5,5))
for i in range(len(c)):
    x[0,:] = c[i] * ratio * np.cos(theta)
    x[1,:] = c[i] * np.sin(theta)
    y = np.dot(U,x)
    plt.plot(y[0,:], y[1,:], "g--")

init_pos = np.array([5*ratio*np.cos(np.pi/6), 5*np.sin(np.pi/6)])
init_rotated = np.dot(U,init_pos)

plt.plot([init_rotated[0]],[init_rotated[1]],"ro")

plt.xlim([-5, 5])
plt.ylim([-5, 5])

plt.show()

・実行結果

当節では下記の事例を改変し、作成を行なった。

勾配法と直線探索を用いた計算例

下記を実行することで直線探索を行うことができる。

A = np.array([[2, 1], [1, 2]])
U = np.array([[1,-1],[1,1]])/np.sqrt(2)
ratio = 1/np.sqrt(3)

plt.figure(figsize=(5,5))
for i in range(len(c)):
    x = np.zeros([2,len(theta)])
    x[0,:] = c[i] * ratio * np.cos(theta)
    x[1,:] = c[i] * np.sin(theta)
    y = np.dot(U,x)
    plt.plot(y[0,:], y[1,:], "g--")
    
init_pos = np.array([5*ratio*np.cos(np.pi/3), 5*np.sin(np.pi/3)])
init_rotated = np.dot(U,init_pos)
alpha = np.arange(-5.0,5.01,0.01)

tangent_l = np.dot(U,np.array([init_pos[0]+alpha, init_pos[1]-alpha*np.sqrt(3)/np.tan(np.pi/3)]))
grad = np.dot(U,np.array([init_pos[0]+alpha, init_pos[1]+alpha*np.tan(np.pi/3)/np.sqrt(3)]))

f_x = np.zeros(len(alpha))
for i in range(len(alpha)):
    f_x[i] = np.dot(grad[:,i],np.dot(A,grad[:,i]))
idx = np.argmin(f_x)
x_new = grad[:,idx]

plt.plot(tangent_l[0,:], tangent_l[1,:], "k--")
plt.plot(grad[0,:], grad[1,:], "b--")

plt.plot([init_rotated[0]],[init_rotated[1]],"ro")
plt.plot([x_new[0]],[x_new[1]],"ro")

plt.xlim([-5, 5])
plt.ylim([-5, 5])

plt.show()

・実行結果

共役勾配法と直線探索を用いた計算例

下記を実行することで共役勾配法における直線探索を行うことができる。

A = np.array([[2,1],[1,2]])
U = np.array([[1,-1],[1,1]])/np.sqrt(2)

x = np.zeros([2,len(theta)])
ratio = 1/np.sqrt(3)

plt.figure(figsize=(5,5))
for i in range(len(c)):
    x[0,:] = c[i] * ratio * np.cos(theta)
    x[1,:] = c[i] * np.sin(theta)
    y = np.dot(U,x)
    plt.plot(y[0,:], y[1,:], "g--")
    
init_pos = np.array([5*ratio*np.cos(np.pi/3), 5*np.sin(np.pi/3)])
init_rotated = np.dot(U,init_pos)
alpha = np.arange(-5.0,5.01,0.01)

tangent_l = np.dot(U,np.array([init_pos[0]+alpha, init_pos[1]-alpha*np.sqrt(3)/np.tan(np.pi/3)]))
tan, grad = np.dot(U,np.array([1, -np.sqrt(3)/np.tan(np.pi/3)])), np.dot(U,np.array([1, np.tan(np.pi/3)/np.sqrt(3)]))
conjugate_grad = grad - np.dot(tan,np.dot(A,grad))/np.dot(tan,np.dot(A,tan))*tan
conjugate_grad_x = init_rotated[0] + alpha*conjugate_grad[0]
conjugate_grad_y = init_rotated[1] + alpha*conjugate_grad[1]

f_x = np.zeros(len(alpha))
for i in range(len(alpha)):
    f_x[i] = 2*(conjugate_grad_x[i]**2 + conjugate_grad_x[i]*conjugate_grad_y[i] + conjugate_grad_y[i]**2)
idx = np.argmin(f_x)
x_new = np.array([conjugate_grad_x[idx],conjugate_grad_y[idx]])

plt.plot(tangent_l[0,:], tangent_l[1,:], "k--")
plt.plot(conjugate_grad_x, conjugate_grad_y, "b--")

plt.plot([init_rotated[0]],[init_rotated[1]],"ro")
plt.plot([x_new[0]],[x_new[1]],"ro")


plt.xlim([-5, 5])
plt.ylim([-5, 5])

plt.show()

・実行結果

Hessian-free optimizationにおける二次形式の計算の軽量化

共役勾配法などにおける行列にヘッセ行列(Hessian Matrix)を用いる場合、ニューラルネットワークのようにパラメータが多い場合はヘッセ行列の要素が多いことで計算が難しくなります。このような際にHessian-free optimizationが有用なので、当記事では概要について取りまとめました。
“Training Deep and Recurrent Networks with Hessian-Free Optimization”の内容を参考に作成を行いました。
https://www.cs.toronto.edu/~jmartens/docs/HF_book_chapter.pdf

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

前提の理解

ヘッセ行列の定義

$p$次元のベクトル$\displaystyle \mathbf{x} = \left(\begin{array}{c} x_1 \\ \vdots \\ x_p \end{array} \right)$に関するスカラー関数$f(\mathbf{x})$のヘッセ行列$H(f)$は下記のように定義される。
$$
\large
\begin{align}
H(f) &= \nabla \otimes \nabla f(\mathbf{x}) \\
&= \left(\begin{array}{c} \displaystyle \frac{\partial}{\partial x_1} \\ \vdots \\ \displaystyle \frac{\partial}{\partial x_p} \end{array} \right) \left(\begin{array}{ccc} \displaystyle \frac{\partial f}{\partial x_1} & \cdots & \displaystyle \frac{\partial f}{\partial x_p} \end{array} \right) \\
&= \left(\begin{array}{ccc} \displaystyle \frac{\partial^2 f}{\partial x_1^2} & \cdots & \displaystyle \frac{\partial^2 f}{\partial x_1 \partial x_p} \\ \vdots & \ddots & \vdots \\ \displaystyle \frac{\partial^2 f}{\partial x_p \partial x_1} & \cdots & \displaystyle \frac{\partial^2 f}{\partial x_p^2} \end{array} \right)
\end{align}
$$

ヘッセ行列について詳しくは下記で取り扱った。

ヘッセ行列の活用例:共役勾配法

共役勾配法では下記のような式に基づいて最適化を行う際にパラメータを移動させる方向を計算する。
$$
\large
\begin{align}
\mathbf{d}_{2} = \mathbf{v}_{2} – \frac{\mathbf{v}_{1}^{\mathrm{T}} A \mathbf{v}_{2}}{\mathbf{v}_{1}^{\mathrm{T}} A \mathbf{v}_{1}} \mathbf{v}_{1} \quad (1)
\end{align}
$$

$(1)$式の詳細は「共役勾配法の原理と具体的な計算例」、「共役勾配法のPythonを用いた計算例」などで取り扱った。ここで$(1)$式の対称行列$A$は数値で与えられることもあるが、多変数関数の最適化などの場合はヘッセ行列$H$が用いられる。

ヘッセ行列と多変数関数のテイラー二次近似

多変数関数$f(\mathbf{x})$の$\mathbf{x}$における近傍$\mathbf{x} + \varepsilon\mathbf{v}$について、下記のような一次のテイラー近似を行うことができる。
$$
\large
\begin{align}
f(\mathbf{x} + \varepsilon \mathbf{v}) & \simeq f(\mathbf{x}) + \frac{\varepsilon \nabla f(\mathbf{x})^{\mathrm{T}} \mathbf{v}}{1!} \quad (2) \\
\mathbf{x} & \in \mathbb{R}^{p}, \, \mathbf{v} \in \mathbb{R}^{p}
\end{align}
$$

同様な表記は下記でも取り扱った。

Hessian-free optimization

問題設定

ニューラルネットワークなどの多変量のパラメータベクトルを$\boldsymbol{\theta} \in \mathbb{R}^{p}$、尤度などに基づく目的関数を$f(\boldsymbol{\theta})$のように定義する。このとき勾配ベクトル$\nabla f(\boldsymbol{\theta})$、ヘッセ行列$H$は下記のように表せる。
$$
\large
\begin{align}
\nabla f(\boldsymbol{\theta}) &= \left(\begin{array}{c} \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_1} \\ \vdots \\ \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_p} \end{array} \right) \\
H &= \left(\begin{array}{ccc} \displaystyle \frac{\partial^2 f}{\partial \theta_1^2} & \cdots & \displaystyle \frac{\partial^2 f}{\partial \theta_1 \partial \theta_p} \\ \vdots & \ddots & \vdots \\ \displaystyle \frac{\partial^2 f}{\partial \theta_p \partial \theta_1} & \cdots & \displaystyle \frac{\partial^2 f}{\partial \theta_p^2} \end{array} \right)
\end{align}
$$

Hessian-vector productの近似

前節の$(2)$式に前項の「問題設定」を適用すると下記のように表せる。
$$
\large
\begin{align}
f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) & \simeq f(\boldsymbol{\theta}) + \frac{\varepsilon \nabla f(\boldsymbol{\theta})^{\mathrm{T}} \mathbf{u}}{1!} \quad (3)
\end{align}
$$

$(3)$式は下記のように変形できる。
$$
\large
\begin{align}
f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) & \simeq f(\boldsymbol{\theta}) + \frac{\varepsilon \nabla f(\boldsymbol{\theta})^{\mathrm{T}} \mathbf{u}}{1!} \quad (3) \\
\nabla f(\boldsymbol{\theta})^{\mathrm{T}} \mathbf{u} & \simeq \frac{f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – f(\boldsymbol{\theta})}{\varepsilon} \\
\left(\begin{array}{c} \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_1} & \cdots & \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_p} \end{array} \right) \mathbf{u} & \simeq \frac{f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – f(\boldsymbol{\theta})}{\varepsilon} \quad (4)
\end{align}
$$

ここで$(4)$式の両辺に対し、$\nabla$を計算すると下記が得られる。
$$
\large
\begin{align}
\left(\begin{array}{ccc} \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_1} & \cdots & \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_p} \end{array} \right) \mathbf{u} & \simeq \frac{f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – f(\boldsymbol{\theta})}{\varepsilon} \quad (4) \\
\nabla \left(\begin{array}{ccc} \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_1} & \cdots & \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_p} \end{array} \right) \mathbf{u} & \simeq \nabla \left( \frac{f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – f(\boldsymbol{\theta})}{\varepsilon} \right) \\
\left(\begin{array}{c} \displaystyle \frac{\partial}{\partial \theta_1} \\ \vdots \\ \displaystyle \frac{\partial}{\partial \theta_p} \end{array} \right) \left(\begin{array}{ccc} \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_1} & \cdots & \displaystyle \frac{\partial f(\boldsymbol{\theta})}{\partial \theta_p} \end{array} \right) & \simeq \frac{\nabla f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – \nabla f(\boldsymbol{\theta})}{\varepsilon} \\
\left(\begin{array}{ccc} \displaystyle \frac{\partial^2 f}{\partial \theta_1^2} & \cdots & \displaystyle \frac{\partial^2 f}{\partial \theta_1 \partial \theta_p} \\ \vdots & \ddots & \vdots \\ \displaystyle \frac{\partial^2 f}{\partial \theta_p \partial \theta_1} & \cdots & \displaystyle \frac{\partial^2 f}{\partial \theta_p^2} \end{array} \right) \mathbf{u} & \simeq \frac{\nabla f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – \nabla f(\boldsymbol{\theta})}{\varepsilon} \\
H \mathbf{u} & \simeq \frac{\nabla f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – \nabla f(\boldsymbol{\theta})}{\varepsilon} \quad (5)
\end{align}
$$

$(5)$式に対し、$\varepsilon \to 0$を考えると下記のように表せる。
$$
\large
\begin{align}
H \mathbf{u} = \lim_{\varepsilon \to 0} \frac{\nabla f(\boldsymbol{\theta} + \varepsilon \mathbf{u}) – \nabla f(\boldsymbol{\theta})}{\varepsilon} \quad (6)
\end{align}
$$

$(6)$式を用いて$(1)$式の$A \mathbf{v}_{1}, \, A \mathbf{v}_{2}$を計算すれば、ヘッセ行列を計算せずにヘッセ行列とベクトルの積を近似することができる。

Ch.27 「ベクトルの微分と条件付き極値問題」の演習問題の解答例〜統計学のための数学入門30講〜

当記事は「統計学のための数学入門$30$講(朝倉書店)」の読解サポートにあたってChapter.$27$の「ベクトルの微分と条件付き極値問題」の章末問題の解答の作成を行いました。
基本的には書籍の購入者向けの解説なので、まだ入手されていない方は購入の上ご確認ください。また、解説はあくまでサイト運営者が独自に作成したものであり、書籍の公式ページではないことにご注意ください。

・書籍解答まとめ
https://www.hello-statisticians.com/answer_textbook_math#math_stat

本章のまとめ

ベクトルによる偏微分

条件付き極値問題

演習問題解答

問題$27.1$

下記のようにラグランジュ関数$L(\mathbf{x},\lambda)$を定義する。
$$
\large
\begin{align}
L(\mathbf{x},\lambda) = \mathbf{x}^{\mathrm{T}} A \mathbf{x} – \lambda(\mathbf{x}^{\mathrm{T}}\mathbf{x}-1)
\end{align}
$$

このとき、極値の必要条件である$\nabla L(\mathbf{x},\lambda)=0$は下記のように変形できる。
$$
\large
\begin{align}
\nabla L(\mathbf{x},\lambda) &= 0 \\
2 A \mathbf{x} – 2 \lambda \mathbf{x} &= 0 \\
A \mathbf{x} &= \lambda \mathbf{x} \quad (1)
\end{align}
$$

上記の$(1)$式に対し、左から$\mathbf{x}^{\mathrm{T}}$をかけると下記が得られる。
$$
\large
\begin{align}
\mathbf{x}^{\mathrm{T}} A \mathbf{x} &= \mathbf{x}^{\mathrm{T}} (\lambda \mathbf{x}) \\
\mathbf{x}^{\mathrm{T}} A \mathbf{x} &= \lambda \mathbf{x}^{\mathrm{T}} \mathbf{x} \\
\mathbf{x}^{\mathrm{T}} A \mathbf{x} &= \lambda \quad (2)
\end{align}
$$

$(2)$式より二次形式$\mathbf{x}^{\mathrm{T}} A \mathbf{x}$の最大値は対称行列$A$の最大固有値、最小値は最小固有値に対応することが確認できる。
$$
\large
\begin{align}
A = \left(\begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{array} \right)
\end{align}
$$

上記の$A$の固有方程式の$\det(A – \lambda I_{3})=0$は下記のように解くことができる。
$$
\large
\begin{align}
\det(A – \lambda I_{3}) &= 0 \\
\left| \begin{array}{ccc} 1-\lambda & 0 & 1 \\ 0 & 1-\lambda & 0 \\ 1 & 0 & 1-\lambda \end{array} \right| &= 0 \\
(-1)^{1+1}(1-\lambda) \left| \begin{array}{cc} 1-\lambda & 0 \\ 0 & 1-\lambda \end{array} \right| + (-1)^{1+3} \left| \begin{array}{cc} 0 & 1 \\ 1-\lambda & 0 \end{array} \right| &= 0 \\
(1-\lambda)^{3} – (1-\lambda) &= 0 \\
(1-\lambda)(\cancel{1} – 2\lambda + \lambda^2 – \cancel{1}) &= 0 \\
-\lambda(\lambda-1)(\lambda-2) &= 0
\end{align}
$$

よって$(1)$式の最大値固有値は$2$、最小値固有値は$0$であることが確認できる。

このとき、$\lambda=2$に対応する固有ベクトル$\displaystyle \frac{1}{\sqrt{2}} \left( \begin{array}{c} 1 \\ 0 \\ 1 \end{array} \right)$、$\lambda=0$に対応する固有ベクトル$\displaystyle \frac{1}{\sqrt{2}} \left( \begin{array}{c} -1 \\ 0 \\ 1 \end{array} \right)$について$\mathbf{x}^{\mathrm{T}}\mathbf{x}=1$が成立し、かつ$\mathbf{x}^{\mathrm{T}} A \mathbf{x}=2, \, \mathbf{x}^{\mathrm{T}} A \mathbf{x}=0$がそれぞれ確認できる。よって得られた結果について十分性も成立する。

問題$27.2$

共役勾配法(Conjugate Gradient Method)のPythonを用いた計算例

共役勾配法(Conjugate Gradient Method)は等高線が同心楕円で表される場合の最適化にあたって有用な手法です。当記事では具体的な二次形式に対して共役勾配法を元に最適化を行う流れをPythonを用いて計算やグラフの描画を行いました。
当記事は「共役勾配法:日本オペレーション・リサーチ学会 $1987$年$6$月号」などを参考に作成を行いました。
https://orsj.org/wp-content/or-archives50/pdf/bul/Vol.32_06_363.pdf

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

二次形式の等高線のグラフ化

二次形式の式定義

$$
\large
\begin{align}
f(\mathbf{x}) &= \mathbf{x}^{\mathrm{T}} A \mathbf{x} \quad (1) \\
A &= \left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \\
\mathbf{x} &= \left(\begin{array}{c} x_1 \\ x_2 \end{array} \right)
\end{align}
$$

$(1)$式のように定義された二次形式$f(\mathbf{x}) = \mathbf{x}^{\mathrm{T}} A \mathbf{x}$の最小点を以下、共役勾配法を用いて計算を行う。

二次形式の図形的解釈

行列$A$の固有値と長さ$1$の固有ベクトルを$\lambda_{1}, \lambda_{2}$、$\mathbf{e}_{1}, \mathbf{e}_{2}$とおくと、それぞれ下記のように得られる。
$$
\large
\begin{align}
\lambda_{1} &= 3, \, \mathbf{e}_{1} = \frac{1}{\sqrt{2}} \left(\begin{array}{c} 1 \\ 1 \end{array} \right) \\
\lambda_{2} &= 1, \, \mathbf{e}_{2} = \frac{1}{\sqrt{2}} \left(\begin{array}{c} -1 \\ 1 \end{array} \right)
\end{align}
$$

上記が正しいことは$A \mathbf{e}_{i} = \lambda_{i} \mathbf{e}_{i}$が成立することから確認できる。ここで任意のベクトル$\mathbf{x}$を係数$a_1, a_2$を用いて下記のように表すことができる。
$$
\large
\begin{align}
\mathbf{x} = a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2} \quad (2)
\end{align}
$$

このとき$(1)$式に$(2)$式を代入すると下記が得られる。
$$
\large
\begin{align}
f(\mathbf{x}) &= \mathbf{x}^{\mathrm{T}} A \mathbf{x} \\
&= (a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2})^{\mathrm{T}} A (a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2}) \\
&= (a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2})^{\mathrm{T}} (\lambda_1 a_1 \mathbf{e}_{1} + \lambda_2 a_2 \mathbf{e}_{2}) \\
&= \lambda_{1} a_{1}^{2} + \lambda_{2} a_{2}^{2} \\
&= 3 a_{1}^{2} + a_{2}^{2} \quad (3)
\end{align}
$$

上記の計算にあたっては$\mathbf{e}_{i}^{\mathrm{T}} \mathbf{e}_{i}=1, \mathbf{e}_{1}^{\mathrm{T}} \mathbf{e}_{2}=0$が成立することを用いた。ここで$f(\mathbf{x})=C$の等高線は$(3)$式より下記のように得られる。
$$
\large
\begin{align}
f(\mathbf{x}) = 3 a_{1}^{2} + a_{2}^{2} &= C \\
\frac{a_{1}^{2}}{\sqrt{C/3}^{2}} + \frac{a_{2}^{2}}{\sqrt{C}^{2}} &= 0 \quad (4)
\end{align}
$$

上記は楕円の方程式に一致する。ここで$(2)$式より、$a_1, a_2$は基底$\mathbf{e}_{1}, \mathbf{e}_{2}$に対応するが、$\mathbf{e}_{1}, \mathbf{e}_{2}$は$\displaystyle \left(\begin{array}{c} 1 \\ 0 \end{array} \right), \, \left(\begin{array}{c} 0 \\ 1 \end{array} \right)$をそれぞれ$\displaystyle \frac{\pi}{4}$回転させたベクトルに一致する。

よって、$(4)$式に基づいて楕円を描き、原点を中心に$\displaystyle \frac{\pi}{4}$回転させることで楕円の等高線の描画を行うことができる。次項でPythonを用いたグラフ化について取り扱った。

楕円の描画

$$
\large
\begin{align}
\frac{x^2}{5^2} + \frac{y^2}{3^2} = 1 \quad (5)
\end{align}
$$

上記の$(5)$式のような楕円の方程式が得られた際に等高線の描画を行うにあたっては媒介変数表示を用いればよい。具体的には$x = 5 \cos{\theta}, \, y = 3 \sin{\theta}$が下記のように$(5)$について成立することを活用する。
$$
\large
\begin{align}
\frac{x^2}{5^2} + \frac{y^2}{3^2} &= \frac{\cancel{5^2} \cos^{2}{\theta}}{\cancel{5^2}} + \frac{\cancel{3^2} \sin^{2}{\theta}}{\cancel{3^2}} \\
&= \cos^{2}{\theta} + \sin^{2}{\theta} = 1
\end{align}
$$

ここまでの内容を元に、$(5)$式は下記を実行することで描画できる。

import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0, 2*np.pi, 100)

x = 5 * np.cos(theta)
y = 3 * np.sin(theta)

plt.figure(figsize=(5,5))
plt.plot(x, y, "g-")

plt.xlim([-7, 7])
plt.ylim([-7, 7])

plt.show()

・実行結果

上記の描画に回転を組み合わせることで$(4)$式が描画できる。$C=1,2,3,4,5$について、下記を実行することで同心楕円を描くことができる。

theta = np.linspace(0, 2*np.pi, 100)
c = np.arange(1,6)

x = np.zeros([2,len(theta)])
U = np.array([[1,-1],[1,1]])/np.sqrt(2)
ratio = 1/np.sqrt(3)

plt.figure(figsize=(5,5))
for i in range(len(c)):
    x[0,:] = c[i] * ratio * np.cos(theta)
    x[1,:] = c[i] * np.sin(theta)
    y = np.dot(U,x)
    plt.plot(y[0,:], y[1,:], "g--")

plt.xlim([-5, 5])
plt.ylim([-5, 5])

plt.show()

・実行結果

回転行列を用いた計算にあたっては下記で詳しく取り扱った。

共役勾配法の実行

同心楕円の接線と勾配

同心楕円上の点は下記のように$\displaystyle \mathbf{x} = \left(\begin{array}{c} \displaystyle \sqrt{\frac{C}{3}} \cos{\theta} \\ \sqrt{C} \sin{\theta} \end{array} \right)$を原点を中心に$\displaystyle \frac{\pi}{4}$回転させることで得られる。
$$
\large
\begin{align}
R \left( \frac{\pi}{4} \right) \mathbf{x} = \left(\begin{array}{cc} 1 & -1 \\ 1 & 1 \end{array} \right) \frac{1}{\sqrt{2}} \left(\begin{array}{c} \displaystyle \sqrt{\frac{C}{3}} \cos{\theta} \\ \sqrt{C} \sin{\theta} \end{array} \right)
\end{align}
$$

上記の式に基づいて、同心楕円上の点は下記を実行することで得ることができる。

x = np.zeros([2,len(theta)])
U = np.array([[1,-1],[1,1]])/np.sqrt(2)
ratio = 1/np.sqrt(3)

plt.figure(figsize=(5,5))
for i in range(len(c)):
    x[0,:] = c[i] * ratio * np.cos(theta)
    x[1,:] = c[i] * np.sin(theta)
    y = np.dot(U,x)
    plt.plot(y[0,:], y[1,:], "g--")
    
init_pos = np.dot(U,x[:,5])

plt.plot([init_pos[0]],[init_pos[1]],"ro")

plt.xlim([-5, 5])
plt.ylim([-5, 5])

plt.show()

・実行結果

また、$\displaystyle \left(\begin{array}{c} \displaystyle \sqrt{\frac{C}{3}} \cos{\theta} \\ \sqrt{C} \sin{\theta} \end{array} \right)$における接線に平行なベクトルを$\mathbf{v}_{1}$、$\mathbf{v}_{1}$の傾きを$\Delta_{1}$、法線ベクトルを$\mathbf{v}_{2}$とおくと、$\Delta_{1}$は下記のように表すことができる。
$$
\large
\begin{align}
\Delta_{1} &= \frac{dx_2}{dx_1} = \frac{\displaystyle \frac{dx_2}{d\theta}}{\displaystyle \frac{dx_1}{d\theta}} \\
&= \frac{\displaystyle \sqrt{C} \cos{\theta}}{\displaystyle -\sqrt{\frac{C}{3}} \sin{\theta}} \\
&= -\frac{\sqrt{3}}{\tan{\theta}}
\end{align}
$$

ここで$\mathbf{v}_{1} \perp \mathbf{v}_{2}$より$\mathbf{v}_{1} \cdot \mathbf{v}_{2}=0$であるので、$\mathbf{v}_{1}$と$\mathbf{v}_{2}$は下記のように得られる。
$$
\large
\begin{align}
\mathbf{v}_{1} &= \left(\begin{array}{c} 1 \\ \displaystyle -\frac{\sqrt{3}}{\tan{\theta}} \end{array} \right) \quad (6) \\
\mathbf{v}_{2} &= \left(\begin{array}{c} 1 \\ \displaystyle \frac{\tan{\theta}}{\sqrt{3}} \end{array} \right) \quad (7)
\end{align}
$$

$(6),(7)$式に基づいて接線と法線は下記のように図示できる。

x = np.zeros([2,len(theta)])
U = np.array([[1,-1],[1,1]])/np.sqrt(2)
ratio = 1/np.sqrt(3)

plt.figure(figsize=(5,5))
for i in range(len(c)):
    x[0,:] = c[i] * ratio * np.cos(theta)
    x[1,:] = c[i] * np.sin(theta)
    y = np.dot(U,x)
    plt.plot(y[0,:], y[1,:], "g--")
    
init_pos = np.dot(U,x[:,5])
alpha = np.arange(-5.0,5.01,0.01)

tangent_l = np.dot(U,np.array([x[0,5]+alpha, x[1,5]-alpha*np.sqrt(3)/np.tan(theta[5])]))
grad = np.dot(U,np.array([x[0,5]+alpha, x[1,5]+alpha*np.tan(theta[5])/np.sqrt(3)]))

plt.plot(tangent_l[0,:], tangent_l[1,:], "b-")
plt.plot(grad[0,:], grad[1,:], "r--")
plt.plot([init_pos[0]],[init_pos[1]],"ro")

plt.xlim([-5, 5])
plt.ylim([-5, 5])

plt.show()

・実行結果

共役勾配法

$\mathbf{d}_{1}=\mathbf{v}_{1}$とおくとき、対称行列$A$について$\mathbf{d}_{1}$と共役であるベクトル$\mathbf{d}_{2}$は、下記の式に基づいて得られる。
$$
\large
\begin{align}
\mathbf{d}_{2} = \mathbf{v}_{2} – \frac{\mathbf{v}_{1}^{\mathrm{T}} A \mathbf{v}_{2}}{\mathbf{v}_{1}^{\mathrm{T}} A \mathbf{v}_{1}} \mathbf{v}_{1} \quad (8)
\end{align}
$$

共役勾配法の式に関しては下記で詳しく取り扱った。

$(8)$式に基づいて下記のような計算を行うことで、$\mathbf{d}_{2}$の描画を行うことができる。

A = np.array([[2,1],[1,2]])
U = np.array([[1,-1],[1,1]])/np.sqrt(2)

x = np.zeros([2,len(theta)])
ratio = 1/np.sqrt(3)

plt.figure(figsize=(5,5))
for i in range(len(c)):
    x[0,:] = c[i] * ratio * np.cos(theta)
    x[1,:] = c[i] * np.sin(theta)
    y = np.dot(U,x)
    plt.plot(y[0,:], y[1,:], "g--")
    
init_pos = np.dot(U,x[:,5])
alpha = np.arange(-5.0,5.01,0.01)

tangent_l = np.dot(U,np.array([x[0,5]+alpha, x[1,5]-alpha*np.sqrt(3)/np.tan(theta[5])]))
tan, grad = np.dot(U,np.array([1, -np.sqrt(3)/np.tan(theta[5])])), np.dot(U,np.array([1, np.tan(theta[5])/np.sqrt(3)]))
conjugate_grad = grad - np.dot(tan,np.dot(A,grad))/np.dot(tan,np.dot(A,tan))*tan
conjugate_grad_x = init_pos[0] + alpha*conjugate_grad[0]
conjugate_grad_y = init_pos[1] + alpha*conjugate_grad[1]

plt.plot(tangent_l[0,:], tangent_l[1,:], "b-")
plt.plot(conjugate_grad_x, conjugate_grad_y, "r--")
plt.plot([init_pos[0]],[init_pos[1]],"ro")

plt.xlim([-5, 5])
plt.ylim([-5, 5])

plt.show()

・実行結果

プログラムの理解にあたっては、19行目で用いたtanが$\mathbf{v}_{1}$、gradが$\mathbf{v}_{2}$、conjugate_gradが$\mathbf{d}_{2}$に対応することに着目すると良い。また、最小点は$\mathbf{d}_{2}$を元にline searchを行うことで得られる。

共役勾配法(Conjugate Gradient Method)の原理と具体的な計算例

勾配に基づく最適化はよく行われる一方で、楕円に対して勾配法を適用する際に収束がなかなか進まない場合があります。このような場合に役立つ手法が共役勾配法(Conjugate Gradient Method)です。当記事では共役勾配法の原理や具体例について取りまとめました。

https://ja.wikipedia.org/wiki/共役勾配法

前提の確認

問題設定

オーソドックスな勾配降下(Gradient Descent)法は最適化にあたって有効である場合が多いが、上記のように同心楕円の等高線を持つ関数を取り扱う場合、初期値によっては楕円の等高線の法線ベクトルや勾配が楕円の中心を向かないことで収束がなかなか進まない。

一方で上記のような同心円の等高線を持つ関数を取り扱う場合は同心円の等高線の法線ベクトルや勾配が円の中心を向くので、最適解をスムーズに得ることができる。

次節では同心楕円における最適化を行うにあたって、共役勾配法(Conjugate Gradient Method)のアルゴリズムの確認を行う。

二次形式を用いた立式

$$
\large
\begin{align}
f(\mathbf{x}) &= \frac{1}{2} \mathbf{x}^{\mathrm{T}} A \mathbf{x} – \mathbf{b}^{\mathrm{T}} \mathbf{x} + c \quad (1.1) \\
\mathbf{x} & \in \mathbb{R}^{2}, \, A \in \mathbb{R}^{2 \times 2}, \, A^{\mathrm{T}}=A, \, \mathbf{b} \in \mathbb{R}^{2}, \, c \in \mathbb{R}
\end{align}
$$

上記のような二次形式$f(\mathbf{x})$の最小化問題を定義する。当項では簡易化にあたって$\mathbf{x}$は$2$次元のベクトルを定義したが、$3$次元以上も同様に取り扱うことができる。

$(1.1)$式は下記のように平方完成を行うことができる。
$$
\large
\begin{align}
f(\mathbf{x}) &= \frac{1}{2} \mathbf{x}^{\mathrm{T}} A \mathbf{x} – \mathbf{b}^{\mathrm{T}} \mathbf{x} + c \quad (1.1) \\
&= \frac{1}{2} \left( \mathbf{x}^{\mathrm{T}} A \mathbf{x} – 2 \mathbf{b}^{\mathrm{T}} \right) + c \\
&= \frac{1}{2} \left( \mathbf{x} – A^{-1} \mathbf{b} \right)^{\mathrm{T}} A \left( \mathbf{x} – A^{-1} \mathbf{b} \right) + c’ \quad (1.2)
\end{align}
$$

上記より、$f(\mathbf{x})$は$A^{-1} b$を中心とする同心楕円を描くことが確認できる。よって基本的には$\mathbf{x}^{\mathrm{T}} A \mathbf{x}$の最適化を取り扱えれば$(1.1)$式も単に平行移動で取り扱うことができる。

二次形式の平方完成については下記で詳しく取り扱った。

二次形式の式と楕円

前項の「二次形式を用いた立式」では$f(\mathbf{x})$と定義したが、関数の等高線の形状を取り扱うにあたっては関数の平行移動である$g(\mathbf{x}) = \mathbf{x}^{\mathrm{T}} A \mathbf{x}$のみを確認すれば十分である。よって以下では$g(\mathbf{x})$の等高線を確認する。

$2$次対称行列$A$の固有値を$\lambda_1, \lambda_2(\lambda_1 \neq \lambda_2)$、長さ$1$に正規化された固有ベクトルを$\mathbf{e}_{1}, \mathbf{e}_{2}$と定義するとき、任意のベクトル$\mathbf{x}$は係数$a_1, a_2$を用いて下記のように表せる。
$$
\large
\begin{align}
\mathbf{x} = a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2}
\end{align}
$$

このとき、$g(\mathbf{x})$は下記のように表せる。
$$
\large
\begin{align}
g(\mathbf{x}) &= (a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2})^{\mathrm{T}} A (a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2}) \\
&= (a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2})^{\mathrm{T}} (\lambda_1 a_1 \mathbf{e}_{1} + \lambda_2 a_2 \mathbf{e}_{2}) \\
&= \lambda_{1} a_1^{2} + \lambda_{2} a_2^{2} \quad (1.3)
\end{align}
$$

上記の計算にあたっては対称行列$A$の固有ベクトルが直交することを用いた。詳しい導出は下記で取り扱った。

ここで等高線$g(\mathbf{x})=C$の方程式は下記のように変形できる。
$$
\large
\begin{align}
g(\mathbf{x}) &= \lambda_{1} a_1^{2} + \lambda_{2} a_2^{2} = C \quad (1.4) \\
\frac{a_1^2}{\sqrt{C/\lambda_{1}}^{2}} + \frac{a_2^2}{\sqrt{C/\lambda_{2}}^{2}} &= 1 \quad (1.5)
\end{align}
$$

上記は$a_1, a_2$に関する楕円の方程式を表すが、ここで$\mathbf{x} = a_1 \mathbf{e}_{1} + a_2 \mathbf{e}_{2}$であることより、$\displaystyle \frac{\sqrt{\mathbf{e}_{1}}}{\lambda_{1}}$と$\displaystyle \frac{\sqrt{\mathbf{e}_{2}}}{\lambda_{2}}$を基底に持つ楕円が式に対応することが確認できる。

対称行列Aに関する直交・共役

$$
\large
\begin{align}
\mathbf{u}^{\mathrm{T}} A \mathbf{v} &= 0 \quad (1.6) \\
\mathbf{u} \in \mathbb{R}^{2}, \, & \mathbf{v} \in \mathbb{R}^{2}, \, A \in \mathbb{R}^{2 \times 2}
\end{align}
$$

上記が成立するとき、ベクトル$\mathbf{u}, \, \mathbf{v}$は対称行列$A$に対して「共役である」または「直交する」という。ここで対称行列$A$が$A = P^{\mathrm{T}} P$のように分解できるとき、$(1.6)$式は下記のように変形できる。
$$
\large
\begin{align}
\mathbf{u}^{\mathrm{T}} A \mathbf{v} &= 0 \quad (1.6) \\
\mathbf{u}^{\mathrm{T}} P^{\mathrm{T}} P \mathbf{v} &= 0 \\
(P \mathbf{u})^{\mathrm{T}} (P \mathbf{v}) &= 0
\end{align}
$$

上記の$P$を用いて$\hat{\mathbf{x}} = P \mathbf{x}$のような変換を行うと、$\hat{\mathbf{x}}$の空間では同心楕円の等高線が同心円の等高線に変換される。ここで$(P \mathbf{u})^{\mathrm{T}} (P \mathbf{v}) = 0$は、二次形式の等高線が同心円で表される空間で$(P \mathbf{u})^{\mathrm{T}} \perp (P \mathbf{v})$が成立することを表す。よって、$\mathbf{x}$の空間では等高線の接点からのベクトルが同心楕円の中心の方向に一致する。

共役勾配法は上記に基づいて、「等高線が同心楕円で表される関数の最適化における収束の効率化」を行う手法である。

グラム・シュミットの正規直交化

グラム・シュミットの正規直交化法は、上図のように線型独立である$\vec{a}_{1}, \vec{a}_{2}$が与えられた際に正規直交基底の$\vec{e}_{1}, \vec{e}_{2}$を得る手法である。以下、簡単に手順の確認を行う。

1) $\vec{a}_{1}$に平行な単位ベクトルを$\vec{e}_{1}$とおく。
2) $\vec{e}_{1}$に垂直なベクトル$\vec{u}_{2}$を$\vec{u}_{2}=\vec{a}_{2}-\vec{x}$を計算することで得る。ここで$\vec{x}$は$\vec{a}_{2}$から$\vec{a}_{1}$への正射影を計算することで得られる。
3) $\vec{u}_{2}$と平行な単位ベクトルを$\vec{e}_{2}$とおく。

上記に基づいて、$\vec{e}_{1}, \vec{e}_{2}$はそれぞれ下記のように計算できる。
$$
\large
\begin{align}
\vec{e}_{1} &= \frac{\vec{a}_{1}}{|\vec{a}_{1}|} \\
\vec{u}_{2} &= \vec{a}_{2} – \vec{x} = \vec{a}_{2} – \frac{\vec{a}_{1} \cdot \vec{a}_{2}}{|\vec{a}_{1}|^2} \vec{a}_{1} \\
\vec{e}_{2} &= \frac{\vec{u}_{2}}{|\vec{u}_{2}|}
\end{align}
$$

グラム・シュミットの正規直交化について詳しくは下記で取り扱った。

共役勾配法

共役勾配法の原理

「共役勾配法」は「グラム・シュミットの正規直交化」における「直交」を「対称行列$A$に関する共役・直交」に基づいて取り扱った手法である。以下、線型独立である$\mathbf{v}_{1}, \mathbf{v}_{2}$を元に対称行列$A$に互いに共役であるベクトル$\mathbf{d}_1, \mathbf{d}_2$をグラム・シュミットの直交化に基づいて生成を行う。$\displaystyle \mathbf{d}_1 = \mathbf{v}_{1}$とおくと、$\mathbf{d}_{2}$は下記のように得られる。
$$
\large
\begin{align}
\mathbf{d}_{2} = \mathbf{v}_{2} – \frac{\mathbf{v}_{1}^{\mathrm{T}} A \mathbf{v}_{2}}{\mathbf{v}_{1}^{\mathrm{T}} A \mathbf{v}_{1}} \mathbf{v}_{1} \quad (2.1)
\end{align}
$$

上記の式は正規化を行わなかったので、前節の式とやや異なることに注意が必要である。$\mathbf{d}_1$の楕円の等高線との接点から$\mathbf{d}_{2}$の方向にline searchなどを用いることで関数の最小点である楕円の中心を得ることができる。また、最小点$\mathbf{x}^{*}$について下記が成立する。
$$
\large
\begin{align}
\nabla f(\mathbf{x}) &= A \mathbf{x} – \mathbf{b} \quad (1.1)’ \\
A \mathbf{x}^{*} – \mathbf{b} &= \mathbf{0} \\
A \mathbf{x}^{*} &= \mathbf{b} \quad (2.2)
\end{align}
$$

$(2.2)$式より、ここで得られた最小点$\mathbf{x}^{*}$は$A \mathbf{x} = \mathbf{b}$の解に一致することも確認できる。

また、$(2.1)$式を一般化することで、線型独立である$\mathbf{v}_{1}, \cdots , \mathbf{v}_{p}$を元に対称行列$A$に互いに共役であるベクトル$\mathbf{d}_{j}, \, (i=2, \cdots, p)$を下記のように得ることができる。
$$
\large
\begin{align}
\mathbf{d}_{j} = \mathbf{v}_{j} – \sum_{i=1}^{j-1} \frac{\mathbf{v}_{i}^{\mathrm{T}} A \mathbf{v}_{j}}{\mathbf{v}_{i}^{\mathrm{T}} A \mathbf{v}_{i}} \mathbf{v}_{i}
\end{align}
$$

共役勾配法の計算例

$$
\large
\begin{align}
f(\mathbf{x}) &= \frac{1}{2} \mathbf{x}^{\mathrm{T}} A \mathbf{x} – \mathbf{b}^{\mathrm{T}} \mathbf{x} + c \quad (2.3) \\
A &= 2 \left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \\
\mathbf{b} &= \left(\begin{array}{c} 0 \\ 0 \end{array} \right), \, c=0
\end{align}
$$

以下、上記のように定義した二次形式$f(\mathbf{x})$の最適化を共役勾配法を用いて行う。$f(\mathbf{x})$は下記のように表すことができる。
$$
\large
\begin{align}
f(\mathbf{x}) &= \left(\begin{array}{cc} x_1 & x_2 \end{array} \right) \left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \left(\begin{array}{c} x_1 \\ x_2 \end{array} \right) \quad (2.3)’ \\
A’ &= \left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right)
\end{align}
$$

ここで行列$A’$の固有値は$\lambda_1=3, \lambda_2=1$、$\lambda_1=3, \lambda_2=1$、固有ベクトルは$\displaystyle \mathbf{u}_{1} = \left(\begin{array}{c} 1 \\ 1 \end{array} \right), \, \mathbf{u}_{2} = \left(\begin{array}{c} -1 \\ 1 \end{array} \right)$が得られる。導出過程はオーソドックスな固有値・固有ベクトルの計算なので省略するが、この固有値と固有ベクトルがそれぞれ適切であることは下記の等式が成立することからそれぞれ確認できる。
$$
\large
\begin{align}
\left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \left(\begin{array}{c} 1 \\ 1 \end{array} \right) &= 3 \left(\begin{array}{c} 1 \\ 1 \end{array} \right) \\
\left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right) \left(\begin{array}{c} -1 \\ 1 \end{array} \right) &= \left(\begin{array}{c} -1 \\ 1 \end{array} \right)
\end{align}
$$

上記に基づいて、二次形式$f(\mathbf{x})$の等高線は下図のように表される。

ここで楕円上の点における勾配は楕円の法線ベクトルに一致するが、勾配ベクトルは下図のように楕円の中心を向かない。

上図では勾配ベクトルが楕円の中心を向かないので、$A’$について共役な方向を計算する必要がある。前項の共役勾配法を用いることで下記のような結果を得ることができる。

共役勾配法のプログラムについては下記で詳しく取り扱った。

対称行列の対角化とスペクトル分解(spectral decomposition)

統計学や機械学習で出てくる行列は対称行列(symmetric matrix)が多いですが、対称行列の取り扱いはやや特殊なので抑えておくと良いです。当記事では対称行列の対角化とスペクトル分解(spectral decomposition)について取りまとめました。
「統計学のための数学入門$30$講」の$23$章の「対称行列の固有値と固有ベクトル」を参考に作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

対称行列の対角化とスペクトル分解

$p$次の対称行列$A$に関して下記がそれぞれ成立する。

$1)$
$p$次の対称行列$A$は直交行列$U$を用いて下記のように対角化が可能である。
$$
\large
\begin{align}
U^{\mathrm{T}} A U &= \Lambda = \left( \begin{array}{ccccc} \lambda_1 & 0 & 0 & \cdots & 0 \\ 0 & \lambda_2 & 0 & \cdots & 0 \\ 0 & 0 & \lambda_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \lambda_{p} \end{array} \right) \quad (1) \\
U &= \left( \begin{array}{ccccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} & \cdots & \mathbf{u}_{p} \end{array} \right) \\
U^{\mathrm{T}} &= \left( \begin{array}{c} \mathbf{u}_{1}^{\mathrm{T}} \\ \mathbf{u}_{2}^{\mathrm{T}} \\ \mathbf{u}_{3}^{\mathrm{T}} \\ \vdots \\ \mathbf{u}_{p}^{\mathrm{T}} \end{array} \right) \\
\mathbf{u}_{i}^{\mathrm{T}} \mathbf{u}_{i} &= 1 \\
\mathbf{u}_{i}^{\mathrm{T}} \mathbf{u}_{j} &= 0, \quad (i \neq j) \\
\mathbf{u}_{i} & \in \mathbb{R}^{p}
\end{align}
$$

$2)$
$p$次の対称行列$A$は固有値$\lambda_1, \cdots , \lambda_{p}$とそれぞれの固有値に対応する長さ$1$の固有ベクトル$\mathbf{u}_{1}, \cdots , \mathbf{u}_{p}$を用いて次のように表せる。
$$
\large
\begin{align}
A &= U \Lambda U^{\mathrm{T}} \quad (2) \\
&= \left(\begin{array}{ccccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} & \cdots & \mathbf{u}_{p} \end{array} \right) \left(\begin{array}{ccccc} \lambda_1 & 0 & 0 & \cdots & 0 \\ 0 & \lambda_2 & 0 & \cdots & 0 \\ 0 & 0 & \lambda_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \lambda_{p} \end{array} \right) \left(\begin{array}{c} \mathbf{u}_{1}^{\mathrm{T}} \\ \mathbf{u}_{2}^{\mathrm{T}} \\ \mathbf{u}_{3}^{\mathrm{T}} \\ \vdots \\ \mathbf{u}_{p}^{\mathrm{T}} \end{array} \right) \\
&= \left(\begin{array}{ccccc} \lambda_{1} \mathbf{u}_{1} & \lambda_{2} \mathbf{u}_{2} & \lambda_{3} \mathbf{u}_{3} & \cdots & \lambda_{p} \mathbf{u}_{p} \end{array} \right) \left(\begin{array}{c} \mathbf{u}_{1}^{\mathrm{T}} \\ \mathbf{u}_{2}^{\mathrm{T}} \\ \mathbf{u}_{3}^{\mathrm{T}} \\ \vdots \\ \mathbf{u}_{p}^{\mathrm{T}} \end{array} \right) \\
&= \lambda_{1} \mathbf{u}_{1} \mathbf{u}_{1}^{\mathrm{T}} + \lambda_{2} \mathbf{u}_{2} \mathbf{u}_{2}^{\mathrm{T}} + \lambda_{3} \mathbf{u}_{3} \mathbf{u}_{3}^{\mathrm{T}} + \cdots + \lambda_{p} \mathbf{u}_{p} \mathbf{u}_{p}^{\mathrm{T}} \\
&= \sum_{i=1}^{p} \lambda_{i} \mathbf{u}_{i} \mathbf{u}_{i}^{\mathrm{T}}
\end{align}
$$

上記の変形をスペクトル分解という。

導出

$p$個の固有値が全て異なる$p$次対称行列$A$の固有値を$\lambda_1, \cdots , \lambda_{p}$、それぞれの固有値に対応する長さ$1$の固有ベクトルを$\mathbf{u}_{1}, \cdots , \mathbf{u}_{p}$とおくと、下記のような式が成立する。
$$
\large
\begin{align}
A \mathbf{u}_{i} = \lambda_{i} \mathbf{u}_{i}, \quad i=1, 2, 3, \cdots , p
\end{align}
$$

上記より、下記が成立する。
$$
\large
\begin{align}
A U &= U \Lambda \quad (3) \\
\Lambda &= \left(\begin{array}{ccccc} \lambda_1 & 0 & 0 & \cdots & 0 \\ 0 & \lambda_2 & 0 & \cdots & 0 \\ 0 & 0 & \lambda_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \lambda_{p} \end{array} \right) \\
U &= \left(\begin{array}{ccccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \mathbf{u}_{3} & \cdots & \mathbf{u}_{p} \end{array} \right)
\end{align}
$$

対称行列の相異なる固有値に対応する固有ベクトルが直交するので、$U$は直交行列である。よって$U^{-1}=U^{\mathrm{T}}$が成立する。したがって$(3)$式の変形に基づいて下記のように$(1), (2)$式の導出を行える。

・$(1)$式の導出
$$
\large
\begin{align}
A U &= U \Lambda \quad (3) \\
U^{-1} A U &= U^{-1} U \Lambda \\
U^{\mathrm{T}} A U &= \Lambda \quad (1)
\end{align}
$$

・$(2)$式の導出
$$
\large
\begin{align}
A U &= U \Lambda \quad (3) \\
A U U^{-1} &= U \Lambda U^{-1} \\
A &= U \Lambda U^{\mathrm{T}} \quad (2)
\end{align}
$$

対称行列(symmetric matrix)の固有値と固有ベクトルの性質とその導出

統計学や機械学習で出てくる行列は対称行列(symmetric matrix)である場合が多いですが、対称行列の取り扱いはやや特殊なので抑えておくと良いです。当記事では対称行列の固有値と固有ベクトルの性質とその導出について取りまとめました。
「統計学のための数学入門$30$講」の$23$章の「対称行列の固有値と固有ベクトル」を参考に作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

対称行列の固有値・固有ベクトルの性質

対称行列の固有値・固有ベクトルでは下記の性質がそれぞれ成立する。

$1. \,$ 対称行列の固有値は全て実数である。
$2. \,$ 対称行列の相異なる固有値に対応する固有ベクトルは直交する

導出

対称行列の固有値は全て実数

対称行列$A$の固有値を$\lambda$、固有値に対応する固有ベクトルを$\mathbf{x}$とおく。このとき、$\lambda$の共役複素数を$\bar{\lambda}$、$\mathbf{x}$の成分の共役複素数を成分に持つベクトルを$\bar{\mathbf{x}}$と定義する。

上記の定義は下記のような数式で表すことができる。
$$
\large
\begin{align}
A \mathbf{x} &= \lambda \mathbf{x} \quad (1) \\
A \bar{\mathbf{x}} &= \bar{\lambda} \bar{\mathbf{x}} \quad (2)
\end{align}
$$

$(1)$式に関して以下の変形が成立する。
$$
\large
\begin{align}
A \mathbf{x} &= \lambda \mathbf{x} \quad (1) \\
\bar{\mathbf{x}}^{\mathrm{T}} A \mathbf{x} &= \bar{\mathbf{x}}^{\mathrm{T}} (\lambda \mathbf{x}) \\
\bar{\mathbf{x}}^{\mathrm{T}} A \mathbf{x} &= \lambda \bar{\mathbf{x}}^{\mathrm{T}} \mathbf{x} \\
(\bar{\mathbf{x}}^{\mathrm{T}} A \mathbf{x})^{\mathrm{T}} &= (\lambda \bar{\mathbf{x}}^{\mathrm{T}} \mathbf{x})^{\mathrm{T}} \\
\mathbf{x}^{\mathrm{T}} A \bar{\mathbf{x}} &= \lambda \mathbf{x}^{\mathrm{T}} \bar{\mathbf{x}} \quad (3)
\end{align}
$$

途中の式変形にあたっては$A^{\mathrm{T}}=A$であることを用いた。

同様に$(2)$式に関して以下の変形が成立する。
$$
\large
\begin{align}
A \bar{\mathbf{x}} &= \bar{\lambda} \bar{\mathbf{x}} \quad (2) \\
\mathbf{x}^{\mathrm{T}} A \bar{\mathbf{x}} &= \mathbf{x}^{\mathrm{T}} (\bar{\lambda} \bar{\mathbf{x}}) \\
\mathbf{x}^{\mathrm{T}} A \bar{\mathbf{x}} &= \bar{\lambda} \mathbf{x}^{\mathrm{T}} \bar{\mathbf{x}} \quad (4)
\end{align}
$$

$(3), (4)$式より下記が成立する。
$$
\large
\begin{align}
\lambda \mathbf{x}^{\mathrm{T}} \bar{\mathbf{x}} &= \bar{\lambda} \mathbf{x}^{\mathrm{T}} \bar{\mathbf{x}} \\
(\lambda – \bar{\lambda}) \mathbf{x}^{\mathrm{T}} \bar{\mathbf{x}} &= 0
\end{align}
$$

ここで$\mathbf{x}^{\mathrm{T}} \bar{\mathbf{x}} > 0$より、$\lambda=\bar{\lambda}$が成立するので$\lambda$は実数である。

対称行列の固有ベクトルが直交

対称行列$A$の相異なる$2$つの固有値を$\lambda_1, \lambda_2$、それぞれの固有値に対応する固有ベクトルを$\mathbf{x}_{1}, \mathbf{x}_{2}$とおく。

上記は下記のような式で表せる。
$$
\large
\begin{align}
A \mathbf{x}_{1} &= \lambda_1 \mathbf{x}_{1} \quad (5) \\
A \mathbf{x}_{2} &= \lambda_2 \mathbf{x}_{2} \quad (6)
\end{align}
$$

このとき$A^{\mathrm{T}}=A$であることに基づいて$(5)$式について下記の変形が成立する。
$$
\large
\begin{align}
A \mathbf{x}_{1} &= \lambda_1 \mathbf{x}_{1} \quad (5) \\
\mathbf{x}_{2}^{\mathrm{T}} A \mathbf{x}_{1} &= \mathbf{x}_{2}^{\mathrm{T}} (\lambda_1 \mathbf{x}_{1}) \\
\mathbf{x}_{2}^{\mathrm{T}} A \mathbf{x}_{1} &= \lambda_1 \mathbf{x}_{2}^{\mathrm{T}} \mathbf{x}_{1} \\
(\mathbf{x}_{2}^{\mathrm{T}} A \mathbf{x}_{1})^{\mathrm{T}} &= (\lambda_1 \mathbf{x}_{2}^{\mathrm{T}} \mathbf{x}_{1})^{\mathrm{T}} \\
\mathbf{x}_{1}^{\mathrm{T}} A \mathbf{x}_{2} &= \lambda_1 \mathbf{x}_{1}^{\mathrm{T}} \mathbf{x}_{2} \quad (7)
\end{align}
$$

また、$(6)$式について下記が成立する。
$$
\large
\begin{align}
A \mathbf{x}_{2} &= \lambda_2 \mathbf{x}_{2} \quad (6) \\
\mathbf{x}_{1}^{\mathrm{T}} A \mathbf{x}_{2} &= \mathbf{x}_{1}^{\mathrm{T}} (\lambda_2 \mathbf{x}_{2}) \\
\mathbf{x}_{1}^{\mathrm{T}} A \mathbf{x}_{2} &= \lambda_2 \mathbf{x}_{1}^{\mathrm{T}} \mathbf{x}_{2} \quad (8)
\end{align}
$$

$(7), (8)$式より下記が成立する。
$$
\large
\begin{align}
\lambda_1 \mathbf{x}_{1}^{\mathrm{T}} \mathbf{x}_{2} &= \lambda_2 \mathbf{x}_{1}^{\mathrm{T}} \mathbf{x}_{2} \\
(\lambda_1 \, – \, \lambda_2) \mathbf{x}_{1}^{\mathrm{T}} \mathbf{x}_{2} &= 0 \quad (9)
\end{align}
$$

ここで$\lambda_1 \neq \lambda_2$を仮定したので、$(9)$式より$\mathbf{x}_{1}^{\mathrm{T}} \mathbf{x}_{2}=0$が成立するが、$\mathbf{x}_{1}, \mathbf{x}_{2}$はどちらも零ベクトルではないので、$\mathbf{x}_{1} \perp \mathbf{x}_{2}$が成立する。

上記を用いることで、「$n \times n$対称行列が$n$個の相異なる固有値を持つとき、大きさ$1$の固有ベクトルを並べた行列は直交行列になる」ということも示すことができる。

統計検定準1級 問題解説 ~2016年6月実施 選択問題及び部分記述問題 問8~

過去問

過去問題は統計検定公式が問題と解答例を公開しています。こちらを参照してください。

[1] 解答

$\boxed{\mathsf{9}}$ : $④$

$S(T) = P(T>t)$ が生存関数なので,累積分布関数 $F(t) = P(T \leq t)$ は,
$$
\begin{equation}
F(t) = P(T\leq t)= 1-S(t)=\begin{cases}
1-\exp(-\lambda t) & \text{($t>0$)} \\
0 & \text{($t<0$)} \ \end{cases} \end{equation} $$ 従って,確率密度関数 $f(t)$ は, $$ \begin{equation} f(t) = F'(T)= \begin{cases} \lambda\exp(-\lambda t) & \text{($t>0$)} \\
0 & \text{($t<0$)}
\end{cases}
\end{equation}
$$
となる.

[2] 解答

$\boxed{\mathsf{10}}$ : $⑤$

平均は,

$$
\begin{align} E[X] &= \int_0^\infty tf(t) dt = \int _0 ^ \infty t\lambda \exp(-\lambda t)dt = \lambda \int_0 ^ \infty t\exp(-\lambda t) dt \\ &= \lambda \left( -\dfrac{1}{\lambda}[ t\exp(-\lambda t)]_0^\infty + \dfrac{1}{\lambda}\int_0 ^\infty \exp(-\lambda t)dt \right) \\ &=\lambda \left( 0 + \dfrac{1}{\lambda}\left[ \dfrac{1}{\lambda}\exp(-\lambda t)\right]_0^\infty \right) \\ &= \dfrac{1}{\lambda} \end{align}
$$


である.


中央値は $1-S(t) = 1/2$ を満たす $t$ である.つまり,$S(t)= 1/2$ を解けばよい.
$S(t)=\exp(-\lambda t) = 1/2$ より,$t = \log 2 /\lambda$ である.

(補足)[1]から確率変数 $T$ の従う分布がパラメータ $\lambda$ の指数分布であることがわかるので,そのことから平均が $1/\lambda$ としてもよい.

[3] 解答

$\boxed{\mathsf{11}}$ : $②$
[2]から平均が $1/\lambda$ がわかっているので,$1/\lambda = 1/2$ より,中央値は $\log2/\lambda = 2\log2 \approx 1.4$ である.

重点サンプリング(Importance Sampling)の数式表記とPythonを用いた計算例

特定の確率分布の期待値を別の確率分布からサンプリングした値に基づいて計算する手法を重点サンプリング(Importance Sampling)といいます。当記事では重点サンプリングの数式表記とPythonを用いた計算例の確認をそれぞれ行いました。
「ゼロから作るDeep Learning④ー強化学習編」の$5.5.2$「重点サンプリング」〜$5.5.3$節「分散を小さくするには」の内容などを参考に当記事の作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

重点サンプリングの基本事項

サンプリングによる期待値の近似

離散確率分布$\pi(x)$に基づいて得られる確率変数$x$の期待値を$E_{\pi}[x]$とおくと、期待値は下記のように定義される。
$$
\large
\begin{align}
\mathbb{E}_{\pi}[x] = \sum x \pi(x)
\end{align}
$$

このとき、確率分布が具体的に判明しなくても確率分布からのサンプリングが行うことができる場合は、確率分布に基づいて得られたサンプルから下記のように期待値$E_{\pi}[x]$を近似することができる。
$$
\large
\begin{align}
& \mathrm{sampling} : \, x^{(i)} \sim \pi \quad i=1,2,\cdots,n \\
& E_{\pi}[x] \simeq \frac{1}{n} \sum_{i=1}^{n} x^{(i)}
\end{align}
$$

ここで上記のインデックス$i$は確率分布$\pi(x)$から得られた$i$番目のサンプルを表すことを抑えておくと良い。

重点サンプリングの仕組み

通常のサンプリングによる近似では前項の『サンプリングによる期待値の近似』のように近似を行う。一方で重点サンプリングは確率分布$\pi(x)$の期待値を別の確率分布$b(x)$から得られたサンプルを元に計算する手法であり、下記のような変形に基づく。
$$
\large
\begin{align}
\mathbb{E}_{\pi}[x] &= \sum x \pi(x) \\
&= \sum x \frac{b(x)}{b(x)} \pi(x) \\
&= \sum x \frac{\pi(x)}{b(x)} b(x) = \mathbb{E}_{b} \left[ \frac{\pi(x)}{b(x)} x \right]
\end{align}
$$

上記より、$\displaystyle \mathbb{E}_{\pi}[x] = \mathbb{E}_{b} \left[ \frac{\pi(x)}{b(x)} x \right]$が成立するので、確率分布$b(x)$に基づいて得られたサンプルを元に確率分布$\pi(x)$の期待値を計算することができる。

Pythonを用いた計算例

確率分布$\pi$に基づくサンプリング

確率分布$\pi(x)$からのサンプリングに基づく期待値$E_{\pi}[x]$の近似値は下記のように計算することができる。

import numpy as np

x = np.array([1, 2, 3])
pi = np.array([0.1, 0.1, 0.8])

# 期待値
e = np.sum(x * pi)
print("E_pi[x]: {}".format(e))

np.random.seed(100)

# モンテカルロ法
n = 100
samples = []
for i in range(n):
    s = np.random.choice(x, p=pi)
    samples.append(s)
    
print("mean: {:.2f}, var: {:.2f}".format(np.mean(samples), np.var(samples)))

・実行結果

E_pi[x]: 2.7
mean: 2.71, var: 0.39

確率分布$b$に基づく重点サンプリング

確率分布$b(x)$からのサンプリングに基づく$\displaystyle \mathbb{E}_{b} \left[ \frac{\pi(x)}{b(x)} x \right]$の近似値は下記のように計算することができる。

np.random.seed(25)

b = np.array([1/3, 1/3, 1/3])
n = 100
samples = []

for i in range(n):
    idx = np.arange(len(b))
    i = np.random.choice(idx, p=b)
    s = x[i]
    rho = pi[i] / b[i]
    samples.append(rho * s)
    
print("mean: {:.2f}, var: {:.2f}".format(np.mean(samples), np.var(samples)))

・実行結果

mean: 2.56, var: 9.70

上記より、重点サンプリングによる近似値はある程度妥当である一方で分散が大きいことが確認できる。次節ではこの分散が大きくなる原因とその対応について確認を行う。

重点サンプリングと分散

分散が大きくなる原因

重点サンプリングを用いた際に分散が大きくなるのは確率分布$\pi(x)$と$b(x)$の差が大きい場合に$\displaystyle \mathbb{E}_{b} \left[ \frac{\pi(x)}{b(x)} x \right]$の$\displaystyle \frac{\pi(x)}{b(x)}$が$1$から大きく離れた値になり、$\displaystyle \frac{\pi(x)}{b(x)} x$の値が安定しないことに起因する。

前節の例では$\pi: (0.1,0.1,0.8)$に対し、$\displaystyle b: \left( \frac{1}{3},\frac{1}{3},\frac{1}{3} \right)$を用いたが、$\displaystyle \frac{\pi(x)}{b(x)}$はそれぞれ下記のように計算される。

・$x=1,2$の場合
$$
\large
\begin{align}
\frac{\pi(x)}{b(x)} &= \frac{0.1}{1/3} \\
&= 0.3
\end{align}
$$

・$x=3$の場合
$$
\large
\begin{align}
\frac{\pi(x)}{b(x)} &= \frac{0.8}{1/3} \\
&= 2.4
\end{align}
$$

上記より、$x=1,2$が得られた場合$\displaystyle \frac{\pi(x)}{b(x)} x$がそれぞれ$0.3, 0.6$であるのに対し、$x=3$が得られた場合$\displaystyle \frac{\pi(x)}{b(x)} x$が$7.2$の値を取る。このように$\displaystyle \frac{\pi(x)}{b(x)}$が$1$から大きく離れた値になることで、$\displaystyle \frac{\pi(x)}{b(x)} x$の値のばらつきが大きくなる。

この解決にあたっては、$\displaystyle \frac{\pi(x)}{b(x)}$を$1$に近づければ良いので、$\pi(x)$に類似の$b(x)$を用いればよい。

また、$\displaystyle \frac{\pi(x)}{b(x)}$を下記のように何らかの文字で表す場合もあるので、下記のような表記も抑えておくと良い。
$$
\large
\begin{align}
\rho(x) = \frac{\pi(x)}{b(x)}
\end{align}
$$

Pythonを用いた計算例

以下では$b:(0.2,0.2,0.6)$を用いて前節と同様に$\displaystyle \mathbb{E}_{b} \left[ \frac{\pi(x)}{b(x)} x \right]$の近似値と分散の計算を行った。

np.random.seed(25)

b = np.array([0.2, 0.2, 0.6])
n = 100
samples = []

for i in range(n):
    idx = np.arange(len(b))
    i = np.random.choice(idx, p=b)
    s = x[i]
    rho = pi[i] / b[i]
    samples.append(rho * s)
    
print("mean: {:.2f}, var: {:.2f}".format(np.mean(samples), np.var(samples)))

・実行結果

mean: 2.82, var: 2.50

実行結果より分散が前節の重点サンプリングの結果より小さくなったことが確認できる。

方策勾配法のアルゴリズムまとめ 〜REINFORCE・ベースライン・Actor-Critic〜

方策勾配法(Policy Gradient Method)を改善させたアルゴリズムには、REINFORCE・ベースライン・Actor-Criticなどのアルゴリズムがあります。当記事ではこれらの$3$つのアルゴリズムについて取りまとめを行いました。
「ゼロから作るDeep Learning④ー強化学習編」の第$9$章の「方策勾配法」や付録Dの「方策勾配法の証明」の内容を参考に当記事の作成を行いました。

・用語/公式解説
https://www.hello-statisticians.com/explain-terms

前提知識

方策勾配法の目的関数と勾配の式

詳しくは上記で取り扱ったが、方策勾配法の目的関数と目的関数の勾配の式はそれぞれ下記のように表される。

・目的関数
$$
\large
\begin{align}
J(\theta) &= \mathbb{E}_{\tau \sim \pi_{\theta}}[G(\tau)] \quad (1.1) \\
G(\tau) &= R_0 + \gamma R_1 + \gamma^{2} R_2 + \cdots + \gamma^{T} R_{T} \quad (1.2)
\end{align}
$$

・目的関数の勾配
$$
\large
\begin{align}
\nabla_{\theta} J(\theta) = \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T} G(\tau) \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)} \right] \quad (1.3)
\end{align}
$$

方策の学習にあたって

方策の学習にあたっては下記のような勾配上昇法を用いることで行える。
$$
\large
\begin{align}
\theta \leftarrow \theta + \alpha \nabla_{\theta} J(\theta) \quad (1.4)
\end{align}
$$

上記はオーソドックスな勾配法の式であり、$\alpha$は学習率を表す。また、実際に$(1.2)$式を用いて学習を行うにあたっては、下記のように実際に方策$\pi_{\theta}$に基づいて軌道$\tau^{(i)}$をサンプリングし、その平均を取れば勾配の近似値が計算できる。
$$
\large
\begin{align}
\tau^{(i)} & \sim \pi_{\theta}, \qquad i=1,2,\cdots,n, \quad \mathrm{i.i.d.} \quad (1.5) \\
j^{(i)} &= \sum_{t=0}^{T} G \left( \tau^{(i)} \right) \nabla_{\theta} \log{ \pi_{\theta} \left( A_t^{(i)} \middle| S_t^{(i)} \right) } \quad (1.6) \\
\nabla_{\theta} J(\theta) &= \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T} G(\tau) \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)} \right] \simeq \frac{1}{n} \sum_{i=1}^{n} j^{(i)} \quad (1.7)
\end{align}
$$

$(1.7)$式が勾配の近似値を表す。上記の$(1.6)$式の$i$は方策$\pi_{\theta}$に基づいて生成された軌道$\tau^{(i)}$のインデックス、$t$は軌道$\tau^{(i)}$における系列$\displaystyle \left(S_0^{(i)},A_{0}^{(i)}, \cdots , S_t^{(i)},A_{t}^{(i)}, \cdots, S_{T+1}^{(i)} \right)$のインデックスを表すことに注意しておくと良い。

$i$番目の軌道$\tau^{(i)}$における収益は$\displaystyle G \left( \tau^{(i)} \right)$によって表されるが、$(1.6)$式はこの$\displaystyle G \left( \tau^{(i)} \right)$による重み付け和であると解釈することができる。

ここで$\displaystyle \nabla_{\theta} \log{ \pi_{\theta} \left( A_t^{(i)} \middle| S_t^{(i)} \right) }$はインデックス$t$を用いて表される一方で、重みが$\displaystyle G \left( \tau^{(i)} \right)$ではこれまでに得た報酬も取り扱うので方策の評価を行う観点では適切ではない。よって、$(1.5)$式の$\displaystyle G \left( \tau^{(i)} \right)$に改善の余地が生じる。次節ではこの改善の余地の視点から方策勾配法の改善の手法についてそれぞれ取り扱う。

方策勾配法の改善

基本方針

前節で取り扱った$(1.3)$式に基づいて勾配は下記のように表される。
$$
\large
\begin{align}
\nabla_{\theta} J(\theta) = \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T} G(\tau) \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)} \right] \quad (1.3)
\end{align}
$$

上記の$G(\tau)$を$\Phi_{t}$に置き換えると下記が得られる。
$$
\large
\begin{align}
\nabla_{\theta} J(\theta) = \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T} \Phi_{t} \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)} \right] \quad (2.1)
\end{align}
$$

以下では$\Phi_{t}$に様々な指標を用いることで方策勾配法の改善について取り扱う。

REINFORCE

REINFORCE(REward Increment $=$ Nonnegative Factor $\times$ Offset Reinforcement $\times$ Characteristic Eligibility)は下記のように$\Phi_{t}$を定義する。
$$
\large
\begin{align}
\nabla_{\theta} J(\theta) &= \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T} \Phi_{t} \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)} \right] \quad (2.1)’ \\
\Phi_{t} &= G_{t} \quad (2.2) \\
G_{t} &= R_t + \gamma R_{t+1} + \gamma^{2} R_{t+2} + \cdots + \gamma^{T-t} R_{T} \quad (2.3)
\end{align}
$$

このようにREINFORCEでは$\displaystyle \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)}$にかかる重みを$S_t$以後に得られる収益$G_t$を用いて表す。

ベースライン付きREINFORCE

ベースライン付きREINFORCEはベースライン$b(S_t)$を導入することで$\Phi_{t}$の分散を小さくする手法であり、下記のように$\Phi_{t}$を定義する。
$$
\large
\begin{align}
\nabla_{\theta} J(\theta) &= \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T} \Phi_{t} \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)} \right] \quad (2.1)’ \\
\Phi_{t} &= G_{t} – b(S_t) \quad (2.4)
\end{align}
$$

ここで上記の収益$G_t$にQ関数、ベースライン$b(S_t)$に価値関数を用いると$\Phi_{t}$はアドバンテージ関数に一致する。詳しくは『ベースラインとアドバンテージ関数』で取り扱う。

Actor-Critic

Actor-Criticでは下記のように$\Phi_{t}$を定義する。
$$
\large
\begin{align}
\nabla_{\theta} J(\theta) &= \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T} \Phi_{t} \nabla_{\theta} \log{\pi_{\theta}(A_t|S_t)} \right] \quad (2.1)’ \\
\Phi_{t} &= R_{t} + \gamma V(S_t) – V(S_t) \quad (2.5)
\end{align}
$$

Actor-CriticはTD法と対応させて理解すると良い。

ベースラインとアドバンテージ関数

ベースライン付きREINFORCE』では$\Phi_{t} = G_{t} – b(S_t)$のようにベースライン$b(S_t)$を導入したが、収益$G_t$にQ関数、ベースライン$b(S_t)$に価値関数を用いると$\Phi_{t}$は下記のように表せる。
$$
\large
\begin{align}
\Phi_{t} &= Q(S_t,A_t) – V(S_t) \\
&= A(S_t,A_t) \quad (2.6)
\end{align}
$$

上記の$A(S_t,A_t)$はアドバンテージ関数であり、PPOのように論文によっては収益ではなくアドバンテージを元に式が表される場合があるので$(2.6)$式は抑えておくと良い。

・PPO論文