LMS Filterの挙動観察中(2)

つまるところ、以下の計算をどうやるか？に尽きる。

\begin{equation*} \mathrm{E}\left[ \delta(\varepsilon(n)) x(n - m) x(n - k) \right] = \lim_{N \to \infty} \frac{1}{N} \sum_{n = 1, \varepsilon(n) = 0}^{N} x(n - m) x(n - k) \end{equation*}

\(\varepsilon(n)\) はi.i.d.（独立に同一の分布）から発生しているので、\(x(n-m), x(n-k)\) には依存しない（予測係数にも依らず）で勝手に揺れると考える。

でもそんな計算は見たことがない。そもそもLADの文脈でこの話は出ているはずで、「Laplace Distribution linear regression」で検索掛けていたら、

実験で試している、勾配ベクトルに分散行列（ヘッセ行列）の逆行列を掛ける行為は、ウィーナーフィルタに等しい。しかしウィーナーフィルタは観測分散行列 \(\ve{XX}^{\mathsf{T}}\) が正則でないと計算できない。そこで、観測分散行列の低ランク近似を行ってその擬似逆行列を使ってフィルタ係数を更新していく手法がある。それを Reduced rank adaptive filters というらしい。

邪念が動いて、パーティクルフィルター（粒子フィルター）でパラメータ決められんか考えてる。でも、LMSは状態空間モデルの範疇に入るのだろうか？（カルマンフィルタの一部だから当てはまったはず）。また、一般の状態空間モデルと違って状態は常に観測できるよな。またパーティクルフィルターもシミュレーションベースなので負荷が高そう。

\(\mathrm{E}\left[ \delta(\varepsilon(n)) x(n - m) x(n - k) \right]\) の解釈

結局 \(\mathrm{E}\left[ \delta(\varepsilon(n)) x(n - m) x(n - k) \right]\) の解釈から逃げている...。もう少し考えていたら、残差 \(\varepsilon\) は入力ベクトル \(x\) と独立であることを思い出した。ここから、次が言える。

\begin{equation*} \mathrm{E} \left[\delta(\varepsilon(n)) x(n - m) x(n - k) \right] = \mathrm{E} \left[ \delta(\varepsilon(n)) \right] \mathrm{E} \left[ x(n - m) x(n - k) \right] \end{equation*}

ここで、 \(\mathrm{E} \left[ \delta(\varepsilon(n)) \right]\) はお察しの通りで、以下の通りに、やはり残差が0となる確率が出てくる。

\begin{align*} \mathrm{E} \left[ \delta(\varepsilon(n)) \right] &= \lim_{N \to \infty} \frac{1}{N} \sum_{n = 1}^{N} \delta(\varepsilon(n)) = \lim_{N \to \infty} \frac{1}{N} \sum_{n = 1, \varepsilon(n) = 0}^{N} 1 \\ &= P(\varepsilon = 0) \end{align*}

よって、

\begin{equation*} \mathrm{E} \left[\delta(\varepsilon(n)) x(n - m) x(n - k) \right] = P(\varepsilon = 0) \mathrm{E} \left[ x(n - m) x(n - k) \right] \end{equation*}

\(P(\varepsilon = 0)\) を考える。まず注意したいのは、連続型確率分布においては一点0をとる確率は0ということ（測度0だから）。近似するしかなく、方針としては、

  1. 残差の閾値を定めて、それ以下の数値を残差0とみなして確率を求める
  2. 離散型確率分布で考える

ラプラス分布の基本的なことをおさらいすると、確率密度関数 \(f(x, \mu, \sigma)\) は、

\begin{equation*} f(x, \mu, \sigma) = \frac{1}{2 \sigma} \exp\left( - \frac{|x - \mu|}{\sigma} \right) \end{equation*}

で、観測 \(x_{1}, ..., x_{N}\) が得られた時の尤度関数 \(L(\mu, \sigma)\) と対数尤度関数は、

\begin{align*} L(\mu, \sigma) &= \prod_{i = 1}^{N} \frac{1}{2 \sigma} \exp\left( - \frac{|x_{i} - \mu|}{\sigma} \right) = \frac{1}{(2 \sigma)^{N}} \prod_{i = 1}^{N} \exp\left( - \frac{|x_{i} - \mu|}{\sigma} \right) \\ \log L(\mu, \sigma) &= -N\log(2\sigma) -\sum_{i = 1}^{N} \frac{|x_{i} - \mu|}{\sigma} = -N\log(2\sigma) - \frac{1}{\sigma}\sum_{i = 1}^{N} |x_{i} - \mu| \end{align*}

\(\mu\) の最尤推定量は標本中央値となる。\(\mu\) が求まったとして、次は \(\sigma\) の最尤推定値を考える。対数尤度関数を \(\sigma\) で偏微分すると、

