雑記: Schur 補行列と Sherman-Morrison-Woodbury の公式とカルマンフィルタとガウス過程回帰の整理

公式の証明はありません(参考文献の1つ目に詳しいのでそちらをご参照ください)。何か問題点がありましたらご指摘いただけますと幸いです。

先にタイトルを回収すると以下です。
  • Schur 補行列とは、「ある行列に対して、その行列の逆行列の左上(または右下)のブロックをとって、その逆行列を取ったもの」。「ある行列に対して、その行列の逆行列の左上(または右下)のブロックをとって、その逆行列を取ったものがほしい」というときにつかう。どういうときにそうなるかというと、「多変量正規分布にしたがう確率ベクトルのうち一部の変数たちを固定して、その条件下で残りの変数たちがしたがう多変量正規分布の分散共分散行列を知りたい」というとき。
  • Sherman-Morrison-Woodbury の公式とは、適当なサイズの4つの行列に対して成り立つ恒等式であって、どんなときに役立つかというと、行列の逐次更新式の形が「時刻 t-1 の行列の逆行列に他の行列を足してから逆行列をとったもの」のようになっているときに、この公式で書き換えた方が計算量的に有利になることがある(逆行列側で更新しろという形だが、諸事情で逆行列側だけでは手順が進められないとき)。別に Sherman-Morrison-Woodbury の公式がないと死ぬわけではない。ただ計算量というのは馬鹿にできないので、やっぱりないと死ぬかもしれない。Schur 補行列との関係は、この公式自体が Schur 補行列による逆行列のブロック表示からきれいに導かれるのと、だから逐次更新したい行列が Schur 補行列の逆行列型になっていたらこの恒等式で書き換えできる。書き換えるとやはり Schur 補行列の形をしている。
  • カルマンフィルタでは、状態ベクトルの分布を、時刻ごとに得られる観測ベクトルにしたがって更新していくが、このうちフィルタ操作(一期先予測分布からフィルタ分布への更新)では、フィルタ後分散共分散行列が Schur 補行列の逆行列型になる。なので、Sherman-Morrison-Woodbury の公式で書き換えて計算量を減らすことが多い(状態ベクトルの次元数よりも、観測ベクトルの次元数の方が小さいときにはこれが有効で、実際カルマンフィルタの適用場面ではそのようなシチュエーションが多い)。フィルタ操作に Schur 補行列をつかうわけではない。
    • ただし、カルマンフィルタにおいて、状態ベクトルと観測ベクトルを連結したベクトルを考えれば、フィルタ操作を「多変量正規分布を一部の変数で条件付けたい」ととらえることはできる。このようにとらえるなら、フィルタ操作に Schur 補行列をつかうことになる。あまりこうとらえないとは思う。このやり方で解くと直接フィルタ後分散共分散行列が Schur 補行列型になる(Sherman-Morrison-Woodbury の公式で書き換えた後の形になる)。
      • 普通(?)にベイズ更新式からカルマンフィルタを解く場合はフィルタ後分散共分散行列は否応なく逆行列型で出てきてしまう=雰囲気としては一期先予測分散共分散行列を  D に関する Schur 補行列ではなく  A に関する Schur 補行列でブロック分割してしまうことになるので、この違いが出る。
  • ガウス過程回帰(PRML下巻 6.4.2 節の)では、目標変数  t_1, \cdots, t_N が観測された下で t_{N+1} のしたがう分布を知りたいが、これは  t_1, \cdots, t_N, t_{N+1} のしたがう事前分布のうち  t_1, \cdots, t_N を固定することで達成される。つまり Schur 補行列をつかう。
    • 逆にこれをカルマンフィルタ的にとらえると、y_{N+1} が状態変数で、 y_1, \cdots, y_N が観測変数であり、その間の共分散がカーネル関数で与えられている。状態の方が1次元で小さいので Sherman-Morrison-Woodbury の公式をつかう必要はないし、というか両辺スカラーになる。
ここからメモ(ガウス過程回帰の話はない)です。
f:id:cookie-box:20180305232608p:plain:w60

以下のサイズの行列 A, \, B, \, C, \, D があるとき、

 A \in \mathbb{R}^{n \times n}, \quad B \in \mathbb{R}^{n \times m}
 C \in \mathbb{R}^{m \times n}, \quad D \in \mathbb{R}^{m \times m}
以下の恒等式逆行列がとられている行列に逆行列が存在するならば、以下の恒等式が成り立ちます。
 (A + BDC)^{-1} = A^{-1}-A^{-1}B(D^{-1}+CA^{-1}B)^{-1}CA^{-1}
これを Sherman-Morrison-Woodbury の公式といいます。

f:id:cookie-box:20180305231302p:plain:w60

え? あ、うん。「そっか、成り立つんだ」って感じなんだけど。なんかうれしいのそれ? あと公式の名前なっが!

f:id:cookie-box:20180305232608p:plain:w60

 A + BDC逆行列が知りたいときに有用です。右辺  A^{-1}-A^{-1}B(D^{-1}+CA^{-1}B)^{-1}CA^{-1} の方を求めればよいということなので。

f:id:cookie-box:20180305231302p:plain:w60

ごめん、いうほど「 A + BDC逆行列が知りたいなー」ってなる? なったとして「右辺の方を求めればいいんだよかったー」ってなる?

f:id:cookie-box:20180305232608p:plain:w60

例えばカルマンフィルタの問題設定(システム・観測は線形、システムノイズ・観測ノイズ・状態の事前分布はガウシアン)の下では、時刻 t におけるフィルタ分布はベイズの定理より以下のようにかけます(各文字の意味は別の記事を参照してください)。

 \begin{split} p(x_{t} \, | \, y_{1:t}) &\propto p(y_t | x_t) p(x_t | y_{1:t-1}) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) \right) \exp \left( -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right)  \\ & \equiv \exp \left( -\frac{1}{2}z \right) \end{split}
ここで z とおいた部分の、x_t に依存する項(x_t のしたがう分布を知りたいので)は以下のようになります。
 \begin{split} z &\propto - {y_t}^{\top} {R_t}^{-1} H_t x_t - {(H_t x_t)}^{\top} {R_t}^{-1} y_t + {(H_t x_t)}^{\top} {R_t}^{-1} H_t x_t \\ & \quad + {x_t}^{\top} {V_{t|t-1}}^{-1} x_t - {x_t}^{\top} {V_{t|t-1}}^{-1} \mu_{t|t-1}  - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1} x_t  \\ &= {x_t}^{\top} ({V_{t|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t) x_t \\ &  \quad - {x_t}^{\top} ({H_t}^{\top} {R_t}^{-1} y_t + {V_{t|t-1}}^{-1} \mu_{t|t-1}) - ({y_t}^{\top} {R_t}^{-1} H_t - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1}) x_t \\&\equiv  {x_t}^{\top} ({V_{t|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t) x_t  - {x_t}^{\top} \alpha - \alpha^{\top} x_t \end{split}
これは x_t の二次形式になので、フィルタ分布もガウシアンであり、2次の項   {x_t}^{\top} ({V_{t|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t) x_t にあらわれる  {V_{t|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t こそが求めたいフィルタ後分散共分散行列  V_{t|t}逆行列に他なりません。上式は  (x_t - \mu_{t|t})^{\top} {V_{t|t}}^{-1}  (x_t - \mu_{t|t}) の形に平方完成できますから(係数比較より  \mu_{t|t} = V_{t|t} \alpha となります。x_t に依存しない定数項の足し引きは x_t の分布の形を変えません)。「 {V_{t|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t逆行列が知りたい(それが求めたい分散共分散行列 V_{t|t} だから)」となったわけです。そしてそれは Sherman-Morrison-Woodbury の公式より  V_{t|t} = V_{t|t-1} - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t|t-1} です。これより、以下の有名なカルマンフィルタのフィルタ分布(の分散共分散行列と平均ベクトル)が導かれます。
 \begin{split} V_{t|t} &= V_{t|t-1} - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t|t-1} \\ \mu_{t|t} &= V_{t|t}\alpha \\ &=  \bigl( V_{t|t-1} - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t|t-1} \bigr) ({H_t}^{\top} {R_t}^{-1} y_t + {V_{t|t-1}}^{-1} \mu_{t|t-1}) \\&= V_{t|t-1} {H_t}^{\top} {R_t}^{-1} y_t + \mu_{t|t-1} \\&\quad - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t|t-1} {H_t}^{\top} {R_t}^{-1} y_t \\&\quad - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t|t-1} {V_{t|t-1}}^{-1} \mu_{t|t-1} \\&= V_{t|t-1} {H_t}^{\top} {R_t}^{-1} y_t + \mu_{t|t-1} \\&\quad - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} (H_t V_{t|t-1} {H_t}^{\top} + R_t ){R_t}^{-1} y_t \\&\quad + V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} R_t {R_t}^{-1} y_t \\&\quad - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t \mu_{t|t-1} \\&= V_{t|t-1} {H_t}^{\top} {R_t}^{-1} y_t + \mu_{t|t-1} - V_{t|t-1} {H_t}^{\top}{R_t}^{-1} y_t \\&\quad + V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} (y_t  - H_t \mu_{t|t-1}) \\&= \mu_{t|t-1} + V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} (y_t  - H_t \mu_{t|t-1})\end{split}
これらは  K_t \equiv  V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} とおくともう少しすっきりかけますね。
 \begin{split} V_{t|t} &= V_{t|t-1} - K_t H_t V_{t|t-1} \\ \mu_{t|t} &= \mu_{t|t-1} + K_t (y_t  - H_t \mu_{t|t-1})\end{split}
