ブログ

【上級】データサイエンス 数学ストラテジスト 公式問題集 解答例まとめ Q.11〜20

「データサイエンス 数学ストラテジスト 上級」はデータサイエンスの基盤である、確率・統計、線形代数、微積分、機械学習、プログラミングなどを取り扱う資格試験です。当記事では「日本数学検定協会」作成の「公式問題集」の演習問題$11$〜$20$の解答例を取り扱いました。

・数学検定まとめ
https://www.hello-statisticians.com/math_certificate

演習問題

Q.11

$\overrightarrow{AB}$と$\overrightarrow{OH}$はそれぞれ下記のように表すことができる。
$$
\large
\begin{align}
\overrightarrow{AB} &= \overrightarrow{OB} – \overrightarrow{OA} \\
\overrightarrow{OH} &= t \overrightarrow{OA} + (1-t) \overrightarrow{OB}
\end{align}
$$

ここで$\overrightarrow{AB} \perp \overrightarrow{OH}$であるので、$\overrightarrow{AB} \cdot \overrightarrow{OH} = 0$が成立するが、下記のように変形できる。
$$
\large
\begin{align}
\overrightarrow{AB} \cdot \overrightarrow{OH} &= 0 \\
(\overrightarrow{OB} – \overrightarrow{OA}) \cdot (t \overrightarrow{OA} + (1-t) \overrightarrow{OB}) &= 0 \\
t &= \frac{70}{127}
\end{align}
$$

よって、$\overrightarrow{OH}$は下記のように表すことができる。
$$
\large
\begin{align}
\overrightarrow{OH} = \frac{70}{127} \vec{a} + \frac{57}{127} \vec{b}
\end{align}
$$

上記より$(4)$が正しい。

Q.12

下記のように式変形を行える。
$$
\large
\begin{align}
\lim_{n \to \infty} \left[ \log_{10}{(n+100)} – \log_{10}{(1000n+1)} \right] &= \lim_{n \to \infty} \log_{10}{\frac{n+100}{1000n+1}} \\
&= \lim_{n \to \infty} \log_{10}{\frac{1+100/n}{1000+1/n}} \\
&= \log_{10}{\frac{1}{1000}} \\
&= \log_{10}{10^{-3}} \\
&= -3
\end{align}
$$

上記より$(5)$が正しい。

Q.13

$$
\large
\begin{align}
y = \frac{6}{(e^{x}+1)(e^{x}+3)}
\end{align}
$$

$e^{x}>0$より、$0 \leq x \leq 1$で$y>0$である。よって、図形の面積を$S$とおくと、$S$は下記のように立式できる。
$$
\large
\begin{align}
S = \int_{0}^{1} \frac{6}{(e^{x}+1)(e^{x}+3)} dx \quad [1]
\end{align}
$$

上記に対し、$t=e^{x}$のような変数変換を行う。このとき$\displaystyle \frac{dt}{dx}=e^{x}$より、$\displaystyle \frac{dt}{t}=dx$が成立する。また、$x$と$t$の区間は下記のように対応する。

$x$ $0 \longrightarrow 1$
$t$ $1 \longrightarrow e$

よって$[1]$式は下記のように値を得ることができる。
$$
\large
\begin{align}
S &= \int_{0}^{1} \frac{6}{(e^{x}+1)(e^{x}+3)} dx \quad [1] \\
&= \int_{1}^{e} \frac{6}{t(t+1)(t+3)} dt \\
&= \int_{1}^{e} \left[ \frac{2}{t} – \frac{3}{t+1} + \frac{1}{t+3} \right] dt \\
&= \left[ 2\log{|t|} – 3\log{|t+1|} + \log{|t+3|} \right]_{1}^{e} \\
&= \left( 2\log{e} – 3\log{(e+1)} + \log{(e+3)} \right) – \left( 2\log{1} – 3\log{(1+1)} + \log{(1+3)} \right) \\
&= 2 – 3\log{(e+1)} + \log{(e+3)} + \log{2}
\end{align}
$$

上記より$(3)$が正しい。

・解説
途中計算では下記のような部分分数分解を行いました。
$$
\large
\begin{align}
\frac{6}{t(t+1)(t+3)} = \frac{2}{t} – \frac{3}{t+1} + \frac{1}{t+3}
\end{align}
$$

上記の分解を発見するにあたっては、$\displaystyle \frac{6}{t(t+1)(t+3)} = \frac{a}{t} + \frac{b}{t+1} + \frac{c}{t+3}$とおき、両辺が一致するように$a, b, c$の値を得れば良いです。
$$
\large
\begin{align}
\frac{a}{t} + \frac{b}{t+1} + \frac{c}{t+3} &= \frac{a(t+1)(t+2) + bt(t+3) + ct(t+1)}{t(t+1)(t+3)} \\
&= \frac{(a+b+c)t^{2} + (4a+3b+c)t + 3a}{t(t+1)(t+3)}
\end{align}
$$

上記より$a+b+c=0, 4a+3b+c=0, 3a=6$が成立するので、$a=2, b=-3, c=1$を得ることができます。

・解説
部分分数分解は下記で詳しく取り扱いました。

Q.14

$$
\large
\begin{align}
f(x) = \int_{0}^{x} f(t) \sin{(x-t)} dt + ax + b \quad [1]
\end{align}
$$

$\sin{(x-t)} = \sin{x}\cos{t} – \cos{x}\sin{t}$より、$[1]$式は下記のように変形できる。
$$
\large
\begin{align}
f(x) &= \int_{0}^{x} f(t) \sin{(x-t)} dt + ax + b \quad [1] \\
&= \int_{0}^{x} f(t) (\sin{x}\cos{t} – \cos{x}\sin{t}) dt + ax + b \\
&= \sin{x} \int_{0}^{x} f(t) \cos{t} dt – \cos{x} \int_{0}^{x} f(t) \sin{t} dt + ax + b \quad [2]
\end{align}
$$

$[2]$式の両辺を$x$で微分すると下記が得られる。
$$
\large
\begin{align}
f'(x) &= \left[ \sin{x} \int_{0}^{x} f(t) \cos{t} dt – \cos{x} \int_{0}^{x} f(t) \sin{t} dt + ax + b \right] \quad [2]’ \\
&= \cos{x} \int_{0}^{x} f(t) \cos{t} dt + \cancel{f(x)\sin{x}\cos{x}} + \sin{x} \int_{0}^{x} f(t) \sin{t} dt – \cancel{f(x)\sin{x}\cos{x}} + a \\
&= \cos{x} \int_{0}^{x} f(t) \cos{t} dt + \sin{x} \int_{0}^{x} f(t) \sin{t} dt + a
\end{align}
$$

上記の両辺を$x$で微分すると下記が得られる。
$$
\large
\begin{align}
f^{”}(x) &= -\sin{x} \int_{0}^{x} f(t) \cos{t} dt + f(x) \cos^{2}{x} + \cos{x} \int_{0}^{x} f(t) \sin{t} dt + f(x) \sin^{2}{x} \\
&= – \left[ \sin{x} \int_{0}^{x} f(t) \cos{t} dt – \cos{x} \int_{0}^{x} f(t) \sin{t} dt \right] + (\cos^{2}{x} + \sin^{2}{x})f(x) \\
&= -(f(x) – ax – b) + f(x) = ax + b
\end{align}
$$

