2019-04-17

雑記

【下のノートについて】多変量正規分布があったとき、 $x$ から $y$ への線形回帰と $y$ から $x$ への線形回帰は一致しないらしいんですが、あまり納得できなかったのでどう説明したら納得できるか考えてみたんですが、やっぱりあまり納得できなかったし、あとグラフに確率密度の等高線を描いた方が一目瞭然だったんじゃないかと思いました。「こう考えれば自明では」というのがあったら教えてください。よろしくお願いします。
gistaeb8a4a5dd292f783c69c28907da96ee

2019-04-10

ベイズ統計の理論と方法：ノート2

本読み機械学習

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。おかしい点がありましたらコメント等でご指摘いただけますと幸いです。

ベイズ統計の理論と方法

作者: 渡辺澄夫
出版社/メーカー: コロナ社
発売日: 2012/03/01
メディア: 単行本
購入: 1人クリック: 4回
この商品を含むブログ (8件) を見る

前回：ノート1 ／次回：まだ

前回のノートでは確率と確率密度がごっちゃになっているところがあります。事象が起きる確率に対しては事象の選択情報量ですが、確率変数の確率密度関数に対しては選択情報量とはよばないと思います。よばないと思いますがよび方に困るので、以下「選択情報量のような量」とよんでいます（自由エネルギーのことなんですが、自由エネルギーって何っていうことばとして）。
- そんな変な呼び方をしなくても「負の対数尤度」などとよべばさしつかえはないですが、個人的にそうよぶとそれって大きいほどどういう意味だっけというのがなんかすぐわかんないので採用しません。

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

前回読んだ1～10ページの内容をふりかえります。この本はベイズ推測について解説する本であり、ベイズ推測とは何かというのは 1.1 節に定義がありました。つまり、「ベイズ推測する」とは

サンプル $x^n$ 内の各点が独立に従っている分布は、確率モデル $p(x|w)$ の事後分布 $p(w|X^n)$ による平均 $\mathbb{E}_w \bigl[p(x|w) \bigr] \equiv p(x|X^n)$ だろうと考える（★）

ことに他なりません。ここで事後分布 $p(w|X^n)$ というのは私たちが何となく知る事後分布ではなく、逆温度 $\beta$ を用いて定義される (1.5) 式であることに注意してください（といっても、 $\beta = 1$ の場合は $p(X^n|w)$ と $\varphi(w)$ からベイズの定理によって導かれる、私たちがよく知る事後分布と同じですが）。

しかし、「～だろうと考える（★）」というだけでは、「じゃあそう考えれば？」という話です。この本で私たちが学ぶのは、「『～だろうと考える（★）』というのは結局どういうことなのか」でしょう。

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

6ページの1.～3.の問題を考えていくことでそれに迫るんだね。きっと「『～だろうと考える（★）』ことで、こんな条件下で、こんな風によい推測を達成できる」という感じの出口になるのかな。そういう出口があることがわかっているから、私たちはベイズ推測を利用することができる。そういう出口を知らずにベイズ推測をするなら、「理由はないけどとりあえず（★）だろうと考えました」ということにしかならない。

しかし「よい推測」とは何かという問題があります。私たちは現実に推測を行う場面で真の分布を知り得ないので、真の分布にどれだけ近づけたかを確認することはできません。しかしそれでも、 $\bigl( \, q(x), \, p(x|w), \, \varphi(w) \, \bigr)$ に拠らない数学的な法則が存在して、推測の「限界」を議論することができるということですが…。1章の続きではその足掛かりとして、「自由エネルギー $F_n(\beta)$ 」と「汎化損失 $G_n$ 」が定義されました。

自由エネルギー $F_n(\beta)$ は事前分布 $\varphi(w)$ と確率モデル $p(x|w)$ と逆温度 $\beta$ とサンプル $X^n$ による量で、分配関数 $Z_n(\beta)$ の対数の温度倍のマイナス1倍です。 $\beta=1$ のとき $F_n(1)$ は周辺尤度（その確率モデルでその事前分布にしたがうあらゆるパラメータの下でサンプル $X^n$ が観測される確率密度） $Z_n(1)$ の対数のマイナス1倍に等しいです。確率密度を確率のようなものと思えば、 $F_n(1)$ は「その確率モデルでその事前分布にしたがうあらゆるパラメータの下でサンプル $X^n$ が観測されるという事象の選択情報量のような量」です。これは、 $F_n(1)$ をあらゆるサンプル $X^n$ の現れ方について平均すると（＝ $q(x^n)$ で平均すると）、 $q(x^n)$ と $Z_n(1)$ の交差エントロピーになります。つまり、「真の分布」と「事前分布で平均した確率モデル」がどれだけ似ているかを測る指標になりそうです。もっとも、知り得ない $q(x^n)$ で平均するということは不可能なので議論が必要ですが…。
汎化損失 $G_n$ は真の分布 $q(x)$ と予測分布 $p^\ast (x)$ の交差エントロピーです。しかし、やはり知り得ない $q(x)$ で平均するということは不可能なので、経験損失 $T_n$ から見積るということですが…。 $T_n$ は「予測分布 $p^\ast (x)$ の下で各サンプル $X_i$ が観測される事象の選択情報量のような量の全サンプル平均」です。 $T_n$ は全サンプルからつくった櫛形の（ $n$ 本のデルタ関数が立った）経験分布と予測分布の交差エントロピーととらえることもできるかもしれません。

そして、 $\beta=1$ のとき $p^\ast (X_{n+1}) = Z_{n+1}(1) /Z_n(1)$ が成り立つので、この両辺の逆数の対数より、「予測分布（ただし $\beta=1$ ）の下で未知データ $X_{n+1}$ が観測される事象の選択情報量のような量」が「事前分布で平均した確率モデルの下でサンプル $X^{n+1}$ が観測される事象の選択情報量のような量（ $F_{n+1}(1)$ ）」から「事前分布で平均した確率モデルの下でサンプル $X^{n}$ が観測される事象の選択情報量のような量（ $F_{n}(1)$ ）」を差し引いたものに等しいことがわかります。前回はここまで読みました。…ここから $X^{n+1}$ で平均してその下の式になりますか？？右辺第1項は $X^{n+1}$ の関数ですが、右辺第2項は $X^{n}$ のみの関数であるようにみえますが…それに左辺は…左辺はこれ何の関数でしたっけ？？

右辺第1項と右辺第2項は丁寧にかくとこうかな。

$\displaystyle \int q(x_1) \cdots q(x_n) q(x_{n+1}) F_{n+1} (1) dx_1 \cdots dx_n dx_{n+1} = \mathbb{E} \bigl[ F_{n+1}(1) \bigr]$
$\displaystyle \int q(x_1) \cdots q(x_n) q(x_{n+1}) F_n (1) dx_1 \cdots dx_n dx_{n+1}$
$\displaystyle = \int q(x_1) \cdots q(x_n) F_n (1) dx_1 \cdots dx_n \cdot \int q(x_{n+1}) dx_{n+1}$
$= \mathbb{E} \bigl[ F_n(1) \bigr]$

左辺は先に $x_{n+1}$ について積分しようか。

$\displaystyle - \int q(x_1) \cdots q(x_n) q(x_{n+1}) \log p^\ast (x_{n+1}) dx_1 \cdots dx_n dx_{n+1}$
$\displaystyle = \int q(x_1) \cdots q(x_n) G_n dx_1 \cdots dx_n$
$\displaystyle = \mathbb{E} \bigl[ G_n \bigr]$

$\mathbb{E} \bigl[ \cdot \bigr]$ はあらゆるサンプル $x^{n+1}$ の出方に対する平均だから、10ページの一番下の式は確率的に変動する項を含まないね。任意の $\bigl( \, q(x), \, p(x|w), \, \varphi(w) \,\bigr)$ に対して、 $\beta=1$ のベイズ推測を実施したときの汎化損失の期待値は、ベイズ推測を実施する前の自由エネルギー（選択情報量のような量）の期待値がどれだけ増加するかに等しいことになる。

