ニュートン法とカルマンフィルタ


## 1. ニュートン法

目的関数 $f: \mathbb{R}^n \to \mathbb{R}$ の最小化を考えます：

```math
\text{minimize} \quad f(x), \quad x \in \mathbb{R}^n
```

ニュートン法による更新式：

```math
dx = -[\nabla^2 f(x)]^{-1} \nabla f(x)
```

ここで、$\nabla f(x) \in \mathbb{R}^n$ は勾配、$\nabla^2 f(x) \in \mathbb{R}^{n \times n}$ はヘッセ行列を表します。

### 2.1 統計多様体とフィッシャー計量

確率分布の空間は統計多様体として捉えることができ、その上でフィッシャー情報行列が自然な計量を与えます。フィッシャー情報計量の各成分は以下で与えられます：

```math
g_{ij}(\theta) = E_p\left[\frac{\partial \log p(x|\theta)}{\partial \theta_i}\frac{\partial \log p(x|\theta)}{\partial \theta_j}\right]
```

これらの成分からなるフィッシャー情報行列 G(θ) は以下のように表されます：

```math
G(\theta) = \begin{bmatrix}
g_{11}(\theta) & g_{12}(\theta) & \cdots & g_{1n}(\theta) \\
g_{21}(\theta) & g_{22}(\theta) & \cdots & g_{2n}(\theta) \\
\vdots & \vdots & \ddots & \vdots \\
g_{n1}(\theta) & g_{n2}(\theta) & \cdots & g_{nn}(\theta)
\end{bmatrix}
```

この計量は、確率分布間の局所的な距離を定義します。

### 2.2 KLダイバージェンスと二次近似

KLダイバージェンスの二次近似は、フィッシャー情報行列と密接に関連します：

```math
D_{KL}(p_{\theta+\delta\theta}||p_{\theta}) = \frac{1}{2}\delta\theta^T G(\theta)\delta\theta + O(||\delta\theta||^3)
```

ここで、G(θ)はフィッシャー情報行列です。この形式は以下の重要な対応関係を示唆します：

1. ニュートン法における目的関数の二次近似：
```math
f(x + \delta x) \approx f(x) + \nabla f(x)^T\delta x + \frac{1}{2}\delta x^T H(x)\delta x
```

2. ガウス分布の対数尤度：
```math
-\log p(x|\mu,\Sigma) = \frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu) + \text{const.}
```

これらの式は全て二次形式で表現され、それぞれG(θ)、H(x)、Σ^{-1}が二次項の係数行列として現れます。特に、フィッシャー情報行列G(θ)はヘッセ行列H(x)と同様の役割を果たし、両者とも最適化における局所的な曲率を表現します。

## 3. ニュートン法の情報幾何学的解釈

### 3.1 古典的なニュートン法と自然勾配法の統一的理解

ニュートン法と自然勾配法は、統計多様体上で本質的に同じ振る舞いをします：

1. **ニュートン法**：ユークリッド空間での二次近似に基づく更新
   ```math
   x_{k+1} = x_k - [H_f(x_k)]^{-1}\nabla f(x_k)
   ```

2. **自然勾配法**：統計多様体上での最適な更新
   ```math
   \theta_{k+1} = \theta_k - [G(\theta_k)]^{-1}\nabla L(\theta_k)
   ```




## 4. カルマンフィルタのKLダイバージェンス最小化としての導出

### 4.1 カルマンフィルタの前提と目的

カルマンフィルタは、真の事後分布 $p(x|y)$ とガウス分布 $q_\theta(x)$ の間のKLダイバージェンスを最小化する自然勾配法として解釈できます。具体的には：

1. **最適化問題としての定式化**
   ```math
   \min_{q \in \mathcal{N}} D_{KL}(q(x) || p(x|y))
   ```

2. **自然勾配法としての解釈**
   - パラメータ更新は統計多様体上での最適な探索方向に従います
   - フィッシャー情報行列 $G(\theta)$ が計量として機能し、更新式は：
   ```math
   \theta_{k+1} = \theta_k - G(\theta_k)^{-1}\nabla_\theta D_{KL}(q_\theta||p)
   ```


### 4.2 ガウス分布の場合の具体的な計算

KLダイバージェンスの二次近似は、より厳密には以下のように表現されます：

```math
D_{KL}(p_{\theta+\delta\theta}||p_{\theta}) = \frac{1}{2}\delta\theta^T G(\theta)\delta\theta + O(||\delta\theta||^3)
```

事後分布の推定問題に対する最適化問題は：

```math
\min_{q \in \mathcal{N}} D_{KL}(q(x) || p(x|y))
```

ここで、$\mathcal{N}$ はガウス分布の族を表します。目的関数は：

```math
f(\theta) = D_{KL}(q_\theta||p) \quad \text{where} \quad \theta = (\mu, \Sigma) \in \mathbb{R}^n \times \mathbb{S}^n_{++}
```