$f^{”}(x) = ax+b$の両辺を$x$で積分すると下記が得られる。
$$
\large
\begin{align}
f'(x) = \frac{a}{2}x^{2} + bx + C_1
\end{align}
$$

ここで$f'(0)=a$より、$C_1=a$が得られるので、$f'(x)$は下記のように表せる。
$$
\large
\begin{align}
f'(x) = \frac{a}{2}x^{2} + bx + a
\end{align}
$$

$f'(x) = \frac{a}{2}x^{2}+bx+a$の両辺を$x$で積分すると下記が得られる。
$$
\large
\begin{align}
f(x) = \frac{a}{6}x^{3} + \frac{b}{2}x^{2} + ax + C_2
\end{align}
$$

ここで$f(0)=b$なので$C_2=b$が得られる。よって、$f(x)$は下記のように表せる。
$$
\large
\begin{align}
f(x) = \frac{a}{6}x^{3} + \frac{b}{2}x^{2} + ax + b
\end{align}
$$

上記より、$(3)$が正しい。

Q.15

$$
\large
\begin{align}
A+B &= \left(\begin{array}{cc} 2 & 2 \\ -3 & 3 \end{array} \right) \quad [1] \\
A-B &= \left(\begin{array}{cc} 4 & 2 \\ 1 & -3 \end{array} \right) \quad [2]
\end{align}
$$

$[1]+[2]$より下記が得られる。
$$
\large
\begin{align}
(A+B) + (A-B) &= \left(\begin{array}{cc} 2 & 2 \\ -3 & 3 \end{array} \right) + \left(\begin{array}{cc} 4 & 2 \\ 1 & -3 \end{array} \right) \\
2A &= \left(\begin{array}{cc} 6 & 4 \\ -2 & 0 \end{array} \right) \\
A &= \left(\begin{array}{cc} 3 & 2 \\ -1 & 0 \end{array} \right)
\end{align}
$$

同様に$[1]-[2]$より下記が得られる。
$$
\large
\begin{align}
(A+B) – (A-B) &= \left(\begin{array}{cc} 2 & 2 \\ -3 & 3 \end{array} \right) – \left(\begin{array}{cc} 4 & 2 \\ 1 & -3 \end{array} \right) \\
2B &= \left(\begin{array}{cc} -2 & 0 \\ -4 & 6 \end{array} \right) \\
B &= \left(\begin{array}{cc} -1 & 0 \\ -2 & 3 \end{array} \right)
\end{align}
$$

よって行列の積$AB$は下記のように得られる。
$$
\large
\begin{align}
AB &= \left(\begin{array}{cc} 3 & 2 \\ -1 & 0 \end{array} \right) \left(\begin{array}{cc} -1 & 0 \\ -2 & 3 \end{array} \right) \\
&= \left(\begin{array}{cc} -7 & 6 \\ 1 & 0 \end{array} \right)
\end{align}
$$

上記より$(2)$が正しい。

Q.16

$$
\large
\begin{align}
A = \left(\begin{array}{ccc} 0 & 1-i & 1 \\ 1+i & 0 & 1-i \\ 1 & 1+i & 0 \end{array} \right)
\end{align}
$$

上記に対し、固有多項式$F_{A}(\lambda)=\det{(\lambda I_{3}-A)}$は下記のように計算できる。
$$
\large
\begin{align}
F_{A}(\lambda) &= \det{(\lambda I_{3}-A)} = \left| \begin{array}{ccc} -\lambda & 1-i & 1 \\ 1+i & -\lambda & 1-i \\ 1 & 1+i & -\lambda \end{array} \right| \\
&= -\lambda^{3} + (1+i)^{2} + (1-i)^{2} – \left[ -\lambda – 2 \lambda(1+i)(1-i) \right] \\
&= -\lambda^{3} + 1+i^2+\cancel{2i} + 1+i^2-\cancel{2i} + \lambda + 2 \lambda(1-i^2) \\
&= -\lambda^{3} + \cancel{1}-\cancel{1}+ + \cancel{1}-\cancel{1} + \lambda + 2 \lambda \cdot (1+1) \\
&= -\lambda^{3} + 5 \lambda
\end{align}
$$

よって、固有方程式$F_{A}(\lambda)=\det{(\lambda I_{3}-A)}=0$は下記のように解ける。
$$
\large
\begin{align}
F_{A}(\lambda) = -\lambda^{3} + 5 \lambda &= 0 \\
\lambda(\lambda^2-5) &= 0 \\
\lambda &= 0, \, \pm \sqrt{5}
\end{align}
$$

ここで$\lambda_{1}<\lambda_{2}<\lambda_{3}$であるので、$\lambda_{1}=-\sqrt{5}, \, \lambda_{2}=0, \, \lambda_{3}=\sqrt{5}$である。よって$\displaystyle \frac{\lambda_{1}}{2} + \frac{\lambda_{2}}{3} + \frac{\lambda_{3}}{4}$は下記のように計算することができる。
$$
\large
\begin{align}
\frac{\lambda_{1}}{2} + \frac{\lambda_{2}}{3} + \frac{\lambda_{3}}{4} &= -\frac{\sqrt{5}}{2} + \frac{\sqrt{5}}{4} \\
&= -\frac{\sqrt{5}}{4}
\end{align}
$$

したがって$(5)$が正しい。

・解説
$3$次以上の正方行列の行列式の計算にあたっては余因子展開を用いることが多い一方で、この問題では余因子展開を用いると$i$の取り扱いにあたって式が複雑になります。一般に$n$次の正方行列の行列式の計算にあたって出てくる項の数は$n!$個であるので、$3$次の場合はうまく使い分けると良いです。

Q.17

$$
\large
\begin{align}
X = \left(\begin{array}{cc} x & y \\ z & u \end{array} \right)
\end{align}
$$

上記のように$X$を表すと、$AX=XA$は下記のように表すことができる。
$$
\large
\begin{align}
AX &= XA \\
\left(\begin{array}{cc} a & b \\ c & d \end{array} \right) \left(\begin{array}{cc} x & y \\ z & u \end{array} \right) &= \left(\begin{array}{cc} x & y \\ z & u \end{array} \right) \left(\begin{array}{cc} a & b \\ c & d \end{array} \right) \\
\left(\begin{array}{cc} ax+bz & ay+by \\ cx+dz & cy+du \end{array} \right) &= \left(\begin{array}{cc} ax+cy & bx+dy \\ az+cu & bz+du \end{array} \right) \quad [1]
\end{align}
$$

$[1]$式より下記の連立方程式が得られる。
$$
\large
\begin{align}
\cancel{ax}+bz &= \cancel{ax}+cy \quad [2] \\
ay+by &= bx+dy \quad [3] \\
cx+dz &= az+cu \quad [4] \\
cy+\cancel{du} &= bz+\cancel{du} \quad [5]
\end{align}
$$

ここで$[2], [5]$式より、$bz=cy$が得られるが、この式が任意の$y,z$で成立する必要十分条件は$b=c=0$である。また、$b=c=0$を$[3], [4]$式に代入すると$a=d$が得られる。したがって$(4)$が正しい。

・解説
$a=d, \, b=c=0$は行列$A$が単位行列の定数倍を表すことも合わせて抑えておくと良いです。

Q.18

