Diary December, 2007

 

02日Sun.

気付いたらいつの間にか12月。師走。

ようやく、修論に筆が入り始めた。 毎日毎日少しずやっつけていく。

今日は調子よく作業が進んだので、 調子にのってk=9\to N_{{\rm pix}} = 3{,}145{,}728でCMBのサンプリングが出来るか実験をした(時間を概算せずに)。 案の定、計算がまったく終わらない。 今ここで計算してみる。 k=6\to N_{{\rm pix}} = 49{,}152のとき10 chain計算するのに約43秒かかるので、

t=43\times 2^{3(9-6)}\sim 6.12 {\rm h}

いくらまっても終わらない訳だ。

やはりfull-rezolutionで計算するのは、現時点では無理がある。

2000CPUだったとしたら、t=43\times 2^{3(9-6)}/500\sim 40{\rm s}で、今の環境のk=6と変わらない。 まったく問題外な時間ではない。 実際にO(N^{3/2})になっているかテストしておくべきであろう。

mac.jpg

03日Mon.

samplingの手法を使えば、template-freeで、 且つ、multi-frequencyなCMB観測の利点を享受することが出来る。 実用されているILC(Internal Linear Combination, 一、二年前に初めてこの手の論文を読んだときは、「Inter」だった気がするのだけれども、勘違いだろうか? 先生は後者で呼ぶ癖がついてしまったようだ。) はfoergroundの差し引き自体はやっていないので、 その精度に限界があり(foregroundを「打ち消している」という感じ)、 また、解像度を最も低いものに合わせなければならない、という制限がある。

その為、WMAPではILCはlarge-scaleでしか使われていない。

samplingの手法であれば、foregorundの推定も行え、解像度は最も高いものまで使うことが出来る。

最低解像度で解析する場合、 先ずはじめに、各周波数の観測結果を最も悪いビームパターンでスムージングすることで、 解像度を合わせてしまう。そうすれば、すべての周波数帯で共通のwindowがかかることになる:

データ回帰の項を実空間で書いた場合、

P( s|d_\nu, C_l)\propto   \exp\left[   -  \sum_\nu (d_\nu-A_\nu s)^tN_\nu^{-1}(d_\nu-A_\nu s)     \right]      \exp\left[ -s^tS^{-1}s\right]

と書ける。実空間で解析する場合、sをサンプルし、A_\nuを畳み込んでデータ回帰を計算するのは効率が悪い。 そこで、A_\nu sをサンプリングし、C_lのサンプリングを行う際にwindowを考慮した(windowで割る; e^{-\frac{1}{2} \sigma_l/(W_l C_l)})サンプルを行えばいい。

後者、つまり最高解像度まで計算する場合、このやり方は使えない。 なぜならそれぞれの周波数で「独立」にCMBを計算することになりcosmic varianceを超えてしまう。 e^{-\frac{1}{2} \sigma_l^\nu/(W_l^\nu C_l)}からC_lを独立にサンプリングしてしまうことになる。

適切な方法は、やはりきちんとsを唯一つだけ推定する方法である。 sをpixelで推定し、


今日はここまで・・・・

http://www.youtube.com/watch?v=HPkAfcw3akg

ヤムチャを名言やことわざの一部にいれると弱々しくなる

ヤムチャも鳴かずば撃たれまい
ヤムチャは登場しなければ死んでいなかった

その通り。

snow01.jpg

04日Tue.

続き
実空間でCMBをサンプリングを行いながら、 周波数ごとに異なる分解能を適切に扱うには、やはりフーリエスペースで扱う必要がある。 今、実空間でのCMBをtとすると、フーリエスペースでは、

a_{lm} = \int d\hat{n}\, Y_{lm}^*(\hat{n}) t(\hat{n})

と書ける。逆変換は、

t(\hat{n})  = \sum_{l,m} Y_{lm}(\hat{n}) a_{lm}

である。

実際観測される量sは、ビームパターンが畳み込まれているので、

s_p  =   \int d\hat{n} B_p(\hat{n})  t(\hat{n})   = \sum_{lm} a_{lm} \int d\hat{n}\, B_p(\hat{n}) Y_{lm}(\hat{n})

となる、今フーリエスペースでのビームパターン

B_{lm}(p) = \int d\hat{n}\, B_p(\hat{n}) Y_{lm}(\hat{n})

を考えると、上式は

s_p = \sum_{l,m} a_{lm} B_{lm}(p)

となる。

実際のビームパターンは複雑な形をしているのが一般的だが、 ここでは空間のスケールlだけによる様なビームパタンを仮定できる(分解能ぎりぎりまで解析能力を期待しない・・・)。 また空のどこを見ても同じビームパターンであると仮定できる。 すると、

B_{lm}(p) = B_{l}Y_{lm}(p)

と書ける。このときB_lはmでのアンサンプル平均から、

B_l = \frac{4\pi}{2l+1} \sum_{m=-l}^l |B_{lm}|^2