…なるほど、もし仮にその事前分布と確率モデルの下で $X^n$ を観測する事象の選択情報量も $X^{n+1}$ を観測する事象の選択情報量も常に変わらないというなら、 $n+1$ 個目のデータには全く「新たな情報」「意外さ」がない、 $X^n$ を観測したらそれがどんな $X^n$ であっても次に観測される $X_{n+1}$ が確実にわかってしまう、そんな状況だね。そんな状況では汎化損失の期待値もゼロだ。確実に $X_{n+1}$ がわかるんだから誤差は生じない。でも、 $n+1$ 個目のデータ $X_{n+1}$ に僅かでも「新たな情報」があれば、その「新たな情報」はベイズ推測に誤差を生じさせる。 $X_{n+1}$ を観測するまで得られない情報がある状況なんだから、完璧な推測はできない。汎化損失の期待値はゼロにならない。…といったけど、正確には汎化損失は「誤差」って感じじゃないね。KL情報量じゃなくて交差エントロピーだから、完璧に予測分布を $q(x)$ にしても汎化損失 $G_n$ はゼロにならない。 $q(x)$ のエントロピー（連続分布なので微分エントロピー）が理論下限だ。

え、えっと？ともかく、何が何の変数で、何が確率変数なのかややこしいですね。改めて整理します。10ページの上・中・下3箇所の数式は、それぞれ「確率密度」「選択情報量のような量」「交差エントロピー」にみえて、2ステップ更新されたようにみえるんです。でも、右辺は2ステップなのですけど、左辺は3ステップあったのですね（中から下への更新に、先に $x_{n+1}$ で積分して、次に残りの変数で積分するという2ステップの更新が含まれています）。つまり、右辺（上側の表）と左辺（下側の表）でそれぞれ主人公が以下のように交代しています。

$Z_n(1)$	事前分布で平均した確率モデルの下でサンプル $X^{n}$ が観測される事象の確率密度 or サンプル $X^{n}$ の下での事前分布で平均した確率モデルの周辺尤度（ $X_n$ により確率的に変動する確率変数）
$F_n(1)$	事前分布で平均した確率モデルの下でサンプル $X^{n}$ が観測される事象の選択情報量のような量 or サンプル $X^{n}$ に対して事前分布で平均した確率モデルを仮定したときの系の自由エネルギー（ $X_n$ により確率的に変動する確率変数）
$\mathbb{E} \bigl[ F_n(1) \bigr]$	「真の分布」と「事前分布で平均した確率モデル」の交差エントロピー or 真の分布の下でのあらゆるサンプル $X_n$ の出方に対する、事前分布で平均した確率モデルを仮定したときの系の自由エネルギーの期待値（確率変数ではない）

$p^\ast(x)$	予測分布の下で点 $x$ が観測される事象の確率密度（ $X_n$ により確率的に変動する $x$ の確率密度関数）
$\displaystyle \log \frac{1}{p^\ast(x)}$	予測分布の下で点 $x$ が観測される事象の選択情報量のような量（ $X_n$ により確率的に変動する $x$ の確率密度関数）
$G_n$	点 $x$ の真の分布と予測分布の交差エントロピー（ $X_n$ により確率的に変動する確率変数）
$\mathbb{E} \bigl[ G_n \bigr]$	真の分布の下でのあらゆるサンプル $X_n$ の出方に対する、点 $x$ の真の分布と予測分布の交差エントロピーの期待値（確率変数ではない）

右辺の表（上側の表）の方のこの色の文字は右辺の主人公を私たちがいくぶんよく知る統計学の言葉でかいたもので、この色の文字はこの本で出てくる言葉にならったものです。意図していなかったんですが自然とこの色の文字の方が「モデルはぶれない、サンプルがぶれる」という、伝統的な頻度論的な統計学に寄った表現になっていますね。比べてこの色の文字は「サンプル所与の下での、モデルの関数」といった感じです。左辺については $G_n$ までは予測分布 $p^\ast (x)$ というものがあるという前提でのこの分布に対する話なので「サンプルとモデルどっちがぶれる？」という要素はありません。ただ、 $G_n$ から $\mathbb{E} \bigl[ G_n \bigr]$ のステップは「モデルがぶれる」という立場ですね。モデル $p^\ast (x)$ がぶれないならば、汎化誤差の「期待値をとる」などという操作は考えられないはずです。

…汎化損失と自由エネルギーにある関係が成り立つのはわかりました。11ページは、なぜ確率密度の対数のマイナスをとるのかという話をしていますね。…これ、「確率密度は $e^{-E}$ の形であることが多いので $E$ を取りたいから」「 $E$ はエネルギーと実感できるから」って、どちらも突拍子もなく感じるんですが。なぜ確率密度がそんな形をしていることが多いなどといえるんです？だいたいどこからエネルギーが出てきたんですか？自由エネルギーとはそのような名前なのだと割り切っていましたが、ここでは紛れもなく物理のエネルギーの話をしていますよね？

まあ確率密度が $e^{-E}$ の形であるかどうかとエネルギーの登場はおいといても、確率 $P$ を対数のマイナス1倍をとって $- \log P$ として「選択情報量」というものさしでみると「大きいほどめずらしい、それが起きたと知ったときの価値が高い出来事だ」って何となくわかりやすかったよね。こっちのものさし方が推測のよさを測るのに感覚に合ってそうだ。もちろん確率 $P$ のままでも推測の誤差を議論することはできると思うけど…でも、元々よく起きる出来事か、レアな出来事かで確率を1%誤る重大さって違う気がするよね。ある年に名古屋では年間に100日くらい雨が降って、10日くらい雪が降ったらしい。真の値より10日多めに110日雨の予報を出してもたぶんあまり怒られないけど、真の値より10日多めに20日雪の予報を出したらたぶんクレームがくるだろうし天気予報を信用してもらえなくなるよね。

いや、雪の日は実際タイヤチェーンなどさまざまな準備が必要ですから、事象自体が誤りの重大さに関係していてその喩えはあまり適切ではないのでは…まあ雰囲気はわかりますが。ただ、それなら確率の誤差ではなく誤差率をみるということもできると思いますが…。

率って扱いにくいし、誤差が小さければ誤差率は対数差分で近似できるしね。 $(P^\ast - P)/P$ がゼロに近ければ $(P^\ast - P)/P \approx \log \bigl( 1 + (P^\ast - P)/P \bigr) = \log P^\ast - \log P$ だよね。

なんと、確かにそうなりますね。

あと確率密度が $e^{-E}$ の形かってのも、何か物理法則にしたがうデータだったらそうなる見込みがあるからね。その辺は統計力学の話になるけど。

やっぱり物理の話じゃないですか…まあ $-\log P$ にさまざまな解釈が与えられるというのはいいです。次節に進みましょう。…事後分布や予測分布を「解析的に計算できない」ことが多い？「解析的に計算できない」って何ですか？ (1.5) 式や (1.8) 式ってそんなに何か困難な要素があるんですか？？

「解析的に解けない」っていうと解が既知の演算や関数でかきあらわせないって意味だね。「5次方程式は一般的に解けない」というのは一般の係数の加減乗除べき根で方程式の解をかきあらわすことができないって意味だし。ただここでいう既知の演算や関数が何かはまだわからないかな。まあこの節に「計算できる例」があるってことなんだからそこから推し量ることはできるんじゃない？