下記のように極限値を得ることができる。
$$
\large
\begin{align}
\lim_{x \to 0} \frac{\sin{bx}}{\sin{ax}} &= \lim_{x \to 0} \frac{\sin{bx}}{bx} \cdot \frac{ax}{\sin{ax}} \cdot \frac{b}{a} \\
&= 1 \cdot 1 \cdot \frac{b}{a} \\
&= \frac{b}{a}
\end{align}
$$

よって$(3)$が正しい。

・解説
この問題では$\displaystyle \lim_{\theta \to 0} \frac{\sin{\theta}}{\theta}=1$を用いました。$\displaystyle \lim_{\theta \to 0} \frac{\sin{\theta}}{\theta}=1$の導出は下記で取り扱ったので合わせて確認しておくと良いです。

Q.19

$$
\large
\begin{align}
f(x,y) = \tan^{1}{(ax+by+1)}, \quad a \neq 0, \, b \neq 0
\end{align}
$$

$2$変数関数$f(x,y)$のマクローリン展開は下記のように表せる。
$$
\begin{align}
f(x,y) = f(0,0) + \frac{1}{1!} \left( x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} \middle) f(x,y) \right|_{x=0,y=0} + \frac{1}{2!} \left( x \frac{\partial}{\partial x} + y \frac{\partial}{\partial y} \middle)^{2} f(x,y) \right|_{x=0,y=0} + \cdots
\end{align}
$$

上記より$xy$の係数は$\displaystyle \frac{\partial^{2}}{\partial x \partial y} f(x,y) \Bigr|_{x=0,y=0}$を計算することで得られる。ここで$\displaystyle (\tan^{-1}{x})’ = \frac{1}{1+x^2}$であるので、$xy$の係数は下記のように得られる。
$$
\large
\begin{align}
\frac{\partial^{2}}{\partial x \partial y} f(x,y) \Bigr|_{x=0,y=0} &= \frac{-2ab(ax+by+1)}{[1+(ax+by+1)]^{2}} \Bigr|_{x=0,y=0} \\
&= \frac{-2ab \cdot (0+0+1)}{[1+(0+0+1)]^{2}} \\
&= -\frac{2ab}{4} = -\frac{ab}{2}
\end{align}
$$

・解説
$\displaystyle (\tan^{-1}{x})’ = \frac{1}{1+x^2}$の導出は下記で取り扱いました。

Q.20

$$
\large
\begin{align}
2020 \times \frac{\displaystyle \int_{0}^{1} (1-x^{20})^{100} dx}{\displaystyle \int_{0}^{1} (1-x^{20})^{101} dx} \quad [1]
\end{align}
$$

$\displaystyle I_{n} = \int_{0}^{1} (1-x^{20})^{n} dx, \, n \geq 1$とおくと、部分積分法により下記のように漸化式が得られる。
$$
\large
\begin{align}
I_{n} &= \int_{0}^{1} (1-x^{20})^{n} dx \\
&= \left[ x(1-x^{20})^{n} \right]_{0}^{1} – \int_{0}^{1} x \cdot n(1-x^{20})^{n-1} \cdot (-20x^{19}) dx \\
&= 20n \int_{0}^{1} x^{20}(1-x^{20})^{n-1} dx \\
&= 20n \int_{0}^{1} (1-x^{20})^{n-1} dx – 20n \int_{0}^{1} (1-x^{20}) (1-x^{20})^{n-1} dx \\
&= 20n \int_{0}^{1} (1-x^{20})^{n-1} dx – 20n \int_{0}^{1} (1-x^{20})^{n} dx \\
&= 20n I_{n-1} – 20n I_{n} \\
(20n+1) I_{n} &= 20n I_{n-1} \\
I_{n} &= \frac{20n}{20n+1} I_{n-1} \quad [2]
\end{align}
$$

$[2]$式に$n=101$を代入すると下記が得られる。
$$
\large
\begin{align}
I_{101} &= \frac{20 \cdot 101}{20 \cdot 101 + 1} I_{101-1} \quad [2]’ \\
I_{101} &= \frac{2020}{2021} I_{100} \quad [3]
\end{align}
$$

$[3]$式より、$[1]$式は下記のように計算できる。
$$
\large
\begin{align}
2020 \times \frac{\displaystyle \int_{0}^{1} (1-x^{20})^{100} dx}{\displaystyle \int_{0}^{1} (1-x^{20})^{101} dx} &= 2020 \times \frac{I_{100}}{I_{101}} \\
&= 2020 \cdot I_{100} \cdot \frac{2021}{2020 I_{100}} \\
&= 2021
\end{align}
$$

上記より$(5)$が正しい。

PPO(Proximal Policy Optimization)まとめ

方策勾配法の学習の安定化にあたっては、TRPO(Trust Region Policy Optimization)やPPO(Proximal Policy Optimization)のようにステップ幅の調整が解決策になります。当記事ではPPOについて詳しく取りまとめを行いました。
「ゼロから作るDeep Learning④ー強化学習編」の$10.2.3$の「TRPO・PPO」やTRPO、PPOの論文などを参考に作成を行いました。

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

・TRPO論文
・PPO論文

前提の確認

方策勾配法まとめ

方策勾配法の目的関数や勾配

軌道$\tau$の収益$G(\tau)$を最大にするような方策$\pi_{\theta}(A_t|S_t)$を得るにあたっては、下記の目的関数を$p$次元ベクトル$\theta$について最大化すればよい。
$$
\large
\begin{align}
J(\theta) &= \mathbb{E}_{\tau \sim \pi_{\theta}}[G(\tau)] \quad (1.1) \\
\theta & \in \mathbb{R}^{p}
\end{align}
$$

上記で定義した$J(\theta)$の勾配ベクトルは下記のように表される。
$$
\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.2)
\end{align}
$$

詳しい導出は下記で取り扱った。

REINFORCEMENT、ベースライン、Actor-Critic

$(1.1)$式では収益$G(\tau)$を元に目的関数$J(\theta)$を定義し、勾配ベクトル$\nabla_{\theta} J(\theta)$を$(1.2)$式で表した。$(1.2)$式は収益$G(\tau)$に基づいて軌道の重み付けを行なったと大まかに解釈できるが、本来的には状態$S_t$後の収益のみを用いるのがより適切である。このように$(1.2)$式の$G(\tau)$については様々な表し方が考えられるので、$(1.2)$式を$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 (1.3)
\end{align}
$$

上記の一例であるベースライン付きREINFORCEは、$\Phi_{t}$を下記のように定義する。
$$
\large
\begin{align}
\Phi_{t} &= G_{t} – b(S_t) \quad (1.4) \\
G_t &= R_t + \gamma R_{t+1} + \gamma^{2} R_{t+2} + \cdots + \gamma^{T-t} R_{T} \quad (1.5)
\end{align}
$$

ここで上記の収益$G_t$にQ関数$Q(S_t,A_t)$、ベースライン$b(S_t)$に価値関数を用いると$\Phi_{t}$はアドバンテージ関数$A(S_t,A_t)$に一致する。
$$
\large
\begin{align}
\Phi_{t} &= Q(S_t,A_t) – V(S_t) \\
&= A(S_t,A_t) \quad (1.6)
\end{align}
$$