と求めることが出来る。

ここでは証明は省くが、ピクセル化(スムージングがかかった)された、データから得られるClと、 実際にClとの間には、

C_l^{\rm pix}   = B_l^2 C_l

の関係が成り立つ。

結局、実空間での観測量は、フーリエスペースでの量を使い、

s(p)   = \sum_{l,m}  Y_{lm}(p)   ( a_{lm} B_l )

と書ける。

Healpix C++ではフーリエスペースで指定したFWHMのガウシアンスムージングを行う関数

template<typename T> void smooth_with_Gauss(Alm<xcomplex<T> > &alm, double fwhm_arcmin);

が存在するため上記の操作は簡単に実装できる。

  1. 実空間でCMBを発生させる。
  2. 球面調和関数でフーリエ変換。
  3. それにガウシアンスムージングを「掛ける(積)」。
  4. 逆変換で実空間に戻す。
  5. データから計算結果をひき、全ピクセルで足し合わせchi-squareを計算する。
Healpix_Map<double> t, s ;
t.Set(order,RING) ;
s.Set(order,RING) ;

arr<double> weight_T ;
weight_T.alloc(2.0*nside) ;
weight_T.fill(1) ;

Alm<xcomplex<double> > alm(lmax,lmax) ;

map2alm_iter( t, alm, i_num, weight_T) ;
smooth_with_Gauss( alm, fwhm_arcmin) ;
alm2map( alm, s) ;

double chisq = 0.0 ;
for( int i=0; i<npix; i++)
 chisq += pow( (d[i]-s[i])/n[i], 2.0)

実空間のCMBをサンプリングしているので、 Clのサンプリングはピクセル化されたClからのサンプリングになり、 これは周波数ごとに共通にすることが可能である。

それぞれの周波数ごとのビームの効果は、s(p)を計算する際に現れている(本当?)。

7.jpg

05日Wed.

今日はAKARIの会議もほっぽっといて(・・・違いますよ、冗談です) 池内了先生の講義・飲み会に参加。 先生の書かれた本は何冊か読んでおります。 うちの指導教官の修士時代の指導教官。

科学と社会、科学者が果たすべき義務とか、 第一線の研究者が真顔でしゃべったら白眼視される(そんなこと考える暇があったら研究でもしてろ)ようなことについて色々お話された。 こんな風に書いているけど、自分はこういうことは邪険にしている人間ではない。 院生の自分がこんなことを書くと、

omae.jpg

となるが、「偉い人」が言うことには意味はある。 その人を「偉い人」と感じるかは人それぞれだが。

世間で「科学」がどのように扱われているかを知るたびに、悲しくなる。 馬鹿な連中には理解できない、と割り切ることも出来る。でも理解されたいとも思っている。

世間をにぎやかしている「○ピリチュアル」「○○イオン」「ゲルマニウム」 「あなたの前世は貴族(貴族は全人口の何パーセントなんだよ!?)」なんかを、 普通の人は本当に信じているのだろうか? ネタなんですよね? あの手の番組は「コメディー」ですよね?

現実を見れば、 院生の自分としては、自分を評価する方々が「業績を稼ぐ」ことが一番だと考えている以上、 こういったことに時間を割くより、蛸壺的に研究することが重要であることは必然である。 社会に対して説明責任を果たしている、果たすべきであるので実践している、 といくら学振の書類に書いても、一遍の論文の前にはそんなものは無意味である。 それが事実。

世の中そんなもの(学振がすべてだとは言わないが・・・)。

  • 院生は苦労しているんだからもっと援助しろとか、
  • 日本は科学技術立国とか言っている割には理系に冷たいとか、
  • そもそも政治を牛耳っている人に科学の専門の教育を受けた人がいない(菅直人、志位和夫は理系出身だが、博士はとっていないので)、

と不満を言っても、職業選択の自由はある訳で、

  • 理系の頭があれば文型理系どちらが得か分かるよね、

といわれれば元も子もない。

科学は文化であるで有るという主張は正しいと思うが、 それが正しいと思っている人が一般社会にどれだけいるのか、正直不安である。 優しい言葉をかければきれいな氷の結晶が出来る、 とかいっている人には、通じないだろうな。

科学をする理由は、 あまり難しいこと考えず、自分が楽しいからやっている、というのが単純で最もいいものだと思う。

と、酔っ払いは思うのであった。

http://jp.youtube.com/watch?v=Pp05Zg3HB8s

らき☆すた見終わったから、次は電脳コイル。 ぼくらのも見たい。ひぐらしもまだ見てないし。


08日Sat.

いつまでも、 「らき☆すた面白かった!」とか、 「電脳コイルは評判通り面白い」とかいっているアレなんで、 今日は修論完成までのロードマップを作成した。

取り敢えず、chapter,section構成は大まかに決め、 全6章+付録2章について筆を入れ始めた。 今まで溜めてきたTeXソースですぐに再利用可能なものを、 文の構成を無視し取り敢えず挿入。 全45頁の修論が完成しました。 これからどこまで増やせるか。