何も逆行列補題を適用せずとも分散共分散行列と平均ベクトルを逐次式としてかくことはできます。ただ、逆行列補題を適用しない場合、毎回のフィルタリングで  V_{t|t} = ({V_{t|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t)^{-1} という逆行列を求めなければなりません。この行列は縦幅も横幅も状態変数の次元数ですが、状態の次元が大きいと大きな行列になりえます。他方、逆行列補題を適用すると  V_{t|t} = V_{t|t-1} - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t|t-1} で、このインバースの中に入っている行列のサイズは縦幅も横幅も観測変数の次元数です(カルマンフィルタの状況では A = V_{t|t-1}^{-1}逆行列が既に求まっているので逆行列を求める必要があるのはこのインバースの中身だけになるわけです)。一般にカルマンフィルタの適用場面では「状態変数の次元数 >> 観測変数の次元数」となることが多いので逆行列補題を適用した方が計算コストが小さくなるんです。状態変数の分散共分散を更新したいだけであれば分散共分散行列の逆行列(精度行列)の方を更新していけば逆行列をとる操作を回避することができそうですが、平均ベクトルの更新では平方完成をするために「精度行列の逆行列」が出てきてしまうので。

f:id:cookie-box:20180305231302p:plain:w60

なっが! えっと、要するに、A逆行列がわかっていて、AD のサイズを比べると A の方が結構大きいときには A + BDC逆行列を直接求めるよりも A^{-1}-A^{-1}B(D^{-1}+CA^{-1}B)^{-1}CA^{-1} を求めた方が楽ってこと? それでそういうシチュエーションになる例がカルマンフィルタ? なんか使用場面が限定的じゃない?

f:id:cookie-box:20180305232608p:plain:w60

適用できる問題が特別だからといって公式の有用性が損なわれるということはないと思いますが。もっとも、この「カルマンフィルタの場合」をもっと抽象化してとらえることができます。逆行列補題に出てくる行列  A, \, B, \, C, \, D \left( \begin{array}{cc} A & B \\ C & D \end{array} \right) の形に配置すると大きな行列になることに気付きましたか?

f:id:cookie-box:20180305231302p:plain:w60

え? あー確かに、全体で  (n+m) \times (n+m) の大きさの行列になるな。それが?

f:id:cookie-box:20180305232608p:plain:w60

カルマンフィルタの問題設定で、時刻 t-1 までの観測 y_{1:t-1} を得た下での  \left( \begin{array}{c} x_t \\ y_t \end{array} \right) なるベクトルのしたがう分布を考えてみてください。x_t は多変量正規分布にしたがうとしていいです。であれば、y_tx_t の線形変換にガウシアンノイズを足したものなのでやはり多変量正規分布にしたがいます。この x_ty_t を縦に連結した  \left( \begin{array}{c} x_t \\ y_t = H_t x_t + w_t \end{array} \right) もプロデューサーさんが何かを勘違いしていなければ多変量正規分布にしたがいます。証明は確率ベクトル  \left( \begin{array}{c} x_t \\ y_t \end{array} \right) のモーメント母関数 \displaystyle E \left[ e^{ \left( \begin{array}{cc} t_1^{\top} & t_2^{\top} \end{array} \right) \left( \begin{array}{c} x_t \\ y_t \end{array} \right)} \right] が多変量正規分布にしたがう確率ベクトルのそれと一致することによりました。あ、証明をここにかくのはもう本当に面倒なので省略します。

f:id:cookie-box:20180305231302p:plain:w60

もう数式打つの疲れてるなこれ。それで  \left( \begin{array}{c} x_t \\ y_t \end{array} \right) が多変量正規分布にしたがったら何なの?

f:id:cookie-box:20180305232608p:plain:w60

モーメント母関数の計算から、 \left( \begin{array}{c} x_t \\ y_t \end{array} \right) のしたがう多変量正規分布は以下になることがわかります。

  \displaystyle \exp \left\{ -\frac{1}{2} \left( \begin{array}{c} x_t - \mu_{t|t-1} \\ y_t - H_t \mu_{t|t-1} \end{array} \right)^{\top} \left( \begin{array}{cc} V_{t|t-1} & V_{t|t-1}H_t^{\top} \\ H_t V_{t|t-1} & H_t V_{t|t-1} H_t^{\top} + R_t\end{array} \right)^{-1} \left( \begin{array}{c} x_t - \mu_{t|t-1} \\ y_t - H_t \mu_{t|t-1} \end{array} \right) \right\}
つまり、分散共分散行列  \Sigma が以下になります( A, \, B, \, C, \, DV_{t|t-1} などとサイズが合うように  \Sigma をブロックに分けた行列とします)。
 \Sigma = \left( \begin{array}{cc} A & B \\ C & D \end{array}\right) = \left( \begin{array}{cc} V_{t|t-1} & V_{t|t-1}H_t^{\top} \\ H_t V_{t|t-1} & H_t V_{t|t-1} H_t^{\top} + R_t \end{array} \right)
さて、ここで時刻 t の観測値 y_t が得られたとしましょう。y_t を観測したもとでの x_t の分布はどうなるでしょうか? つまり、上の式から「 x_t の分散共分散行列」だけを切り出したいということです。もちろん V_{t|t-1} は不正解です。分散共分散行列は逆行列の形で挟まっているので、分散共分散行列から左上の n \times n のサイズのブロックを切り出せばいいということにはなりません(精度行列であればそれでいいんですが)。
  • 元々 x_t の分散共分散行列が V_{t|t-1} であったはずなのになぜ V_{t|t-1} にならないのか疑問に思う方もいるかもしれません。 V_{t|t-1}y_t を知らない下での x_t の分散共分散行列です。y_t を知ってしまうと分散共分散行列は変わってきます。斜めに傾いた(独立でない)2変量正規分布を思い浮かべて、それの横軸を x_t、縦軸を y_t と考えてみてください(この記事のグラフの真ん中のようなイメージです)。もし y_t について何も知らないなら、x_t の分布はこの2変量正規分布y_t 方向に積分することになります(先ほどの記事の Histgram of x)。しかし、ある y_t を手に入れてしまった場合は、その値に引いた水平線で2変量正規分布を切った断面になります。グラフの上の方の値なのか、下の方の値なのかで x_t の広がりは異なりますよね。y_t を観測したもとでの x_t の分布を得るには、「2変量正規分布からの切り出し」が必要なんです。

f:id:cookie-box:20180305231302p:plain:w60

確かに逆行列だから単純に左上のブロックを取るわけにはいかないな。って、あれ? y_t を観測した下での x_t の分布って、それフィルタ後分布じゃん。じゃあフィルタ後分散共分散行列はさっき上で求めたこれだろ?

V_{t|t} = ({V_{t|t-1}}^{-1} + {H_t}^{\top} {R_t}^{-1} H_t)^{-1} = V_{t|t-1} - V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} H_t V_{t|t-1}

f:id:cookie-box:20180305232608p:plain:w60

さすがにばれてしまいましたか。その  V_{t|t} \Sigma の部分行列  D に関するシューア補行列(Schur complement matrix)といいます。 A, \, B, \, C, \, D でかくと  A - B D^{-1} C ですね。

f:id:cookie-box:20180305231302p:plain:w60

へ? シュ、シューア補行列? 何それ? 確かに V_{t|t} の式と A, \, B, \, C, \, D を見比べると  A - B D^{-1} C になってるけど。

f:id:cookie-box:20180305232608p:plain:w60

シューア補行列とは「ある行列に対して、その行列の逆行列の左上のブロックをとって、その逆行列を取ったもの」ですよ(あるいは右下のブロックでもよく、その場合は A に関するシューア補行列  D - C A^{-1} B となる)。上の多変量正規分布からの切り出しは、「一部の変数を固定した残りの変数のみの分散共分散行列を知りたいので、もとの分散共分散行列の逆行列をとって、その左上のブロックをとって、それの逆行列をとって分散共分散行列に戻したい」ということなので、これはシューア補行列そのものです。ちなみにシューア補行列を用いた逆行列のブロック表示からの Sherman-Morrison-Woodbury の公式の証明が参考文献にあるので参照してください。

中途半端ですがつづかない

パターン認識と機械学習 下(第6章:その4)

以下の本を読みます。上巻を読むこともあります。何か問題点がありましたらご指摘いただけますと幸いです。

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

前回:その3 /次回:まだ
f:id:cookie-box:20180305232608p:plain:w60

結局等価カーネルというのは、線形基底関数モデルの重みパラメータをベイズ的に解いたときの、重みパラメータの事後分布の平均ベクトルによる、ある対象の点における予測値が、訓練データ中の各被説明変数の線形結合になっているという見方をすることができて、そのような見方をしたときに各被説明変数にかかる重みが、対象の点とその訓練データ点の「等価カーネル」なんですね。…この等価カーネルと 6.1 節のカーネル関数は異なるものですね。等価カーネル k(x, x')=\beta \phi(x)^{\top} S_N \phi(x') は特徴の内積ではなく間に重みパラメータの事後分散共分散行列 S_N を挟んでいますし。しかし  \psi(x) \equiv \beta^{1/2} S_N^{1/2} \phi(x) と新しい特徴を定義すれば等価カーネルはこの新しい特徴の内積にはなっていますね。

f:id:cookie-box:20180305231302p:plain:w60

あれ? ジュン今ドイツじゃないの?

f:id:cookie-box:20180305232608p:plain:w60

それは…これはテレビ電話です。

f:id:cookie-box:20180305231302p:plain:w60

そこまでして勉強会してる設定なの俺たち!?

f:id:cookie-box:20180305232608p:plain:w60

6.3.1 節の Nadaraya-Watson モデル(カーネル回帰)のシチュエーションは、点 x における値が t である同時分布 p(x,t) を、各データ点から x-x_n 離れた点における値が t-t_n である同時分布 f(x-x_n,t-t_n) の混合分布で表そうとしていますね。これも Parzen 推定法なんですね。Parzen 推定法は密度推定の方法であって回帰というイメージではありませんでしたが…しかし、被説明変数 t を説明変数 x 側にくっつければ回帰問題を密度推定問題ととらえることもできそうですね。

f:id:cookie-box:20180305231302p:plain:w60

f(x,t)t についての平均は零であるとする」っていうのは、ここでは x をある点に固定したとき tt_n を中心にした1次元の正規分布(例えば)みたいな分布になっているのを想定してるってことだよな。…じゃあ結局1つ1つの f(x-x_n, t-t_n) については t=t_n のところでの確率密度が一番大きくて、ある x における y(x) がどうなっているかは f(x-x_n, t-t_n) の形が x 方向にどうなっているかによるけど、x_n から離れるほど小さくなるような密度だったら、結局その x に一番近い x_n に対応する t_n の影響を一番色濃く受けることになって…それって等価カーネルと似てるな。それに (6.45) 式も (3.61) 式と同じ形だ。でも今は線形基底関数モデルをベイズ的に解いたわけじゃないのになんで答えの形が同じになるの?

f:id:cookie-box:20180305232608p:plain:w60