アドバンテージ関数は状態$S_t$における行動$A_t$の評価を単体で行うことができるように定義される関数である。また、オーソドックスなREINFORCEやActor-Criticなどベースライン付きREINFORCE以外の$\Phi_{t}$の定義は下記で詳しく取りまとめた。

TRPO

TRPOの最適化問題は下記のように表される。
$$
\large
\begin{align}
\mathrm{maximize}_{\theta} &: \, \mathbb{E}_{t} \left[ \frac{\pi_{\theta}(a_t|s_t)}{\pi_{\theta_{\mathrm{old}}}(a_t|s_t)} A_{t} \right] \quad (1.7) \\
\mathrm{subject \, to} &: \, \mathbb{E}_{t} \left[ KL[\pi_{\theta_{\mathrm{old}}}(\cdot|s_t),\pi_{\theta}(\cdot|s_t)] \right] \leq \delta \quad (1.8) \\
\theta & \in \mathbb{R}^{p}, \, \theta_{\mathrm{old}} \in \mathbb{R}^{p}
\end{align}
$$

$(1.7)$式を$(1.8)$式の制約を元に最適化するにあたって、TRPOでは共役勾配法と直線探索を用いる。詳細は下記で取り扱った。

$(1.8)$式の制約により、繰り返し演算における変化の幅の上限値を設定することが可能で学習の安定化が実現できる一方で、計算コストがやや大きい。PPOはこの計算コストを解決した手法である。

PPO

PPOの目的関数と学習の流れ

Clipped Surrogate Objective

TRPOにおける目的関数(Objective function)を”conservative policy iteration”に基づいて$L^{CPI}(\theta)$とおくと、$L^{CPI}(\theta)$は確率の比率である$r_{t}(\theta)$を用いて、下記のように表すこともできる。
$$
\large
\begin{align}
L^{CPI}(\theta) &= \mathbb{E}_{t} \left[ \frac{\pi_{\theta}(a_t|s_t)}{\pi_{\theta_{\mathrm{old}}}(a_t|s_t)} A_{t} \right] = \mathbb{E}_{t} \left[ r_{t}(\theta) A_t \right] \quad (2.1) \\
r_{t}(\theta) &= \frac{\pi_{\theta}(a_t|s_t)}{\pi_{\theta_{\mathrm{old}}}(a_t|s_t)} \quad (2.2)
\end{align}
$$

$(2.1)$式に対して、PPOでは下記のような目的関数の$L^{CPI}(\theta)$を用いる。
$$
\large
\begin{align}
L^{CLIP}(\theta) &= \mathbb{E}_{t} \left[ \min \left( r_{t}(\theta)A_t, \, \mathrm{clip}(r_{t}(\theta), 1-\varepsilon, 1+\varepsilon)A_t \right) \right] \quad (2.3)
\end{align}
$$

$(2.3)$式は$\min$の内側に下記の二つの項を持つ。
$$
\large
\begin{align}
& r_{t}(\theta)A_t \\
& \mathrm{clip}(r_{t}(\theta), 1-\varepsilon, 1+\varepsilon)A_t
\end{align}
$$

第$1$項を選ぶ場合、$(2.3)$式は$(2.1)$式に一致するので$L^{CLIP}(\theta)$は$L^{CPI}(\theta)$に一致する。

第$2$項は$r_{t}$の変化を$1-\varepsilon$以上$1+\varepsilon$以下に抑える関数である。第$2$項があることで、一度のパラメータのUpdateにおけるパラメータの変化が大きくなり過ぎないように制約を課すことができる。

$(2.3)$式は下記の図を元に理解すると良い。

PPO論文 Figure$1$

図の理解にあたっては、「$1-\varepsilon$と$1+\varepsilon$のどちらを用いるかはアドバンテージ関数$A_t$の符号による」ことに注意する必要がある。アドバンテージ関数が正であるときは$1+\varepsilon$を用い、負の場合は$1-\varepsilon$を用いる。

Adaptive KL Penalty Coefficient

PPO論文では目的関数$L^{CPI}(\theta)$(conservative policy iteration)に加えて、下記のような目的関数$L^{KLPEN}(\theta)$に基づいた実験の結果の記載がある。
$$
\large
\begin{align}
L^{KLPEN}(\theta) &= \mathbb{E}_{t} \left[ \frac{\pi_{\theta}(a_t|s_t)}{\pi_{\theta_{\mathrm{old}}}(a_t|s_t)} A_{t} \, – \, \beta \mathrm{KL}\left[ \pi_{\theta_{\mathrm{old}}}(\cdot|s_t), \pi_{\theta}(\cdot|s_t) \right] \right] \quad (2.4) \\
d &= \mathbb{E}_{t} \left[ \mathrm{KL}\left[ \pi_{\theta_{\mathrm{old}}}(\cdot|s_t), \pi_{\theta}(\cdot|s_t) \right] \right] \quad (2.5)
\end{align}
$$

目的関数$L^{KLPEN}(\theta)$の”KLPEN”は”KL penalty”の略である。また、$\beta$の値は固定(fixed)する場合と、$(2.5)$式の$d$の値に基づいて適応的に(Adaptive)変化させる方法の二通りが取り扱われる。

適応的な手法を用いる場合はハイパーパラメータのd_targの値を元に、下記のようなプログラムに基づいて$\beta$の値を変化させる。

beta = 
d_targ = 

for i in range(n_epoch):
    d = calc_KL()
    if d < d_targ/1.5:
        beta = beta/2
    elif d > d_targ*1.5:
        beta = beta*2

上記の処理は$(2.5)$式に基づいて計算を行ったKLダイバージェンスの値が小さければペナルティの係数$\beta$を小さくし、KLダイバージェンスの値が大きければペナルティの係数$\beta$を大きくすることに対応する。よって、「確率分布の変動が大きい場合にペナルティを大きくする」と解釈することができる。また、betaの初期値もハイパーパラメータであるが、アルゴリズムがすぐに適応するので値は重要でないとPPO論文に記載がある。

PPOの学習の流れ

$(2.3), (2.4)$式の目的関数を元に、PPOの学習の流れは下記のように表される。

PPO論文 Algorithm$1$

PPO論文における実験

PPO論文では下記の数式に基づいてそれぞれ実験が行われる。

PPO論文 $6.1$ Comparison of Surrogate Objectives より

上記のClippingは「Clipped Surrogate Objective」、KL penaltyは「Adaptive KL Penalty Coefficient」で取り扱った内容にそれぞれ対応する。このそれぞれの目的関数を用いて、OpenAI Gymで実装されているHalfCheetah、Hopper、InvertedDoublePendulum、InvertedPendulum、Reacher、Swimmer、Walker2dの$7$つのsimulated robotic tasksについて学習を行った結果が下記である。

PPO論文 Table$1$

上記は$7$つのタスクのスコアの正規化された平均であり、$\epsilon=0.2$のClippingが最も良い結果であることが確認できる。また、Adaptive KLの$\beta$の初期値はここでは$\beta=1$に設定される。

PPO論文 Figure$3$

また、論文のFigure$3$ではPPOに加えてA$2$CやTRPOなど、その他アルゴリズムの学習の結果についてまとめられているので合わせて確認しておくと良い。紫がPPOであり、相対的に良い結果であることが確認できる。