文章を書くのに苦は感じないのだが、 如何せん私の文章は読みにくく、意味分からんのが多いので、 無駄に頁を増やすことは避けます。

宇宙論グループの修士論文には力作もそろっていまして、 特に小松さん・荒井さんの修論は、 そのまま参考書・辞書(勿論研究もすばらしいですが)として使えるすばらしいものです。 小松さんは200頁越え(投稿論文もついてる!?)荒井さんは300頁越え(文字も小さい!?)です。

修士論文は研究だけではなく、 学部・修士前期までに自分が勉強してきたことの集大成でもある。 お二人の真似をすることはできないが、 将来後輩の役に立つよな修論を書きたいものである。

気合も入ったところで、明日も休日出勤だ。

book.jpg

そういえば、読もう読もうと思ってまだ手をつけてない論文があったな。 明日読むか。


10日Mon.

サンプリングの問題を考える。 ある周波数nuを完全にサンプリングするには、 その倍以上の周波数でサンプリングする必要がある。 これを「ナイキストのサンプリング定理」という。

一般に人間の可聴周波数領域は、 20〜20kHzと言われているので、 CDのサンプリング周波数はは44.1kHzである。 しかしながら上記の範囲はあくまでも「一般的な可聴領域」であるので、 個人によって違ったり、そもそも聞こえる音がすべての音とは限らないので、 聞こえないからと言って、消してしまうのは不十分である(高周波は実際は聞く事が出来ていて、 音色の違いとして感知している?)。

そういった意味で、スーパーオーディオCDやDVDオーディオはCDの倍以上のサンプリング周波数を採用している。 スパーオーディオCDは確かに普通のCDとは違う。聞けば分かる。

音楽の話はここまで。

画像に関してほぼ同様なことが考えられる。

ある空間スケールの情報を記録するには、 その倍以上のサンプリングが必要である。

つまり、x degスケールの情報を記録するためには、 pixelのサイズはx/2 deg以下でなければならない。

以下のような画像を考える:

  • 6.87 arcmin / pixl
  • smoothed 1 deg

画像自体は、FWHM = 1 deg のガウシアンでスムージングがかかっているので、 これより小さい角度スケールの情報は手に入らない。 pixelサイズはスムージングのスケールに比べて十分小さいため、 この画像は1 degの角度スケールの情報を完全に保持しているといる。

これ以上pixelのサイズを小さくするのは無意味である。

この画像をFWHM = 3 deg でさらにスムージングすることを考える。 このとき、できるだけピクセルのサイズは大きくしたい。 サンプリングの定理から少し余裕を持たせて、pixelの大きさはFWHMの三分の一、1 degとする。 このとき

  • FWHM = 3 deg でスムージングをかけてから、pixelの大きさを 1 degにする。
  • pixelの大きさを1 degにしてから、FWHM = 3 degのスムージングをかける。

の二通りの変換手順がある。 手順によって結果がどのように変わるか?

結果:

  • smoothing -> downgrade
    plot_smoothing_downgrade.png
  • downgrade -> smoothing
    plot_downgrade_smoothing.png

ぱっと見変化していない。定量的に調べてもほとんど変化していないことがわかる。

因みに、smoothing,downgradeともに値は保存している。

前者の場合、先ずはじめにFWHM = 3 degでスムージングをかけてしまっているので、 この角度スケール以下の情報はこの段階で失われてしまっている(pixelのサイズは十分。 でもスムージングなんだから、必ずしも「失われている」と言えるのだろうか?)。 そしてこれを、pixelサイズ 1 degでリサンプリングしている。 このpixelの大きさはナイキストのサンプリング定理をみたしている。

後者の場合、FWHM = 1 degの角度スケールの情報を持っている画像を、 pixelサイズ 1 degでリサンプリングしているので、 この時点で2 deg以上の角度スケールしか再現できない。 この画像に対して、FWHM = 3 degのスムージングをかける。 最終的にはpixelの大きさはナイキストのサンプリング定理をみたしている。

普通に考えれば、前者が正しく後者がも違っている気がするが、 結果はどちらも(ほぼ)正しい。

確かに、最終的にほしい角度スケールが3 degなので、 途中でいくらそれ以下の角度スケールに変換したところで(値が保存している限り)問題ないように思える。

smoothingは、FWHMの角度スケールに再分配、 downgradeはpixelサイズのカットオフという微妙な違いが、 直感的な理解の妨げになっていると思うのだが。

  • (l,m)空間で考えれば、前者はalmにb_l (beam)をかけてからp_l (pixel window)をかける、つまりalm->(alm*b_l)*p_lのに対し、後者はその逆でalm->(alm*p_l)*b_l。もちろん積は交換するから、両者はほぼ等しいですね。「ほぼ」とする理由はHealpixのharmonic transformが近似的であるためと、pixel spaceのdowngradeと、(l,m)空間でp_lをかける作業が完全に等しくはないためです。 -- コマツ 2007-12-12 (水) 22:01:52