そうですね、この節の例を計算してみましょう。 $\displaystyle p(x|w) = v(x) e^{f(w) \cdot g(x)}$ で $\displaystyle \varphi(w|\phi) = \frac{1}{z(\phi)} e^{\phi \cdot f(w)}$ ということです。 $f( \cdot)$ と $g( \cdot)$ の制約はわかりませんが確率モデルが積分できることは確かですね。このとき、分配関数及び事後分布を計算すると、
$\displaystyle Z_n(\beta) = \int_W \frac{1}{z(\phi)} e^{\phi \cdot f(w)} \prod_{i=1}^n \bigl( v(X_i) e^{f(w) \cdot g(X_i)} \bigr)^{\beta} dw$
$\displaystyle \qquad \; \; \, = \frac{1}{z(\phi)} \left( \prod_{i=1}^n v(X_i)^\beta \right) \int_W e^{\phi \cdot f(w) + \sum_{i=1}^n \beta f(w) \cdot g(X_i)} dw$
$\displaystyle \qquad \; \; \, = \frac{z \bigl( \phi + \sum_{i=1}^n \beta g(X_i) \bigr)}{z(\phi)} \prod_{i=1}^n v(X_i)^\beta$
$\displaystyle p(w|X^n) = \frac{1}{Z_n} \frac{1}{z(\phi)} e^{\phi \cdot f(w)} \prod_{i=1}^n \bigl( v(X_i) e^{f(w) \cdot g(X_i)} \bigr)^{\beta}$
$\displaystyle \qquad \quad \; \; \; = \frac{1}{Z_n} \frac{1}{z(\phi)} \left( \prod_{i=1}^n v(X_i)^\beta \right) e^{\phi \cdot f(w) + \sum_{i=1}^n \beta f(w) \cdot g(X_i)}$
$\displaystyle \qquad \quad \; \; \; = \frac{1}{z \bigl( \phi + \sum_{i=1}^n \beta \cdot g(X_i) \bigr)} e^{\bigl( \phi + \sum_{i=1}^n \beta g(X_i) \bigr) \cdot f(w)}$
$\displaystyle \qquad \quad \; \; \; = \varphi \left(w \middle| \phi + \sum_{i=1}^n \beta g(X_i) \right) \equiv \varphi(w | \hat{\phi})$
これは…ハイパーパラメータ $\phi$ が更新された形になりましたね。

…確率モデル $\displaystyle \prod_{i=1}^n p(X_i|w)^\beta$ と事前分布 $\varphi(w|\phi)$ の積が $a \varphi(w|\hat{\phi})$ の形にかける（ $a$ は $w$ に依存しない係数）なら $Z_n(\beta) = a$ になって事後分布が $\varphi(w|\hat{\phi})$ になるってことか。逆に事後分布をこの形にしたいなら…やっぱり $\displaystyle \prod_{i=1}^n p(X_i|w)^\beta$ と $\varphi(w|\phi)$ の積が $a \varphi(w|\hat{\phi})$ の形にかけないといけない。 $\displaystyle \prod_{i=1}^n p(X_i|w)^\beta$ と $\varphi(w|\phi)$ の積の $w$ への依存性はこの形でないといけない。「事前分布を事後分布にする」ことを「ハイパーパラメータ $\phi$ を更新する」形で達成することは、以下を満たす確率モデルと事前分布を選ぶことと同じように感じる。

$\varphi(w|\phi) \displaystyle \prod_{i=1}^n p(X_i|w)^\beta = a \varphi(w|\hat{\phi})$ 指数型分布と共役な事前分布はばっちりこれを満たす。指数型分布でない $p(x|w)$ であったら絶対にこれを満たさないのかはすぐわからないな…でも、14ページをみるとどうもそうみたいだね。

「解析的に解ける」というのは結局「パラメータの更新で済む」という意味だったんでしょうか….？もちろん、「パラメータの更新で済む」ような確率モデルと事前分布を選ばなければならないという意味ではないと14ページにありますね。あくまで人間の都合です。…1.3節はさまざまな推測方法ということですが…最尤推測って $\beta=\infty$ のベイズ推測なんですか？？

最尤推測では $w$ の分布でなくただ1つの $w$ を求めることになる。 $p(w|X^n)=\delta(w - w_{ML})$ って感じかな。15ページの $\varphi(m|\hat{\phi_1}, \hat{\phi_2})$ は $\beta \to \infty$ でデルタ関数に近づくと思う。近づく先は $\sum_{i=1}^n X_i / n$ だからパラメータ $m$ の最尤推定値だね。一般の場合にも最尤推定値にたつデルタ関数っていうのを示すのはどうやるのかすぐ思いつかないけど…。

2019-05-07 追記
＞一般の場合にも最尤推定値にたつデルタ関数っていうのを示すのはどうやるのかすぐ思いつかない
以下のような感じな気がします。
https://twitter.com/CookieBox26/status/1125552313572020224

事後確率最大化というのは、ベイズ推測をした上で事後分布の最頻値を採用する方式なのですかね。他方、平均値を採用するのが平均プラグイン？なぜプラグインというのか…。他の方式としては…平均場近似？？

148ページからかいてあるよ。これはどちらかというと、ベイズ推測において事後分布を平均場近似によって得ることを変分ベイズ法というのかな。

17ページの (3) にある場合にベイズ推測と最尤推測が同じといえるのかいえないのかといった気になることがかいてありますね…まあいまはとばしますか。18ページからは解析的に解ける場合で数値実験しているようですね。…節変わって23ページに、「確率モデルが仮のものである場合」とありますね。もちろん確率モデルがわかっている場合とわかっていない場合とあるのはわかるんですが、だから何を言いたかったんでしょうか…。

「1000人の試験結果のときと100万人の試験結果のときで推測されるパラメータがずれていてもいい」って感じに読めるね。仮に100万人の試験結果が正規分布を $J$ 個混合した分布で上手くフィッティングできても「中学生全体が $J$ 個のグループに分けられると結論されたのではない」というのは、 $J$ は「100万人の試験結果の解析に有用なパラメータ」に過ぎないってことだね。

じゃあ、中学生全体の試験結果が何個にクラスタリングされるかを研究することはできないってことですか？

やっちゃいけないのは、特に根拠なく混合正規分布を仮定して、100万人の試験結果から推測した $J$ でもって「中学生というものは $J$ 個のグループに分けられます」と結論付けることだ。「この100万人の中学生は」と付けるなら間違ってはないだろうけどね。本当に「中学生は $J$ 個のグループに分けられるのではないか」という仮説を検証したいなら、きっと人数を増やして収束していくかとか確認すべきで、現実にはそこで頓挫するんじゃないかな。

（ノート3があれば）つづく

2019-03-28

カーネル法入門：ノート2

本読み機械学習

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらコメント等でご指摘いただけますと幸いです。

カーネル法入門―正定値カーネルによるデータ解析 (シリーズ多変量データの統計科学)

作者: 福水健次
出版社/メーカー: 朝倉書店
発売日: 2010/11/01
メディア: 単行本
クリック: 19回
この商品を含むブログ (11件) を見る

以前カーネル法を勉強したときの記事：パターン認識と機械学習下（第6章：その6〈終〉） - クッキーの日記

前回：ノート1 ／次回：まだ

前回のあらすじをまとめると、

データを分析するには（つまり、ある変数とある変数が依存し合っているかを調べたり、それに基づいて回帰したり分類したりするには）、線形な依存性のみを調べるのでは不十分なこともあるので、非線形な依存性も調べたい。
言い換えると、データを $\Phi : \Omega \to H$ なる特徴写像で特徴空間 $H$ に写像することで非線形な特徴も引き出して、この特徴空間で分析をしたい。
しかし、あらゆる高次項を特徴に追加したのでは、組み合わせが爆発してしまい、主成分分析のために標本共分散行列の固有値分解をしたり、回帰分析のためにデータ行列の擬似逆行列を求めたりすることができない。
が、少なくとも特徴空間における主成分分析は、行列の各成分が $k(X_i, X_j) = \langle \Phi(X_i), \Phi(X_j) \rangle$ であるような行数と列数がデータサイズの（有限の大きさの）行列の2乗の固有値分解をすることで達成できる。ので、 $k(X_i, X_j)$ の値さえ評価できればよいことになる。ここでポイントとなったのは、

