情報科学と人工知能のノート

初等的な知識から最新論文の解説まで色々集めていきます.備忘録兼用.

有向グラフの強連結成分分解

ナイーブなアルゴリズム

ある頂点  v からグラフ探索(深さ優先探索など.なんでもよい)を始めて訪れることができた頂点の集合  V_1 と,グラフの全ての辺を逆向きにしてから  v からグラフ探索を始めて訪れることのができた頂点の集合  V_2 に対し, V_1 \cap V_2 によって誘導される部分グラフが  v を含む強連結成分です.
よってせいぜい  O(|V|(|V| + |E|)) くらいの時間計算量のアルゴリズムで全ての強連結成分を調べ上げることができます.
これを線形時間で達成するのが以降の話です.

アルゴリズム

強連結成分分解は深さ優先探索を合計で二回分行うことによって達成できます.

  1. トポロジカルソートと全く同じ方法で頂点に番号を付ける.
  2. グラフの全ての辺を逆向きにする.
  3. 番号の最も大きい頂点から深さ優先探索を行う.訪れた頂点が一つの強連結成分.
  4. 得られた強連結成分をグラフから取り除き,まだグラフが空でなければ 3 に戻る.

ステップ 1 が深さ優先探索一回分,ステップ 3,4 が合計で一回分です.

アルゴリズムの正しさ

頂点に番号を付ける以外に特にフラグをたてたりしていないのですが,この方法できちんと全ての強連結成分を調べ上げることができます.

まず,ステップ 1 ではグラフに対してトポロジカルソートのような番号付けを行っています.もし入力グラフが DAG であれば,当然トポロジカルソートによる順序付けそのものが各頂点の番号になります.また入力グラフが有向サイクルを含むグラフであったとしても「グラフの各強連結成分を一頂点に縮約し,それぞれの強連結成分の頂点の番号の中で最も大きな値を縮約後の頂点の番号とする」という処理を行うと,縮約後の番号は縮約後の DAG のトポロジカルソートになっています.

すなわち番号の最も大きな頂点を含む強連結成分に入る辺は存在しないので,グラフの辺を逆向きにすると,この強連結成分から出る辺が存在しないことになります.よって,最も大きな番号を持つ頂点からグラフ探索を行うと,この頂点を含む強連結成分の各頂点だけをかき集めることができます.
この強連結成分を取り除いた後も同様に,縮約後の DAG のトポロジカルソートで次に来る強連結成分に含まれる頂点が次に選ばれますので,以下同様に正しく強連結成分のみを集めることができます.

主成分分析ってたぶんこういうもの

例によってさも分かってるかのように説明していますので間違ってたら教えてください。

たたき台

以下のように四次元の点が 10 個並んだ表があったとします。

> round(a, 3)
[,1] [,2] [,3] [,4]
[1,] -0.043 0.026 -0.005 -0.013
[2,] -1.314 0.785 -0.150 -0.414
[3,] 1.757 -1.049 0.201 0.554
[4,] -0.119 0.071 -0.014 -0.037
[5,] -0.963 0.575 -0.110 -0.303
[6,] 1.059 -0.633 0.121 0.334
[7,] 1.248 -0.745 0.143 0.393
[8,] -1.608 0.961 -0.184 -0.507
[9,] -1.149 0.686 -0.132 -0.362
[10,] 1.132 -0.676 0.130 0.357

これだけ見てもどういうデータなのかはよくわかりません。
しかし、これらのデータは実は以下のようにして作られたものだったのです。

> d1 <- rnorm(4) # 適当に方向ベクトルを定める
> d1 <- d1 / sqrt(d1 %*% d1) # 後の都合のために単位ベクトルにしておく
> c1 <- rnorm(10) # 適当に 10 個の倍率を定める
> c1 <- c1 - mean(c1) # 後の都合のために平均を 0 にしておく
> a <- t(mapply(function(x) return (d1 * x), c1)) # 方向を倍率倍

言葉で言うと、方向ベクトルを適当に定めて、適当に何倍かしただけです。
上記のデータがこのように作られたものであったならば、それぞれの点を読んでもよくわからない4次元の絶対座標で表すよりも、方向ベクトルの何倍なのかで表した方がデータを把握しやすくなりますよね。
このデータに対して主成分分析を行うと、今まさに欲しかった、各点が方向ベクトルの何倍なのかを求めることができます。
R で実行するのは簡単で、行列を princomp 関数に与えるだけです。

> res1 <- princomp(a)
> summary(res1)
Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4
Standard deviation 1.414748 2.980232e-08 2.634178e-09 0
Proportion of Variance 1.000000 4.437535e-16 3.466824e-18 0
Cumulative Proportion 1.000000 1.000000e+00 1.000000e+00 1