11日Tue.

今日は朝から晩までがんばってしまった。 実働時間は、0730-2230の15時間。 この時間が精神的・肉体的にも限界だと思われる。

今日と同じ生活をすると、修論仮提出までに500時間超ほど時間がある。 今日は修論を文章だけで7頁ほど進めたので、この調子で書くと200頁超は硬い。 こんな風に皮算用していると悲しくなるので、一刻も早く研究を進める。

今日は、FWHMとsigmaの違いで時間を無駄にしてしまった。 window functionをはじめに書いたとき、 引数FWHMを何変換もせず、sigmaに代入してしまっていたため、 Clの計算が上手くいかない、という結果になった。

sigmaとFWHMの関係は(Full Width Half Maximumが示す通り)、

\Large  \exp\left[ -\left({\rm FWHM}\right)^2/(2\sigma^2)      \right]  =0.5


\Large  \sigma={\rm FWHM}/\sqrt{8\log 2} = 0.4245{\rm FWHM}  \quad\Longleftrightarrow\quad {\rm FWHM} = 2.35\sigma

となる。結構違う。

p_Cl_02_50.png
p_Cl_02_150.png
FWHM 3degなのでl=50まではと思ったが、微妙に小さい。線を太くしただけではごまかせないカットオフであるl>60からは一致しない

様々プログラムを書いているが、 ソースの分割、Makefile、ライブラリ化等が場当たり的なものになってしまっている。 一先ずソースの管理・ヴァージョン管理には「subversion」を使うことにした。 subversionを使えば以前のヴァジョンに戻る事が簡単であり、 ヴァージョン自体はsubversionで管理してくれるため、いちいちファイルにナンバリングする必要がない。 ローカル、勿論ネット越しにも同期が出来るので、 (ライブラリ当最低限の設定が出来ていれば)他のマシーンで同時に開発・実行することが可能である。 現在はlocalなマシーンで開発・テストを行い、PC cluster上で同期を行い実際に実行するといった形で開発を進めている。

後の問題は、ファイルの分割・ライブラリ、モジュール化だが、 こればかりはやはり経験がものを言うのだろうか。 一定のやり方があるのは知っているが、 それを今覚えるのは得策ではないな。


12日Wed.

今日は08000-2500。 ここ一週間AKARIの作業が遅々として進んでいないので、 集中的に取り掛かる。 今は、diffuseにマップ作成に向け、 beam patternの決定が自分の仕事。 NGC6543のデータを使い、これを行う。

AKARI内では、二次元ビームパターンの解析をする予定はないらしい。 zero sheet法と呼ばれる、他の方法とは「格」が違う奴に取り組んでいるのだが、 その難易度も「格」が違う。

国内外で文献がとても少ない。 参考になるのは唯一以下の本だけ:

先生と議論した末、 参考文献書いている人に直接聞こう、ということに。 これで解決に向かうといいのだが、 無理だった場合は、広く使われているiterativeな方法に移行する必要がある。

作業としては、zero sheet法を見越し、画像のフーリエ変換に関する実装を行った。 世間にFFTのライブラリは数多く存在するが、 C++と相性がいいのがなかなか見つからなかったので、取り敢えず自前で実装。 vector + GSL FFTで画像の離散フーリエ変換を行う簡単なライブラリを作成した。

これを使って遊んでいる最中に、会議開始・・・・。

ta.jpg

修論が進まない・・・。


13日Thu.

今日は0830-2330の15時間。 CMBを中心に、途中COEシンポに顔を出したりした。

CMB+Foregorundを分離するMCMCの実装は研究の要なので、 自分で開発する必要がある(汎用ツールがあり、一部を自分流に変更すればいい、という状況ではない)。

MCMCを使い、CMB and other dataから宇宙論パラメータの制限をつけるのには、 CosmoMCがよく使われる(godzillaみると、近頃I来さんがバリバリ走らせている)。 これは、宇宙論パラメータのMCMCなので、 自分の研究の後の段階で関わってくる(私は正確なClを計算したい。このClを宇宙論パラメータの関数として計算し、 宇宙論パラメータを推定するのがCosmoMC)。 なので、今の段階では関係ないと思い、触ってこなかった。

ちょっとしたきっかけ(前述のI来さん。CosmoMCが走っているのを見て、そういえば・・・、という感じ)でWeb頁を見たのだが、 どうやらCosmoMCには、宇宙論的計算をするプログラム「cosmomc」以外に、 Markov Chainを解析する「getdist」があることが分かった。 これを使えば、

weight like param1 param2 param3 ...

の様なデータを解析し、1d,2d,3dのmarginalizedされた分布や平均、分散、likelihoodなんかを計算してくれるらしい。 さらに、その結果をプロットするための命令を、以下に示すツールについて出力してくれる:

  • supermongo 1D plots
  • matLab '.m' files to do 1D, 2D and 3D plots
    • custom matlab plots