特徴空間における第1主軸 $f$ は $\displaystyle f = \sum_{i=1}^N a_i \Phi(X_i)$ の形でかける（もしかけない成分があったら、かけない成分をゼロにしてその分かける成分を引き伸ばした方が各データの $f$ との内積の分散が大きくなる）。
そうすると解くべき最大化問題の中にデータは $k(X_i, X_j)$ の形でしか出てこない。

…普通に考えれば $\Phi(X_i)$ の標本分散共分散行列を固有値分解しなければならないところを、 $f$ は $H$ の部分空間 $H_0$ から探せば十分ということに着目して $f$ を表現し、最適化問題から $\Phi(X_i)$ を排除してしまうとは、なかなかトリッキーですね…。

カーネルトリックというくらいだからね。

しかし、 $f$ は $H_0$ から探せば十分、というのは通常の主成分分析でも成り立つ話ですよね？通常の主成分分析の言葉でいうと、第1主軸 $u$ を求めるのに、 $u$ をあらゆる向きの $m$ 次元単位ベクトルから探すのではなくて、 $\displaystyle u = \sum_{i=1}^N a_i X_i$ でかけるものを探せば十分、ということになりますよね。

まさしくそうだけど、通常の主成分分析では固有値分解すべき標本分散共分散行列が $m \times m$ の大きさだったのが、カーネルPCAとして解くと固有値分解すべき行列が $N \times N$ の大きさになるから、どちらがありがたいのかという話だね。

あー、元データの次元数 $m$ よりデータサイズ $N$ の方が大きければかえって計算量が大きくなってしまいますね…。8ページからは主成分分析からまた変わって、リッジ回帰という例がありますね。リッジ回帰とは、線形モデルを最小2乗フィッティングするのにパラメータベクトルのL2ノルムの2乗に比例する正則化項を付けて、「誤差2乗和を小さくするパラメータベクトルを求めよ。ただし、パラメータベクトルのL2ノルムの2乗は小さく抑えること」といったもののようですね。その解 $\hat{a}$ は…「容易に求められるように」？何気ない「容易に」が読者を傷つけるのでは？

本だと尺の都合があるからね。最小化したいのは以下だよね。

$\displaystyle \sum_{i=1}^N || Y_i - a^\top X_i ||^2 + \lambda || a ||^2$ これを $a$ で偏微分しようか。ベクトルで微分するときの公式をつかうよ。 $\begin{split} \displaystyle \frac{\partial}{\partial a} \left( \sum_{i=1}^N || Y_i - a^\top X_i ||^2 + \lambda || a ||^2 \right) &= 2 \sum_{i=1}^N ( Y_i - a^\top X_i ) X_i + 2 \lambda a \\ &= 2 \sum_{i=1}^N Y_i X_i - 2 \sum_{i=1}^N (a^\top X_i) X_i + 2 \lambda a \\ &= 2 X^\top Y - 2 X^\top Xa + 2 \lambda a \end{split}$ あとはこれが $0$ になるように $a$ を解けばいい。8ページ中ほどの式になるのはわかるんじゃないかな。

まあそこまで説明してもらえばわかります。その式の下の、リッジ回帰はどんなときに利用されるかというので、「 $X^\top X$ が退化していたり」というのは何ですか？

退化というのは $X^\top X$ が正則でない、逆行列をもたないという意味だね。もし今の回帰問題で正則化項がなかったら、つまり、 $\lambda = 0$ だったら、回帰係数の解は $\hat{a} = (X^\top X)^{-1} X^\top Y$ になるけど、 $X^\top X$ が逆行列をもたなかったらこれは求まらない。

確かに $X^\top X$ が逆行列をもっていなかったら $\hat{a}$ は求まりませんね。いえでも、解けないというのはどういうことなんでしょう？それに、解けないから正則化項 $\lambda || a||^2$ を導入して解決、というのもどういうことなのか…。

$X^\top X$ が正則でないというのは、 $X^\top X$ の階数が $m$ に満たない、さらに言い換えると $X$ の階数が $m$ に満たないって状況だね。そうなっちゃうのは以下のようなときかな。

データ数 $N$ を上回るほど、データの次元 $m$ が巨大であるとき。
データセットが線型従属なとき（例えば、同じ点が重複して含まれているとか）。

こんな状況だと、誤差2乗和を最小にする $\hat{a}$ は無数に出てきちゃうんだよね。そうだな…以下の PRML の 146 ページに Figure 3.4 があるよね。

Christopher Bishop at Microsoft Research

Figure 3.4 の左側の図を例にすると、これはちゃんと $X^\top X$ が退化していないケースで、誤差2乗和（青い等高線）が最小になる点があるんだよね（青い点）。誤差2乗和はこの青い点を底にしたお椀になってる。でも、 $X^\top X$ が退化しているとこれみたいなお椀にならない。スノボのハーフパイプみたいにずーっと谷底が続いてたり、最悪まっさらな平面で「全体的に底」って感じになっちゃう。そうなると最小点を選べなくて困る。その解決策として、「じゃあ別のお椀をもってくればいいじゃん」というのが正則化項だ。Figure 3.4 の左側の図の赤い等高線が原点を底にした「L2ノルムのお椀」だね。このお椀を重ね合わせれば最小点が求まるって寸法だね（※）。もちろん、元々の $X^\top X$ が退化していない場合でも正則化のお椀を重ねていい。Figure 3.4 は誤差2乗和の青いお椀と正則化の赤いお椀を重ねて、重ねてできた新しいお椀の底になるのが $w^\ast$ ってことだね。

※	厳密には、もし $X^\top X$ が固有値 $-\lambda$ をもっていたら、正則化項を導入しても $X^\top X + \lambda I_N$ は正則になりません。が、ピンポイントでそうなることはあまりないのではないかと思っています。

なるほど…誤差2乗和のみでは回帰係数が選びきれないときでも選ぶ基準を新たに与えたんですね。原点を底にしたお椀ということは「選べなかったら原点になるべく近い解にしといて」って感じですよね…？

そうだね。ちなみに、L2ノルムのお椀は Figure 3.4 の左側の赤い等高線のように丸いけど、L1ノルムのお椀はその右側のようにひし形をしてるね。お椀っていうかひっくり返した四角錐だね。だから、L1ノルムで正則化するラッソ回帰では一部の成分がゼロになるような解が求まりやすい。この図でも $w_1 = 0$ の点が解になってるね。

じゃあラッソ回帰は「なるべく平面 $w_1 = 0, \; w_2 = 0, \; \cdots$ に貼り付く解にしといて」って感じですかね。…まあラッソ回帰は置いておいて、今回はリッジ回帰を特徴空間でやりたいわけです。特徴空間でのリッジ回帰の目的関数は式 (1.9) ですね。また特徴空間 $H$ での内積とノルムが出てきていますね…。

データ空間でやっていた「 $X_i$ の各成分を回帰係数 $a$ で足し合わせてそれを予測値としようとする」というのは、特徴空間 $H$ でやろうとすると「 $\Phi(X_i)$ とベクトル $f$ との内積をとってそれを予測値をしようとする」ということだからね。

しかし $\Phi(X_i)$ というベクトルはいま得体の知れないもので、カーネルPCAのときは最終的に $\Phi(X_i)$ は排除されたのですよね。リッジ回帰ではどうなるのでしょう。リッジ回帰でも解は $\displaystyle f = \sum_{j=1}^N c_j \Phi(X_j)$ の形でかけるとありますね。

もし $f$ に直交補空間の成分 $f_\bot$ があったらそれを丸々削った方が正則化項の $|| f||^2$ が抑えられるからね。

