独立成分分析
2025/03/09 解析 行列 数学
独立成分分析について
どんな時に使うか
例えばこんな問題を考えます。2人の人がしゃべっていて、それを2つのマイクで録音する。収録した音源を解析することで、元の二人の声に分離できないか?
他にも多数の電極で記録した脳波を成分ごとに分離したいとかに使えます。
問題の定式化
必要な仮定
ここでいくつかの仮定を置くことにします。
まず信号源はn個あり、それぞれ独立しています。ある信号源は他とは無関係で、どの一つをみてもその情報からは別の信号源について何もわかりません。
n個の信号源に対し受信機はn個あり、それぞれ信号源からの情報が決まった割合でMIXされたものを受け取ります。例えば、マイク1は人1の声を1/2の音量で受け取り、さらに人2の声を1/3の音量で受け取り、それの和が記録されるといった具合です。
また、それぞれの信号は時間差なく届くものとします。時間差がある場合に対応した方法もあるようですが、解説しません。
ノイズはないものとします。ノイズに対応した方法もあるようですが、これも解説しません。
数学的には
仮定を数学的にとらえると
\boldsymbol{s}をある時点でn個の信号源が出す信号を、 \boldsymbol{x}を同じ時点でのn個の観測値を表すベクトルとすると、ある行列𝑨が存在して \boldsymbol{x}=𝑨·\boldsymbol{s}と表される (※rank(A)=nとする)
また、それぞれの信号源から出る信号は平均が0で、E[\boldsymbol{s}]=\boldsymbol{0}である。(固定値を足したり引いたりしても独立性とは関係ない)
信号が\boldsymbol{s}となる確率は
\boldsymbol{s}=\{𝑠_1,𝑠_2, \dots 𝑠_𝑛\}^Tのそれぞれの確率の積
𝑝(\boldsymbol{s})=\prod_{i=1}^{n} p(s_i) である
分析の目標
観測値\boldsymbol{x}の確率分布を基にして、
\boldsymbol{y}=W·\boldsymbol{x}としたときに \boldsymbol{y}=\{y_1,y_2, \dots y_𝑛\}^Tが 𝑝(\boldsymbol{y})=\prod_{i=1}^{n} p(y_i) を満たすような分離行列𝑾を求めたい。このような行列は必ず存在して、例えばAの逆行列がそうである。
白色化+ICAという方法を紹介する。
前処理 白色化
まず、観測値\boldsymbol{x}=\{x_1,x_2, \dots x_𝑛\}^Tを白色化する。
そのため主成分分析を行い、各成分の分散を1に正規化しておく。その計算法は省略するが、結果としては白色化行列を𝑽とおくと
白色化された観測値 \boldsymbol{z}=V·\boldsymbol{x} は無相関かつ分散が1になる。(\boldsymbol{z}の共分散行列が単位行列になる)
※\boldsymbol{z}=\{z_1,z_2, \dots z_𝑛\}^Tに関して 𝑧_1と𝑧_2は無相関だし、𝑧_1の分散は1
ちなみに、無相関と独立は別の概念で無相関は線形の関係(一方が増えるほどもう一方が増えるみたいな)がないことなのに対し、独立は非線形な関係を含めた概念である。
ここで、白色化された観測値\boldsymbol{z}=V·\boldsymbol{x} を分離された観測値\boldsymbol{y}にさらに変換する、つまり\boldsymbol{y}=U·\boldsymbol{z}を満たすような行列𝑼はすべての行ベクトルが直交する
※そうでないとせっかく\boldsymbol{z}を無相関にしたのに\boldsymbol{y}が相関を持つようになってしまう。
Uのある行だけを定数倍しても\boldsymbol{y}が満たすべき条件は満たされるので、 𝑼は直交行列だけを考えても存在が保証される。
非ガウス性の評価
独立な確率分布同士を足すと、ガウス分布に近づいていく。逆に言えば、ガウス分布から離れるほど独立に近づくということである。 この「ガウス分布にどれだけ近いか」を表す指標には様々あるが、ここで尖度を計算する。
平均0の確率変数Xの確率密度関数をp(X)とすると尖度は
Kurt(X)=E[X^4]/E[X^2]^2-3=\int{x^4p(x)}dx/(\int{x^2p(x)}dx)^2-3
で計算される。 正規分布の尖度は0であり、それと比べて尖度が大きければ鋭いピークと長く太い裾をもった分布であり、 尖度が小さければより丸みがかったピークと短く細い尾をもつ分布である。
\boldsymbol{y}は \boldsymbol{y}=U·\boldsymbol{z}を満たすため、各成分が平均0、分散1かつ、仮定より独立である。 U=(\boldsymbol{u_1},\boldsymbol{u_2}, \dots \boldsymbol{u_n})^Tとすると、 \boldsymbol{u_i}と \boldsymbol{u_j} を入れ替えてもUの条件を満たす。
したがってKurt(y_1)\geq Kurt(y_2)\geq \dots \geq Kurt(y_n)を仮定してよい。
よって\boldsymbol{\alpha}をノルム1のベクトルとすると
Kurt(\boldsymbol{\alpha}^T\boldsymbol{y})
=E[(\boldsymbol{\alpha}^T\boldsymbol{y})^4]/E[(\boldsymbol{\alpha}^T\boldsymbol{y})^2]-3
=E[(\sum_{i=1}^{n}\alpha_i y_i)^4]/(E[(\sum_{i=1}^{n}\alpha_i y_i)^2])^2-3
=\sum_{i=1}^{n} E[\alpha_i^4 y_i^4]/(\sum_{i=1}^{n} E[\alpha_i^2 y_i^2])^2-3 (\because \boldsymbol{y}が独立)
=\sum_{i=1}^{n} ({\alpha_i}^4 E[y_i^4])-3
\leq\sum_{i=1}^{n} ({\alpha_i}^2 E[y_i^4])-3 (\because E[y_i^4]>0)
=\sum_{i=1}^{n} {\alpha_i}^2 Kurt(y_i)
\leq Kurt(y_1)
が成り立ち、 等号成立は |\alpha_1|=1の時である。
よって\boldsymbol{\beta}をノルム1のベクトルとすると Kurt(\boldsymbol{\beta}^T\boldsymbol{z})=Kurt(\boldsymbol{\beta}^TU^{-1}\boldsymbol{y})\leq Kurt(y_1)である。 等号成立は\boldsymbol{\beta}^T\boldsymbol{z}=\pm y_1の時である。 すなわちKurt(\boldsymbol{\beta}^T\boldsymbol{z})を最大化するように \boldsymbol{\beta}を選べば等号成立条件を満たしy_1が求まることになる。 このときの\boldsymbol{\beta}を \boldsymbol{u_1}とする。
\boldsymbol{\alpha}をノルム1、 \alpha_1=0, \dots,\alpha_{i-1}=0を満たすベクトルとすると
Kurt(\boldsymbol{\alpha}^T\boldsymbol{y})\leq Kurt(y_i) が成り立つ。
このことからy_2以降の計算は、Uが直交行列である、つまり \boldsymbol{\beta}とこれまでに決定した \boldsymbol{u_i}のすべてが直交するという条件つきで Kurt(\boldsymbol{\beta}^T\boldsymbol{z})を最大化するような \betaを次々に選び \boldsymbol{u_2}以降を決定すればよい。
ニュートン法によるICAの計算(Fast ICA)
問題はKurt(\boldsymbol{u_i}^T\boldsymbol{z})を最大化するような\boldsymbol{\boldsymbol{u_i}}を計算することである。 𝑖の小さい順にこれを求める。
\boldsymbol{\boldsymbol{u_i}}の満たす条件は、 ノルムが1かつ\boldsymbol{\boldsymbol{u_j}} (𝑗=1,2,\dots 𝑖−1)と直交
G(u)=u^4-3とおくと、E[G(\boldsymbol{u_i}^T\boldsymbol{z})]の最大化をすればよい。 つまりこの関数の極値を与えるu_iを求めたい
まず\boldsymbol{\boldsymbol{u_i}}の初期値を何か決めておいて
① y_i=\boldsymbol{u_i}・\boldsymbol{z}
② \boldsymbol{u_i}\leftarrow \boldsymbol{u_i}-E[G'(y_i)\boldsymbol{z}]/ E[G''(y_i)] ニュートン法
③ \boldsymbol{u_i}\leftarrow \boldsymbol{u_i}-\sum_{j=1}^{i-1}(\boldsymbol{u_j}^T\boldsymbol{u_i})\boldsymbol{u_j} グラムシュミットの直交化
④ \boldsymbol{u_i}\leftarrow \boldsymbol{u_i}/||\boldsymbol{u_i}|| ノルムの正規化
①から④を収束するまで繰り返す。 実際に計算する場合には、n個の観測値からなるデータをt回分得られたとすると、②で出てくる期待値はこれらt回分の平均値ということになる。
その他の評価関数
尖度を使って説明したが、尖度は外れ値に弱いという特徴がある。 実用的な評価関数としては
G(u)=\log(\cos(a_1h))/a_1, G'(u)=\tanh(u), G''(u)=1/\cosh^{2}(u)
G(u)=-\exp(-a_2u^2/2)/a_2, G'(u)=u\exp(u^2/2), G''(u)=u^2\exp(u^2/2)
といった関数が用いられる。