とても便利そうだ。 chainデータを解析するためにいちいち解析プログラムを組むのが、面倒だと感じ始めたところだった。

最もきれいな(ベストな)解は、「MATLAB」を買うことだが、 今は金がない。当分はGnuplotでがんばるしかないのだろうか。 MATLABのクローンが各種あるとI藤さんに伺ったのだが、それも試して見る価値はありそうだ。

ba.jpg

「今日」は修論の題名提出期限。


14日Fri.

0730-2300。途中少し飲み会。

AKARI、CMBともコツコツと進める。 CMB+Foregroundの分離のために、C++のソースを書いているが、

・・・・

今日は眠いので寝ます


15日Sat.

今日は0900-2000。 本とはもう少し早く大学にいけるはずが、 ○イ○がとんでもないことになっていたので、 時間を食ってしまった。

修士論文発表会の日時が決定した。

  • 2月5日(火曜日)
  • 12:00-14:00
  • A6SM3006茅根裕司
  • 階層ベイズ法を使った宇宙マイクロ波背景放射・前景放射の分離法の開発

MATLABのクローンについて調べてみた。 といってもグラフィックス中心に。

MATLABはとても高いソフト。 でも、研究の分野ではよく使われるツールなので、 金がない人にも使えるように、クローンのソフトが有志によって開発されている。 その代表が「GNU Octave」。 MATLABの基本的な機能は網羅している。 しかしグラフィックスがGNUPLOTである。 (GNUPLOTもそれ自体でがんばればきれいだが)プッロターとしては貧弱である。 そもそも、自分の目的は、CosmoMCが吐くファイルから、 MCMCの解析・プロットを綺麗・簡単にすることである。

次が「FreeMat」。

  • Q. Is FreeMat? 100% compatible with MATLAB? What about IDL?
  • No. FreeMat? supports roughly 95% (a made up statistic) of the features in MATLAB.

なかなかの互換性をもつ。グラフィックスに関しても、contour plotをサポートしている。 実際、CsomoMCが吐くファイルを正確に読めるかは不明。 インストールを試みたが、「QT」のver4が必要といわれた。 現在手元のver.は3。 アップデートすればいいのだが、如何せんこれを上げることはDebianのver.をsargeからetchにアップすることを意味する。 X及び基本ライブラリ群がアップデートしてしまうので、何らかの不具合がでないとはいえない。

普段であれば、強行するところだが、何せ今は修士論文執筆中、 何かあったらいやなので、今のところ見送っている。 やるんであれば、手元のcoLinuxに入れてみるつもり。

Octaviz」はOctaveのグラフィックスの「拡張」なので意味なし。

最後は「OctPlot」。

OctPlot is a handle graphics package for Octave, the free alternative to MATLAB

Octaveのグラフィックス部分を、GNUPLOTとは別に開発したもの。 contour、subplot等に対応している。 もしかしたら上手くいくかもしれない。 ところがこれも、Debian向けはetch以上に上げなければならない。

godzilla上で出来ると一番いいのだが、どうやら今の段階では入っていないらしい。

何がベストな解なんだろうか・・・。 こんなことしてる間に時間が過ぎていく・・・・・。

esya.jpg

16日Sun.

MATLABは、評価版を使うことで、取り敢えずの解決をしようかと思います。

I藤さんにさらにscilabを教えていただいた。 PC Clusterにユーザー権限で入れてみようとしたが、 ライブラリが不足している。 これもやはり手間がかかる。 又、

コマンドは、MATLABに似ているが互換性があるわけではない。

だそうだ。

ここ数日間MATLABのクローン捜索を行ったが、 自分の目的に合致したものは見つからなかった。

しかし、

  • 研究を止める訳にはいかない
  • お金がない
  • 今後も解析ツールを自分で作っていくには、「車輪の再発見」に過ぎない。折角getdist -> MATLABで綺麗な絵(結構よく論文でも見る)をかける(解析が本質ですが・・・)使わない手はない。

ので、今後MATLABを導入することを考え、 現段階では取り敢えず「評価版」使うことに決めた。

評価版は、ユーザー登録をすれば入手可能だそうなので、 早速登録した。 期間は30日間。

絵を描ききるには十分だと思われる。

自分のプログラムの実行時間が遅い。 今まではchain 15 secだったのが、 1 min近くかかる。 前と今での違いは、バンドを1 -> 4に増やしたこと。 プログラム自体は位置から書き直しているので、 バグが混在しているかも知れない。

バンドが4つになったのだから、 時間が4倍になってもいいかとも思うが、 このプログラムで一番時間がかかるのは、 球面調和関数関連の箇所である「はず」であるので、 バンドの計算中球面調和関数が計算されるのは一度だけであり、 時間が四倍になることは無い。