全くですね。ということは…今回もまた最小化したい式の中にデータは $k(X_i, X_j)$ の形でしか出てきませんね。 $\Phi(X_i)$ は出てきません。

最小化したい式の中に出てこないだけじゃなく、予測値を計算するときにも $\Phi(X_i)$ が出てこないことが肝心だね。 $\Phi(X_i)$ を直接取り扱いたくないから。でも未知の点 $x$ に対して「 $\Phi(x)$ とベクトル $f$ との内積をとるということをするとやっぱりデータは内積 $k(X_j, x) = \langle \Phi(X_j), \Phi(x) \rangle$ の形でしか出てこない。

えっと、まとめると、特徴空間でそのままリッジ回帰する（予測したい未知データの特徴 $\Phi(x)$ と内積を取るべき特徴空間内のベクトル $f$ そのものをもとめる）のではなく、 $f$ が $\displaystyle f = \sum_{j=1}^N c_j \Phi(X_j)$ の形でかけることに注目して $c_j$ を求める問題にしたのですよね。うーん、カーネルPCAでもカーネルリッジ回帰でも「解はこの形でかけるじゃん」というところがどうも降って湧いたようなんですが…。

最適化問題中の $\Phi(X_j)$ を意地でも $\langle \Phi(X_i), \Phi(X_j) \rangle$ の形にするために、 $f$ から $\Phi(X_i)$ を絞り出したって考えればそんなに降って湧いてもないんじゃないかな。

絶対に $\Phi(X_i)$ を直接取り扱いたくないのだという強い意思といったところですか…。ところで、カーネルPCAとカーネルリッジ回帰ではデータを特徴写像 $\Phi$ で特徴空間 $H$ にとばしましたが、 $\Phi$ や $H$ の詳細には立ち入っていませんでした。 $H$ についてわかっているのは、ベクトル空間であることと、内積 $\langle \cdot, \cdot \rangle$ をもつということだけですね。 $\Phi$ もどんな写像なのかまるでわかりません。直接取り扱っていないので。常にカーネル関数 $k$ を通して取り扱っています。…この $H$ と $\Phi$ は適当に決めていいんでしょうか。

ではないね。まずカーネル関数 $k(X_i, X_j) = \langle \Phi(X_i), \Phi(X_j) \rangle$ の値が評価できないといけないからこれが発散するような $\Phi$ だと困ると思うし…そもそも私たちは直接 $\Phi$ を決めないし…まあ [tex\Phi] とは何かは次章だね。

そもそも非線形な依存性を調べたくて $\Phi$ などという特徴写像を持ち出したのですから、どのような $\Phi$ が可能なのかが肝心な気がしますが…。ところで今回カーネルリッジ回帰を学びましたが、python の scikit-learn でカーネルリッジ回帰ができるようなので実行してみます。iris データの全品種ごちゃ混ぜで、がく片の長さからがく片の幅を予測するという何がしたいのかわからない例ですが…カーネル関数 $k$ には色々なものが指定できるし、自分で実装したものを渡すこともできそうです。組み込みのものを利用するときはおそらく以下の pairwise.kernel_metrics にあるものが指定できると思います。以下のスクリプトではいじっていませんが、パラメータも変更できます。

sklearn.kernel_ridge.KernelRidge ‐ scikit-learn 0.20.3 documentation

sklearn.metrics.pairwise.kernel_metrics ‐ scikit-learn 0.20.3 documentation

f:id:cookie-box:20190328234405p:plain:w480

もともとが線形関数でも非線形関数でもフィッティングできなさそうなのでそもそも例として適当ではありませんが、線形フィッティングは駄目そうですね。 $k$ にpolyカーネル $k(x, y) = (x^\top y)^3$ を選んだときよりもrbfカーネル $k(x, y) = \exp(-\gamma || x-y ||)^2$ を選んだときの方が点の密度が大きいところを通れているといえばいるような。

（ノート3があれば）つづく

スクリプト

import numpy as np
from sklearn.kernel_ridge import KernelRidge
from sklearn.datasets import load_iris

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 8, 6
rcParams['font.size'] = 16

iris = load_iris()
#x = np.array(iris['data'][:,0][iris['target']==0]) # setosa の sepal length (cm)
#y = np.array(iris['data'][:,1][iris['target']==0]) # setosa の sepal width (cm)
x = np.array(iris['data'][:,0]) # sepal length (cm)
y = np.array(iris['data'][:,1]) # sepal width (cm)

print(x.shape)
print(y.shape)

x = x.reshape(-1, 1)
y = y.reshape(-1, 1)

print(x.shape)
print(y.shape)

# ----- オリジナル -----
plt.scatter(x, y, label='original')
# ----- linearカーネル -----
kr = KernelRidge(kernel='linear')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(linear)')
# ----- polyカーネル -----
kr = KernelRidge(kernel='poly')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(poly)')
# ----- rbfカーネル -----
kr = KernelRidge(kernel='rbf')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(rbf)')

plt.legend()
plt.show()

(150,)
(150,)
(150, 1)
(150, 1)

2019-03-22

カーネル法入門：ノート1

カーネル法

カーネル法入門―正定値カーネルによるデータ解析 (シリーズ多変量データの統計科学)

作者:健次, 福水
朝倉書店

Amazon

以前カーネル法を勉強したときの記事：パターン認識と機械学習下（第6章：その6〈終〉） - クッキーの日記

次回：まだ

まえがきから読み始めると冒頭からカーネル法とは何かきっぱり書いてありますね。つまり、「正定値カーネルの定める再生核ヒルベルト空間を用いたデータ解析の方法論」とのことです。…データ解析に「空間を用いる」とは何でしょう？「ほにゃららモデルを用いる」、ならわかるんですが…まえがきの続きによるとこの再生核ヒルベルト空間とは関数空間であるらしいですが、何が何やら…まあ置いておきましょう。冒頭にはカーネル法とは「データの非線形性や高次モーメントを扱うための方法論」ともありますね。データの非線形性や高次モーメントを扱うとは？

おそらく1章の2～3ページで挙げられているのが最も簡単な例だね。「データの線形性しか扱わない」というのは、「データの空間のどこかにデータを仕切るような超平面を置いて、その超平面から表側に垂直に離れるほどプラスのスコア、裏側に垂直に離れるほどマイナスのスコア、のようなスコア付けでもってデータを分類したり回帰したりしようとする」ということだろう。でもまっすぐな仕切りで仕切れないのが図1.1の例だね。「高次モーメントを扱う」というのは「2つの変数 $X, \, Y$ に対して $X$ の $n$ 次かつ $Y$ の $m$ 次のモーメントからわかる依存関係を捉える」って感じなんだろうか。図1.2は $X$ の方を2乗しないと相関係数が小さくなる例だね。相関係数がゼロであることは独立であることを意味しないというのは有名だね（ウィキペディアの「相関係数」の頁の絵）。

相関係数 - Wikipedia

そこまでごり押されなくても「線形の解析では不十分」であることに異存はないですよ…。確かに、図1.1で分類しようとしたり、図1.2で $X$ がどうなると $Y$ がどうなるかという関係を見出したりするには線形関数では足りませんね…もっと柔軟なモデルでなければ…解決策を思い付きました。決定木を使いましょう。

解決できるのはできると思うけど、本のタイトルが「決定木入門」になっちゃうかな。

図1.1の方は図の2次元空間の上に $f(x_1, x_2)$ という「その場所では＋か ○ のどちらか」という関数を見つけたいのですよね。それを $f(x_1, x_2) = \hat{f}(\phi(x_1, x_2)) = \hat{f}(z_1, z_2, z_3)$ とすることで $\hat{f}$ を線形関数で足りるようにしたんですね。図1.2の方は、 $f_1(x)$ と $f_2(y)$ の相関を最大化するのに、 $f_1(x) = \hat{f_1}(\phi(x)) = \hat{f_1}(x^2)$ とすることで $\hat{f_1}$ を線形関数で足りるようにしたという感じでしょうか。…なんというか、これらをみるとカーネル法とはむしろ「何が何でも線形関数で解く。線形関数で解けるようにデータの方が動くべき」というものにみえますね。ニューラルネットワークの方が人としての器が大きいのではないですか？あちらは「データはそのままでどんな関数でも表現する」でしょう？

