2013-11-07

multicore (R)

Rのライブラリ multicore を試した結果です．

mclappy()という関数を使いました．

この関数は lapply() を並列でできるようにしたもので，
独立に繰り返し計算する必要がある場合（乱数や初期値を変えた結果が欲しいときなど）
に簡単に並列で計算できます．

例えば k-means 法を R で行う場合に kmeans() という関数があります．
この関数は

kmeans(x, centers, iter.max = 10, nstart = 1,
       algorithm = c("Hartigan-Wong", "Lloyd", "Forgy",
                     "MacQueen"))

という形になっていて，x としてクラスタリングしたいデータ，centers に
クラスタ数をとります．このときクラスターの初期値を nstart 回ランダムに選び，
その結果もっともクラスター内平方和（tot.withinss）を小さくした結果を返します．

k-meansの各初期値からの反復は，ほかの初期値を選んだ場合の結果に依存しないので
簡単に並列化できると考えられます．

そこで並列で計算させた場合にどれくらい早くなるか実験しました．
まず mclapply() は

mclapply(X, FUN, ..., mc.preschedule = TRUE, mc.set.seed = TRUE,
         mc.silent = FALSE, mc.cores = 1L,
         mc.cleanup = TRUE, mc.allow.recursive = TRUE)

となっています．

lapply() と同じで，X の要素に関数 FUN を作用させた結果をリストで返します．
mc.cores が（最大限）使うコアの数です．

これを踏まえ，ライブラリMASSにあるデータ Boston を4つのクラスタに k-means 法を使ってクラスタリングしました．
初期値は 5000*n 個 (n=1,2,...,10) 与え，並列したときとかかった時間の差を調べました．

並列させないときは

t <- proc.time()
invisible(kmeans(Boston, 4, nstart=5000*n))
duration <- (proc.time() - t)[3]

で時間を計り，並列させたときは

t <- proc.time()
results <- mclapply(rep(5000*n/8, 8), function(nstart) kmeans(Boston, 4, nstart=nstart))
#collect results
objective.values <- sapply(results, function(result) result$tot.withinss)
result.best <- results[[which.min(objective.values)]]
duration <- (proc.time() - t)[3]

で時間を計りました．
#collect results の後の二行はクラスター内平方和（tot.withinss）を小さくした結果を返すための部分です．

mclapply() は mc.cores を指定しないと，自動で使えるコアの数を調べます．
自分が実験した環境だとコアの数が8だったので，5000*n個の初期値を
8個に分割しています．分割数に差があるかどうかを調べるために
9個に分割した場合も調べました:

t <- proc.time()
results <- mclapply(rep(5000*n/9, 9), function(nstart) kmeans(Boston, 4, nstart=nstart))
objective.values <- sapply(results, function(result) result$tot.withinss)
result.best <- results[[which.min(objective.values)]]
duration <- (proc.time() - t)[3]

細かいですが kmeans() は nstart の値が小数でもエラーがでず，nstartの小数以下を切り捨てた回数だけ初期値をランダムに発生させて計算しています（1:4.5と1:4が同じため）．

結果としては次のようになりました:

f:id:mkprob:20131106234450j:plain

横軸は初期値の数(5000*n)，縦軸はかかった時間です．
凡例の cores の値は使ったコアの数で，list の値は 5000*n 個の初期値の分割数です．

コア数が多い方が当然早く計算ができることが分かりました．
また初期値の分割数はコア数の倍数にした方がいいと思われます．
実際 multicore のマニュアルを読むと

"By default (mc.preschedule=TRUE) the input vector/list X is split into as many parts as there are cores (currently the values are spread across the cores sequentially, i.e. first value to core 1, second to core 2, ... (core + 1)-th value to core 1 etc.) and then one process is spawned to each core and the results are collected."

とあります．そこでコア数の倍数である16, 24分割した場合は8分割の場合とほとんど同じ結果になりました（凡例の一番上はcores=8の間違いです）:

f:id:mkprob:20131106234907j:plain

以上をまとめると

・並列計算したほうが早く，multicoreのmclapply()で簡単にできる
・初期値や乱数の分割はコア数の倍数にした方が良い（一つ一つの計算量が同じなら）

snowとかも調べてみたいと思います．

2013-07-15

Contiguity

van der Vaart

統計にでてくるContiguityの概念についての問題(van der Vaart: Asymptotic Statistics. Chap 6.)を解いたので，そのまとめです．

問題

(1) を標準正規分布，を平均分散1の正規分布とする．このときとが互いにcontiguousであることと，は同値であることを示せ．

(2) total variation を

とすると　

ならばとは互いにcontiguousである．

(3) 任意のに対して，とは互いにcontiguousであるが，となる例を見つけよ．

(4) はに対してcontiguousであるが，はに対してcontiguousでない例を挙げよ．

(5) 二つの確率測度, が互いに絶対連続であることと，確率測度の列, が互いにcontiguousであることは同値である．

証明

(1) 互いにcontiguousなら，が有界であることを対偶を使って示す．

とする．このときの部分列でまたはとなるものが存在する．この部分列をとし，

とすると，

より互いにcontiguousでない．よって互いにcontiguousなら，が有界である．

逆を示す．が有界であるとする．とする．が有界であるから，任意のに対して，あるがあって，

となる（tight）．このとき

となるから，

となる．は任意であったから

よってはに対してcontiguous．同様にはに対してcontiguousであることが分かる．よって同値性が示された．

(2) とする．total variationが0に収束するという仮定より，任意のに対して，あるがあって，

となる．よってに対して，

．

よってはに対してcontiguousであることが示された．条件の対称性から同様にして，互いにcontiguousであることが分かる．

(3) をに台を持つディラック測度とする．

