多変数関数の定義とテイラー展開
連続最適化では、目的関数が多変数関数であり、出力はスカラー値です。この多変数関数を $ f(\mathbf{x}) $ と表し、ベクトル $ \mathbf{x} $ によって定義されます。
多変数関数のテイラー展開は、特定の点での関数の値とその周辺での動きを表す有用なツールです。展開式は次のようになります：二次形式で近似した場合
f(\mathbf{x} + \Delta \mathbf{x}) \approx f(\mathbf{x}) + \nabla f(\mathbf{x})^T \Delta \mathbf{x} + \frac{1}{2} \Delta \mathbf{x}^T \nabla^2 f(\mathbf{x}) \Delta \mathbf{x}
ここで、$ \nabla f(\mathbf{x}) $ は勾配ベクトル、$ \nabla^2 f(\mathbf{x}) $ はヘッシアン（二次の微分情報を含む行列）を表します。
この形式は、二次形式（quadratic forms）と非常に類似しています。二次形式は、ベクトルと行列を用いて表される数学的表現で、一般的には以下の形式を持ちます。
q(\mathbf{x}) = \mathbf{x}^T A \mathbf{x}
ここで、$ A $ は対称行列であり、$ \mathbf{x} $ は任意のベクトルです。
テイラー展開の二次の項 $\frac{1}{2} \Delta \mathbf{x}^T \nabla^2 f(\mathbf{x}) \Delta \mathbf{x}$ は、二次形式の特殊なケースと見なすことができます。ここで、$ A = \nabla^2 f(\mathbf{x}) $ となり、この項が関数の曲率を表すのに対し、二次形式はベクトルの長さや方向に依存した値を生成します。この類似性により、関数の局所的な挙動を解析しやすくなり、特に最適化の文脈で有用です。
ちょっとまって、関数ってベクトルで微分できるの？
勾配 (Gradient)
勾配は、多変数スカラー関数（つまり、複数の変数を持ち、一つの数値を出力する関数）の最も基本的な微分形式です。
\nabla f(\mathbf{x}) = \begin{bmatrix}
\frac{\partial f}{\partial x_1} \\
\frac{\partial f}{\partial x_2} \\
\vdots \\
\frac{\partial f}{\partial x_n}
\end{bmatrix}
この勾配ベクトルは、関数が最も急激に増加する方向を指し示します。つまり、勾配が示す方向に進むことで、関数の値は最も速く増加します。
ヘッセ行列 (Hessian)
ヘッセ行列は、スカラー関数の二次導関数を含む行列です。関数 $ f(\mathbf{x}) $ のヘッセ行列は以下のように表されます。
\nabla^2 f(\mathbf{x}) = \begin{bmatrix}
\frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\
\frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n^2}
\end{bmatrix}
ラプラシアン (Laplacian)
ラプラシアンはスカラー関数の二次導関数のトレースであり、すべての二次偏微分の和です。次元を一般化して表すと、ラプラシアンは以下のようになります。
\Delta f = \sum_{i=1}^{n} \frac{\partial^2 f}{\partial x_i^2}
ここで、$ n $ は次元数を表します。ラプラシアンは関数がどのように平坦化するか、またはどのように極値を持つかを理解するのに役立ちます。
ヘッセ行列とラプラシアンの関係性は次のとおりです：
\Delta f = \text{trace}(\nabla^2 f)
つまり、ラプラシアンはヘッセ行列の対角成分の和に等しいです。
ヤコビアン (Jacobian)
ヤコビアンは、ベクトル値関数（複数の変数から複数の出力を持つ関数）の微分を一般化したものです。関数 $ \mathbf{f}(\mathbf{x}) = (f_1(\mathbf{x}), f_2(\mathbf{x}), ..., f_n(\mathbf{x})) $ に対して、そのヤコビアンは以下のような行列です。
J = \begin{bmatrix}
\frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_m} \\
\vdots & \ddots & \vdots \\
\frac{\partial f_n}{\partial x_1} & \cdots & \frac{\partial f_n}{\partial x_m}
\end{bmatrix}
この行列は、ベクトル値関数が各変数に対してどのように変化するかを表します。
ヘッシアン行列（ヘッセ行列）は、勾配ベクトルのヤコビアンとして理解することができます。
ニュートン法の適用
ニュートン法では、上記のテイラー展開を用いて、$ f $ の局所的な最小値を見つけるための反復式を導出します。具体的には、二次の項までを考慮して、次のように反復式を設定します。
\mathbf{x}_{\text{new}} = \mathbf{x} - [\nabla^2 f(\mathbf{x})]^{-1} \nabla f(\mathbf{x})
この式は、$ \nabla f(\mathbf{x}) = 0 $ となる $ \mathbf{x} $ を見つけることを目的としています。ヘッシアンが正定値行列であれば、この点での $ f $ は局所的な最小値を持つことになります。
ニュートン法ってあの関数の根を求めるやつ？
やっていることはほとんど同じですが、文脈が違うと内容も変わります。連続最適化におけるニュートン法は、$ f $ そのものではなく $ f $ の勾配 $ \nabla f(\mathbf{x}) = 0 $ になる点を求めることで、その積分値である多変数関数の極値とその位置を見つける作業になります。
最小二乗法問題をニュートン法で解く
最小二乗法は統計学やデータサイエンスで広く用いられる方法ですが、この問題を連続最適化の観点からニュートン法を用いて解く方法を探ります。特に、線形最小二乗問題の場合、ニュートン法は非常に効果的です。
最小二乗問題の目的関数 $f(\mathbf{x})$ は次のように定義されます：
f(\mathbf{x}) = \|A\mathbf{x} - \mathbf{b}\|^2
ここで、$A$ はデザイン行列、$\mathbf{b}$ は観測された出力ベクトルです。この関数を展開すると、以下の二次形式が得られます：
f(\mathbf{x}) = (A\mathbf{x} - \mathbf{b})^T(A\mathbf{x} - \mathbf{b}) = \mathbf{x}^T A^T A \mathbf{x} - 2\mathbf{b}^T A \mathbf{x} + \mathbf{b}^T \mathbf{b}
勾配とヘッセ行列
この関数の勾配とヘッセ行列は以下のように計算されます：
- 勾配（グラディエント）
\nabla f(\mathbf{x}) = 2A^T A \mathbf{x} - 2A^T \mathbf{b}
- ヘッセ行列:
H = 2A^T A
ニュートン法の更新式は以下のようになります：
\mathbf{x}_{\text{new}} = \mathbf{x} - H^{-1} \nabla f(\mathbf{x})
解の導出
ヘッセ行列 $H$ は定数行列 $2A^T A$ であり、これが逆行列を持つと仮定すると（$A$ がフルランクの場合）、更新式は次のように単純化されます：
\mathbf{x}_{\text{new}} = \mathbf{x} - (2A^T A)^{-1} (2A^T A \mathbf{x} - 2A^T \mathbf{b})
\mathbf{x}_{\text{new}} = \mathbf{x} - \mathbf{x} + (A^T A)^{-1} A^T \mathbf{b}
\mathbf{x}_{\text{new}} = (A^T A)^{-1} A^T \mathbf{b}
以下に、擬似逆行列を用いた解とニュートン法による解の比較を追加したQiita記事のセクションを加筆しました。
ヘッセ行列の逆行列と分散共分散行列
楕円体の形状は $ \nabla^2 f(\mathbf{x}) $ の固有値と固有ベクトルによって決定されます。ヘッセ行列の逆行列 $ (\nabla^2 f(\mathbf{x}))^{-1} $ は、楕円体の軸の長さを制御し、分散共分散行列の役割を果たします。
多変量ガウス分布との関係
1. 関数の二次形式での近似
関数 $ f(\mathbf{x}) $ のある点 $ \mathbf{a} $ でのテイラー展開は以下のように表されます（ここでは2次まで展開）:
f(\mathbf{x}) \approx f(\mathbf{a}) + (\mathbf{x} - \mathbf{a})^T \nabla f(\mathbf{a}) + \frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a})
ここで、$ \mathbf{H} $ は $ f $ のヘッセ行列（二階偏微分の行列）です。最適化問題などでよく使われるのは $ \mathbf{a} $ が極値（例えば最小値）である場合で、この時 $ \nabla f(\mathbf{a}) = \mathbf{0} $ となります。したがって、
f(\mathbf{x}) \approx f(\mathbf{a}) + \frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a})
2. 負の対数尤度関数としての表現
尤度関数 $ L(\mathbf{x}) $ が $ \exp(-f(\mathbf{x})) $ であると仮定しましょう。このとき、$ f(\mathbf{x}) $ が上記のように近似される場合、負の対数尤度は
-f(\mathbf{x}) \approx -f(\mathbf{a}) - \frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a})
3. 正規化
正規化定数 $ C $ を使って、尤度の積分が1になるように調整します。$ C $ は次の積分で求められます:
C = \int \exp\left(-\frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a})\right) d\mathbf{x}
ガウス積分の性質から、$ C $ は $ \mathbf{H} $ の行列式 $ |\mathbf{H}| $ と次元に依存します。結果として、
C = (2\pi)^{n/2} |\mathbf{H}|^{-1/2}
4. 多変量ガウス分布
結局、次の式で示される $ L(\mathbf{x}) $ は、平均 $ \mathbf{a} $ と共分散行列 $ \mathbf{H}^{-1} $ を持つ多変量ガウス分布です：
$$ L(\mathbf{x}) = \frac{|\mathbf{H}|^{1/2}}{(2\pi)^{n/2}} \exp\left(-\frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a})\right). $$
多次元ガウス分布の一般的な形式は次のようになります：
$$ p(\mathbf{x}) = \frac{1}{(2\pi)^{n/2} |\Sigma|^{1/2}} \exp\left(-\frac{1}{2} (\mathbf{x} - \mu)^T \Sigma^{-1} (\mathbf{x} - \mu)\right). $$
ここで、平均 $\mu = \mathbf{a}$ とし、共分散行列を $\Sigma = \mathbf{H}^{-1}$ と置き換えると、上記の導出した式と一致します。このため、ヘッセ行列の逆行列が共分散行列と一致し、見慣れた多次元ガウス分布の形になります。
多次元ガウス分布は、簡単に $ N(\mu, \Sigma) $ と表すことができます。
1次元の場合
1次元の場合、関数 $ f(x) $ を点 $ a $ でテイラー展開し、負の対数尤度を考えると以下のようになります：
$$ f(x) \approx f(a) + \frac{1}{2} f''(a)(x - a)^2. $$
ここで、$ f''(a) $ は $ f $ の2階導関数です。負の対数尤度は次のようになります：
$$ -f(x) \approx -f(a) - \frac{1}{2} f''(a)(x - a)^2. $$
正規化定数 $ C $ を用いて尤度の積分が1になるように調整すると、以下のように積分の結果を得ます：
$$ C = \int \exp\left(-\frac{1}{2} f''(a)(x - a)^2\right) dx = \sqrt{\frac{2\pi}{f''(a)}}. $$
この定数を用いて、尤度関数 $ L(x) $ は次のように記述できます：
$$ L(x) = \frac{\sqrt{f''(a)}}{\sqrt{2\pi}} \exp\left(-\frac{1}{2} f''(a)(x - a)^2\right). $$
これは、平均 $ a $ と分散 $ (f''(a))^{-1} $ を持つ正規分布の形です。正規分布の一般的な形式は次のようになります：
$$ p(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right). $$
ここで、平均 $\mu = a$ とし、分散を $ \sigma^2 = (f''(a))^{-1} $ と置き換えると、先ほど導出した式と一致することがわかります。
2次元の場合
2次元の場合、関数 $ f(\mathbf{x}) $ は、点 $\mathbf{a} = (a_1, a_2)$ で次のように展開されます：
$$ f(\mathbf{x}) \approx f(\mathbf{a}) + \frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a}), $$
ここで、$\mathbf{H}$ は $ f $ のヘッセ行列（2階微分行列）です。この場合の正規化定数 $ C $ は次の通りです：
$$ C = \int \exp\left(-\frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a})\right) d\mathbf{x} = (2\pi)^{n/2} |\mathbf{H}|^{-1/2}. $$
この定数を使うことで、多変量ガウス分布（多変数正規分布） $ L(\mathbf{x}) $ は次のように表されます：
$$ L(\mathbf{x}) = \frac{|\mathbf{H}|^{1/2}}{(2\pi)^{n/2}} \exp\left(-\frac{1}{2} (\mathbf{x} - \mathbf{a})^T \mathbf{H} (\mathbf{x} - \mathbf{a})\right). $$
この式は、平均 $\mathbf{a}$ と共分散行列 $\mathbf{H}^{-1}$ を持つ多次元ガウス分布です。
多次元ガウス分布の一般的な形式は次のようになります：
$$ p(\mathbf{x}) = \frac{1}{(2\pi)^{n/2} |\Sigma|^{1/2}} \exp\left(-\frac{1}{2} (\mathbf{x} - \mu)^T \Sigma^{-1} (\mathbf{x} - \mu)\right). $$
ここで、平均 $\mu = \mathbf{a}$ とし、共分散行列を $\Sigma = \mathbf{H}^{-1}$ と置き換えると、上記の導出した式と一致します。このため、ヘッセ行列の逆行列が共分散行列と一致し、2次元でも見慣れた多次元ガウス分布の形が現れます。
最尤推定とニュートン法の関係
最尤推定は、与えられたデータに対してパラメータを最適化するための一般的な手法であり、負の対数尤度を最小化することを目的とします。これは、最適化アルゴリズムであるニュートン法と密接な関係があります。ニュートン法では、関数の2階微分情報（ヘッセ行列）を利用して、次の推定値を効率的に計算します。以下では、最尤推定の導出をニュートン法との関係に着目して解説します。
平均ベクトルの推定
多変量正規分布のサンプル集合 $\mathbf{X} = {\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_m}$ に対する対数尤度関数は以下のように表されます：
$$ L(\boldsymbol{\mu}, \Sigma) = -\frac{mn}{2}\log(2\pi) - \frac{m}{2}\log |\Sigma| - \frac{1}{2}\sum_{i=1}^m (\mathbf{x}_i - \boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{x}_i - \boldsymbol{\mu}). $$
この対数尤度関数を $\boldsymbol{\mu}$ について最大化するため、導関数を0に設定します。すなわち、
$$ \nabla_{\boldsymbol{\mu}} L = \Sigma^{-1}\sum_{i=1}^m (\mathbf{x}_i - \boldsymbol{\mu}) = 0. $$
この方程式は、ニュートン法の1回の更新ステップのように見えます。式を $\boldsymbol{\mu}$ について解くと、次のような算術平均で表される推定値が得られます：
$$ \hat{\boldsymbol{\mu}} = \frac{1}{m}\sum_{i=1}^m \mathbf{x}_i. $$
ここで、ニュートン法を適用するイメージでは、1回のステップでこの解に到達します。これは、平均ベクトルについての最尤推定が線形で簡単なためです。
共分散行列の推定
次に、共分散行列 $\Sigma$ について同様に対数尤度関数の導関数を0に設定します。
$$ \nabla_{\Sigma} L = -\frac{m}{2} \Sigma^{-1} + \frac{1}{2} \Sigma^{-1}\left(\sum_{i=1}^m (\mathbf{x}_i - \boldsymbol{\mu})(\mathbf{x}_i - \boldsymbol{\mu})^T\right)\Sigma^{-1} = 0. $$
この方程式を $\Sigma$ について解くと、共分散行列の推定値 $\hat{\Sigma}$ は以下のように表されます：
$$ \hat{\Sigma} = \frac{1}{m}\sum_{i=1}^m (\mathbf{x}_i - \hat{\boldsymbol{\mu}})(\mathbf{x}_i - \hat{\boldsymbol{\mu}})^T. $$
共分散行列の最尤推定もまたニュートン法と類似した考え方で、2階導関数（ヘッセ行列）に基づく情報を活用して推定されます。最尤推定の結果は、ニュートン法が収束するポイントと同じになります。
https://seetheworld1992.hatenablog.com/entry/2019/07/01/103515
https://qiita.com/AnchorBlues/items/8fe2483a3a72676eb96d
混合ガウス分布と最尤推定の複雑性
混合ガウス分布は、複数のガウス分布の加重和でモデル化される確率分布です。これは、データが異なるガウス分布から生成されるという仮定に基づいています。数学的には、$K$個のガウス分布の混合で表され、各ガウス分布は特定の平均 $\boldsymbol{\mu}_k$、共分散行列 $\Sigma_k$、そして混合係数（混合比率を表す重み）$\pi_k$ を持ちます。混合ガウス分布の確率密度関数は次のように表されます：
$$
p(\mathbf{x}) = \sum_{k=1}^K \pi_k \mathcal{N}(\mathbf{x}; \boldsymbol{\mu}_k, \Sigma_k)
$$
ここで、$\mathcal{N}(\mathbf{x}; \boldsymbol{\mu}_k, \Sigma_k)$は平均 $\boldsymbol{\mu}_k$ と共分散 $\Sigma_k$ を持つ多変量正規分布の確率密度関数です。
最尤推定の問題点
混合ガウス分布のパラメータを最尤推定する際には、対数尤度関数を最大化する必要があります。この対数尤度は次のように表されます：
$$
L(\pi, \boldsymbol{\mu}, \Sigma) = \sum_{i=1}^m \log \left( \sum_{k=1}^K \pi_k \mathcal{N}(\mathbf{x}_i; \boldsymbol{\mu}_k, \Sigma_k) \right)
$$
この式での主な計算上の問題は、$\log$の中の和（log-sum-exp）です。これは以下の理由で計算が難しくなります：
- 数値的安定性: $\mathcal{N}(\mathbf{x}_i; \boldsymbol{\mu}_k, \Sigma_k)$が非常に小さくなる場合、数値的に不安定になることがあります。
- 導関数の計算: 各パラメータに関する導関数を求める際に、内部の和に対して外部の$\log$関数を微分しなければならず、これが複雑で計算コストが高くなります。
Log-Sum-Exp関数の有用性
一方で、log-sum-expは計算上の問題を引き起こす一方で、その数学的特性から非常に便利なツールとしても使用されます。特に、
-
Smooth maximum approximation: Log-sum-expは複数の値の中での最大値を「滑らかに」近似することができます。つまり、複数の項の最大値を取る代わりに、それらを柔軟に統合して一つの値にすることが可能です。これにより、最大値関数の微分不可能な点を避けることができ、最適化アルゴリズムで扱いやすくなります。
-
多峰性のモデリング: 混合ガウス分布のように、複数の異なる統計的集団を一つのモデルで表現する場合、log-sum-expはこれらの異なる分布を効果的に統合し、データの多峰性を自然に表現することができます。
このように、log-sum-exp関数は計算の複雑性を増す一方で、モデルの表現力と最適化の柔軟性を大きく向上させるため、統計学や機械学習で広く用いられています。