なにより、入力バンド数パラメータをN_NU=1にしても変化しないのがその証拠。

明日は、時間食うのは「map2alm_iter」なので、 この周りを調べる。

ne.jpg

この一週間は睡眠時間が減った。


17日Mon.

1000-2330

MCMCの問題点と研究

MCMCは静的モンテカルロ(円とそれに外接する正方形からπを求めるやつ)で起こる困難を克服するために使われる。 このとき、 円の中の一様分布p(x,y)を、外接する正方形ないの一様密度q(x,y)で近似し、 規格化定数Zからπが求まる。 これを拡張すればN次元の球の体積を求めることが可能である。 しかし次元が上がるにつれ、 計算の効率は著しく落ちる。 この原因はqがpの(大局的な)よい近似になっていないからである。

大局的に未知の分布を、 大局的に適当に仮定してしまっていることがこの問題の本質である。

この様に静的モンテカルロは、 対象が多次元・多変量になると、 「次元の呪い」により使い物にならない。

そこで、「動的」なモンテカルロである、Markov Chain Monte Carlo法が開発された。 MCMCは大局的に分布を近似するのではなく、局所的に分布を近似している。 大局的には複雑な形状をしている多次元・多変量分布も、 局所的に見ると、単純な分布(例えばガウシアン)で十分近似しうるものとなる。

今研究では自由度十万単位で最適化を行いたいと考えている。 この様な場合、局所的に見ても複雑な分布をする場合がある。

今の場合、更新すべきピクセルの数は四万である。 もし一度にすべてのピクセルを更新するならば、この四万次元空間からのサンプリングとなる。 これは現実的ではない。

こういった場合、 ピクセルを適当に分割(10〜1000?)程度に分割し、 グループごとにサンプリングをすることが考えられる。 四万に比べれば、この程度の次元からのサンプリングは(割合)ラクである(いずれにせよ更新の刻み幅で左右される。 妥当な刻み幅を設定するのが問題)。

ここで書いたことは、Metropolis-Hastings法の場合である。 大局的な分布が既知の分布でかけるとき、そのサンプリング効率は100%である。 CMB+Foregroundの場合、 その事後分布はガウシアンxガウシアンなので、 一つのガウシアンでかける。 しかしながら、このガウシアンに変形する計算は、 サンプリングごと実行する必要があり、計算時間はピクセル数の三分の二乗のオーダーである。

先にあげたMetropolis-Hastingの場合、 この様な変形は必要なく、計算時間を食うのは唯一球面調和展開のみである(Gibbs samplerでもこれは必要)。 問題は先に挙げたサンプリング効率でる。

あるグループの提案分布を計算するには、球面調和展開が必要である。 これをグループの数だけ実行する必要がある。

Gibbs samplerは球面調和関数展開を一度したあとに(本当に「最初」の一度だけ?)、 サンプルごとにガウシアン化を計算する。

Metropolis-Hastingは効率が悪いかもしれないが、1chainの計算回数が、 Gibbs samplerに比べ早ければ、数で対抗できる。

実際、この二つの計算時間がどれだけかかる見積もりが出来ていない。 オーダーの計算が今ひとつ分からない(最初は分かってたと思って、 計算オーダーがコンパラになるように分割するライブラリを作ったりした)。

ガウシアン化には、

For white noise, Nij = σiδij ,
this is the spherical harmonic transformrequired between
pixel (for noise covariance matrix multiplication) and
harmonic (for beam convolution and signal covariance
matrix multiplication) space, with a scaling of O(N^3/2pix ).
For correlated noise, it is the multiplication with a dense
Npix ×Npix inverse noise covariance matrix, with a scaling of O(N^2pix).

とあるが・・・(球面調和展開はN^3/2のオーダー)。

Eriksenの論文では自分と同じパラメータで(自分はdustの項は入れてない)

The wall-clock
time for generating one single sample is 50 s,
parallelized over five 2.6GHz AMD Opteron 2218 processors,
one for each frequency band. We generate five chains
with 1000 samples each, for a total wall clock time of
14 hr. The total computational cost is 350CPU hours.

とあるが、自分は4cpu(3.00GHz Xeon)を使い、全バンドが1chainが20secで計算できるのだが・・・・。

08.jpg
注:絵です

修論の審査員は公表されるのでしょうか?

昼までに天文学会にアブスト出さなきゃ。 大学のなんかもあったような、大学院GP? やったほうがいい?


18日Tue.

0730-2300

dame.jpg

今帰宅してこれを書いている訳ですが、わかってしまった、すべて。 馬鹿みたいだ。やばい。

「The wall-clock time for generating one single sample is 50 s, parallelized over five 2.6GHz AMD Opteron 2218 processors, one for each frequency band. We generate five chains with 1000 samples each, for a total wall clock time of 14 hr. The total computational cost is 350CPU hours.」