ここで、$\mathbb{S}^n_{++}$ は n×n の正定値対称行列の空間を表します。

目的関数 f(θ) のヘッセ行列は以下の形を取ります：

```math
H_k = \begin{bmatrix} 
\Sigma_k^{-1} + H_k^TR_k^{-1}H_k & 0 \\
0 & -\frac{1}{2}\Sigma_k^{-1} \otimes \Sigma_k^{-1}
\end{bmatrix}
```

ここで $H_k$ は観測モデル $h(x)$ の時刻kにおけるヤコビアン行列です。

ヘッセ行列の逆行列を明示的に書き下すと、ブロック行列の逆行列の公式から：

```math
H_k^{-1} = \begin{bmatrix} 
(\Sigma_k^{-1} + H_k^TR_k^{-1}H_k)^{-1} & 0 \\
0 & -2(\Sigma_k \otimes \Sigma_k)
\end{bmatrix}
```

そして、勾配ベクトル $g_k$ は以下の形を取ります：

```math
g_k = \begin{bmatrix} 
-\Sigma_k^{-1}(\mu_k-\mu_{k|k-1}) - H_k^TR_k^{-1}(y_k-h(\mu_k)) \\
-\frac{1}{2}\Sigma_k^{-1} + \frac{1}{2}(\Sigma_k^{-1} + H_k^TR_k^{-1}H_k)
\end{bmatrix}
```

ニュートン法の更新式 $-H_k^{-1}g_k$ を計算することで、以下の更新式が同時に得られます：

```math
\begin{bmatrix} 
\mu_{k+1} - \mu_k \\
\Sigma_{k+1} - \Sigma_k
\end{bmatrix} = -\begin{bmatrix} 
(\Sigma_k^{-1} + H_k^TR_k^{-1}H_k)^{-1}[-\Sigma_k^{-1}(\mu_k-\mu_{k|k-1}) - H_k^TR_k^{-1}(y_k-h(\mu_k))] \\
(\Sigma_k \otimes \Sigma_k)[-\Sigma_k^{-1} + (\Sigma_k^{-1} + H_k^TR_k^{-1}H_k)]
\end{bmatrix}
```

ここで、μ に関する更新式の変形にWoodburyの公式を適用します：

```math
(\Sigma^{-1} + H^TR^{-1}H)^{-1}H^TR^{-1} = \Sigma H^T(H\Sigma H^T + R)^{-1}
```

この変形により、計算効率が大幅に改善されます。なぜなら：
- 変形前：状態空間の次元（n×n行列）の逆行列計算が必要
- 変形後：観測空間の次元（m×m行列）の逆行列計算で済む（通常 m < n）

この結果を用いて、時系列更新式として以下が得られます：

```math
\begin{align}
\mu_{k+1} &= \mu_k + \Sigma_kH^T(H\Sigma_kH^T + R)^{-1}(y - h(\mu_k)) \\
\Sigma_{k+1} &= (I - \Sigma_kH^T(H\Sigma_kH^T + R)^{-1}H)\Sigma_k
\end{align}
```

ここで、$\Sigma_kH^T(H\Sigma_kH^T + R)^{-1}$ をカルマンゲイン $K$ と定義すると：

```math
K = \Sigma_kH^T(H\Sigma_kH^T + R)^{-1}
```

最終的に、カルマンフィルタの標準的な更新式が得られます：

```math
\begin{align}
\mu_{k+1} &= \mu_k + K(y - h(\mu_k)) \\
\Sigma_{k+1} &= (I - KH)\Sigma_k
\end{align}
```

ここで $K = \Sigma_kH^T(H\Sigma_kH^T + R)^{-1}$ はカルマンゲインです。

> **補足：Woodburyの公式**
> 
> Woodburyの公式（または Sherman-Morrison-Woodburyの公式）は、以下の行列恒等式です：
> ```math
> (A + UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}
> ```
> 
> ただし、全ての逆行列が存在すると仮定します。

> **補足：Woodburyの公式と擬似逆行列の関係**
> 
> カルマンゲイン $K = \Sigma_{-}H^T(H\Sigma_{-}H^T + R)^{-1}$ の形式は、
> 重み付き最小二乗問題における擬似逆行列と密接な関係があります。
> 
> 観測方程式 $y = Hx + v$ において、重み行列 $W = R^{-1}$ を用いた
> 重み付き最小二乗解は以下の形をとります：
> ```math
> x = (H^TWH)^{-1}H^TWy
> ```
> 
> これは $H$ の重み付き擬似逆行列を用いた解に相当します。
> カルマンフィルタの更新式は、この最小二乗解に事前分布の情報 $\Sigma_{-}$ を
> 組み込んだ形になっています。
> 
> つまり、カルマンゲインは「事前分布による重み」と「観測ノイズによる重み」の
> 両方を考慮した一種の擬似逆行列と解釈することができます。