人じゃないからね。カーネル法はカーネル法で（あるクラスの最適化問題に対する正定値カーネルの定める再生核ヒルベルト空間内での）最適解が線形関数の形でかけるという理論的裏付けと、さまざまな拡張があるんだよ。「線形関数で解けるようにデータの方が動くべき」というようにデータの方が一旦動いてくれたら、その上に線形な依存性に基づいたさまざまな分析手法を応用することができるんだから。

はあ。まあ何にせよカーネル法ではデータを都合よく動かさなければならないということですよね。しかし4ページをみると、あらゆる高次のモーメントを考慮するには組合せが爆発してしまうようですが。前途多難では？

にもかかわらず効率的に高次モーメントを抽出できますよというのがカーネル法の売りだからね。

どうやってそんなことを…まあそれがこれからの話題ですか。1.1.2節に移りましょう。データの空間が $\Omega$ で、それを移す実ベクトル空間が $H$ ですか。図1.1や図1.2の左側が $\Omega$ 、右側が $H$ ということですよね？ …実ベクトル空間というのは？

集合であって、和やスカラー倍が入ってて、和やスカラー倍が結合律やら分配律やらを満たすやつ。「実」ベクトル空間といったときはスカラーは $\mathbb{R}$ の元だね。

それは結局、ベクトルを集めた集合ではないんですか？ 3次元実数ベクトルが入っている $H$ は結局 $\mathbb{R}^3$ ということではないんですか？

どちらかというと実ベクトル空間の条件を満たす集合があったとき、その元をベクトルというかな。あと、3次元実数ベクトルを含む $H$ でも、 $\mathbb{R}^3$ 内のある平面だけとかある直線だけとかの実ベクトル空間はありうるね。

あー平面や直線の中で和やスカラー倍はできますね…まあ、データを和やスカラー倍がちゃんとできる空間にもってくるのだと理解しました。そうですね、元のデータが単語や文章などであったら、和やスカラー倍をしようとしても和やスカラー倍とは何か定かではありませんね。あれでも、 $H$ は内積ももたなければならないと続きにありますね。

データ分析では点がどのように密集（分布）しているかに興味があるからね。「ここに仕切りを入れられそうだ」「これが主成分だ」とかいうには点どうしの近さが決まっていないと始まらない。

それって決めないと決まらないものですか？内積って $x_1 y_1 + x_2 y_2 + x_3 y_3$ のように授業で習った気がするのですけど。

それは実数ベクトルの標準内積だね。それに決めてもいいよ。でもそれとは違う内積に決めてもいい。例えば、 $100 x_1 y_1 + x_2 y_2 + x_3 y_3$ みたいに、1つ目の軸方向だけやたら重視する内積でもいいよ。そうしたかったらね。内積が満たすべき性質は、付録の193～194ページだね。

なんと、そういう内積もありなんですか？一旦193ページにとびますか。…って、内積より前にまずノルムというものが出てきますね。その後に内積はこういうものというのも出てきて…内積が定義できればその $1/2$ 乗をノルムにできるってことですね。なるほど、副部長がいう定義も内積の性質を満たしますね。

対称行列 $A$ について、 $x^{\top}Ay$ は内積になるからね。

まあ何にせよ何かしら内積を決めたとして、でも $H$ が高次元だと内積計算にはコストがかかると…それはそうかもしれませんが。でもカーネル法ではその内積計算がカーネル関数 $k$ の評価で済む？えっと、話を整理しましょう。いま、手元のデータを分析するのに、線形の依存性を調べるのでは不十分なので、高次元の依存性を調べたくて、つまり高次元の項も含めた内積を計算したくて、それにはカーネル値 $k(X_i, X_j)$ を計算すればよい。…ちょっと最後が飛躍しすぎですね。

それはおいおいかな。そこがカーネル法の基本原理だからね。5ページからは実用例みたいな話が始まってるね。

カーネル主成分分析、ですか。通常の主成分分析もわからないのですが…python の scikit-learn でも使えるようなのでまずは実行してみますか。

sklearn.decomposition.KernelPCA — scikit-learn 1.0 documentation

データはそうですね、非線形性があるものがいいのですよね。以下にスイスロールのデータがあったのでこれをお借りすることにします。3次元で、y軸に垂直に切ると渦巻きになっているようなデータですね。

a Swiss Roll Dataset

上のデータに sklearn.decomposition の PCA と KernelPCA を適用したのが下図です（スクリプトは記事の最下部）。左上が元のデータをxy平面に落とした図、右上が元のデータをxz平面に落とした図ですね。左下が通常の主成分分析の第1主成分-第2主成分、右下がカーネル主成分分析の第1主成分-第2主成分です。上段の実線は通常の主成分分析の第1主軸の方向、点線が第2主軸の方向ですね。データ点は上のリンク先の gif と同じ4色の塗分けをしています。よくみると、スイスロールはこの4つのクラスタに切り分けられるのかもしれませんね。色が同じ部分がやや密集しているようです。
f:id:cookie-box:20190321180703p:plain:w540

f:id:cookie-box:20190321180703p:plain:w540

f:id:cookie-box:20190321180731p:plain:w540

…通常の主成分分析はほとんどxy平面に落としただけだね。データがよく広がっている方向をとらえられたという感じはしない。元々がスイスロールなので無理はあるんだけど。それに比べてカーネル主成分分析は、最初の2つの主成分でデータが広がっている方向をとらえることができたのかもしれない。非線形性を扱えるから、丸めたカーペットのように丸まったスイスロールを少し開くことができたって感じかな。

えっと、5ページの通常の主成分分析の復習をみると、通常の主成分分析というのは、データの空間の中で、すべてのデータをその単位ベクトルに正射影したときに正射影した後のデータたちの分散が最大になるような単位ベクトルの向きを求めて、それを第1主軸とするのですね。続けて第2主軸以降も求めるときは、第1主軸と直交する空間の中でそのようなベクトルを探していくと…それが「標本分散共分散行列 $V$ の大きい方から $d$ 個の固有値に対応する固有ベクトル」？ $u^{\top}Vu$ を最大化する $u$ を求めたいというのはいいですが、それが固有ベクトルというのはどういうことですか？？

落ち着いて追いかけようか。 $u^{\top}Vu$ を $u^\top u = 1$ の下で最大化したいわけだ。じゃあラグランジュ関数を $\displaystyle L(u, \lambda) = u^{\top}Vu - \lambda \left(u^\top u - 1\right)$ とおこう。これを $u$ で偏微分する。微分公式 $\displaystyle \frac{\partial}{\partial u} u^\top A u = (A + A^\top) u$ をつかえば $\displaystyle \frac{\partial}{\partial u} L(u, \lambda) = 2Vu - 2\lambda u$ だ。つまり、 $\displaystyle \frac{\partial}{\partial u} L(u, \lambda) = 0 \; \Leftrightarrow \; Vu = \lambda u$ になる。 $u$ は $V$ の固有ベクトルでなければならないね。それも、 $u^{\top}Vu = \lambda u^{\top} u =\lambda$ を最大化したいんだから、最大の固有値 $\lambda$ に対応する $u$ を第1主軸に選ぶべきだ。第2主軸以降はそれまでの軸と直交するように選びたいわけだけど、いま $V$ は対称行列だから固有ベクトルどうしが直交している（証明は以下などを参照してね）。

対称行列の固有値と固有ベクトルの性質の証明 | 高校数学の美しい物語