\begin{equation*} \frac{\partial}{\partial \sigma} \log L(\mu, \sigma) = -2 \frac{N}{2\sigma} + \frac{1}{\sigma^{2}} \sum_{i = 1}^{N} |x_{i} - \mu| = -\frac{N}{\sigma} + \frac{1}{\sigma^{2}} \sum_{i = 1}^{N} |x_{i} - \mu| \end{equation*}

\(\frac{\partial}{\partial \sigma} \log L(\mu, \sigma) = 0\) とおいて \(\sigma\) について解くと、

\begin{equation*} \frac{N}{\sigma} = \frac{1}{\sigma^{2}} \sum_{i = 1}^{N} |x_{i} - \mu| \Rightarrow \sigma = \frac{1}{N} \sum_{i = 1}^{N} |x_{i} - \mu| \end{equation*}

\(\sigma\) の最尤推定値は偏差の絶対値の標本平均となる。次に離散ラプラス分布を考える（ ここ を参考にしている）。離散ラプラス分布は次の確率（質量）関数 \(P\) を持つ:

\begin{align*} P(X = k) &= \frac{f(k, 0, \sigma)}{\sum_{j = -\infty}^{\infty} f(j, 0, \sigma)} = \frac{\exp\left( -\frac{|k|}{\sigma} \right)}{\sum_{j = -\infty}^{\infty} \exp\left( -\frac{|j|}{\sigma} \right)} \\ &= \frac{\exp\left( -\frac{|k|}{\sigma} \right)}{1 + 2 \sum_{j = 1}^{\infty} \exp\left( -\frac{j}{\sigma} \right)} \end{align*}

ここで、

\begin{equation*} \sum_{j = 1}^{\infty} \exp\left( -\frac{j}{\sigma} \right) = \lim_{n \to \infty} \frac{\exp(-1/\sigma)(1 - \exp(-n/\sigma))}{1 - \exp(-1/\sigma)} = \frac{\exp(-1/\sigma)}{1 - \exp(-1/\sigma)} \end{equation*}

よって、

\begin{align*} P(X = k) &= \frac{\exp\left( -\frac{|k|}{\sigma} \right)}{1 + 2 \frac{\exp(-1/\sigma)}{1 - \exp(-1/\sigma)}} = \frac{1 - \exp(-1/\sigma)}{1 + \exp(-1/\sigma)} \exp\left(-\frac{|k|}{\sigma}\right) \\ &= \frac{1 - p}{1 + p} p^{|k|}, \quad p = \exp(-1/\sigma) \end{align*}

これは離散型確率分布であることに注意。

連続版かつ \(\mu=0\) で、 \(|x|\) がある閾値 \(\delta > 0\) 以下となる確率は次のように計算できる:

\begin{align*} P(|x| \leq \delta) &= \int^{\delta}_{-\delta} f(x, \mu, \sigma) dx = \frac{1}{2 \sigma} \int^{\delta}_{-\delta} \exp\left(-\frac{|x|}{\sigma} \right) dx \\ &= \frac{2}{2\sigma} \int^{\delta}_{0} \exp\left(-\frac{x}{\sigma} \right) dx = \frac{1}{\sigma} (-\sigma) \int^{\delta}_{0} \left\{ \exp\left(-\frac{x}{\sigma} \right) \right\}^{\prime} dx \\ &= -\left[ \exp\left(-\frac{x}{\sigma} \right) \right]^{\delta}_{0} = \exp(0) - \exp\left( - \frac{\delta}{\sigma} \right) \\ &= 1 - \exp\left( - \frac{\delta}{\sigma} \right) \end{align*}

この式により分散行列にかける係数を決めることを考えると、次が考察される。

  • \(\delta\) が大きい（分散 \(\sigma\) が小さい）と確率が1に近づき、分散行列はLMSのそれと近くなる。
  • 逆に \(\delta\) が小さい（分散 \(\sigma\) が大きい）と分散行列に小さいスカラーを乗じる。分散行列の逆行列をとると、大きいスカラーを乗じることになり、勾配ベクトルのノルムが大きくなりそう。
  • ノイズレベル（\(\approx\) 分散）が小さいときは勾配が小さくなり極値付近を精密に調べ、大きい場合は勾配が大きくなりダイナミックに探索空間を動き回りそう。

LMSの性能解析に関する文献

色々さまよっているうちに出てきた。

TODO

  • 評価を続ける。評価がまとまったら結果共有に入りたい。
    • LMSの適応動作は、単層パーセプトロンの学習にも該当する。NNの観点からも引き続き論文調査を行うべし。
  • OMPを使う。
  • メッセージパッシング使えない？
    • 何らかの確率モデル化をせよ、というふうに受け取った。
    • AMP, Survay-Propagation（三村さん、樺島さん）がありえる。
    • → AMP, Survay-Propagationについて調査すべし。
  • いろんな論文で自然勾配をどうやって定義しているか要観察。

優先度低

  • パーティクルフィルター使えない？
    • 今日検討した結果、ちょっと今は保留。大量のサンプルが必要そうに見える。計算負荷を気にした結果、優先度を低くした。