とすれば互いにcontiguousであるが，total variation は集合{0}のときにとなる．

(4)

とすれば条件を満たす．

(5) 二つの確率測度, が互いに絶対連続であるとする．このときラドン--ニコディム微分

はそれぞれの確率測度の下で可積分である．

であるから，

よってはに対してcontiguousである．全く同様にしてはに対してcontiguousである．

次に確率測度の列, が互いにcontiguousであるとする．とする．とすれば，contiguityより．よってQはPに関して絶対連続．同様にPはQに関して絶対連続．[証明終]

2013-05-10

Ostrowski-Taussky Inequality

行列

Ostrowski-Taussky Inequality
Aを正方行列とし，が正定値行列であるとする．このとき以下の不等式が成立する:

証明
は正定値対称行列であるから，ある直交行列Pがあって

となる．とすれば (P:直交)であり，

である．これよりBの(i,j)成分は

となっている．このような行列に関して次の不等式が成立する(下の補第参照):

これが成り立てば

これより証明が終わる．あとは次の補題1を示せばよい．ただし帰納法の成立のため主張を拡張する．

補第1

証明
Bの次元に関する帰納法で示す．n=1のときは明らかに成立する．n=k-1のとき成立するとする．まず次のように行列Bを分割する．

ここでCは(k-1)×(k-1)行列である．まず任意のに対して，として，

となることがまずわかる．

次に，よりブロック行列の行列式を考えれば

帰納法の仮定からCは正則であるから，Sherman-Morrison-Woodburyの公式より

ここでであることがわかる．実際帰納法の仮定より

Cは正則だから任意のxに対してあるyがあって

よって任意のxに対して

となるためである．これより

[証明終]

主張が言いたいことを考えると，これは複素数zに対する不等式

の行列への拡張ではないかという指摘を受けました．というのも

と分けます．右辺の第一項はAを対称化したもので，第二項は反対称化したものです．実対称行列の固有値は実数で，反対称行列の固有値は純虚数です．つまり固有値の視点から見ると，上の分解は行列の実部と虚部への分解に対応するものではないかと考えられます．そうすると主張の左辺はAの固有値の実部の積，右辺は固有値の絶対値の積なので，不等号が成り立つと思われます．これを正当化するには，右辺の第一項と第二項が同時三角化可能で，行列の和の固有値が，各行列の固有値の和になることを言わないとダメだと思いますが...

ちなみに主張の右辺の絶対値は要らないことが証明から分かります．

2013-04-27

ゼータ分布に従う独立な確率変数が互いに素になる確率[D.Williams]

確率

主張

X,Yは独立にゼータ分布に従うものとする．つまり

とする．このとき

(gcdは最大公約数)とすると，

となる．[D.Williams: Probability with martingales, p226]

証明

とすると，達は独立．何故ならば

であるため．これより

となる．

次につまりX,Yが互いに疎になる確率を計算する．，とすれば達は独立である．実際

よりわかる．従って

\begin{eqnarray*}
P(H=1)
&=&P\bigg(\bigcap_{p:{\rm prime}}(E_p\cap F_p)^c\bigg)\\
&=&\prod_{p:{\rm prime}} P((E_p\cap F_p)^c)\\
&=&\prod_{p:{\rm prime}} (1-1/p^{2s})\\
&=&\frac{1}{\zeta(2s)}
\end{eqnarray*}

Xがpの倍数であるときのXの分布は，またゼータ分布になる:

nが素数でないときも

と定義すれば，条件付き期待値の性質から

となって示された．（包除原理を使って計算する方が正確だと思われます）[証明終]

ちなみにgcd(X,Y)の可測性は

から分かります．また結局互いに素な確率は

ということもわかり，面白いと思いました．

2013-02-20

Hilbert--Schmidt作用素はコンパクト作用素

Kuo

主張

を可分なヒルベルト空間とし，をからへのHilbert-Schmidt作用素とする．このときはコンパクト作用素である．

証明

事実1

をコンパクト作用素とし，

であるならば，はコンパクトである(コンパクト作用素の全体はバナッハ空間なので)．

事実2

有界線形作用素の値域が有限次元空間であればコンパクトである(有限次元と局所コンパクトは同値)．

この二つの事実を用います．

をの内積とし，を完全正規直交系とする．有界線形作用素を

と定義する．事実2よりはコンパクト．がに強収束することを示せば，事実1よりはコンパクトとなり，主張が示されたことになる．

は完全正規直交系であるから，

したがって任意のに対して

となる．はHilbert-Schmidt作用素であるから，

以上よりはに強収束する．[証明終]

Hui-Hsiung Kuoの Gaussian Measures in Banach Spacesという本を読み始めました．上の主張はこの本のExerciseになっています．抽象ウィナー空間が分かるようになりたいです．

Gaussian Measures in Banach Spaces
Hui-Hsiung Kuo

2006-06-29
売り上げランキング : 1234523

Amazonで詳しく見る by G-Tools

2012-11-14

Problem5.4.4

KaratzasAndShreve

Problem5.4.4
連続で適合している確率過程がd次元のブラウン運動であることの必要十分条件は

が任意のに対して連続な局所マルチンゲールになることである．ただし

とする．

証明
がd次元のブラウン運動であれば，が連続な局所マルチンゲールであることは伊藤の公式から明らか(cf. Karatzas and Shreve, Proposition5.4.2, p312)．逆に任意のに対してが連続な局所マルチンゲールになるとする．とくにについて考えれば，は連続な適合した局所マルチンゲールとなる．さらにであることが，Karatzas and Shreve, Proposition5.4.2から分かる．よってLevy's characterization of Brownian motion よりこれらがd次元のブラウン運動であることがわかる．[証明終]