細かい説明は省きますが Proportion of Variance の左端が 1 になっていますよね。
これはおおざっぱに言うと、データが1次元に並んでいましたという意味です。
また、res1$score[,1] で各点が方向ベクトルの何倍なのかを得ることができます。

> res1$score[,1]
[1] -0.05176209 -1.59257228 2.12906814 -0.14387618 -1.16726426 1.28359077
[7] 1.51268231 -1.94956546 -1.39275710 1.37245617

c1 の値を確認してみると(向きは逆ですが)正しい答えが得られていることがわかります。

> c1
[1] 0.05176209 1.59257228 -2.12906814 0.14387618 1.16726426 -1.28359077
[7] -1.51268231 1.94956546 1.39275710 -1.37245617

これを一般化して、n 次元空間の点を m 本の方向ベクトルの線形和で生成したときも、点の座標の表を主成分分析すれば、各点について m 個の成分に情報を削減することができます。
ただし、方向ベクトルが単位ベクトルでなかったり、方向ベクトルが直交していなかったり、成分ベクトルの平均が0でない場合には、始めに与えた成分とは異なる m 次元のデータが得られます。

主成分分析

ここまでは説明のために方向ベクトルの線形和でデータを生成してきましたが、主成分分析は、データが実際にはどう作られたかはさておき、与えられた成分をうまいこと少数の方向ベクトルの線形和に変換するための手法です。

座標変換(といっていいのかわかりませんが)の中心はデータの標本平均です。
最初の一本目の方向は、その方向に全てのデータを射影したときに、その方向上での標本分散が最大になるものを選びます。
分散がなるべく大きい方が、その方向ではデータがばらついている=その方向でデータの差異を説明し易い、ということになります。
二本目は、一本目に直交する方向の中で、同じく全てのデータをその方向に射影したときに方向上での標本分散が最大になるものを選びます。
以下同様に分散を用いて必要な本数だけ方向を選んでいきます。
なお、i 番目の方向(の各係数)を第 i 主成分といい、データが主成分の線形和の中でその方向に何倍されているかを第 i 主成分得点といいます。

princomp の結果を summary に与えたときの各列の値は各方向に関する情報を表しています。
左から第 1 主成分、第 2 主成分、となっています。
Standard deviation はその方向にデータを射影したときの標準偏差です。
Proportion of Variance (寄与率)は全ての主成分に関する分散の和のうち、それぞれの主成分の分散の大きさがどれくらいかを表す指標です。
Cumulative Proportion (累積寄与率)は第 1 主成分からその主成分までの分散の和が全体に対してどれくらいの大きさかを表しています。
もし第 1 主成分の寄与率が 1 であれば、データは直線上に載っていることになります。
もし第 m 主成分の累積寄与率が 1 であれば、わざわざ元の n 次元の全ての成分を見なくても、m 次元の情報でデータを表せることになります。
実際に使用するときには、先頭の少ない主成分(第 1, 2, 3 くらい)の累積寄与率が十分高ければそれ以降は無視してしまいます。
よくある使い方としては、第 1 主成分得点と第 2 主成分得点を使って二次元散布図にプロットして可視化したりします。

ちなみに、主成分を求める問題はデータの分散共分散行列の固有ベクトルを求める問題と同等だそうです。

R と MATLAB の実行時間計測

乗り換えを検討するために R と MATLAB で単純な処理の実行時間を計測してみました。

なお、以下の実験ではそもそもソフトウェアを動かしているマシンが違います。
自宅PC の R と大学PC の MATLAB と自分にとってどっちが都合がいいかを計っただけです(どこかに書いておかないと忘れますので虚空に向かって載せてるだけです)。
よって、以下の結果から一般的にどっちがどうと言えるものではありませんのでご注意ください。
でも MATLAB 走らせてるマシンの方がハヤイはずなんですけどねー。

R

> times <- rep(0, times=20)
> for (i in 1:20) {
+ t <- proc.time()
+ for (j in 1:100)
+ diag(1000) * diag(1000)
+ times[i] <- as.vector((proc.time() - t)[3])
+ }
> times
[1] 2.607 2.596 2.644 2.625 2.677 2.677 2.603 2.612 2.659 2.659 2.603 2.608
[13] 2.598 2.598 2.595 2.594 2.596 2.619 2.628 2.679

MATLAB

