多次元テイラー展開・多次元ニュートン法と連立方程式の近似解の計算

一様最強力不偏検定で両側検定の両側の棄却域を計算のような連立方程式はそのまま解くことができないので、多次元ニュートン法が適用されます。当記事では多次元テイラー展開を用いた多次元ニュートン法の導出や、多次元ニュートン法を用いた連立方程式の近似解の計算について取り扱いました。

・$1$次元のニュートン法
https://www.hello-statisticians.com/optimization/newton_taylor1.html

テイラー展開

$1$変数関数のテイラー展開

関数$f(x)$に関して点$x=a$を中心とするテイラー展開は下記のように表せる。
$$
\large
\begin{align}
f(x) &= \sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x-a)^{n} \\
&= \frac{f(a)}{0!}(x-a)^{0} + \frac{f'(a)}{1!} (x-a)^{1} + \frac{f^{”}(a)}{2!} (x-a)^{2} + \cdots \\
&= f(a) + f'(a)(x-a) + \frac{f^{”}(a)}{2} (x-a)^{2} + \cdots
\end{align}
$$

多変数関数のテイラー展開

以下具体的に表記を行うにあたって$2$変数関数$f_1(x_1,x_2)$と$f_2(x_1,x_2)$を元にテイラー展開の式を確認する。$(x_1,x_2)=(a_1,a_2)$を中心とするテイラー展開の式は下記のように表せる。
$$
\large
\begin{align}
\left(\begin{array}{c} f_1(x_1,x_2) \\ f_2(x_1,x_2) \end{array} \right) &= \left(\begin{array}{c} f_1(a_1,a_2) \\ f_2(a_1,a_2) \end{array} \right) + \left(\begin{array}{cc} \displaystyle \frac{\partial f_1}{\partial x_1} & \displaystyle \frac{\partial f_1}{\partial x_2} \\ \displaystyle \frac{\partial f_2}{\partial x_1} & \displaystyle \frac{\partial f_2}{\partial x_2} \end{array} \right) \left(\begin{array}{c} x_1-a_1 \\ x_2-a_2 \end{array} \right) + \cdots \quad (1) \\
&= \left(\begin{array}{c} f_1(a_1,a_2) \\ f_2(a_1,a_2) \end{array} \right) + \left(\begin{array}{c} \displaystyle \frac{\partial f_1}{\partial x_1}(x_1-a_1) + \frac{\partial f_1}{\partial x_2}(x_2-a_2) \\ \displaystyle \frac{\partial f_2}{\partial x_1}(x_1-a_1) + \frac{\partial f_2}{\partial x_2}(x_2-a_2) \end{array} \right) + \cdots
\end{align}
$$

ここで偏微分が含まれる行列を下記のように定義する。
$$
\large
\begin{align}
J(x_1,x_2) = \left(\begin{array}{cc} \displaystyle \frac{\partial f_1}{\partial x_1} & \displaystyle \frac{\partial f_1}{\partial x_2} \\ \displaystyle \frac{\partial f_2}{\partial x_1} & \displaystyle \frac{\partial f_2}{\partial x_2} \end{array} \right)
\end{align}
$$

上記の行列をヤコビ行列という。ここでは$2$変数関数を取り扱ったが、多変数関数でも同様な式で表される。

多次元ニュートン法の導出

$1$変数関数のテイラー近似とニュートン法の導出

$$
\large
\begin{align}
f(x_{n+1}) \simeq f(x_n) + f'(x_n)(x_{n+1}-x_{n}) = 0 \quad (2)
\end{align}
$$

$1$変数関数のニュートン法は上記の式を$x_{n+1}$ついて解くことで下記のように得られる。
$$
\large
\begin{align}
x_{n+1} = x_{n} – \frac{f(x_{n})}{f'(x_{n})}
\end{align}
$$

詳しい式変形は下記で取り扱ったので省略を行なった。

多次元ニュートン法の導出

前項の$(2)$式では$x_{n+1}-x_{n}$について式を解くことで$\displaystyle x_{n+1} = x_{n} – \frac{f(x_{n})}{f'(x_{n})}$を得ることができた。前節の$(1)$式についても同様な式変形を下記のように行うことができる。
$$
\begin{align}
\left(\begin{array}{c} f_1(x_1^{(n+1)},x_2^{(n+1)}) \\ f_2(x_1^{(n+1)},x_2^{(n+1)}) \end{array} \right) \simeq \left(\begin{array}{c} f_1(x_1^{(n)},x_2^{(n)}) \\ f_2(x_1^{(n)},x_2^{(n)}) \end{array} \right) + & \left(\begin{array}{cc} \displaystyle \frac{\partial f_1}{\partial x_1} & \displaystyle \frac{\partial f_1}{\partial x_2} \\ \displaystyle \frac{\partial f_2}{\partial x_1} & \displaystyle \frac{\partial f_2}{\partial x_2} \end{array} \middle) \right|_{x_1=x_1^{(n)},x_2=x_2^{(n)}} \left(\begin{array}{c} x_1^{(n+1)}-x_1^{(n)} \\ x_2^{(n+1)}-x_2^{(n)} \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \end{array} \right) \\
\left(\begin{array}{c} f_1(x_1^{(n)},x_2^{(n)}) \\ f_2(x_1^{(n)},x_2^{(n)}) \end{array} \right) + J(x_1^{(n)},x_2^{(n)}) \left(\begin{array}{c} x_1^{(n+1)}-x_1^{(n)} \\ x_2^{(n+1)}-x_2^{(n)} \end{array} \right) &= \left(\begin{array}{c} 0 \\ 0 \end{array} \right) \\
\left(\begin{array}{c} x_1^{(n+1)} \\ x_2^{(n+1)} \end{array} \right) &= \left(\begin{array}{c} x_1^{(n)} \\ x_2^{(n)} \end{array} \right) – J(x_1^{(n)},x_2^{(n)})^{-1} \left(\begin{array}{c} f_1(x_1^{(n)},x_2^{(n)}) \\ f_2(x_1^{(n)},x_2^{(n)}) \end{array} \right) \quad (3)
\end{align}
$$

上記は$2$変数の場合を表したが、$3$変数以上の場合も同様に導出することができる。

連立方程式の解

$$
\large
\begin{align}
\left(\begin{array}{c} f_1(x_1,x_2) \\ f_2(x_1,x_2) \end{array} \right) = \left(\begin{array}{c} 0.9 – \exp(-x_1) + \exp(-x_2) \\ x_1\exp(-x_1) – x_2\exp(-x_2) \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \end{array} \right)
\end{align}
$$

上記に対し、ヤコビ行列$J(x_1,x_2)$は下記のように計算できる。
$$
\large
\begin{align}
J(x_1,x_2) &= \left(\begin{array}{cc} \displaystyle \frac{\partial f_1}{\partial x_1} & \displaystyle \frac{\partial f_1}{\partial x_2} \\ \displaystyle \frac{\partial f_2}{\partial x_1} & \displaystyle \frac{\partial f_2}{\partial x_2} \end{array} \right) \\
&= \left(\begin{array}{cc} \exp{(-x_1)} & -\exp{(-x_2)} \\ (1-x_1)\exp{(-x_1)} & -(1-x_2)\exp{(-x_2)} \end{array} \right)
\end{align}
$$

上記の連立方程式の解は$(3)$式を元に下記のように計算を行うことができる。

import numpy as np

def calc_f1(x_1, x_2):
    return 0.9 - np.exp(-x_1) + np.exp(-x_2)

def calc_f2(x_1, x_2):
    return x_1*np.exp(-x_1) - x_2*np.exp(-x_2)

x = np.array([[0.7],[1.5]])

for i in range(10):
    J = np.array([[np.exp(-x[0,0]), -np.exp(-x[1,0])], [(1.-x[0,0])*np.exp(-x[0,0]), -(1.-x[1,0])*np.exp(-x[1,0])]])
    x = x - np.linalg.solve(J, np.array([[calc_f1(x[0,0],x[1,0])], [calc_f2(x[0,0],x[1,0])]]))

print(x)

・実行結果

array([[ 0.08381479],
       [ 3.93214595]])

上記の結果が適切であることは、下記を実行することで確認することができる。

print(calc_f1(0.08381479, 3.93214595))
print(calc_f2(0.08381479, 3.93214595))

実行結果がほぼ$0$であるので、連立方程式の近似解の$x_1 \simeq 0.083, x_2 \simeq 3.932$が得られたと解釈できる。