…6.3.1 節では最初から t 方向にも x 方向にもばらつきをもつ分布を混合しようとして自然に式 (6.45) の形が出てきたようにみえます。どちらかというと 3.3 節の結果が式 (3.61) の形になる方が直感的ではない気がしますが…。ただ式の形はともかく、やりたいこととしては、3.3 節で重みパラメータを分布として考えたのは結局 t 方向のばらつきを考慮した回帰をしたい、ということになると思います。6.3.1 節で密度推定をしているのも同じだと思います。結局知りたいのは各点での p(t|x) で、これを解くのに重みパラメータ w を主役にするか各訓練データ点が張る密度 f(x-x_n, t-t_n) を主役にするかというだけの違いな気もします。ただこちらの 6.3.1 節の枠組みでは分布の形 f を等方的なガウシアンではなく柔軟に選択できると。いきなり分布の形を顕わにもつ (6.42) 式を仮定として出発しているのでそれはそうですね。…また 6.3 節冒頭に戻ると、(3.61) 式、(6.45) 式に加えてもう一つ同様の形の (6.40) 式が登場していて、これは (6.39) 式から出発して導かれたものですね。まとめると、以下のケースではすべて結果の回帰モデルが「各訓練データ点における被説明変数の値を、各訓練データ点との距離のRBFで混合したもの」になるんですね。

  • 線形基底関数モデルを、重みパラメータの分布を事前分布から事後分布にベイズ的に更新するやり方で解いたときの、事後分布の平均値による解(3.6.1節)。
  • 入力変数側にノイズが含まれる回帰モデルを変分法で解いたときの解(6.3節冒頭)。
  • 同時分布を Parzen 推定したときの、各点における被説明変数の期待値(6.3.1節)。
RBF で混合する解になったということは「x_n に近い点での値は t_n に近い値であるべき」という結論になったということですが、上の3つのケースのそれぞれ「重みパラメータの分布」「入力変数のノイズ」「説明変数と被説明変数の同時分布への仮定」がその由来になっていると思います。逆に、訓練データとして説明変数と被説明変数の組を与えられても、それだけでRBFを混合した解にはなりませんね。線形基底関数モデルを最小2乗法で解くことができる場合なら (6.9) 式の形にはなりますが、これはそもそも t_nカーネル関数の重みで混合した形をしていませんし。

f:id:cookie-box:20180305231302p:plain:w60

ふーん…まあそれで13ページの下の方にさ、「より一般的には、同時分布 p(t,x) を混合ガウス分布で表すことも可能であり」ってあるけど、そもそも (6.42) 式って f がガウシアンなら混合ガウス分布じゃん。何言ってんのこれ? あと9章ぱらぱら見ても「9章で紹介するテクニック (Ghahramani and Jordan, 1994)」がどれのことかわかんなかった(見落としたかもだけど)…。

f:id:cookie-box:20180305232608p:plain:w60

一瞬わかりにくいですよねここ…次のページまで読むと、ここで「混合ガウス分布で表す」といっているのは各訓練データの周りに密度を張ることではなくて、訓練データをEMアルゴリズムである程度混合ガウス分布クラスタリングして、クラスタ単位で密度を張るという意味ですね。「9章で紹介するテクニック」のくだりは「9章で紹介するEMアルゴリズムを用いてクラスタリングしてからカーネル回帰する」ということだと思うんですが…Ghahramani and Jordan, 1994 自体は 9.3 節に関係していそうですが…。

f:id:cookie-box:20180305231302p:plain:w60

まあいっか。結局 6.3 節の話って、この節に紹介されているようなケースだと RBF 型の等価カーネルの重みで訓練データ中の被説明変数の値を混合する形の答えになる、ってだけだよな…。ただここにきて初めて、カーネル関数の形に理由がついたってことなのかな。6.2 節までの時点ではどのカーネル関数を選ぶといいのかの理由って「このカーネル関数なら上手く学習できるかも?」くらいしかなかった気がするし…。といっても等価カーネルって 6.1 節で最初に紹介されたカーネル関数とはなんか違う(素直な特徴の内積じゃない)からちょっと引っかかるけど…。それで 6.4 節の「ガウス過程」もこれたぶんカーネル関数への理由づけだよな。6.4.1 節は「線形回帰再訪」? なんかさっきちらって見た9章にも「混合ガウス分布再訪」ってあったし、再訪しまくってんな…。

f:id:cookie-box:20180305232608p:plain:w60

いいじゃないですか何度訪ねても。そうでなくても訪ねなければ忘れてしまいますし。

f:id:cookie-box:20180305231302p:plain:w60

いや悪くないけど、「カーネル関数はあなたがすでに知っている問題にも現れる、こんなにも正当なものなんですよ」って話より、早くカーネル関数特有のすごい話みたいなのを聞きたいんだけど…。それで 6.4.1 節の問題設定は 3.3 節と同じだな。でも  y(x_1), \cdots, y(x_N) の同時分布を考えようとしてる…え、これらを並べたベクトルってガウス分布にしたがうの?

f:id:cookie-box:20180305232608p:plain:w60

(6.51) 式ってこうですよね。それぞれ独立に多変量正規分布にしたがうベクトルを足し合わせた形をしていますから、y も多変量正規分布にしたがうんじゃないですか。

 \displaystyle y = \left( \begin{array}{c} w_1 \phi_1(x_1) + w_2 \phi_2(x_1) + \cdots + w_M \phi_M(x_1) \\ w_1 \phi_1(x_2) + w_2 \phi_2(x_2) + \cdots + w_M \phi_M(x_2)  \\ \vdots \\ w_1 \phi_1(x_N) + w_2 \phi_2(x_N) + \cdots + w_M \phi_M(x_N)  \end{array} \right)
ほら、以下などに互いに独立な多変量正規分布にしたがう確率ベクトルの和は多変量正規分布にしたがうとかいてありますよ。

f:id:cookie-box:20180305231302p:plain:w60

なるほど。それで…え、E[y]=0 になっちゃうの? 回帰モデルで「どこの点でも答えの期待値はゼロです」っておかしくない? …って、w がまだ事前分布 p(w)=N(w|0, \alpha^{-1}I) だからか。重みの平均がゼロならそりゃ回帰される値はゼロになるよな。それで、これがガウス過程? ガウス過程って確率過程の一種だよな? 確率過程ってなんか時間変化していくような確率変数のことじゃないの? ウィキペディアにもそうかいてあるし。いまどこに時間変化要素あったの??

f:id:cookie-box:20180305232608p:plain:w60

「時間変化していく」というのは確率過程の定義に含まれていません。だいたい確率過程の定義は数学の話ですし、数学の世界に時間という概念などないでしょう? 確率変数の添え字は添え字であって、時間を意味している必要はありません。もちろん実際の応用では時間変化する対象を記述しようとして、時間の経過を添え字に背負わせることが圧倒的に多いわけですが…。ともかくここでは「 x_1, \cdots, x_N をどのように選んでも  y(x_1), \cdots, y(x_N) の同時分布が多変量正規分布にしたがう」ことをガウス過程といっているにすぎません。なお、多変量正規分布の平均ベクトルと分散共分散行列は常に同じである必要はありません。ここでは平均ベクトルは常にゼロベクトルになりますが、分散共分散行列は (6.53) ですから  x_1, \cdots, x_N によって変わりますね。そしてこの分散共分散行列は、各成分が訓練データ点 x_nx_mカーネル関数になっています(正確にはそれに重みパラメータの事前分布の分散をかけたものですが)。「これらの訓練データ点上での値がどうなるか特徴で回帰したいが、まだ何も値を観測していないというとき、値はこれくらいの不確かさ(分散共分散行列)で広がっている」というのを出すと、カーネル関数が現れたということなんですね。他方、特徴は直接現れていません。

f:id:cookie-box:20180305231302p:plain:w60

ふーん…まあガウスカーネルだったら近い点どうしは大きい値に(共分散大)、遠い点どうしは小さい値に(共分散小)なるわけで、じゃあこの分散共分散行列は、「まだ全然どんな値になるかわからなくて不確かだけど、近い点どうしだったら不確かさも同じ方向に広がっているだろ」「遠い点どうしだったら全然関係ないだろ」って感じなのかな…だから何って感じだけど…。で、この図 6.4 は何? 本文には「ガウス過程からのサンプル」ってあるけど、説明足りなさすぎない? 本文でいうどの式をプロットしてるの?? グラフにはちゃんと縦軸と横軸は何かかきましょうって習わなかったの??

f:id:cookie-box:20180305232608p:plain:w60

いきなり図だけだとわかりにくいですね…あ、描けました。

f:id:cookie-box:20181205224243p:plain:w320f:id:cookie-box:20181205232147p:plain:w320

f:id:cookie-box:20180305231302p:plain:w60

うおおいどうやって描いたんだよそれ!?

f:id:cookie-box:20180305232608p:plain:w60

すみません、ぶっちゃけ以下にあったスクリプトを丸々コピーして実行しました…刊行に先立って意見を募集されているということなので、何か意見をお送りできるように善処します…。

f:id:cookie-box:20180305231302p:plain:w60

何その本!? めっちゃわかりやすそう!!

f:id:cookie-box:20180305232608p:plain:w60

ともかくこれでわかりましたが、図 6.4 を描く方法は、例えば x_1 = -1, \cdots, x_{100}=1 と-1~1の点を均等に100点でもとって、それらの点から (6.54) 式にしたがって 100×100 次元のグラム行列を生成、それを分散共分散行列とした多変量正規分布から1点をサンプリングしたもの y_1, \cdots, y_{100} をプロットしているだけです(1点といっても100変量正規分布からサンプリングした点なので100次元ですね)。100変量正規分布にしたがう乱数ベクトルを1つ生成してプロットしたものが赤色のグラフで、同じ分布からもう1回乱数ベクトルを生成してプロットしたものが青色のグラフで、…といった具合です。もし100変量正規分布の異なる成分間が独立であったらホワイトノイズのようなグラフになってしまいますが(分散の大きさが各点で同じなら実際ホワイトノイズですが)、そうではなく各成分がカーネル関数であるような分散共分散構造をもっているので、近い点どうしは近い値に=滑らかになるわけです。指数カーネルの場合はちょっとギザギザしますけど。

f:id:cookie-box:20180305231302p:plain:w60

ああこれ折れ線グラフだけど横軸方向に点がめっちゃいっぱいあって滑らかに見える感じなのか…てっきり何か関数の重ね合わせ的なのかと…。

f:id:cookie-box:20180305232608p:plain:w60