>> for i = 1:20 tic(); for j = 1:100 eye(1000) * eye(1000); end; toc(); end
Elapsed time is 4.100130 seconds.
Elapsed time is 4.216998 seconds.
Elapsed time is 4.228933 seconds.
Elapsed time is 2.800470 seconds.
Elapsed time is 2.593125 seconds.
Elapsed time is 2.561670 seconds.
Elapsed time is 2.556924 seconds.
% 中略
Elapsed time is 2.626316 seconds.
Elapsed time is 2.570071 seconds.
Elapsed time is 2.614598 seconds.
Elapsed time is 2.581647 seconds.

とりあえずこれだけ見る分には R のままで良さそうです。
もっと複雑な処理だと値渡しが云々といったことなどが聞いてくるんでしょうかね。

余談ですが、昨日 Revolution R の導入を試してみたところ途中でこけてしまいました。
インストールできました。
ちなみに、うちの環境は Mac OS X 10.6 で、configure が FORTRAN がどうこうで止まっていたのですが、MacPorts から g95 を入れたら通りました。
で結果はこれ。

> times <- rep(0, times=20)
> for (i in 1:20) {
+ t <- proc.time()
+ for (j in 1:100)
+ diag(1000) * diag(1000)
+ times[i] <- as.vector((proc.time() - t)[3])
+ }
> times
[1] 6.400 6.277 6.290 6.416 6.627 6.608 6.619 6.571 6.499 6.580 6.633 6.623
[13] 6.599 6.476 6.450 6.481 6.501 6.463 6.623 6.591

あれ?

大学マシンに R が入ってたのでついでに。

> times <- rep(0, times=20)
> for (i in 1:20) {
+ t <- proc.time()
+ for (j in 1:100)
+ diag(1000) * diag(1000)
+ times[i] <- as.vector((proc.time() - t)[3])
+ }
> times
[1] 5.425 6.004 6.634 6.683 6.655 6.681 6.656 6.680 6.655 6.703 6.666 6.682
[13] 6.645 6.672 6.689 6.625 6.671 6.644 6.672 6.642

んんんん??
計り間違いではないはずなのですが。

更に別の、自宅のより良い Mac Book Pro。

> times <- rep(0, times=20)
> for (i in 1:20) {

  1. t <- proc.time()
  2. for (j in 1:100)
  3. diag(1000) * diag(1000)
  4. times[i] <- as.vector((proc.time() - t)[3])
  5. }

> times
[1] 0.949 0.849 0.848 0.864 0.836 0.808 0.823 0.821 0.821 0.821 0.808 0.822
[13] 0.826 0.821 0.822 0.810 0.821 0.821 0.819 0.822

R 備忘録 2013/05/11

今日は短めです。

ベクトルの要素のうち条件を満たすもののインデクスを集める

ということが which 関数を使うと実現できます。

> a <- seq(1, 30, 3)
> a
[1] 1 4 7 10 13 16 19 22 25 28
> which(a %% 2 == 0) # 偶数
[1] 2 4 6 8 10

> which(max(a) == a)
[1] 10

この which(max(a) == a) が初見時にとても不思議な文に見えたのですが、which の中身だけを評価してみると仕組みがよくわかりました。

> max(a) == a
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE

要するに、TRUE のインデクスを集めてきて返すだけみたいですね。

なお、行列に対して使用しても

> b <- matrix(1:8, 2, 4)
> b
[,1] [,2] [,3] [,4]
[1,] 1 3 5 7
[2,] 2 4 6 8
> b %% 2 == 0
[,1] [,2] [,3] [,4]
[1,] FALSE FALSE FALSE FALSE
[2,] TRUE TRUE TRUE TRUE
> which(b %% 2 == 0)
[1] 2 4 6 8

となり、ベクトルとしてのインデクスを返すようです。

要素のソート・ランキング

sort 関数を使うとベクトルの要素を(デフォルトでは)昇順に並べ替えることができます。
また、order 関数を使うとベクトルの要素が昇順となるインデクスの順序を求めることができます。

> c <- runif(10)
> c
[1] 0.82190691 0.98909512 0.05540597 0.23543929 0.63394228 0.87780395
[7] 0.55227800 0.20187754 0.05960969 0.98804746
> sort(c)
[1] 0.05540597 0.05960969 0.20187754 0.23543929 0.55227800 0.63394228
[7] 0.82190691 0.87780395 0.98804746 0.98909512
> order(c)
[1] 3 9 8 4 7 5 1 6 10 2

なお、降順にするには decreasing を TRUE にします。

> sort(c, decreasing=TRUE)
[1] 0.98909512 0.98804746 0.87780395 0.82190691 0.63394228 0.55227800
[7] 0.23543929 0.20187754 0.05960969 0.05540597
> order(c, decreasing=TRUE)
[1] 2 10 6 1 5 7 4 8 9 3