なるほど、ありがとうございます。6ページからは主成分分析を今度は特徴空間 $H$ で実行しようとしているんですね。先ほどは各データ $X_i$ を $u$ 方向に正射影した $u^\top X_i$ たちの分散を最大化しようとしていましたが、今度は $\langle f, \, \Phi(X_i)\rangle$ なるものの分散を最大化しようとしていますね。 $||f||=1$ の下で。

内積 $\langle \cdot, \cdot \rangle$ とノルム $||\cdot||$ さえ決まっている空間であれば、主成分分析はできるということだね。

それで、解は $\displaystyle f = \sum_{i=1}^N a_i \tilde{\Phi} (X_i)$ の形で探せば十分？？え、えっと、なぜこの形で十分かもわからないですが、 $f$ とは何なんですか？先ほど通常の主成分分析でデータの分散が最大となる向きの単位ベクトル $u$ を探したように、 $f$ も $H$ 内の単位ベクトルというわけではないんですか？これ、ベクトルなんですか？？

無論ベクトルだよ。でも部長はまだベクトルを「成分を縦に並べたもの」とイメージしているかな。ベクトルというのはあくまで「和やスカラー倍ができる集合の元」だ。それでしかない。そして「データの分散が最大となる向き」というのは、「各データとその単位ベクトルとの内積をとったものの分散が最大になるような単位ベクトル（ノルムが1のベクトル）」だ。この「向き」は、成分を縦に並べたものである必要はないよね。

そういうことですか…ではイメージを捨てて、内積で考えるようにしましょう。…すみません、部分空間と直交補空間とは何でしょう。

うーん、よく知っているベクトルのイメージに戻っちゃうけど、 $\mathbb{R}^3$ というベクトル空間を考えて、この空間の中のある平面だけみると、この平面は和とスカラー倍で閉じているから $\mathbb{R}^3$ の部分空間だ。直交補空間というのは、その平面に垂直な直線の中で閉じている空間だね。 $\displaystyle \left( \begin{array}{c} x \\ y \\ 0 \end{array}\right)$ とかけるベクトルたちの空間と $\displaystyle \left( \begin{array}{c} 0 \\ 0 \\ z \end{array}\right)$ とかけるベクトルたちの空間って感じかな。それでいま考えている話に戻ると、 $H$ という特徴空間の中には $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ の和やスカラー倍ではつくれない「ベクトル」がある。 $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ の和やスカラー倍でつくれない、というのは、 $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ たちがつくる部分空間に射影したベクトルを差し引いてもゼロでないベクトルが残っているってことだね。ただ、もし $f$ というベクトルがそのようなゼロでないベクトルを含むようなものであったとしても、その差し引き残ったベクトルってあってもなくても $\tilde{\Phi}(X_i)$ との内積の値に影響しない。いま「データの分散が最大となる向き」 $f$ を探したいわけだけど、もし $f = f_0 \oplus f_{\bot}$ の $f_{\bot}$ がゼロでなかったとしよう。このとき $||f_0||^2 = ||f|| - || f_\bot ||^2 =1-|| f_\bot ||^2 < 1$ だけど、 $f' = 1/(1-|| f_\bot ||^2)^{1/2} f_0 \oplus 0$ とすれば、 $f'$ への射影の方が $f$ への射影よりデータの分散が単純に $1/(1-|| f_\bot ||^2)^{1/2}$ 倍大きくなるんだよね。いま $f$ のノルムは1という制約があるんだから、 $\tilde{\Phi}(X_i)$ との内積の値に影響しない成分に割り振ってる余裕はないってとこかな。

だから $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ たちがつくる部分空間内で探せば十分、つまり、(1.6) 式の形で探せば十分ということですか。次元がよくわからない、というか次元という概念があるのかどうかもよくわからないベクトル空間でも、有限個の係数を探せばよいということになるんですね。そうなると、さっきとそっくりな解き方ができるんですね。しかし、先ほどと違って $a^\top \tilde{K}^2 a$ の $a$ というのはもはや「ベクトルの向き」を意味しませんね。そして、7ページの下部の式のように第 $p$ 主成分まで求まるんですね。これがカーネルPCA…。

（ノート2があれば）つづく

スクリプト

import numpy as np
import pandas as pd

df = pd.read_csv('swissroll.dat', sep='\s+', header=None)
df = df.apply(lambda x: x / x.std(), axis=0) # 標準化
df[1] = df[1] - df[1].mean() # y軸だけ中央がゼロでないのでシフト
print(df.head())
c = np.array(['tomato', 'dodgerblue', 'gold', 'mediumseagreen']).repeat(400) # 4色に色分け

# PCA
from sklearn.decomposition import PCA
pca = PCA()
pca.fit(df)
df_pca = pca.transform(df)

# KernelPCA
from sklearn.decomposition import KernelPCA
kpca = KernelPCA(kernel='rbf', gamma=1.5)
kpca.fit(df)
df_kpca = kpca.transform(df)

# プロット
%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['font.size'] = 16

x = np.array([-10, 10])

fig, axes = plt.subplots(ncols=2, figsize=(12, 5), sharex=True, sharey=True)
axes[0].scatter(df.iloc[:,0], df.iloc[:,1], c=c)
axes[0].plot(x, x * pca.components_[0][1] / pca.components_[0][0], c='black', linestyle='solid', linewidth=2) # 第1主成分
axes[0].plot(x, x * pca.components_[1][1] / pca.components_[1][0], c='black', linestyle='dotted', linewidth=2) # 第2主成分
axes[0].set_xlim([-3, 3])
axes[0].set_ylim([-3, 3])

axes[1].scatter(df.iloc[:,0], df.iloc[:,2], c=c)
axes[1].plot(x, x * pca.components_[0][2] / pca.components_[0][0], c='black', linestyle='solid', linewidth=2) # 第1主成分
axes[1].plot(x, x * pca.components_[1][2] / pca.components_[1][0], c='black', linestyle='dotted', linewidth=2) # 第2主成分

plt.subplots_adjust(wspace=0.15)
plt.show()

fig, (axL, axR) = plt.subplots(ncols=2, figsize=(12, 5))
axL.scatter(df_pca[:,0], df_pca[:,1], c=c)
axR.scatter(df_kpca[:,0], df_kpca[:,1], c=c)
plt.show()

2019-03-15

機械学習のための特徴量エンジニアリング：ノート2

本読み機械学習

機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)

作者: Alice Zheng,Amanda Casari,株式会社ホクソエム
出版社/メーカー: オライリージャパン
発売日: 2019/02/23
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

正誤表リンク： https://www.oreilly.co.jp/books/9784873118680/

前回：ノート1 ／次回：まだ

5章は「カテゴリ変数の取り扱い」というタイトルですね。以前参加した Kaggle のコンペティションでも「職業」というカテゴリ特徴があって「会社役員、会社員、…」といったカテゴリ値があった気がします。83ページの Effect コーディングというのは初めて聞きました。これ、サンフランシスコのデータの予測値が $w_1 + b$ 、ニューヨークのデータの予測値が $- w_1 - w_2 + b$ 、シアトルのデータの予測値が $w_2 + b$ になるということですか。だから切片 $b$ が全体平均になると。

正確にいうと、「サンフランシスコのデータとニューヨークのデータとシアトルのデータを等しい重みで平均した値」だと思うな。元データ中にカテゴリ値の偏りがあったら $b$ はそのデータ全体の平均ということにはならないはず。