・注意事項
PPO論文では「$\epsilon=0.2$のClipping」が最も良い結果とされるが、InstructGPT/ChatGPTでは「Clipped Surrogate Objective」ではなく、「Adaptive KL Penalty Coefficient」に基づく目的関数が用いられることに注意が必要である。

TRPO(Trust Region Policy Optimization)まとめ

方策勾配法の学習の安定化にあたっては、TRPO(Trust Region Policy Optimization)やPPO(Proximal Policy Optimization)のようにステップ幅の調整が解決策になります。当記事ではTRPOとPPOについて取りまとめました。
「ゼロから作るDeep Learning④ー強化学習編」の$10.2.3$の「TRPO・PPO」やTRPO、PPOの論文などを参考に作成を行いました。

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

・TRPO論文
・PPO論文

前提の確認

方策勾配法まとめ

方策勾配法の目的関数や勾配

軌道$\tau$の収益$G(\tau)$を最大にするような方策$\pi_{\theta}(A_t|S_t)$を得るにあたっては、下記の目的関数を$p$次元ベクトル$\theta$について最大化すればよい。
$$
\large
\begin{align}
J(\theta) &= \mathbb{E}_{\tau \sim \pi_{\theta}}[G(\tau)] \quad (1.1) \\
\theta & \in \mathbb{R}^{p}
\end{align}
$$

上記で定義した$J(\theta)$の勾配ベクトルは下記のように表される。
$$
\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.2)
\end{align}
$$

詳しい導出は下記で取り扱った。

REINFORCEMENT、ベースライン、Actor-Critic

$(1.1)$式では収益$G(\tau)$を元に目的関数$J(\theta)$を定義し、勾配ベクトル$\nabla_{\theta} J(\theta)$を$(1.2)$式で表した。$(1.2)$式は収益$G(\tau)$に基づいて軌道の重み付けを行なったと大まかに解釈できるが、本来的には状態$S_t$後の収益のみを用いるのがより適切である。このように$(1.2)$式の$G(\tau)$については様々な表し方が考えられるので、$(1.2)$式を$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 (1.3)
\end{align}
$$

上記の一例であるベースライン付きREINFORCEは、$\Phi_{t}$を下記のように定義する。
$$
\large
\begin{align}
\Phi_{t} &= G_{t} – b(S_t) \quad (1.4) \\
G_t &= R_t + \gamma R_{t+1} + \gamma^{2} R_{t+2} + \cdots + \gamma^{T-t} R_{T} \quad (1.5)
\end{align}
$$

ここで上記の収益$G_t$にQ関数$Q(S_t,A_t)$、ベースライン$b(S_t)$に価値関数を用いると$\Phi_{t}$はアドバンテージ関数$A(S_t,A_t)$に一致する。
$$
\large
\begin{align}
\Phi_{t} &= Q(S_t,A_t) – V(S_t) \\
&= A(S_t,A_t) \quad (1.6)
\end{align}
$$

アドバンテージ関数は状態$S_t$における行動$A_t$の評価を単体で行うことができるように定義される関数である。また、オーソドックスなREINFORCEやActor-Criticなどベースライン付きREINFORCE以外の$\Phi_{t}$の定義は下記で詳しく取りまとめた。

共役勾配法

対称行列$A$と$A$の固有ベクトルを$\mathbf{v}_{1}, \, \mathbf{v}_{2}$と定義する。ここで$\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 (1.7)
\end{align}
$$

TRPOでは最適化における各ステップでの修正の方向の計算にあたって、通常の勾配ベクトルではなく上記に基づいて計算を行なった$\mathbf{d}_{2}$を用いる。共役勾配法の式の詳しい解釈に関しては下記で詳しく取り扱った。

直線探索

$$
\large
\begin{align}
\mathbf{x}_{k+1} &= \mathbf{x}_{k} – \alpha_{k} \nabla f(\mathbf{x}_{k}) \quad (1.8) \\
\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.8)$式は勾配降下法の式であるが、このような式の$\alpha_{k}$を$f(\mathbf{x}_{k} – \alpha_{k} \nabla f(\mathbf{x}_{k}))$が最大値を取るように設定するのが直線探索(line search)である。

直線探索のプログラム作成にあたってはnp.argmaxなどを用いるとプログラムをシンプルに作成できる。詳しくは下記で取り扱った。

フィッシャー情報行列

$p$次元パラメータベクトル$\theta \in \mathbb{R}^{p}$を用いて尤度関数を$L(\theta|x)$のように定義するとき、フィッシャー情報行列(FIM; Fisher Information Matrix)の$I(\theta)$は下記のように定義される。
$$
\large
\begin{align}
I(\theta) &= E \left[ \frac{\partial}{\partial \theta} \log{L(\theta|x)} \frac{\partial}{\partial \theta^{\mathrm{T}}} \log{L(\theta|x)} \right] \\
&= E \left[ \left(\begin{array}{c} \displaystyle \frac{\partial \log{L(\theta|x)}}{\partial \theta_{1}} \\ \vdots \\ \displaystyle \frac{\partial \log{L(\theta|x)}}{\partial \theta_{p}} \end{array} \right) \left(\begin{array}{ccc} \displaystyle \frac{\partial \log{L(\theta|x)}}{\partial \theta_{1}} & \cdots & \displaystyle \frac{\partial \log{L(\theta|x)}}{\partial \theta_{p}} \end{array} \right) \right]
\end{align}
$$

詳しくは下記で取り扱った。

Hessian-Free Optimization

$$
\large
\begin{align}
\nabla = \frac{\partial}{\partial \theta} = \left(\begin{array}{c} \displaystyle \frac{\partial}{\partial \theta_1} \\ \vdots \\ \displaystyle \frac{\partial}{\partial \theta_p} \end{array} \right)
\end{align}
$$

上記のように方向微分の演算子$\nabla$を定義するとき、ヘッセ行列$H$とベクトル$\mathbf{u}$の積$H \mathbf{u}$について下記のような近似が成立する。
$$
\large
\begin{align}
H \mathbf{u} & \simeq \frac{\nabla f(\theta + \varepsilon \mathbf{u}) – \nabla f(\theta)}{\varepsilon} \quad (1.10)
\end{align}
$$

上記の式の導出については下記で詳しく取り扱った。

TRPO

TRPOの目的関数と制約

TRPO(Trust Region Policy Optimization)では下記のように目的関数(objective function)と制約を定義する。
$$
\large
\begin{align}
\mathrm{maximize}_{\theta} &: \, \mathbb{E}_{t} \left[ \frac{\pi_{\theta}(a_t|s_t)}{\pi_{\theta_{\mathrm{old}}}(a_t|s_t)} A_{t} \right] \quad (2.1) \\
\mathrm{subject \, to} &: \, \mathbb{E}_{t} \left[ KL[\pi_{\theta_{\mathrm{old}}}(\cdot|s_t),\pi_{\theta}(\cdot|s_t)] \right] \leq \delta \quad (2.2) \\
\theta & \in \mathbb{R}^{p}, \, \theta_{\mathrm{old}} \in \mathbb{R}^{p}
\end{align}
$$

$(2.1)$式は$(1.1)$式の$G(\tau)$をアドバンテージ関数で置き換え、さらに重点サンプリング(importance sampling)を行うにあたって式変形を行なったものである。$(2.1)$式のように定義する目的関数をTRPOやPPOの論文ではsurrogate objectiveということも抑えておくと良い。