共有型メモリで二千個CPUの並列化すげーと思っていたが、そんなものはない!? 「five」って書いてあるじゃん!?、5だよ伍!? 「2218」ってOpteronの型番だよ!?

・・・・・・・(´・ω・`)

me.jpg

これで計算のオーダーについて再考してみる。

自分の場合一番計算時間がかかるのは、 グループ化しそれをグループごとに提案分布を計算する箇所である。 球面調和展開がNの二分の三乗のオーダー、 採択率を上げつつ計算回数をリーズナブルに抑えるために、 分割するグループの数をNの二分の一乗オーダーにとる。 結局一回のチェインでの計算回数はNの二乗となる。

CG法を使う前に、 Gibbs sampler(ガウシアン化)の場合もまず球面調和展開が必要である。 (CG法は反復ごとに球面調和展開が必要?) このオーダーはNの二分の三乗だが、 もしノイズに相関がある場合、これとの積のオーダーNの自乗が卓越する。

そしてCG法だが、この方法は反復法なので必ずしもいつも同じ計算回数で、 解が求まるとは限らない。解の安定性の問題もある。 最低でもオーダー一乗かける反復回数が必要である。 結局ガウシアン化の方法では、オーダーはいくら?

眠いので寝る


19日Wed.

0900-2400

計算オーダー

  • Eriksenはcorrelated noiseを考えているので、オーダーはNの自乗
  • 自前のはuncorrelated noseでやっているが、反復計算でNの自乗オーダー

なので、前者50sec/chain, 15sec/chain、オーダーであっている、でFA?

ここまで高次元で1 chainあたりの計算時間が膨大だと、 サンプリングの効率の問題は避けては通れない。 自由度自体は変わらないが、 CMBを実空間とフーリエ空間でサンプリングするのはどちらが効率がいいか。 普通効率を上げるためには、 一回で複数個をまとめる。 このとき、互いに依存している場合はその効果が現れやすい。

CMBを実空間で適当にグルーピングした場合、それらは独立である。

CMBは一様等方であるから、ある角度スケールlについては同じような振る舞いをするはずである。 つまりあるlについての温度揺らぎ(-l<m<l)は同じような振る舞いをする(本当?)。 よってある角度スケールlごとに2l+1個の値を更新することは効率のよいグルーピングだと考えられる。

Foregroundは一様等方ではないので、この考えは使えないと思われる。

今日もがんばった、眠い。

nya.jpg


  • N川さんに資料を送る
  • d井さんに位置決めのために二つ山のFITSを送る

20日Thu.

1200-2500

検証しなければならないことがでてきた。 時間が・・・。

a_{lm}を実数で振っても問題ない?

syo.GIF
  • ん?どういう事でしょう? -- コマツ 2007-12-22 (土) 23:33:33
  • 何時もありがとうございます。 疑問は、実空間ではサイコロを振ってある空のCMBの値をサンプリングしますが、フーリエスペースでa_{lm}に対してサイコロを振るとき、 位相の情報をどのように扱うのが適切なのかということです。P(a_{lm}|d)\propto \exp [ -(d_{lm}-a_{lm})^tN_{lm}^{-1}(d_{lm}-a_{lm})/2]からa_{lm}をサンプリングするときは、位相はd_{lm}と同じで、絶対値についてサイコロを振るのが正しいのか(実部・虚部、絶対値・位相に対してサイコロを振るのではなくて)?P(a_{lm}|d)\propto \exp [ -(a_{lm}-\hat{a}_{lm})^t(N_{lm}^{-1}+S_{lm}^{-1})(a_{lm}-\hat{a}_{lm})/2] からa_{lm}をサンプリングするとき、位相も含めて\hat{a}_{lm}を解き、位相は\hat{a}_{lm}と同じで、絶対値をgaussianからサンプリングするのが正しいのか(複素ガウシアンからのサンプリングではなくて)? -- 茅根 2007-12-24 (月) 03:22:15
  • もしalmをCMBと思っているのであれば、複素ガウシアンで実部、虚部に対してさいころを振るのが一番簡単だと思いますよ。絶対値はガウシアンでなくレイリー分布、位相は一様分布なので。あと細かい事ですが、指数の中身は(N^{-1}+S^{-1})ではなくて、(N+S)^{-1}ですよね? -- コマツ 2007-12-25 (火) 02:32:32
  • あ、(N^{-1}+S^{-1})で良いですね。Eriksenの(10)式を書いていたのですね。 -- コマツ 2007-12-25 (火) 02:49:37

23日Sun.

もう23日。昨日はあまり研究できず、そのまま飲み会へ。 28時帰宅。今日は12時起床。それから大学へ。

飲み会はとても楽しかったのですが、後に及ぼす影響は甚大です。

色々考えた結果+現在の結果から、 実空間でCMBのサイコロを振るのは効率が悪いということが確信できた。 foregroundを発散させないようにすれば、 長いchainを作ることで、実空間でもおそらく収束すると思われるが、 効率が悪い。二、三日計算まわしていればなんとかなるが。

koumei.jpg

やはりフーリエ空間で角度スケールごとにサンプリングした方がいい。 この際、事後分布はガウシアン型にかけるのでサンプリングの効率は上がる。 このとき逆行列を求める必要があるが、 これはCG法等によって解決する。

Erikseの論文では、 線形パラメータを全部この方法でやっている。 それが出来るのは、foregroundに複雑なpriorを仮定していないからである。 foregroundの線形なパラメータでも、 MEMや隣接ピクセルの情報を使うような形のpriorを導入した場合、 単純なガウシアンの形にはかけない。

そこで、 CMBシグナルだけをGibbs samplerからsamplingし、 foregroundはMetropolis-Hastingsでサンプリングすることで、 Foregroundについて様々なpriorを仮定することが出来る。 foregroundはgibbs sampler(条件付確率の意味での)の形で書けば、 基本的にpixelごとに値を更新していくことで、 最大事後分布に収束する。 (CMBはpriorに全天の情報が入っている。このためpixel毎に更新しようとすると、 そのたびに全天でのフーリエ変換をしなければならず、時間がかかりすぎる。 foregroundは注目しているピクセルだけの計算で済むので、ピクセルの数が増えても最大でNの自乗のオーダー(correlated noise)にしかならない)。

問題はCG法と、多次元ガウシアンからのサンプリングをする方法。 取り敢えず調べたことをメモ:

  • BFL :期待値、共分散をセットすると、そこからboxmuller法でサンプリングしてくれる関数がある。
  • Intel Math kernel
  • IML++ :複素数無理?
  • SparseLib++ :複素数無理?

逆行列は位相を含めて解く。 そしてそのサンプリングは「絶対値」で行い、 サンプリングされた「絶対値」と「位相」を組み合わせて、 実空間に戻す。

  • COCG :COCG法は係数行列が対称であるがエルミートでない線形方程式に対して有効であり、そのアルゴリズムはCG法と極めて酷似している。
  • ScaLAPACK
  • SILC

どれから片付けるべきか・・・。


24日Mon.

1000-2000

今朝バスを待っていたら、全然くる気配がなかった。 十数分待ったと後で(このぐらい遅れることはある)、 今日は天皇誕生日の振り替え休日それ以上でもそれ以下でもないということに気付いた。 休日ダイアだった。

20071224.png

どうやら、サンプリングについては勘違いをしていたようだ。 取り敢えず、uncorrelatedな場合であれば、 フーリエスペースでもサンプリングできる(自分でプログラム組める)確信がもてた。

correlated noiseの場合、

  • やはりサンプリングの際の位相の取り扱い。
    • 位相も含めて方程式を解くのか?
      • ⇒多分Yes
      • よってCGではなくCOCG法が必要。
    • 昨日の意味ではない。昨日のは間違い。
  • {\bf N}^{-1}を求める手間。これは一度求めればいい。皆使うときは厳密に求めているのか?
    • {\bf N}^{-1/2}といった行列の平方根。コレスキー分解、特異分解等で求めるが、やはりこれも位相込み。
    • これらもやはり位相込みで計算が必要なので、複素数までかばしたライブラリが必要。

当面は、correlatedノイズは扱う予定はないので問題ないが、 根本的なことが理解できていないので、気分が悪い。

  • Nがピクセル空間のものであれば、N^{-1}は厳密に求まります。というか、N^{-1}が最初に求まるもので、チャレンジはNを計算する事にあります。大抵の場合はNを計算する事はないですけど。N^{-1/2}は、ちのねくんの言うようにCholeskyやSVDを使って求めますが、あまりピクセル数が大きいと計算時間がかかりすぎますよね。なので、僕らはnside=16までは厳密に扱い、それ以降はノイズはwhite noise (uncorrelated)としています。ところで「位相込み」というところが良くわかりませんが、Nに位相ってないですよね? -- コマツ 2007-12-25 (火) 02:59:08

25日Tue.

1200-2400

今日は薬飲んだので眠いです。 寝ます。


添付ファイル: file20071224.png 753件 [詳細] filekoumei.jpg 748件 [詳細] filesyo.GIF 747件 [詳細] filenya.jpg 897件 [詳細] filedame.jpg 733件 [詳細] fileme.jpg 692件 [詳細] file08.jpg 715件 [詳細] filene.jpg 709件 [詳細] fileesya.jpg 729件 [詳細] fileba.jpg 715件 [詳細] fileta.jpg 674件 [詳細] filep_Cl_02_50.png 725件 [詳細] filep_Cl_02_150.png 719件 [詳細] fileplot_smoothing_downgrade.png 745件 [詳細] fileplot_downgrade_smoothing.png 725件 [詳細] filebook.jpg 729件 [詳細] fileomae.jpg 753件 [詳細] file7.jpg 710件 [詳細] filesnow01.jpg 704件 [詳細] filemac.jpg 883件 [詳細]