2021-01-28

グラフラプラシアンとPCAと固有値分解

点集合の座標が行列Xとして与えられているとき、そのペアワイズ内積行列 $P =X^TX$ は対称行列であり、それを固有値分解し
回転行列Vと非負固有値を対角成分とする対角行列 $\Sigma$ とを使って
$P = V \Sigma V^T$ と分解できる
これは、点をVで回転して $P = X_{new}^T X_{new}$ となるような $X_{new} = \sqrt{\Sigma} V^T$ なる、新しい座標を点に与え、その配置が、分散の観点で整然とした状態になるようにすること、と説明される。
一方、単純無向グラフには、ラプラシアン行列Lなるものがさだまる
Lはその対角成分がグラフ頂点の次数であるような行列Dと、グラフの隣接行列Aとを使って $L = D - A$ として与えられる
Lも対称行列であり、固有値分解すると、回転行列と固有値を対角成分とする対角行列の積に分解できて、行列Xの場合と同様に、グラフの各頂点に座標を持たせたとみなすことができる
いわゆるPCAでは、大きな固有値に対する固有ベクトルが与える座標軸は、点の存在状態をおおまかにとらえる力が大きい
他方、グラフラプラシアンの場合には、小さな固有値に対する固有ベクトルが与える座標軸が、グラフをおおまかに分ける力が大きい
ちなみに、連結グラフのラプラシアンの固有値の１つは０で、それ以外はすべて正であることが知られている
PCAの固有値の働きとグラフラプラシアンの固有値の働きが、大小に関して逆転していることを解消することにする。それによりグラフにとっても何かしらの対称行列があり、その固有値分解が、頂点に座標のようなものを与え、値の大きい固有値に対応する軸がそのグラフ頂点のばらつきを大まかに説明するような、そんな構成にしたい
グラフラプラシアンを離れて、正定値対称正方行列　Q　を考える
Qには逆行列が存在するとして、それを $Q^{-1}$ とする
今、Qを固有値分解すると、 $Q = V \Sigma V^T$ となるとする。ここでVは回転行列であり、 $\Sigma$ は対角行列ですべての対角成分が正である
Vは回転行列であるから $V^T = V^{-1}$ である
実は、 $Q^{-1} = V \Gamma V^{-1}$ となる。ただし、Vは回転行列でQを固有値分解したものと同じであり、 $\Gamma$ は対角行列であり、その値はすべて正で、かつ、 $\Sigma$ の大きい順にi番目の固有値と、 $\Gamma$ の小さい順にi番目の固有値との積は、i=1,2,....のすべてのiについて１になる

# 適当に正定値対称行列を作る
library(GPArotation)
n <- 5
# 回転行列
V <- Random.Start(n)
# 正固有値となるように正成分の対角行列
Sigma <- diag(runif(n)*5)
# 正定値対称行列
Q <- V %*% Sigma %*% t(V)
# その逆行列
Qinv <- solve(Q)
# 固有値分解する
eout.Q <- eigen(Q)
eout.Qinv <- eigen(Qinv)
# 固有値のセットを昇順・降順を入れ替えて掛け合わせると、すべて１
eout.Q[[1]] * eout.Qinv[[1]][n:1]
# 固有値ベクトルは、符号の入れ替えが起きるので、すべての固有ベクトルが、一致しているか、符号違いで一致しているかを以下の式で確かめる
# 全成分がゼロになるので、固有値ベクトルが同じになっていることがわかる
(eout.Q[[2]] - eout.Qinv[[2]][,n:1]) * (eout.Q[[2]] + eout.Qinv[[2]][,n:1])

つまり、正定値対称行列の固有ベクトルと、その行列の逆行列の固有ベクトルは等しいことがわかる
ただし、片方で大きい固有値に対応する固有ベクトルが、逆行列では小さい固有値に対応する固有ベクトルである、という対応関係にある
したがって、ラプラシアンの0でない最小固有値に対する固有ベクトルは、ラプラシアンの逆行列の、2番目に大きな固有値に対応する固有ベクトルに一致する
ちなみに、連結グラフのラプラシアンは１個の固有値を持つ。それに対応する「ラプラシアンの逆行列」の固有値は無限大となる
このことから、実際には、ラプラシアンの逆行列を計算することはできないが、ラプラシアンの逆行列の大きな固有値に対応する固有ベクトルを求める代わりに、ラプラシアンの小さな固有値に対応する固有ベクトルを求め、その固有ベクトル軸が重要である、と言っている
さらにちなみに、ラプラシアンの固有値ゼロに対応する固有ベクトルというのは、すべての成分の値が等しいベクトルのことである。逆に言うと、ラプラシアンの逆行列が取り扱えるとして、最大固有値に対応する固有ベクトルも、すべての成分の値が等しいベクトルになる、と言っている
今、PCAでは点をうまくばらつかせることを目指しており、ラプラシアンの固有値分解で得る固有ベクトルにもその役割を持たせようとしているわけなので、すべての頂点の座標が同じになるような固有値に興味はないから、ラプラシアンの逆行列の2番目に大きな固有値＝ラプラシアンの2番目に小さな固有値が、頂点座標を決めるうえで最もおおまかな情報を持っていることになる
グラフラプラシアンの逆行列は「無限大の固有値」を持たせる必要があるので、作れないが、「十分大きな固有値」を持たせることで近似することにすれば、以下のようになる

# 適当に単純無向グラフを作り
library(igraph)
n <- 7
A <- matrix(sample(0:1,n^2,replace=TRUE,prob=c(0.7,0.3)),n,n)
A <- A + t(A)
diag(A) <- 0
A <- sign(A)
g <- graph.adjacency(A,mode="undirected")
plot(g)
# ラプラシアン行列を作る
L <- diag(degree(g)) - A
# ラプラシアン行列を固有値分解する
eout <- eigen(L)
# 固有値ベクトルの束の行列と固有値を対角成分とする対角行列とでラプラシアン行列を再計算してみる
eout[[2]] %*% diag(eout[[1]]) %*% t(eout[[2]]) - L # essentially zero
# ラプラシアン行列の逆行列は０なる固有値があるのでうまく行かない
#Linv <- solve(L)
# 0の固有値に対応させて大きな正の値を、非0の固有値に対応させて、その逆数を固有値とし
# 固有値ベクトルはラプラシアンのそれをそのまま用いて
# ラプラシアンの逆行列の近似行列を作る
Linv. <- (eout[[2]]) %*% diag(c(1/eout[[1]][1:(n-1)],10^8)) %*% t(eout[[2]])

# ラプラシアンの逆行列を固有値分解してやる
einvout <- eigen(Linv.)
# ラプラシアンの固有値と
# ラプラシアンの逆行列の固有値とは相互に逆数になっている
eout[[1]] * einvout[[1]][n:1]
# ラプラシアンの固有値ベクトルと
# ラプラシアンの逆行列の固有値ベクトルとは
# 符号の入れ替わりはあるが
# それ以外は同一であることを確認する
eout[[2]] - einvout[[2]][,n:1]
eout[[2]] + einvout[[2]][,n:1]

(eout[[2]] - einvout[[2]][,n:1]) * (eout[[2]] + einvout[[2]][,n:1])