上記の最適化にあたっては、重点サンプリングを元に目的関数の近似値を計算し、近似値の計算結果を元に共役勾配法直線探索(line search)に基づく繰り返し演算を行う。詳細は次項以降で取り扱う。

最適化問題の定式化

前項の「TRPOの目的関数と制約」の$(2.1), \, (2.2)$式を下記のように式の簡略化を行う。
$$
\large
\begin{align}
\mathrm{Maximize} & : \, L(\theta) \quad (2.1)’ \\
\mathrm{subject \, to} & : \, \overline{D}_{KL}(\theta_{\mathrm{old}}, \theta) \leq \delta \quad (2.2)’ \\
\theta & \in \mathbb{R}^{p}, \, \theta_{\mathrm{old}} \in \mathbb{R}^{p}
\end{align}
$$

TRPOでは上記の最適化問題を下記のように制約式ではなくペナルティを用いて再定義する。
$$
\large
\begin{align}
\mathrm{Maximize} & : \, L(\theta) – \lambda \overline{D}_{KL}(\theta_{\mathrm{old}}, \theta) \quad (2.3)
\end{align}
$$

制約やペナルティに用いられるKLダイバージェンスは下記のように二次近似(quadratic approximation)を行う。
$$
\large
\begin{align}
\overline{D}_{KL}(\theta_{\mathrm{old}}, \theta) & \simeq \frac{1}{2} (\theta_{\mathrm{old}}, \theta)^{\mathrm{T}} A (\theta_{\mathrm{old}} – \theta) \quad (2.4) \\
A_{ij} &= \frac{\partial^{2}}{\partial \theta_{i} \partial \theta_{j}} \overline{D}_{KL}(\theta_{\mathrm{old}}, \theta) \quad (2.5) \\
\theta & \in \mathbb{R}^{p}, \, \theta_{\mathrm{old}} \in \mathbb{R}^{p}
\end{align}
$$

フィッシャー情報行列は前節の$(1.9)$式のように定義されることが一般的であるが、TRPO論文では$(2.5)$式によって$A$の要素の推定が行われる。同様に目的関数$L(\theta)$について一次近似を行い方向微分を計算すると、$Ax=g$形式の方程式が得られるので、共役勾配法を用いて解けばよい。

共役勾配法と線形探索

「共役勾配法」「Hessian-Free Optimization」で確認した計算に基づいてパラメータ$\theta \in \mathbb{R}^{p}$をUpdateする方向の$\mathbf{d}_{t+1}$が得られた際に、$\theta_{t+1} = \theta_{t} + \beta \mathbf{d}_{t+1}$によって$\theta_{t+1}$の計算を行う。この際に制約式における$\delta$は下記に対応する。
$$
\large
\begin{align}
\delta = \overline{D}_{KL}(\theta_{\mathrm{old}}, \theta) \simeq \frac{1}{2}(\beta \mathbf{d}_{t+1})^{\mathrm{T}} A (\beta \mathbf{d}_{t+1}) = \frac{1}{2} \beta^{2} \mathbf{d}_{t+1}^{\mathrm{T}} A \mathbf{d}_{t+1} \quad (2.6)
\end{align}
$$

$(2.6)$式を$\beta$について解くと下記が得られる。
$$
\large
\begin{align}
\delta & \simeq \frac{1}{2} \beta^{2} \mathbf{d}_{t+1}^{\mathrm{T}} A \mathbf{d}_{t+1} \quad (2.6) \\
\beta^{2} & \simeq \frac{2 \delta}{\mathbf{d}_{t+1}^{\mathrm{T}} A \mathbf{d}_{t+1}} \\
\beta & \simeq \sqrt{\frac{2 \delta}{\mathbf{d}_{t+1}^{\mathrm{T}} A \mathbf{d}_{t+1}}} \quad (2.7)
\end{align}
$$

したがって、線形探索にあたっての$\beta$による探索幅の最大値を$\displaystyle \sqrt{\frac{2 \delta}{\mathbf{d}_{t+1}^{\mathrm{T}} A \mathbf{d}_{t+1}}}$に設定することで、$(2.2)$式の制約の範囲で線形探索を行うことができる。

フィッシャー情報行列(FIM; Fisher Information Matrix)の式定義

フィッシャー情報行列(FIM; Fisher Information Matrix)は多変数スカラー関数の二次近似(quadratic approximation)を行う際に計算を行う行列です。当記事ではフィッシャー情報行列の定義式について取りまとめました。
当記事は「Wikipedia:フィッシャー情報量」などを参考に作成を行いました。
https://ja.wikipedia.org/wiki/フィッシャー情報量

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

フィッシャー情報量の定義

確率関数や確率密度関数を$f(x|\theta)=P(X=x|\theta)$のように定義するとき、パラメータ$\theta$の尤度$L(\theta|x)$は下記のように表せる。
$$
\large
\begin{align}
L(\theta|x) = f(x|\theta)
\end{align}
$$

このとき、対数尤度$\log{L(\theta|x)}$を元にスコア関数$V(x|\theta)$を下記のように定義する。
$$
\large
\begin{align}
V(x|\theta) = \frac{\partial}{\partial \theta} \log{L(\theta|x)}
\end{align}
$$

上記で定義した対数尤度やスコア関数を元に、フィッシャー情報量$I(\theta)$は下記のように定義される。
$$
\large
\begin{align}
I(\theta) &= E \left[ V(x|\theta)^{2} \right] \\
&= E \left[ \left( \frac{\partial}{\partial \theta} \log{L(\theta|x)} \right)^{2} \right]
\end{align}
$$

フィッシャー情報行列の定義

下記のように$p$次元パラメータベクトル$\boldsymbol{\theta} \in \mathbb{R}^{p}$と、方向微分の演算子$\nabla$を定義する。
$$
\large
\begin{align}
\boldsymbol{\theta} &= \left(\begin{array}{c} \theta_{1} \\ \vdots \\ \theta_{p} \end{array} \right) \\
\nabla &= \left(\begin{array}{c} \displaystyle \frac{\partial}{\partial \theta_{1}} \\ \vdots \\ \displaystyle \frac{\partial}{\partial \theta_{p}} \end{array} \right)
\end{align}
$$

このときフィッシャー情報行列を$I(\boldsymbol{\theta})$と定義すると、尤度関数$L(\boldsymbol{\theta}|x)$に基づいて$I(\boldsymbol{\theta})$は下記のように表される。
$$
\large
\begin{align}
I(\boldsymbol{\theta}) &= E \left[ \nabla \log{L(\boldsymbol{\theta}|x)} \nabla^{\mathrm{T}} \log{L(\boldsymbol{\theta}|x)} \right] \\
&= E \left[ \frac{\partial}{\partial \boldsymbol{\theta}} \log{L(\boldsymbol{\theta}|x)} \frac{\partial}{\partial \boldsymbol{\theta}^{\mathrm{T}}} \log{L(\boldsymbol{\theta}|x)} \right] \\
&= E \left[ \left(\begin{array}{c} \displaystyle \frac{\partial \log{L(\boldsymbol{\theta}|x)}}{\partial \theta_{1}} \\ \vdots \\ \displaystyle \frac{\partial \log{L(\boldsymbol{\theta}|x)}}{\partial \theta_{p}} \end{array} \right) \left(\begin{array}{ccc} \displaystyle \frac{\partial \log{L(\boldsymbol{\theta}|x)}}{\partial \theta_{1}} & \cdots & \displaystyle \frac{\partial \log{L(\boldsymbol{\theta}|x)}}{\partial \theta_{p}} \end{array} \right) \right]
\end{align}
$$