本当はいくらでも密になるんだと思うんですが…。それはさておき、14ページの 6.4 節の冒頭にあった、「関数 y(x,w) に関する事前分布」というのがこの色々な色のグラフたちの上の分布なんですね。100変量正規分布からもっとたくさんの点をサンプリングすれば「関数の分布」らしさが出るでしょうか。

f:id:cookie-box:20181205233338p:plain:w320

f:id:cookie-box:20180305231302p:plain:w60

ごちゃごちゃしすぎ!!

(次回があれば)つづく

パターン認識と機械学習 下(第6章:その3)

以下の本を読みます。上巻を読むこともあります。何か問題点がありましたらご指摘いただけますと幸いです。

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

前回:その2 /次回:まだ
f:id:cookie-box:20180305232608p:plain:w60

特徴の内積として定義されたカーネル関数でしたが、現実には特徴を意識せずにいきなりカーネル関数を構成したい。でも勝手な関数にするわけにはいかないので、前回は、「k(x,x') が有効なカーネル関数である」、つまり「任意の訓練データの組合せに対して k(x_n, x_m)k(x_n, x_m) = \phi(x_n)^{\top} \phi(x_m) の形にかけるような、入力空間から特徴空間への写像 \phi(x) が存在する」ことの必要十分条件は、「nm 列の成分が k(x_n, x_m) であるような行列(グラム行列)が半正定値である」ことを確かめました。といっても、適当な関数を選んで「カーネル関数たる資格があるかどうか」を調べるのも面倒なので、既知のカーネル関数から6ページの公式 (6.13)~(6.22) を利用して新しいカーネル関数を構成するのが便利なようです。前回はここの証明はとばしたので証明の概略を考えてみます。

  • (6.13) は、新しい特徴を \phi'(x)=\sqrt{c}\phi_1(x) とすればよいことから明らかですね。
  • (6.14) も、新しい特徴を \phi'(x)=f(x)\phi_1(x) とすればよいです。
  • (6.17) ですが、半正定値行列は「全ての固有値が非負」の他に「2次形式が非負」という定義づけをすることもできます(証明略)。ので、k_1(x_n, x_m)+k_2(x_n, x_m) を成分とする行列も半正定値です。
  • (6.18) は、新しい特徴が \phi'(x)=\phi_1(x)\phi_2(x) になるはずです。
  • (6.15) は、(6.13)、(6.17)、(6.18) を繰り返せば導出できるはずです。
  • (6.16) は、(6.15) に  \exp (x)マクローリン展開を適用すればよいと思います。
  • (6.19) は、\phi'(x)=\phi_3(\phi(x)) ですかね。
  • (6.20) は、A は半正定値なので、A固有値の非負の平方根を対角成分に並べた対角行列を \hat{\Lambda} と、適当な直交行列 P によって A=P\hat{\Lambda}(P\hat{\Lambda})^{\top} とかけるので、k(x,x')=((P\hat{\Lambda})^{\top}x)^{\top}(P\hat{\Lambda})^{\top}x' ですね。
  • (6.21) (6.22) は \phi(x_a)=x_b と考えればここまでの公式から導かれる気がします。
いずれの公式でも、新しいカーネル関数もまた特徴の内積になっていることがわかりますね。

f:id:cookie-box:20180305231302p:plain:w60

最後雑だなー。それで、前回の続きには、カーネル法の利点として入力が数値ベクトルでなくてもカーネル関数が定義できればよい、というのが挙げられているな。確かに、テキスト分類なんかの場合、各単語の特徴ベクトルを考えるより、各単語間のカーネル関数を考える方がやりやすそうかも…。あと、カーネル関数を確率的生成モデルにしたり、隠れマルコフモデルにすることもできるのか。

f:id:cookie-box:20180305232608p:plain:w60

(6.29) 式は xx' が同じ i 番目の分布に所属しているかどうかといった感じがしますね。重み p(i) が乗じられますが。(6.30) 式も、配列 X と配列 X' が同じ隠れ状態 Z から生成されたかどうかをみているようです。なるほどこれらのモデルは xx'(配列 X と配列 X')の類似度を直接測るよりも訓練データへの依存性が低いモデルになりそうですが、データの性質をよく表現する生成モデルや隠れマルコフモデルが用意できるかというのはまた別の問題ですよね。

f:id:cookie-box:20180305231302p:plain:w60

それでさらにフィッシャーカーネルというのもある、と。…なんかこの 6.2 節って要するに「あんなカーネル関数もこんなカーネル関数もつかえるよ」って話で、確かにそんな関数をパターン認識につかえたら便利なのかもって感じはするけど、具体例が伴わないからなんかふわふわしてるんだよな…。それでフィッシャースコアとフィッシャー情報行列って何?

f:id:cookie-box:20180305232608p:plain:w60

フィッシャースコアは式 (6.32) の通り対数密度のパラメータベクトルに関するナブラですよ。そしてフィッシャー情報行列はその分散共分散行列です。あるパラメータベクトルにおける。式 (6.34) ですね。

f:id:cookie-box:20180305231302p:plain:w60

え、式 (6.34) ってフィッシャースコアの分散共分散行列になってる?

f:id:cookie-box:20180305232608p:plain:w60

フィッシャースコアの期待値はゼロベクトルになるので。

f:id:cookie-box:20180305231302p:plain:w60

ふーん…それにしても、カーネル関数って、一方の点から他方の点につくる密度みたいな感じだったけど、これだとそれぞれの点での対数密度をパラメータで偏微分してから掛け合わせてるから密度ではないよな。パラメータが1次元だと考えると、(6.33) 式ってこうなる?

 \displaystyle k(x, x';\theta) = \frac{g(\theta, x) \cdot g(\theta, x')}{E_x \left[ g(\theta, x)^2 \right]} = \frac{ \displaystyle \frac{\partial}{\partial \theta} \ln p(x|\theta) \cdot \frac{\partial}{\partial  \theta}  \ln p(x'|\theta) }{E_x \left[ \left( \displaystyle \frac{\partial}{\partial \theta} \ln p(x|\theta) \right)^2 \right] }
分母は x,\,x' によらないから差し当たり無視すると、このカーネル関数ってこういうことだよな。
  • パラメータに関する偏微分がゼロになるような点(例えば正規分布の平均をパラメータと考えているとしたら、その平均にあたる点)については、どこの点とのカーネル関数を計算してもゼロになる。
  • ある点と同じ点とのカーネル関数は、その点でのパラメータに関する偏微分の絶対値が大きいほど大きい(ゼロ以上)。ガウスカーネルのようにどの地点でも一定値ではない。
  • 例えば正規分布の平均をパラメータと考えているとしたら、このカーネル関数が最小になる地点の組み合わせは、正規分布の左側の傾きの絶対値が最大になる点と正規分布の右側の傾きの絶対値が最大になる点で、最大になる地点の組み合わせはそれらの点の同じ点どうし。
…これが「類似度」になるの? 正規分布の中心の点はどこの点とも類似度がゼロっておかしくない?

f:id:cookie-box:20180305232608p:plain:w60

「生成モデルのパラメータを少し動かしたときに密度がどれだけ急に変わるか」を主役にした類似度ですね…確かに密度のパラメータに関する偏微分の絶対値が大きい点ではパラメータを少し動かすと全体の尤度への影響が大きいですが…。ある点とある点のカーネル関数が大きなプラスだったら、その点どうしは「自分たちの尤度を大きくするためにパラメータを同じ方向に動かすことに強く同意し合っている」、大きなマイナスだったら「パラメータを動かす方向について意見が全然合わない」、ゼロだったら「どちらでもない」といった感じですね。と考えればある意味「点どうし気が合うか」を測っているような気はしますが…。フィッシャーカーネルは文書や画像の認識で威力を発揮するようですが、テキストの9ページからはその心を読み取りづらいですね。そこに言及がある情報幾何を知っていればもっと理解できるのでしょうが…。

f:id:cookie-box:20180305231302p:plain:w60

でもなんかフィッシャーカーネルって聞くとケンタッキーフライドフィッシュ思い出さない? そんなの一時期売ってなかった?

f:id:cookie-box:20180305232608p:plain:w60

え? すいません僕ファストフードに詳しくないんで…。

f:id:cookie-box:20180305231302p:plain:w60

で最後のカーネル関数の例がシグモイドカーネルか。え、必ずしも半正定値にはならないって、「カーネル関数は特徴の内積」っていう大前提ぶん投げてるじゃん…。それで、カーネル法ニューラルネットワークは実は似てるんですよって話が出てきてるな。というかここにくるまでカーネル法とは何かが定義されてないんだけど、6.1節の双対表現を解くってことでいいのか?

f:id:cookie-box:20180305232608p:plain:w60

ニューラルネットワークとの関連の話は後の節でまた出てくるようなので、6.3節に進みましょうか。

f:id:cookie-box:20180305231302p:plain:w60

6.3節は「RBFネットワーク」か。基底関数、カーネル法の文脈で言い換えると特徴の1つ1つの成分、としてどのような形のものをとるとよいかって話で、RBF(動径基底関数;radial basis function)っていう  \phi_j(x) = h(|| x - \mu_j || ) の形のものが多いってかいてあるな。

f:id:cookie-box:20180305232608p:plain:w60

ただこの6.3節の導入部ではまずはカーネル法は意識されていませんね。単にモデルを構成するパーツとしての RBF です。なぜ RBF が導入されたのかとありますが、例えば…  x=10 において  y = 1 x=20 において  y = 2 x=30 において  y = 1 という値が観測されたときに、では x \in \mathbb{R} ではどうなっているんだろうと考えたいわけですが、x=10, \, 20, \, 30 を中心にガウス関数でも置いとけばいいだろという精神ですよね(頂点での値はちゃんとその地点での観測値になるようにして)。

f:id:cookie-box:20180305231302p:plain:w60

いやいや雑すぎるだろ! そんな山3個置いただけのがあらゆる x に対するモデルになってる気がしないんだけど! いやそもそも3点しかないから全体の姿はなんかもう全然わかんないけど!

f:id:cookie-box:20180305232608p:plain:w60

ええ、通常このようにはしないというような記述がありますね。

f:id:cookie-box:20180305231302p:plain:w60

次の段落の正則化理論と、グリーン関数って何?

f:id:cookie-box:20180305232608p:plain:w60

正則化理論はそういう機械学習の理論なのでしょうね。よくわかりませんが。グリーン関数というのはその…電磁気学ポアソン方程式のような微分方程式を解くときに…求めたい静電ポテンシャルが電荷分布と何らかの関数の畳み込み積分なのではと考えて…それをもとの微分方程式に代入して電荷分布にデルタ関数を代入することで何らかの関数がわかって…それを電荷分布と畳み込み積分すればめでたく求めたかった静電ポテンシャルになって…その何らかの関数がグリーン関数とよばれて点電荷がつくるポテンシャルに相当するというような…(以下などを参照してください)。

f:id:cookie-box:20180305231302p:plain:w60

なんでそんな歯切れ悪いの?

f:id:cookie-box:20180305232608p:plain:w60

プロデューサーさんは大学で散々グリーン関数を扱ったはずなのに何も説明できないんですよね。

f:id:cookie-box:20180305231302p:plain:w60

残念すぎるな。

f:id:cookie-box:20180305232608p:plain:w60

しかしグリーン関数はまるでカーネル関数ですね、というより、カーネル関数がグリーン関数ですね。訓練データという点電荷たちが空間の各点につくる静電ポテンシャルを求めるための。

f:id:cookie-box:20180305231302p:plain:w60

でも訓練データは点電荷じゃないじゃん。

f:id:cookie-box:20180305232608p:plain:w60

はい。

f:id:cookie-box:20180305231302p:plain:w60

式 (6.39) は確率的なノイズがある場合の回帰モデルの最適化ってことだよな。Nadaraya-Watson モデルっていうのか。それでなんか 6.3 節導入部の最後に、カーネル密度推定を用いた回帰では各データ点に基底関数が関連付けられているから計算コストがかかるけど、基底関数(ここでは各データ点を中心にした RBF を全部特徴の成分にするつもりだったんだよな、それは計算コストがやばそうだな)の数を削減する方法もあるってかいてあるな。削減方法の例には直交最小2乗法(誤差の2乗和を最も削減する点から順に選んでいく)や K-means クラスタリングして重心を採用するとかがあるのか。…で 6.3.1 節がいきなり「3.3.3 節では(中略)等価カーネルによって与えられた」からはじまってるんだけど、俺等価カーネルって何か知らないんだけど、というか他にも 6.3 節ってちょいちょい3章を回顧してない? 上巻の3章ってどんな内容だったんだろ?

f:id:cookie-box:20180305232608p:plain:w60

3章は「線形回帰モデル」というタイトルですね。3.1節ではまず特徴の線形和でかけるモデルが導入されていて…なんか図 3.1 に図 6.1 と全く同じ絵が…。そのパラメータ推定ですが、加法的な正規ノイズの場合は最尤推定と最小2乗推定が等価だと…それはそうですよね。図3.2 は面白いですね。特徴が2次元の場合の絵にみえます。\varphi_1 は全データの特徴の1つめの成分を並べたベクトルで、\varphi_2 は2つめの成分を並べたベクトルで、これらの線形和で各データに対応する被説明変数を並べたベクトル t をつくりたいわけですが、線形和は面 \mathcal{S} 内しか動けないので t\mathcal{S} への正射影を最適解とするしかないという状況ですね。あとは正則化の話もあって…3.2節にはバイアス-バリアンス分解による推定のやり方のよさの議論なんかもあるんですね。節の最後で「ここからはベイズ的にいく」って感じで決別されてる感じですけど。3.3節はベイズ線形回帰というタイトルで、パラメータ w の分布を更新しようという話になってますね。このとき事後分布の平均は (3.53) 式  m_N = \beta S_N \Phi^{\top} t となり、これによる予測は(3.61) 式  y(x, m_N) = \sum_n \beta \phi(x)^{\top} S_N \phi(x_n) t_n \equiv \sum_n k(x,x_n) t_n となってこの k(x,x') を等価カーネルとよぶと。これは、未知の点 x での値の予測において、訓練データ中の被説明変数をどれだけの割合ずつ混ぜるべきかと解釈できるんですね(\sum_n k(x,x_n)=1)。図 3.10 のように、ガウス基底関数を選んだときのカーネル関数を x' についてプロットすると、相手の点 x の周りに局在するというのは僕が前回プロットで示したのと同じですね(もっとも、前回プロットしたカーネル関数は単なるガウス基底関数の内積であって、事後分散共分散行列 S_N を挟んでいませんでしたが)。つまり、未知の点に対する値を予測したいときは、その点に近い説明変数に対する被説明変数をより多く混ぜよ、という直感的な結果になりますね。

(次回があれば)つづく

パターン認識と機械学習 下(第6章:その2)

以下の本を読みます。上巻を読むこともあります。何か問題点がありましたらご指摘いただけますと幸いです。

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

前回:その1 /次回:まだ
f:id:cookie-box:20180305232608p:plain:w60

前回の 6.1 節では、「各訓練データを特徴空間に埋め込んで線形回帰する」という問題は、「入力空間内の各点に各訓練データが張る密度を線形回帰する」という問題に読み替えることができることを確認しましたね。この場合は D 次元(D は特徴空間の次元数)の係数ベクトル w を求めるのではなく N 次元(N は訓練データの個数)の係数ベクトル a を求めることになります。位置 x に位置 x' から張られる密度がカーネル関数 k(x,x') = \phi(x)^{\top} \phi(x') でした(もっとも、カーネル関数は対称なので、位置 x から位置 x' に張る密度といっても構いませんが)。

f:id:cookie-box:20180305231302p:plain:w60

なあジュン、その「密度」って言い方引っかかるんだけど、特徴空間への写像  \phi(\cdot) って別に内積を取ったら確率密度っぽくなるように選んでるわけじゃないだろ? 適当に選んだ  \phi(\cdot)k(x,x_0) = \phi(x)^{\top} \phi(x_0) を入力空間全体で積分したら発散しちゃったりしない? ていうかだいたいいつも発散するんじゃない?

f:id:cookie-box:20180305232608p:plain:w60

発散したらダメなんですか?

f:id:cookie-box:20180305231302p:plain:w60

え?

f:id:cookie-box:20180305232608p:plain:w60

僕は「確率」密度とはいってませんよ。確かに 2.5.1節「カーネル密度推定法」では目的が確率密度推定でした。でも今の目的は回帰や分類です。それぞれのデータ点が周囲に張り巡らす「場」は空間全体で積分したら発散するようなものであっても構いませんし、それを各点で足し上げる重みも負であっても構いません。2.5.1 節では推定する確率密度の構成要素として明示的に確率密度の条件を満たすカーネル関数 k(x - x_n) を選びましたが、本節におけるカーネル関数 k(x, x') はハヤトの言う通りこんなカーネル関数がいいと選んだものではなく「特徴の内積」として結果的に出てきたものです。といっても、実際にカーネル法で解こうとするような場合はやはり特徴ではなくカーネル関数の方を決めるようですね。それを 6.2 節でみていきましょう。

f:id:cookie-box:20180305231302p:plain:w60

こんがらがってきたな…いま特徴を決めて線形回帰しようとしてたら、「特徴と各訓練データの特徴との内積カーネル関数)」を線形回帰してもいいじゃんって言われたのが 6.1 節で、さらに特徴じゃなくてカーネル関数を決めろって言われてる状況? どんだけカーネル関数ごり押しだよ…まあ先に進むか。6.2 節には、「カーネル関数を構成するには、特徴から構成するアプローチと、カーネル関数を直接定義するアプローチとがある」って書いてあるな。え、だから特徴の内積がほしいカーネル関数になるように狙って特徴つくるの難しくない? 基底関数って? あと5ページの上の図 6.1 何これ?

f:id:cookie-box:20180305232608p:plain:w60

ちょっと 6.2 節の冒頭趣旨がわかりづらいですね…先に図 6.1 ですが、要するにこういうことですね(※ スクリプトは記事の最下部です)。パラメータは目算で合わせました。

f:id:cookie-box:20181114200048p:plain:w660

f:id:cookie-box:20180305231302p:plain:w60

あ、図 6.1 と同じ絵…だけど、何も要するにになってないし。あと、図 6.1 では上段と下段しかないけどその中段は何? まあ上段と下段も何なのかわかんないけど。

f:id:cookie-box:20180305232608p:plain:w60

真ん中の列が一番わかりやすいんじゃないかと思うんです。手書き数字のグレースケール画像を分類するような場合を想像してください。もっともその場合は入力空間は1次元でなく画像のピクセル数の次元があるわけですが、上のグラフは1次元ですし無理やり1次元だと考えてください。同じ正解ラベルをもつ画像は入力空間での距離が近いと考えるのが自然でしょう。未知データに対する予測も、その未知データに近い訓練データがより大きな影響を及ぼすようなものであってほしいです。なので、個々の訓練データは入力空間に、そのデータの位置で最も濃く、データから離れるほど薄くなるようなカーネル関数の煙幕を張りたいんです。中央列の最下段のように。

f:id:cookie-box:20180305231302p:plain:w60

あー、この図の最下段は「x 軸上の各点とバツ印の点(にある訓練データ)とのカーネル関数」なのか。

f:id:cookie-box:20180305232608p:plain:w60

しかし問題はこのようなカーネル関数を得るのにどのような特徴をとればいいのかということですが…天下りになってしまいますが、このようなカーネル関数を得るには「入力空間で少しずつ中心をずらしたガウス関数たちを特徴とすればよい」んです。中央列の最上段のカラフルな関数たちが特徴  \phi(x) の各成分  \phi_i(x) です。そして中段は  \phi_i(x) \phi_i(x') を各 i についてプロットしたものです。中段を全て足し合わせたのが最下段 k(x,x') = \phi(x)^{\top} \phi(x') です。最上段のプロットには訓練データの位置に破線を入れていますが、最上段の例えばピンク色の曲線を、「破線との交点の値」倍すると中段のピンク色の曲線に移ります。もちろん他の色の曲線(特徴)もそうです。そして、スクリプトをコピーして x_dash の値を変えてみてもらえると確かにわかると思いますが、バツ印の位置を横にずらせば中央列最下段の分布も横にずれます。色々な位置の訓練データ点に対してこのように訓練データを中心にした密度を張り巡らせることができるので理想的なわけです。

f:id:cookie-box:20180305231302p:plain:w60

えー Jupyter Notebook 起動すんのめんどい。バツ印をずらした場合もわかるようにアニメーションとかにしといてよ。

f:id:cookie-box:20180305232608p:plain:w60

忙しいから無理。

f:id:cookie-box:20180305231302p:plain:w60

じゃあ左の列と右の列は何なの? これらは訓練データを中心に広がる分布、って感じじゃないけど…。

f:id:cookie-box:20180305232608p:plain:w60

右の列は2値分類のようなイメージですね。迷惑メールかどうかの分類で、「このメールよりあちら側なら迷惑メールっぽいし、こちら側なら迷惑メールっぽくない」という煙幕でしょうか。この最上段の各特徴は中心を少しずつずらしたシグモイド関数です。左の列は解釈が難しいので説明を割愛します。

f:id:cookie-box:20180305231302p:plain:w60

ぶん投げた!

f:id:cookie-box:20180305232608p:plain:w60

ただ、左の列の最上段の各特徴は \phi_i(x) = x^i です。あ、それに関連して、基底関数というのは「あるクラスに含まれる関数たちの全てがその関数の線形結合であらわせる」ようにとった関数たちのことです。例えば  f(x)=a_0 + a_1 x + a_2 x^2 + \cdots の形でかけるような関数全体(多項式関数というクラス)の基底関数は \phi_i(x) = x^i \; (i=0, 1, 2, \cdots) ととることができます。そのままですね。なので、多項式関数状のカーネル関数の煙幕を張りたいならばその基底関数である \phi_i(x) = x^i を特徴たちとして採用する、ということもあり得るのかもしれませんが…ただ、特徴たちをカーネル関数に統合するときに否応なく個々の訓練データ点での \phi_i(x') が線形結合の重みとなりますから、煙幕の形をそう自由にデザインできるわけでもないと思うんですよね。「こんな形状のカーネル関数にしたいから、その基底関数を特徴にしよう」という流れは僕にはどうにも想像しにくいんです。なので、ここで「基底関数」というのもどうなのかなと。しかし、僕が不勉強なだけかもしれません。あるいは、ことさら関数解析を意識しているわけではなく、この節の用語としての基底関数(basis function)の可能性もありますけど。…6.2 節の1段落目はこれくらいにしておいて、カーネル関数を直接定義するアプローチの方に移りましょう。

f:id:cookie-box:20180305231302p:plain:w60

2段落目からはそのアプローチだな。カーネル関数を直接決めていいならそりゃ好きな形にできそうだけど、でも全く好きにカーネル関数を選ぶわけにはいかないよな。特徴の内積の形でかけなきゃだから。で、k(x,x') がそういう形になっているかは (6.12) 式みたいに個々にまともに調べる方法もあるけど、それじゃ大変だよな。そこで k(x,x') が有効なカーネルであるための必要十分条件があるのか。それは、成分が k(x_n, x_m) で与えられるグラム行列 K が半正定値であること…半正定値って何?

f:id:cookie-box:20180305232608p:plain:w60

実対称行列の固有値が全て非負である場合、半正定値行列であるといいますね。他の同値な条件で表現されることもありますが。

本には証明がありませんが、グラム行列 K(サイズ  N \times N)が半正定値行列ならば k(x,x') が有効なカーネルであることを軽く確認しておきましょうか。
  • K は実対称行列なので、適当な直交行列 P によって P^{\top} KP = \LambdaK の各固有値(重複含む)を対角成分とする対角行列とすることができます。証明をかくのは面倒なので以下を参照してください。
  • K固有値は全て非負なのでルートをとることができます。\Lambda の対角成分に全てルートをかぶせた行列を \hat{\Lambda} とします。P^{\top} KP = \Lambda = \hat{\Lambda} \hat{\Lambda} とかけます。さらに  K = P \hat{\Lambda} \hat{\Lambda} P^{\top} = P \hat{\Lambda} ( P \hat{\Lambda} )^{\top} \equiv \Phi \Phi^{\top} とできます。 \Phi の各行をとれば k(x_n, x_m) = \phi(x_n)^{\top} \phi(x_m) とかけます。これだけ満たされれば 6.1 節の議論は問題なく成り立ちますね。つまり、k(x,x') は有効なカーネル関数であるといえるでしょう。

f:id:cookie-box:20180305231302p:plain:w60

ふーん…でも、K が半正定値行列かどうかの確認もそれはそれで大変なんじゃない? K って縦幅も横幅もデータサイズなんだろ?

f:id:cookie-box:20180305232608p:plain:w60

正定値性判定のアルゴリズムは色々あるんじゃないですか。詳しくないですが。まあでも、本にも、新しいカーネルをつくるには既知のカーネルから (6.13)~(6.22) の式変形を用いてつくるのが便利みたいにかいてありますね。というかこれらの式の証明全部演習問題なんですね…。

f:id:cookie-box:20180305231302p:plain:w60

この6ページの画像のピクセルがどうって何?

f:id:cookie-box:20180305232608p:plain:w60

(6.12) 式の例だと  \phi(x) = (x_1^2, \sqrt{2} x_1 x_2, x_2^2)^{\top} のように特徴に入力変数の2次の組み合わせが全てあらわれていますよね。それと同様に、特徴にあらゆる M 次の組み合わせがあらわれるということですね。まあそれを画像のピクセルに喩えられても、M 個のピクセルを掛け合わせるときに個々に重みを設定するわけでもないので畳み込みともちょっと違う気がしますが。しいていうなら、未知画像のあらゆる M 個のピクセルを掛け合わせた特徴たちを、訓練データたちが与える重みで畳み込むという感じはするかもしれませんが。

f:id:cookie-box:20180305231302p:plain:w60

カーネル法は画像認識にもつかえるって言いたいのかな? わかんないけど。その次のガウスカーネルはジュンが前回やっちゃったな。

f:id:cookie-box:20180305232608p:plain:w60

7ページには、ガウスカーネルの2乗ノルムの部分を好きな非線形カーネルに置き換えてもよいとありますね。

(次回があれば)つづく

スクリプト
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 13,11
rcParams['font.size'] = 15

# カーネル関数 k(x, x') = Σ_{i=1}^M phi_i(x)*phi_i(x')
def kernel(phi, x, x_dash, M):
    k_ = 0
    for i in range(1, M + 1):
        k_ = k_ + phi(i, x) * phi(i, x_dash)
    return k_

# プロット用
def my_plot(phi, x, x_dash, M, axes):
    # (上段) 基底関数をそれぞれプロット
    for i in range(1, M + 1):
        axes[0].plot(x, phi(i, x))
    xmin, xmax, ymin, ymax = axes[0].axis()
    axes[0].plot(x_dash, ymin, 'kx', markersize=10)
    axes[0].vlines([x_dash], ymin, ymax, 'black', linestyles='dashed')
    # (中段) 足し合わせる前のカーネル関数をプロット
    for i in range(1, M + 1):
        axes[1].plot(x, phi(i, x) * phi(i, x_dash))
    axes[1].plot(x_dash, axes[1].axis()[2], 'kx', markersize=10)
    # (下段) カーネル関数をプロット
    axes[2].plot(x, kernel(phi, x, x_dash, M), 'k-')
    axes[2].plot(x_dash, axes[2].axis()[2], 'kx', markersize=10)

# =========================
# 色々な基底関数を定義
# =========================
# 多項式型の基底関数 phi_i(x) = x^i
def phi_poly(i, x):
    phi_ = 1
    for i_ in range(i):
        phi_ = phi_ * x
    return phi_

# ガウス型の基底関数 phi_i(x) = exp(-0.5 * (x-mu)^2 / s^2)
mu = np.arange(-1.2, 1.01, 0.2) # 平均の位置は -1.2(0番目は無視), -1.0, -0.8, ..., 1.0
s = 0.2 # 標準偏差
def phi_gaussian(i, x):
    return np.exp(-0.5 * (x - mu[i]) * (x - mu[i]) / (s * s))

# シグモイド型の基底関数 phi_i(x) = 1 / (1 + exp(-(x-mu)))
a = 5.0 # シグモイド関数のゲイン
def phi_sigmoid(i, x):
    return 1.0 / (1.0 + np.exp(-a * (x - mu[i])))

# =========================
# プロット
# =========================
fig = plt.figure()
M = 11 # 基底関数をいくつとるか=特徴の次元数
x = np.arange(-1.0, 1.01, 0.01) # プロットする x の範囲
x_dash = -0.5 # ある訓練データ点
my_plot(phi_poly, x, x_dash, M, [fig.add_subplot(3, 3, 1), fig.add_subplot(3, 3, 4), fig.add_subplot(3, 3, 7)])
x_dash = 0.0 # ある訓練データ点
my_plot(phi_gaussian, x, x_dash, M, [fig.add_subplot(3, 3, 2), fig.add_subplot(3, 3, 5), fig.add_subplot(3, 3, 8)])
my_plot(phi_sigmoid, x, x_dash, M, [fig.add_subplot(3, 3, 3), fig.add_subplot(3, 3, 6), fig.add_subplot(3, 3, 9)])
plt.show()

パターン認識と機械学習 下(第6章:その1)

以下の本を読みます。上巻を読むこともあります。何か問題点がありましたらご指摘いただけますと幸いです。

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

次回:まだ
f:id:cookie-box:20180305232608p:plain:w60

上巻読んでないですよね?

f:id:cookie-box:20180305231302p:plain:w60

プロデューサーが面倒だから上巻はとばしてもいいんじゃないかって。

f:id:cookie-box:20180305232608p:plain:w60

何がいいんですか!?

f:id:cookie-box:20180305231302p:plain:w60

だから下巻の最初の「第6章 カーネル法」から読み始めてみたんだけど、最初の段落は「ふつう回帰/分類モデルを学習するときにはパラメータ  w を学習して、学習しちゃったら訓練データ自体はもう使わないよね」って話で、それはそうだよねって感じなんだけど、次の段落では「でも予測時にも訓練データを使う手法も中にはあったよね」って言ってるんだ。上巻でそういう手法をいくつか扱ったみたいなんだけど、特に 2.5.1 節では「カーネル関数」について触れたみたいだから、上巻に戻って読んでおいた方がいい気がするんだ。だってこの章のタイトル「カーネル法」だし…。

f:id:cookie-box:20180305232608p:plain:w60

上巻に戻るの早くないですか!?

f:id:cookie-box:20180305231302p:plain:w60

だからここから上巻の「2.5.1節 カーネル密度推定法」にとぶけど、ここではデータは  {\mathbb R}^D の要素みたいだな。観測データから真の密度  p(x) を推定したいって状況で…なぜかわかんないけど  {\mathbb R}^D のある部分領域  \mathcal{R} に注目すると、データが  \mathcal{R} 内に入る確率は  P=\int_{\mathcal{R}}p(x)dx のはずで、この  P をつかうと「 N 個のサンプルを得たときに  \mathcal{R} 内に入っているサンプルの割合」の期待値と分散が出るのか。…なんで  PP(1-P)/N になるんだっけ?

f:id:cookie-box:20180305232608p:plain:w60

いやそこに思いっきり二項分布書いてありますし、この本の読者はそこで引っかからない想定だと思うんですが…本題から逸れますけど復習しておきますか。いまの状況で  N 個のサンプルのうち  K 個が  \mathcal{R} 内に入る確率は、表が出る確率が  P のコインを  N 回投げるという試行をしたときに内  K 回が表になる確率に等しいですよね。この試行で例えば最初の  K 回のみが表になる確率は  P^K (1-P)^{N-K} ですが、 N 回のうち  K 回が表になるパターンは他にもあります。つまり、全部で  _{N}{\rm C}_{K} パターンありますから、 K 回表になる確率は結局  _{N}{\rm C}_{K} P^K (1-P)^{N-K} ですね。であれば、表が出る回数の期待値と分散は計算するだけです。

   \displaystyle E[K] = \sum_{K = 0}^N K \, _{N}{\rm C}_{K} P^K (1-P)^{N-K} = \sum_{K = 1}^N \frac{K \cdot N!}{K! (N-K)!} P^K (1-P)^{N-K}
   \displaystyle \qquad \,  \, = NP \sum_{K = 1}^N \frac{(N - 1)!}{(K - 1)! (N-K)!} P^{K-1} (1-P)^{N-K} = NP
   \displaystyle V[K] = E[K^2] - E[K]^2 = \sum_{K = 0}^N K^2 \, _{N}{\rm C}_{K} P^K (1-P)^{N-K} - N^2 P^2
   \displaystyle \qquad \,  \, = \sum_{K = 1}^N \frac{K(K-1) \cdot N!}{K! (N-K)!} P^K (1-P)^{N-K} + \sum_{K = 1}^N \frac{K \cdot N!}{K! (N-K)!} P^K (1-P)^{N-K} - N^2 P^2
   \displaystyle \qquad \,  \, = N(N - 1)P^2 \sum_{K = 2}^N \frac{(N - 2)!}{(K - 2)! (N-K)!} P^{K-2} (1-P)^{N-K} + NP - N^2 P^2 = NP(1-P)
したがって、 E[K/N]=E[K]/N=P および  V[K/N]=E[K]/N^2=P(1-P)/N になるでしょう。

f:id:cookie-box:20180305231302p:plain:w60

ありがとジュン。1個1個のサンプルが  \mathcal{R} の中にあるかどうかがコインの表か裏かって感じなのか。続く (2.244) 式は「すごくたくさんの N 個のサンプルを得たら \mathcal{R} の中に入ってるサンプル数  K のばらつきは小さくなってもうぴったり NP も同然だろ」ってことだよな、たぶん。で、(2.245) 式は「領域 \mathcal{R} が小さかったらこの領域内で確率密度は一様だろ」ってことか。(2.246) 式はそれらを合わせて「じゃあもう領域 \mathcal{R} 内の確率密度はサンプルがこの中に入った割合 K/N を体積  V で割った値  \displaystyle \frac{K}{NV} でいいだろ」って感じかな。

f:id:cookie-box:20180305232608p:plain:w60

1次元でイメージするとヒストグラム状の経験分布のようなものでしょうね。

f:id:cookie-box:20180305231302p:plain:w60

それで、K を固定し V を推定するのが K 近傍法…? K 近傍法って前にプロデューサーがやってたよな。確か、未知データをクラス分類するときに、その未知データに近い既知データを K 個取ってきて、その K 個の中に多く含まれるクラスに分類するやつだっけ。

f:id:cookie-box:20180305232608p:plain:w60

まあそんな感じですね。なるほど、K 近傍法をつかうのであれば訓練データを捨ててはいけませんね。しかし、 K 近傍法はクラス分類の手法であって密度推定の手法では…ああ、122~123ページを読むと密度推定手法としての K 近傍法が紹介されていますね。K=5 なら 5 と決めて、その場所場所でサンプルを 5 個含む領域の大きさを測って、その大きさほどの広がりをもつ滑らかな関数を重ね合わせて推定密度とするということですか。この節では観測データから真の密度を推定するのに、空間全体の密度の形を考えようというのではなく、空間を少しずつ切り分けて、その場所場所の密度はどうなっているんだろうというアプローチを取るんですね。だから部分領域  \mathcal{R} などというものを持ち出したんですね。

f:id:cookie-box:20180305231302p:plain:w60

まだ120~121ページ読んでないのに勝手に先進むなよ! …120ページに戻ると、 V の方を固定するのがカーネル密度推定法なのか。(2.247) 式の k(u) は、 u を中心とする超立方体内では 1、その外では 0 を取る関数だな。これがカーネル関数の一例なんだ。ただの超立方体じゃん…。で、(2.248) 式と (2.249) 式は何やってんのこれ?

f:id:cookie-box:20180305232608p:plain:w60

(2.246) 式の主張は「ある点 x における密度は x を含む小さな領域  \mathcal{R} 内に含まれるサンプルの割合を、その領域の体積で割った値に等しい」でしょう。これを具体的に各サンプルの座標の式で書き表したのが (2.249) 式というだけですよ。小さな領域に具体的な形は必要ですから、ここでは x を中心とする一辺の長さが h の超立方体に決めてしまって、するとその中に含まれるサンプルの個数は「各サンプルの周りに超立方体の領域を張り巡らせたときに位置 x に何重にその領域が重なるか」と定式化できるので、結局 (2.249) 式が成り立つわけですね。

f:id:cookie-box:20180305231302p:plain:w60

ふーん…あれ、てことは、もうこれで密度推定終わったんだ。なんかすぐ終わったな。ん? でも、「この場合、立方体の淵で、人為的な不連続が生じてしまう(上巻121ページ)」?…そりゃそうだよな。ガッタガタだよこれ…だから、ふつうは超立方体じゃなくガウス関数をつかうのか。…なあジュン、(2.250) 式って結局「N 個の観測データからの推定密度 p(x) を、各サンプルの周りに張り巡らせたガウス分布 1/N ずつの比率で混合したものとする」ってことだよな。そういわれたらあーなんかそういう分布を考えるのってアリそうだなって何となくわかるしこの式からでよかったんじゃない?

f:id:cookie-box:20180305232608p:plain:w60

それを x における推定密度とすることをどうやって正当化するんです? どのような条件下でそのような近似ができるか明らかですか? 何となくじゃダメでしょう?

f:id:cookie-box:20180305231302p:plain:w60

はい。

f:id:cookie-box:20180305232608p:plain:w60

(2.250) 式以降は、h の大きさのバランスを取ることが肝要なことや、カーネル関数 k(u) にはもっと別の関数を選んでもよいことなどが書かれていますね。

f:id:cookie-box:20180305231302p:plain:w60

それで話は最近傍法に移って…K 近傍密度推定法だと空間全体上の積分が発散するってなんで?  V の方を固定してたのを K の方を固定するように変えただけでなんでそんなことになっちゃったの?

f:id:cookie-box:20180305232608p:plain:w60

僕も厳密には追えてませんが、空間内の各点について K 番目に近いサンプルまでの距離を h として (2.246) 式を適用すると空間全体上での積分が発散するんでしょうね。本節における密度推定法は「領域  \mathcal{R} はその内部で密度一定とみなせるほどじゅうぶん小さい」という仮定から出発していますが、体積 V の方を可変にしてしまうともはや領域を小さいままに保てません。そこが正規化できなくなった敗因なのでは。

f:id:cookie-box:20180305231302p:plain:w60

あー確かに。既知データたちからめっちゃ遠い地点から既知データを K 個囲むように領域 \mathcal{R} を取ろうと思ったらめっちゃでかくなるよな、たぶん。

f:id:cookie-box:20180305232608p:plain:w60

123ページは、クラス分類としての K 近傍法の話に触れられていますね。未知データがクラス  c に帰属する事後確率が K_c / K になることが示されています(K_c は未知データに近い既知データを K 個取ってきたときにその中に含まれているクラス  c のデータの個数)。「未知データの近所に多いクラスに分類する」というのが結局どのような分類をすることを意味するのかこれでわかりましたね。…これくらいで下巻の1ページに戻りましょうか。

f:id:cookie-box:20180305231302p:plain:w60

そうだ、元々下巻読んでたんだったよな。1ページの第2段落に戻ると、2.5.1節でみてきたようなカーネル密度推定法や最近傍法では、入力空間中の任意の2点の類似度を測れなければならないっていってるな。そりゃ空間内の各点に既知データからの距離に応じた密度が割り当てられるしな。それに、「予測には時間がかかる(1ページ)」…ってのは、カーネル密度推定法だったら全ての既知データとの間の距離を測って密度を足し上げなきゃいけないからしんどいってことだよなたぶん。で次の段落は…双対表現って何??

f:id:cookie-box:20180305232608p:plain:w60

いや、それをこの後 6.1 節で扱うので…ただ少し調べると双対表現(dual representation)という用語は数学や心理学にもあるようですが、ここでの双対表現とはカーネル法用語というか、PRML用語といった感じがしますね…そもそも英語では「対となる表現」というごく一般的な意味の言葉でしょうし。

f:id:cookie-box:20180305231302p:plain:w60

ふーん? まあ今は気にしなくていいってことか。この第3段落によると、「目的の回帰/分類を達成するようなパラメータ  w を学習する」って問題を対となる別の問題に書き換えることができて、そっち側で解いたときは予測時にカーネル関数をつかうことになるってことかな? でそっち側で解くときは何か x を特徴空間へ飛ばす写像 \phi(\cdot) を決めておいて、カーネル関数は k(x,x')=\phi(x)^{\top}\phi(x') という式になるってことか。え、カーネル関数ってガウス関数とかだったよな? こんな内積みたいな感じじゃなかったと思うんだけど。

f:id:cookie-box:20180305232608p:plain:w60

ちょうど章末の演習問題の 6.11 が、ガウスカーネルが無限次元の特徴ベクトルの内積であることを示せ、となっていますよ。ちょっとやってみましょうか。(6.23)~(6.25) 式よりはじめて、

   k(x, x') = \exp(-|| x - x' ||^2 / 2\sigma^2)=\exp(- x^{\top}x / 2\sigma^2 - x'^{\top}x' / 2\sigma^2 + x^{\top}x' / \sigma^2)
   \qquad \quad \, \, = \exp(- x^{\top}x / 2\sigma^2) \exp(- x'^{\top}x' / 2\sigma^2) \exp(x^{\top}x' / \sigma^2)
   \qquad \quad \, \, \displaystyle = \exp(- x^{\top}x / 2\sigma^2) \exp(- x'^{\top}x' / 2\sigma^2)\left( 1 + \frac{x^{\top}x'}{ \sigma^2} + \frac{(x^{\top}x')^2}{ 2\sigma^4} + \cdots \right)
これを眺めると、例えば  x が2次元の場合、 \phi(\cdot) を以下のような無限次元への非線形変換とすればよいのではないでしょうか。これでガウス関数k(x,x')=\phi(x)^{\top}\phi(x') という内積みたいな感じになりますよね。
   \displaystyle \phi(x) = \exp(- x^{\top}x / 2\sigma^2) \left( \begin{array}{c} 1 \\ x_1 / \sigma \\ x_2 / \sigma \\ x_1^2 / \sqrt{2} \sigma^2 \\ x_2^2 / \sqrt{2} \sigma^2  \\ x_1 x_2 / \sigma^2 \\ \vdots \end{array} \right)

f:id:cookie-box:20180305231302p:plain:w60

うおおいだから勝手に先進むなよ! あと無限次元て! そんなホイホイ次元を無限にしていいのかよ!?

f:id:cookie-box:20180305232608p:plain:w60

まだ予測時にカーネル関数をつかうような解き方の中で特徴ベクトル  \phi(x) がどのような制約を受けるのかわかりませんね。

f:id:cookie-box:20180305231302p:plain:w60

序文に戻るぞ。といっても、序文の続きにはカーネル関数の例はこうですとか、カーネル関数を導入するとカーネルトリックというまだよくわからないテクニックが使えるようになりますってあるだけで、カーネル関数によって何がよくなるのかとかはよくわかんないな。

f:id:cookie-box:20180305232608p:plain:w60

第2段落で「訓練データを記憶しておくような(=ノンパラメトリックな、ですね)方法は訓練データ間の距離の定義も要るし時間もかかる」といっていたので、カーネル関数がそれを解決してくれるのではないかと思うんですがね。よくわからないので 6.1 節に進みましょうか。

f:id:cookie-box:20180305231302p:plain:w60

ここでは (6.2) 式を最小化したいっていうシチュエーションなのか。線形回帰モデルの重みを求めるときの目的関数は確かにこんな感じだよな。じゃあ t_n は被説明変数か。

f:id:cookie-box:20180305232608p:plain:w60

そして右辺第2項は、重みベクトルが密になりすぎないようにするための正則化項でしょうね。

f:id:cookie-box:20180305231302p:plain:w60

上巻ぶっとばしてるから色々察する必要があるな…それで、(6.3) 式は  \displaystyle \frac{\partial J}{\partial w}=0 を式変形してるんだよな。え、行列  \Phi って? どっから出てきたの? あと計画行列って何?

f:id:cookie-box:20180305232608p:plain:w60

N=3 で特徴空間が2次元の場合だとこうですよ。ちゃんと行列ですよね。計画行列は上巻の第3章139ページで出てきた言葉のようですが、見たまま各行が各データの特徴量ベクトルになっているような行列のことのようですね。なので、行数はデータサイズ、列数は特徴空間の次元数になりますね。

 \displaystyle \sum_{n=1}^3 a_n \phi(x_n) = \begin{pmatrix} a_1 \phi(x_1)_1 + a_2 \phi(x_2)_1 + a_3 \phi(x_3)_1 \\ a_1 \phi(x_1)_2 + a_2 \phi(x_2)_2 + a_3 \phi(x_3)_2 \end{pmatrix} = \begin{pmatrix} \phi(x_1)_1 & \phi(x_2)_1 & \phi(x_3)_1 \\ \phi(x_1)_2 & \phi(x_2)_2 & \phi(x_3)_2 \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \\ a_3 \end{pmatrix}
 \displaystyle  \qquad \qquad \quad \, \, = \begin{pmatrix} \phi(x_1)_1 & \phi(x_1)_2 \\ \phi(x_2)_1 & \phi(x_2)_2 \\ \phi(x_3)_1 & \phi(x_3)_2 \end{pmatrix} ^{\top} \begin{pmatrix} a_1 \\ a_2 \\ a_3 \end{pmatrix} = \Phi ^{\top} a

f:id:cookie-box:20180305231302p:plain:w60

結局また上巻読みにいったのか。

f:id:cookie-box:20180305232608p:plain:w60

上巻 3.1.1 節の最小2乗法はいま下巻 6.1 節でみている最適化問題正則化項がないバージョンなんですね。しかし、こちらの下巻では正則化項があるから同じ形にはなりません。というか、(6.3) 式で w から a に主役を交代させてしまっていますね。こうして w最適化問題a最適化問題にすりかえるのが双対表現ですか。

f:id:cookie-box:20180305231302p:plain:w60

aw の関係は (6.4) 式だから、変数変換みたいなもんなんだな。でもさジュン、変数変換ってそうすることで簡単になるときとかにするんじゃないの? (6.5) 式とかかえって見た目おぞましくなってない?

f:id:cookie-box:20180305232608p:plain:w60

見た目で判断するのはよくないのでは…ほら、この目的関数は (6.7) 式のように各要素がデータの各組合せのカーネル関数になっているような行列=グラム行列でかくことができますよ。(6.6) 式は丁寧にかくと以下ですね。各要素が m 番目のデータが n 番目のデータの位置につくる密度になっています。

 \displaystyle K = \Phi \Phi ^{\top} = \begin{pmatrix} \phi(x_1)_1 & \phi(x_1)_2 \\ \phi(x_2)_1 & \phi(x_2)_2 \\ \phi(x_3)_1 & \phi(x_3)_2 \end{pmatrix} \begin{pmatrix} \phi(x_1)_1 & \phi(x_2)_1 & \phi(x_3)_1 \\ \phi(x_1)_2 & \phi(x_2)_2 & \phi(x_3)_2 \end{pmatrix}
 \displaystyle  \qquad \qquad \, \, \, =  \begin{pmatrix} \phi(x_1)^{\top} \phi(x_1) & \phi(x_1)^{\top} \phi(x_2) & \phi(x_1)^{\top} \phi(x_3) \\ \phi(x_2)^{\top} \phi(x_1) & \phi(x_2)^{\top} \phi(x_2) & \phi(x_2)^{\top} \phi(x_3) \\ \phi(x_3)^{\top} \phi(x_1) & \phi(x_3)^{\top} \phi(x_2) & \phi(x_3)^{\top} \phi(x_3) \end{pmatrix}

f:id:cookie-box:20180305231302p:plain:w60

あっカーネル関数…って別に、カーネル関数だったら何がうれしいのかわからないし…。とりあえず先に行くか。(6.8) 式って、(6.3) 式を (6.4) 式に代入するとほんとにこうなる?

f:id:cookie-box:20180305232608p:plain:w60

こうですね。 K^{\top}=K であることに注意してください。この後は自分でやってください。

 \displaystyle a = - \frac{1}{\lambda} \begin{pmatrix} w^{\top} \phi(x_1) - t_1 \\ w^{\top} \phi(x_2) - t_2 \\ w^{\top} \phi(x_3) - t_3 \end{pmatrix} = - \frac{1}{\lambda} \begin{pmatrix} a^{\top} \Phi \phi(x_1) \\ a^{\top} \Phi \phi(x_2)  \\ a^{\top} \Phi \phi(x_3) \end{pmatrix} + \frac{t}{\lambda}= - \frac{1}{\lambda} (a^{\top} \Phi \Phi^{\top})^{\top} + \frac{t}{\lambda}= - \frac{1}{\lambda} K a + \frac{t}{\lambda}

f:id:cookie-box:20180305231302p:plain:w60

あーほんとだ…でもやっぱりこの答えってさ、上巻 3.1.4 節の正則化最小二乗法の (3.28) 式とそっくりじゃない? w じゃなくて a を主役にした意味とかあったの? 双対表現とかいう名前まで付けてさ。

f:id:cookie-box:20180305232608p:plain:w60

wa の見た目がそうまるで違うというわけではないですね。ただある位置における予測を a の方をつかってかくと「各訓練データ点を相手としたカーネル関数」のみの式(  \phi(x) は出てこない )でかける点が違うんじゃないでしょうか。この式 (6.9) に関連してメモしておくと  a^{\top} \Phi \phi(x) = \phi(x)^{\top} \Phi ^{\top} a なので  k(x) = \Phi \phi(x) であってこれはデータサイズの次元をもつ縦ベクトルですね。各要素は各データが位置 x につくる密度ですね。この k(x) の係数が a です。w が特徴空間における最適化問題の答えだったのなら、a は「実空間の場所場所に各データ点が張る密度」の世界における最適化問題の答えだったんですね。

f:id:cookie-box:20180305231302p:plain:w60

ふーん、そういわれると係数の意味がちょっとわかりやすくなるのかな…あれでも、4ページに、双対表現にすると求める逆行列のサイズが大きくなってしまうってかいてある?

f:id:cookie-box:20180305232608p:plain:w60

上の方で書き下した例でも、w は2次元で a は3次元でしたよね。それに現実では訓練データ点が3点などということはなく、いくらでも大きくなりますからね。特徴空間の次元数はたかがしれているかもしれませんが。そしていまハヤトがいったように、わかりにくい特徴空間を明示的に扱うことを避けることができ、無限次元の特徴すら取り扱うことができるとありますね。

f:id:cookie-box:20180305231302p:plain:w60

特徴を無限次元とれるっていわれると確かにすごそうだけど…特徴の次元って多ければ多いほどいいってもん? それに、無限次元っていっても特徴量どうしの内積が入力空間でまともな密度になるような特徴量を取らないとダメだよな…? そんなすごい自由って感じでもないような…。なんかいまいちカーネル法のよさが見えてこないな…この先を読むとわかるのかな…?

(次回があれば)つづく