84ページの最下部の「参照カテゴリに対する各カテゴリの相対的な影響をエンコードすることは」という箇所を読んでイメージが湧きました。ダミーコーディングにおける $(e_1, e_2)$ は、ニューヨークを原点にとったときの各カテゴリ値の特徴量なんですね。確かに、「なぜニューヨークが原点なのか」という感じはしますね。決定木のような原点の場所がどこかということは何も関係ないモデルでは関係なさそうですが。5.2.1節の特徴量ハッシングというのは、ランダムにカテゴリ値をまとめてしまうということですよね…例 5-3 のコードでは、おそらくレビュー文章か何かの単語列（word_list）を、m 次元の数値ベクトルに変換していますね。単語毎にベクトルのどの成分をインクリメントするかは、その単語の生のハッシュ値を m で割った余りで決めています。例 5-4 の方は似ていますが、インクリメントするかデクリメントするかもまたハッシュ値で決めています。こうすると大きなバイアスが発生しない？？どういうことでしょう。

ちょっと例を考えてみようか。元々単語のユニーク数が9個だったとする。この時点でじゅうぶん少ないけどあくまで例だからね。いま手元の2つの文章を、どの単語IDが何回現れるかで9次元にエンコードしたとする。仮に以下のような感じになったとする。

文章X： $(0, 1, 0, 1, 1, 0, 2, 3, 0)$
文章Y： $(1, 1, 0, 0, 1, 0, 2, 0, 1)$

これらのベクトルの内積は 6 だね。ここで、9 次元だと多すぎるから 3 次元に削減したいとなったとする。例 5-3 に基づく方法なら圧縮後の特徴ベクトルは以下だ。ここで、ハッシュ値は元々の単語IDの番号そのものとするよ。だから、上のベクトルを 3 つずつぱたぱたと折りたたむだけだね。

文章X： $(0, 1, 0) + (1, 1, 0) + (2, 3, 0) = (3, 5, 0)$
文章Y： $(1, 1, 0) + (0, 1, 0) + (2, 0, 1) = (3, 2, 1)$

他方、例 5-4 に基づく方法で圧縮するなら以下。元のベクトルの奇数番目にマイナスをかける。

文章X： $(-0, 1, -0) + (1, -1, 0) + (-2, 3, -0) = (-1, 3, 0)$
文章Y： $(-1, 1, -0) + (0, -1, 0) + (-2, 0, -1) = (-3, 0, -1)$

例 5-3 に基づく方法でも例 5-4 に基づく方法でも 9 次元のベクトルを 3 次元に圧縮できるけど、圧縮後の文章Xと文章Yの内積が前者の方法では 19、後者の方法では 3 になっている。まあどっちの方法でも圧縮前の 6 を保つことはできていないけど、前者の方法では明らかにどんな文章間も内積がインフレしそうだ。単に元のベクトルを折りたたむだけだからね。でも、文章を示す特徴量ベクトルの間の内積というのは、文章どうしがどれだけ類似しているかを示す肝心な量だから、次元を圧縮した表現にしただけで文章間がどんどん類似してしまうというのは好ましくないはずだ。後者の方法ではハッシュ値によって符号を変えることで内積が一方的に増えてしまうという事態を抑えている。もちろん次元を削減しているから元の内積を完全に保つことはできないけど。「内積の期待値が変わらない」の意味は原論文を読まないとわからないけど、おそらく上の思考実験をもっとたくさんの文章でやってみたり、あらゆるハッシュ関数でやってみたりしたら内積の平均値が保たれるんじゃないかな。

確かに bag-of-words ではコサイン類似度などで類似度が測られるのでしたっけ。次元を圧縮したいからといって内積が保たれない表現にしてしまっては台無しですね。元々の bag-of-words は「意味が近い文章どうしは距離が近くなっている」がゆえに文章を表現する特徴量たりえたのですから。意味が近くない文章どうしでなくても距離が近い表現など、適切な特徴量とはいえません。5.2.2 節は、カテゴリ毎の何かの最小値や最大値などでもよいのでしょうかね。94ページの最小カウントスケッチというのは？これはレアではないカテゴリも含めて $d$ 種類の $m$ 値へのマッピングを用意するということでしょうか。そして最小値を正式に採用する？うーん、やり方はわかるんですが、ハッシュ関数を $d$ 個にすると結局何がよかったのかとかなぜ最大値などではなく最小値をとるのかとかよくわかりません…。

Count-Min Sketch の原論文の [Cormode & Muthukrishnan, 2005] というのはおそらく以下の記事にリンクがあるものだね。

Count-Min Sketch

http://dimacs.rutgers.edu/~graham/pubs/papers/cmsoft.pdf（これは2012年の文書）

上の文書の3ページに書いてある手続きは本の94ページと全く一緒だ（絵も似ているね）。 $d \times w$ のテーブルに合計 $N$ のカウントを加えるなら、この手続きによるアイテム $i$ のカウントの推定値は $1 - (1/2)^{d}$ 以上の確率で誤差 $2N/w$ 以内になるらしい。理由は簡単だね。1つの行にのみ着目すると、アイテム $i$ が入っているマスに他のアイテムのカウントがどれだけ混入するか（誤差）の期待値は $N/w$ だ。となると、マルコフの不等式より、誤差が $2N /w$ 以上になる確率は $1/2$ 以下だ。これが $d$ 行あるから、全ての行で誤差が $2N /w$ 以上になる確率は $(1/2)^d$ になる。

そうか、「余計なアイテムのカウントが一番小さいマスを選びたいのだ」という気持ちであれば最小値を選ばなければなりませんね。統計量が何かのカウントであるとは限らないと思いますが。

$d \times w$ のマスがあるならそれを一行に展開して1つのハッシュ関数のみつかうのでは駄目なのか、って気もするけど、それだとあるアイテムたちはとても誤差が大きく、あるアイテムたちは誤差がないという偏りが出そう。 $d$ 種類のハッシュ関数をつかうことで最悪のアイテムでも誤差が少ないというようにできる。

95ページの一番下の段落はどういう意味でしょうか。「任意のデータ点の有無によって統計量の分布がほぼ変わらない」？

音楽を推薦するモデルをつくるのに、アーティストを特徴量にしたいけど、アーティストはきっと多いからそのままカテゴリ値として扱いたくない。だから、「レディー・ガガ」というカテゴリ値の変わりに「レディー・ガガの曲の再生回数の全ユーザ合計」のような連続値にしたい。けど、1人だけ異様にレディー・ガガの曲を再生しているユーザがいたらよくない。任意のユーザを抜いたとしても、あらゆるアーティストの再生回数合計の分布が変わらない必要がある、ということかな…いや、任意のアーティストを抜いても分布が変わらない、かもしれないかな…もしレディー・ガガだけ再生回数が断トツで多かったら、どのユーザにもレディー・ガガばかり推薦されることになっちゃいそうだし…。

そんなに再生回数が多かったらもう万人にレディー・ガガを推薦しておけばよくないですか？

個々のユーザの嗜好を予測しようとして？？

6章に入りますね。特異値分解ですか…以前に特異スペクトル変換法で扱いましたね。

異常検知と変化検知： 9章メモ（部分空間法による変化検知） - クッキーの日記

式 6-6 から式 6-7 はこうですね。 $\displaystyle (x_1^\top w)^2 + \cdots + (x_N^\top w)^2 = \left( \begin{array}{c} x_1^\top w \\ \vdots \\ x_N^\top w \end{array} \right)^\top \left( \begin{array}{c} x_1^\top w \\ \vdots \\ x_N^\top w \end{array} \right) = \left( X w \right)^\top \left( X w \right) = w^\top X^\top X w$ $w^\top w = 1$ の制約のもとで $Xw$ の長さを最大にするには、 $w$ は $X$ の最大の特異値によって引き延ばされる向きを向いていなければなりませんね。…せっかく次元削減できたのに、110～111ページで計算コストが高いだの色々言われていますね…。111ページに、この手法の実用場面として「時系列の異常検出」と言及されていますね。112ページで、ZCA は「相関関係を取り除くことができ」「画像のより面白い構造を見つけ出すことに集中」ってどんな画像になるんでしょう？

以下に ZCA をやっている記事があったよ。一番下の方に画像があるね。

CIFAR-10のPCA白色化/ZCA白色化

つづきは後で