【上級】データサイエンス 数学ストラテジスト 公式問題集 解答例まとめ Q.1〜10

「データサイエンス 数学ストラテジスト 上級」はデータサイエンスの基盤である、確率・統計、線形代数、微積分、機械学習、プログラミングなどを取り扱う資格試験です。当記事では「日本数学検定協会」作成の「公式問題集」の演習問題$1$〜$10$の解答例を取り扱いました。

・数学検定まとめ
https://www.hello-statisticians.com/math_certificate

演習問題

Q.1

$$
\large
\begin{align}
x^{2} + \frac{1}{x^2} = 10
\end{align}
$$

上記の式は下記のように変形できる。
$$
\large
\begin{align}
x^{2} + \frac{1}{x^2} &= 10 \\
\left( x + \frac{1}{x} \right)^{2} – 2 &= 10 \\
\left( x + \frac{1}{x} \right)^{2} &= 12
\end{align}
$$

ここで$x>1$より、$\displaystyle x + \frac{1}{x} = 2 \sqrt{3}$が成立するが、下記のように変形を行える。
$$
\large
\begin{align}
x + \frac{1}{x} &= 2 \sqrt{3} \\
x^{2} – 2 \sqrt{3} x + 1 &= 0 \\
x &= \sqrt{3} \pm \sqrt{3-1} = \sqrt{3} \pm \sqrt{2}
\end{align}
$$

$x>1$なので$x=\sqrt{2}+\sqrt{3}$が成立する。

・解説
$x>1$を元に$\pm$を外す際に注意が必要です。$\displaystyle x + \frac{1}{x} > 1$や$\sqrt{3}-\sqrt{2}<1$を途中計算で用いました。

Q.2

$n$角形の$1$つの点に対し、線分が対角線である点は$n-3$、辺である点は$2$つ存在する。ここで$n-3=2$であれば$p=q$である。よって$n=5$であれば良い。

・解説
問題集の解答では$2$点選ぶことを${}_{n} C_{2}$で表現されますが、このように表すと式がやや複雑です。$1$点はどの点を選んでも他の点に関する取り扱いは同じであることに着目することで計算を簡略化することが可能です。

Q.3

$$
\large
\begin{align}
y &= \log_{2}{\frac{1}{x}} \\
&= \log_{2}{x^{-1}} \\
&= -\log_{2}{x}
\end{align}
$$

上記より$(3)$のグラフが正しい。

Q.4

$$
\large
\begin{align}
\sqrt{10}x^2 – 2x + k = 0
\end{align}
$$

上記の方程式の解が$\sin{\theta}, \, \cos{\theta}$であるとき、二次方程式の解と係数の関係より下記が成立する。
$$
\large
\begin{align}
\sin{\theta} + \cos{\theta} &= \frac{2}{\sqrt{10}} \quad [1] \\
\sin{\theta} \cos{\theta} &= \frac{k}{\sqrt{10}} \quad [2]
\end{align}
$$

ここで$(1)$式の両辺の$2$乗は下記のように変形できる。
$$
\large
\begin{align}
(\sin{\theta} + \cos{\theta})^{2} &= \frac{2^2}{\sqrt{10}^2} \quad [1]’ \\
\sin^{2}{\theta} + \cos^{2}{\theta} + 2 \sin{\theta} \cos{\theta} &= \frac{2}{5} \\
1 + 2 \sin{\theta} \cos{\theta} &= \frac{2}{5} \\
2 \sin{\theta} \cos{\theta} &= -\frac{3}{5} \\
\sin{\theta} \cos{\theta} &= -\frac{3}{10} \quad [3]
\end{align}
$$

$[2], \, [3]$式より下記のように$k$の値が得られる。
$$
\large
\begin{align}
\frac{k}{\sqrt{10}} &= -\frac{3}{10} \\
k &= -\frac{3}{\sqrt{10}}
\end{align}
$$

上記より、元の二次方程式は下記のように表せる。
$$
\large
\begin{align}
\sqrt{10}x^2 – 2x + -\frac{3}{\sqrt{10}} &= 0 \\
10x^2 – 2 \sqrt{10} x + -3 &= 0
\end{align}
$$

二次方程式の解の公式より、上記の方程式の解は下記のように得られる。
$$
\large
\begin{align}
x &= \frac{\sqrt{10} \pm \sqrt{\sqrt{10}^{2} + 10 \cdot 3}}{10} \\
&= \frac{\sqrt{10} \pm 2 \sqrt{10}}{10} \\
&= -\frac{\sqrt{10}}{10}, \, \frac{3 \sqrt{10}}{10}
\end{align}
$$

上記より$(1)$が正しい。

Q.5

$$
\large
\begin{align}
\left( \frac{2}{x} \right)^{\log_{e}{2}} &= \left( \frac{3}{y} \right)^{\log_{e}{3}} \quad [1] \\
3^{-\log_{e}{x}} &= 2^{-\log_{e}{y}} \quad [2]
\end{align}
$$

上記に対し、$\displaystyle X=\frac{1}{x}, \, Y=\frac{1}{y}$のようにおくと、$[1], \, [2]$式は下記のように表せる。
$$
\large
\begin{align}
(2X)^{\log_{e}{2}} &= (3Y)^{\log_{e}{3}} \quad [1]’ \\
3^{\log_{e}{X}} &= 2^{\log_{e}{Y}} \quad [2]’
\end{align}
$$

上記の$[1]’$式の両辺の自然対数を取ると、下記が得られる。
$$
\large
\begin{align}
(2X)^{\log_{e}{2}} &= (3Y)^{\log_{e}{3}} \quad [1]’ \\
\log_{e}{2} \log_{e}{(2X)} &= \log_{e}{3} \log_{e}{(3Y)} \\
\log_{e}{2} (\log_{e}{2} + \log_{e}{X}) &= \log_{e}{3} (\log_{e}{3} + \log_{e}{Y}) \quad [1]^{”}
\end{align}
$$

同様に$[2]’$式の両辺の自然対数を取ると、下記が得られる。
$$
\large
\begin{align}
3^{\log_{e}{X}} &= 2^{\log_{e}{Y}} \quad [2]’ \\
\log_{e}{X} \log_{e}{3} &= \log_{e}{Y} \log_{e}{2} \\
\log_{e}{Y} &= \frac{\log_{e}{3}}{\log_{e}{2}} \log_{e}{X} \quad [2]^{”}
\end{align}
$$

ここで$[2]^{”}$式を$[1]^{”}$式に代入することで下記のように変形できる。
$$
\large
\begin{align}
\log_{e}{2} (\log_{e}{2} + \log_{e}{X}) &= \log_{e}{3} (\log_{e}{3} + \log_{e}{Y}) \quad [1]^{”} \\
\log_{e}{2} (\log_{e}{2} + \log_{e}{X}) &= \log_{e}{3} \left( \log_{e}{3} + \frac{\log_{e}{3}}{\log_{e}{2}} \log_{e}{X} \right) \\
\frac{1}{\log_{e}{2}} \cancel{\left( (\log_{e}{2})^{2} – (\log_{e}{3})^{2} \right)} \log_{e}{X} &= -\cancel{\left( (\log_{e}{2})^{2} – (\log_{e}{3})^{2} \right)} \\
\log_{e}{X} &= -\log_{e}{2} \\
\log_{e}{X} &= \log_{e}{2^{-1}} \\
X &= \frac{1}{2}
\end{align}
$$

上記を$[2]^{”}$式に代入することで下記が得られる。
$$
\large
\begin{align}
\log_{e}{Y} &= \frac{\log_{e}{3}}{\log_{e}{2}} \log_{e}{X} \quad (2)^{”} \\
\log_{e}{Y} &= \frac{\log_{e}{3}}{\log_{e}{2}} \cdot -\log_{e}{2} \\
\log_{e}{Y} &= -\log_{e}{3} \\
\log_{e}{Y} &= \log_{e}{3^{-1}} \\
Y &= \frac{1}{3}
\end{align}
$$

ここで$\displaystyle X=\frac{1}{x}, \, Y=\frac{1}{y}$であるので、$x=2, y=3$である。よって$x^2+y^2=13$であるので$(1)$が正しい。

Q.6

$$
\large
\begin{align}
f(x) = x|x-3|
\end{align}
$$

上記より、$x<3$のとき$f(x)=-x(x-3)=-x^2+3x, \, f'(x)=-2x+3$、$3 \leq x$のとき$f(x)=x(x-3)=x^2-3x, \, f'(x)=2x-3$である。よって$f'(2)-f'(4)$は下記のように計算できる。
$$
\large
\begin{align}
f'(2) – f'(4) &= (-2 \cdot 2 + 3) – (2 \cdot 4 – 3) \\
&= -1 – 5 \\
&= -6
\end{align}
$$

上記より$(1)$が正しい。

Q.7

$$
\large
\begin{align}
f(x) = \int_{-1}^{x} (t^2-t-2) dt
\end{align}
$$

$g(t)=t^2-t-2$の原始関数を$G(t)$とおくと、$f'(x)$は下記のように得ることができる。
$$
\large
\begin{align}
f'(x) &= \frac{d}{dx} \int_{-1}^{x} g(t) dt \\
&= \frac{d}{dx} (G(x)-G(-1)) \\
&= g(x) \\
&= x^2 – x – 2 \\
&= (x+1)(x-2)
\end{align}
$$

上記より、関数$f(x)$は$x=-1$で極大値$M$、$x=2$で極小値$m$を取る。$M$と$m$はそれぞれ下記のように計算できる。
$$
\large
\begin{align}
M &= f(-1) = \int_{-1}^{-1} (t^2-t-2) dt \\
&= 0 \\
m &= f(2) = \int_{-1}^{2} (t^2-t-2) dt \\
&= \left[ \frac{1}{3}t^{3} – \frac{1}{2}t^{2} – 2t \right]_{-1}^{2} \\
&= \left( \frac{8}{3} – 2 – 4 \right) – \left( -\frac{1}{3} – \frac{1}{2} + 2 \right) \\
&= \frac{16-36+2+3-12}{6} \\
&= -\frac{27}{6} \\
&= -\frac{9}{2}
\end{align}
$$

上記より$(4)$が正しい。

Q.8

$\displaystyle X \sim \mathrm{Bin} \left( 90, \frac{1}{3} \right)$であるので、期待値$E[X]$と標準偏差$SD[X]=\sqrt{V[X]}$はそれぞれ下記のように計算できる。
$$
\large
\begin{align}
E[X] &= 90 \times \frac{1}{3} \\
&= 30 \\
SD[X] &= \sqrt{V[X]} \\
&= \sqrt{90 \times \frac{1}{3} \times \frac{2}{3}} \\
&= \sqrt{20} = 2 \sqrt{5}
\end{align}
$$

上記より、$(5)$が正しい。

Q.9

$$
\large
\begin{align}
I_{A} &= -\log_{2}{p_{A}} \\
p_{A} &= \frac{2}{10} \times \frac{1}{9} \\
&= \frac{1}{45}
\end{align}
$$

上記より選択情報量$I_{A}$は下記のように計算できる。
$$
\large
\begin{align}
I_{A} &= -\log_{2}{\frac{1}{45}} \\
&= \log_{2}{(3^{2} \times 5)} \\
&= 2 \log_{2}{3} + \log_{2}{5} \\
&= 2 \times 1.585 + 2.322 \\
&= 5.492
\end{align}
$$

よって$(3)$が正しい。

Q.10

$$
\large
\begin{align}
\vec{a} = \left(\begin{array}{c} 3 \\ -2 \\ 1 \end{array} \right), \quad \vec{b} = \left(\begin{array}{c} -6 \\ k \\ -2 \end{array} \right)
\end{align}
$$

$\vec{a}$と$\vec{b}$が平行のとき、$-2\vec{a} = \vec{b}$が成立する必要がある。このとき$k=4$である。逆に$k=4$のとき$\vec{a}$と$\vec{b}$が平行であるので、$k_1=4$が必要十分である。

$\vec{a} \perp \vec{b}$のとき、$\vec{a} \cdot \vec{b} = 0$より下記が成立する。
$$
\large
\begin{align}
\vec{a} \cdot \vec{b} &= 0 \\
-18 – 2k -2 &= 0 \\
k &= -10
\end{align}
$$

よって$k_2=-10$である。したがって$(3)$が正しい。

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

過去問

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

[1]解答

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

$(A)$ 触媒,$(B)$ 温度,$(A\times B)$ 触媒と温度の交互作用 と表すことにする.交互作用を誤差にプールすると,$S’_E = S_E + S_{A\times B} = 346 + 168 = 514$ .$\phi’ E = \phi _E + \phi_{A\times B} = 13$ であるから,$V’_E = 514/13 = 39.54$ .これは$V_E = 38.4$ より大きいので,触媒,温度それぞれの $F$ 値はいずれも小さくなる.

$①$. : それぞれの $F$ 値は小さくなるため,$F$ 値の和も変化する.誤り.
$②$:$F$ 値と残差の自由度が変化するため $P$ 値も変化する.誤り.
$③$:$②$ と同様.誤り.
$④$:正しい.
$⑤$:残差の自由度は変化する.誤り.

[2]解答

$\boxed{\mathsf{13}}$ : $④$
触媒と温度の交互作用は有意でないため,$① \sim ③$ は誤り.
触媒について,$F_A = 0.624 < F_{0.05}(2,9) = 4.26$ であるから,有意水準 $5\%$ で触媒の効果は有意でない.温度については,$F_B = 4.370 > F_{0.05}(2,9) = 4.26$ であり,有意水準 $5\%$ で有意である.温度についての表をみると,最適温度は $450°c$ ,点推定値は $88\%$ であることがわかる.
したがって,正しいのは $④$ .

注意
[1] に交互作用は有意でないと判断したと記述があるため,$① \sim ③$ は誤り.
$④$,$⑤$ については,有意水準 $1\%$ で考えると触媒,温度のどちらも有意でない.選択肢の中から適切なものを選ぶ必要があるため,有意水準$5\%$ で考えた結果を採用している.

直線探索(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を行うことで得られる。