Fire Engine

化学の修士号→消防士→ITエンジニア(2016年11月〜)

企業研究者の立場からKaggleに取り組む意義を考えた

先日,Kaggleで初めてコンペに挑戦し,その振り返りをブログに書きました. 現在,企業で研究者として働いている私は,Kaggleのコンペに取り組むことは非常に学びが多く,自身の研究活動にも良い貢献をするだろうと確信しました. 私自身Kaggleに取り組むまでKaggleと研究に繋がりを見いだせておらず,実際にコンペに取り組むことでその繋がりが見えてきました. まだ1コンペしか参加経験がないビギナーではありますが,私が考えている研究者としてKaggleに取り組む意義について現時点の考えをまとめたいと思います.

自分の立場について

Kaggleへの取り組みが研究に良い貢献をするかどうかは,研究の分野や内容に依存すると思うので,私の立場をはっきりさせておきます. 私は現在,インターネットインフラサービスを提供するさくらインターネット株式会社の組織内研究所であるさくらインターネット研究所で研究員として働いています. 私は,自分たちのドメインにおける課題の中で人間の手や簡単なルールベースのアルゴリズムでは解決が難しい課題に対して,機械学習を適用して解決する研究に取り組んでいます. したがって,私は機械学習の理論の研究をしているわけではなく,機械学習を応用して研究を行っている立場です.

機械学習を応用する研究

機械学習を課題解決の手段として用いる研究は,多くの分野で増えていると思います. 私が参加した2020年11月開催のIBIS2020の発表の中でも,農業・土木工学・機械設計など様々な分野に機械学習が応用されていました. 私の分野においても例外でなく,2020年に開催されたクラウドコンピューティング技術に関する国際会議であるIEEE CLOUD 2020において,機械学習をはじめとした数理的アプローチで課題解決を行っている研究が多数発表されていました.

我々の研究所ではこういった流れに取り残されず研究を行っていくためにも,研究の問題設定を行うドメイン知識や経験が豊富な研究員と,その課題を機械学習統計学を駆使して解決する研究員で役割分担をしてチームとして研究を行う取り組みを積極的に行っています.以下はその役割分担について,同僚の松本さん発表資料から抜粋したものです.

f:id:hirotsuru314:20210114111528p:plain

Kaggleを知って危機感を覚えた

前述のような自分たちのドメインにおける課題に対して機械学習を適用する研究において,機械学習を用いた問題解決力が研究に直結することは明らかです. 一方,私がKaggleの存在を知ってから,公開されている過去コンペの解法についての発表資料をいくつか読んでみると,その中の多く解法が私が知らない,もしくは名前だけ知っているけど理論や使い方もわからない手法を用いていました.(例を挙げると,LightGBM・BERT・Metric Learningなどなど) 私はこのことに非常に危機感を覚えました.というのも,もし研究における問題設定がうまくできていたとしても,問題解決力が不足していて,実用上有効な精度が出せないと研究が進まないといったことが起きると思います. 例えば,我々の分野においてWebシステムの異常検知という問題設定をした場合,その検知精度が60%であったとしたら,実用性があまりにも低いので,どれだけ素晴らしい問題設定や提案をしたとしても研究として成立しないと思います. 世界中のデータサイエンティストが精度を競い合うKaggleにおいて頻出のモデルやアルゴリズムをキャッチアップできていない今の自分の力ではそういった問題が頻繁に起きるのではないかと思いました.

Kaggleを通じて問題解決力を高める

そこで,Kaggleを通じて機械学習を用いた問題解決力を高めることは,自分の研究活動において有効だと考えました. Kaggleで得られる知識や経験はいろいろありますが,私が特に研究に直結すると思ったのは以下の3点です.

1.最新かつ実用的な手法のキャッチアップ

Kaggleでは枯れた技術のみが使われているわけではなく,最新の研究や論文の成果がどんどん実践で使われています.逆に言うと,state-of-the-artを追わないと勝てないのかもしれません.このため,コンペに参戦することや参戦後に上位解法を学ぶことで,数年以内に提案されたような新しい手法をキャッチアップできると思いました. さらに,新しいだけでなく,Kaggleのコンペで実データに適用されて,精度が検証されていることが非常に有用だと思いました. 手法の提案論文では,(査読の段階で取り除かれることも多いと思いますが)提案に有利な評価であったり,実用的でない擬似データを用いた評価であることが多いです.それがKaggleで実戦投入されて,論文とは異なる実データで精度や実用性が検証されていることには非常に価値があると思いました.

2.実装力

Kaggleのコンペの中でも「Code Competition」と呼ばれるコンペでは,与えられたNotebook環境で成果物をSubmitする必要があるため,メモリの制限やコードの実行時間に制約が課されます.これにより,省メモリかつ高速な実行コードを書く力要求されます. また,前述のような最新の手法の中にはライブラリが公開されていないものもあるため,時には論文を読みながら自分で実装する力も必要になると思います.

3.データの捉え方

Kaggleでは様々な種類のデータ・問題設定でコンペが開催されるため,継続的に参加することで幅広いデータの捉え方を学べると思いました. これは自分たちの手持ちのデータだけを分析していてはなかなか得られない経験なので,非常に有用だと思います.

また,これらの習得できる力を踏まえた問題解決力が,コンペの順位,メダル,称号などを通して客観的な指標として評価できることも魅力の一つだと思います.Kaggleで学べることはまだあると思いますが,上記の3つだけでも,どれも自分の研究にも必要な力であるので,私は研究に必要な機械学習を用いた問題解決力の習得・維持のためにKaggleに継続的に取り組んでいきたいと考えました.
ここまでが企業研究者の立場である私がKaggleに取り組んでいこうと思う主な理由なのですが,Kaggleでコンペに取り組むことが直接研究の業績に繋がるような事例もあるようなのでまとめました.

Kaggleの成果が直接的に研究業績に繋がる事例

ここからは,Kaggleのコンペで上位入賞するような場合に限った話ですのでハードルは高いです.事例を以下の3点にまとめました.

1.Kaggleの解法を論文としてまとめる

  • 「PLAsTiCC Astronomical Classification」というコンペの優勝者が書いた論文がジャーナルに採択

    iopscience.iop.org

  • Kaggleの過去コンペの結果をまとめたようなジャーナル論文もある

    www.sciencedirect.com

2.コンペで入賞して学会で発表する

  • 私が前回参加した「Riiid! Answer Correctness Prediction」では上位入賞者は「AAAI-2021 Workshop on AI Education」で発表できる.

  • 国際学会が主催しているコンペもある.例えば,IEEEが主催した「IEEE-CIS Fraud Detection」で上位入賞者はIEEE CIS Conferenceで発表できる.

3.コンペで入賞してプレスリリースを出す

多くの企業がKaggleのコンペで上位入賞した場合,以下のようにプレスリリースを出しているようです.

www.nssol.nipponsteel.com

www.jri.co.jp

さいごに

私は自分の研究を進めていく上で必要な機械学習を用いた問題解決力の習得・維持のためにKaggleに継続的に取り組んでいきたいと考えています.いずれはコンペで上位入賞できるような力をつけ,Kaggleの成果自体が研究業績になるようにしていきたいと思います.このブログが私と同じような立場や危機感を持っている方々の参考になれば嬉しいです.Kaggleやっていきましょう!

Kaggle初コンペの振り返り〜Riiidコンペで銀メダル獲得〜

Kaggleで開催されていたRiiid! Answer Correctness Predictionに参加しました.結果を簡単にまとめると以下の通りです.

  • 順位:139位(3406チーム中)
  • メダル:銀メダル(上位5%以内)
  • 解法:LightGBMとSAKTのアンサンブル
  • チーム:1人で参加

自身初となるデータ解析コンペでしたが,なんとか銀メダルを獲得することができました.今回はその振り返りを時系列で書いていきたいと思います. 先に感想を述べると,Kaggleは本や論文ではなかなか得られない学びに溢れている上に,世界中の人々と順位を競い合うゲームのような感じでとにかく楽しいので最高すぎました. この記事がこれからKaggleを始める方々の参考になれば嬉しいです.

目次

Kaggleとの出会い

Kaggleの存在を知ったのは4年くらい前ですが,Kaggleを始めるきっかけとなったのは2020年11月に開催された第23回情報論的学習理論ワークショップ (IBIS2020)という学会でした.IBIS2020のチュートリアルセッションにて「Kaggle・実践データ解析入門」という発表がありました.この発表では,Kaggleの称号の最高位であるGrandmasterを持つ登壇者の方が,Kaggleの概要や,「こういうデータのコンペでは初手としてこういう手法が有効である」といった実践的なテクニック,研究や論文の成果をどのようにKaggleに使っているか,など非常に興味深い内容が盛り沢山でした.動画や資料は公開されていませんが,Kaggleの概要の部分に関しては,発表中で紹介されていた以下のスライドがまとまっています.

結局、Kagglerは何を必死にやっているのか? / What is Kaggle? - Speaker Deck

この発表を聞いて,単純に「Kaggle面白そう!」と思いました.さらに,Kaggleでは最新の研究や論文の成果がどんどん実戦で使われていることを知り,Kaggleと研究に繋がりが見えたことが始める大きなきっかけとなりました.実際に登壇者の方は常に最新の論文から新しく提案された機械学習モデルをチェックしていて,精度が良さそうなものがあれば,論文を読みながら実装してKaggleで実戦投入するといったことをやっているようです.現在私は,さくらインターネット研究所で,機械学習をソルバーとして用いて自分たちのドメインクラウドコンピューティングなど)における課題を解決する研究に取り組んでいるため,機械学習を用いた問題解決能力は研究に直結します.そのため,Kaggleを通して最新のモデル,かつ提案論文内での評価とは異なるデータで実用性が検証されているモデルを把握しておくこと,また実践を通じて自身の問題解決能力を高めることは,研究活動に大きく貢献するだろうと確信しました.そんなこんなで興味を持ってすぐにアカウントを作ってKaggleを始めました.

コンペ参戦準備

Kaggleを始めるといっても,コンペへの参加の仕方さえ分からなかったので,まずはKaggleの全体像を掴むところから始めました.具体的には,以下の本を一通り手を動かしながら読みました.

bookclub.kodansha.co.jp

この本は,Kaggleをゼロから始める人にとっては非常に良い本で,アカウント作成からコンペへの参加,予測結果のSubmitまでの流れを一通り学ぶことができるため,初めてのコンペ参加へのハードルをかなり下げてくれます.

コンペに登録

Kaggleスタートブックを読み終えた時点ですぐにコンペに参加登録をしようと思い,良さそうなコンペを探しました. Kaggleスタートブックの中で初めてのコンペを選ぶときのコツとして,「すでに開催から一定時間が経過していてコンペ終了まで一ヶ月程度のコンペがおすすめ」と書いてありました.その意図としては,開催から一定時間が経過していると,Discussionや公開Notebookなどの充実した情報が集まっているので,それを参考に自分のモデルを作ることができるためです.これを参考にRiiid! Answer Correctness Predictionに登録しました.コンペの概要は以下の通りです.

  • 開催期間:2020/10/06〜2021/01/08
  • データ:TOEICの学習アプリで各ユーザが過去に取り組んだ問題とその問題に正解したかどうかのデータ
  • 問題:将来あるユーザはある問題に正解できるかを予測する二値分類問題

このような人間の知識の習得をモデリングする問題はKnowledge Trackingというそうで,教育分野では重要な問題なようです.
私がコンペに登録したのが2020/11/25で締切まで約1ヶ月半程度全力で取り組みました.ちなみに,取り組んだ時間としては,毎日寝る前2~3時間程度で,大晦日は格闘技を見るためにやりませんでしたが,それ以外は休まず毎日やりました.

初Submit

まず簡単な特徴量とモデルでSubmitまでやってみようと思ったのですが,このSubmitまでが鬼門でした. 本コンペは,KaggleのNotebook環境を用いてSubmitする必要がある「Code Competition」と呼ばれるコンペでした.このコンペの特性上,以下の2点がSubmit時の大きな制約となりました.

  1. Notebook環境のメモリ上限は16GB

  2. Notebook環境上での予測時間は9時間以内(テストデータは250万行)

この2点の制約により,機械学習やデータサイエンスの力だけでなくエンジニアリング力も試されました. まず,メモリ上限についてですが,学習データのCSVファイルを普段通りpandasでread_csvしたところ,読み込みに数十分かかる上にメモリが上限の16GBを超えてしまいエラーになりました.というのも学習データのCSVは,行数が1億を超えていました.さっそくデータの読み込みさえできずに焦りましたが,公開Notebookに「Tutorial on reading large datasets」という素晴らしいNotebookがあり,データ型をきちんと指定することや,読み書きが早いファイル形式を利用することなど知見が溜められていました.

次に,予測時間の制限についてですが,私は最初,学習データで作った独自の特徴量をPandasのDataFrameで保持し,テストデータに対してPandasのmerge関数でLeft joinすることで特徴量を追加していましたが,どうやらこのDataFrameのjoinの処理が遅いらしく,推論時間を9時間以内おさめることができませんでした.これもまた素晴らしいNotebookに助けられて「LGBM with Loop Feature Engineering」に特徴量を辞書型で保持し,テストデータのDataFrameの行に対してfor文を回し,特徴量を追加していくといった手法でした.これを真似ると推論時間は大幅に改善されました.このように悩んだことがあっても多くのKagglerが有益な情報を公開してくれていて学びながら取り組んで行けることもKaggleの魅力だと実感しました.(補足:Pandasのleft joinを357倍高速化したというNotebookも公開されていてこちらも参考になりそうですがコンペ中は追いきれませんでした)

これらの問題を解決して,やっとSubmitまで行き着きました.モデルには,LightGBMを選択しました.LightGBMは,ここ数年のテーブルデータのコンペの初手として多くのKagglerが用いているようです.なぜ初手にLightGBMを用いるかについては,以下の記事にわかりやすくまとまっていて非常に参考になります.

upura.hatenablog.com

初Submit後にLeaderboard(コンペの順位表)に自分の名前が載っているのを見てとにかくモチベーションが上がりました.

特徴量エンジニアリング

コンペ期間において一番時間を割いたのは特徴量エンジニアリングです.まずは簡単に思いつくようなユーザの正解率や,問題の正解率など,直感的に有用そうな特徴量を追加していきました.あとから知ったのですが,これはTarget encodingと呼ばれる特徴量生成のテクニックでした.Target encodingとは,説明変数に含まれるカテゴリ変数と目的変数を元にして特徴量を作り出すことです.
ここでTarget encodingにより明らかに有効そうな特徴量を追加して手元のスコアが上がっているのに,Submitしたスコアが全然上がらないという問題に悩まされました.これはTarget encodingにおける定番的なミスで「リーク」という問題を起こしていました.例えば,ユーザの正解率という特徴量を作るときに,学習データ全てを使って算出したユーザの正解率を特徴量として使うと,ある時点では将来の情報が漏れていることとなってしまい,精度が過剰に高くなってしまいます.そこで,ある時点のユーザの正解率はそれより過去のユーザの問題の回答から正解率を算出するように変更したところ,手元とSubmit時のスコアに相関が取れるようになりました.私はこれまで数年ほど機械学習の勉強をしてきたのにいざ実戦でやってみると知らないことが多いなと辛い気持ちになりましたが,これが成長というやつです.
自分がコンペに取り組む上で重要な知識が欠落していること認識したので,知識を補うために以下の本を読みました.

gihyo.jp

この本は,Kaggleのコンペに取り組む上で必要な情報が網羅的に書かれていてめちゃくちゃ参考になりました.前述のTarget encodingのリークの話もしっかり書いてありました.この本を参考に思いつく特徴量をどんどん追加していきました.最終的な特徴量の数は30~40程度でした.

モデルの構築

モデルは,最初のSubmitの時に採用したLightGBMを主に使っていました.特徴量エンジニアリングの成果が出て,早い段階でLightGBMのシングルモデルで銅メダル圏内に届きました.しかし,特徴量エンジニアリングが行き詰まってきてスコアの伸びが小さくなってきたので,モデルのパラメータチューニングに着手しました.パラメータチューニングにはOptunaを用いました.

github.com

Optunaには,LightGBMに特化したLightGBM Tunerという拡張機能があり,Pyhtonコードのimport文を1行変更するだけで簡単にチューニングができて便利でした.しかし,このチューニングはスコアにはあまり貢献しませんでした.
この点については,コンペのDiscussionの「LGBM Parameter Tuning」で議論が展開されており,パラメータチューニングは大きく精度に貢献しないことなどが書かれていました.このようにコンペ開催中にも惜しみなく知識や経験を共有する文化は素晴らしいと思いました.
締切まで残り3週間程度の時期にLightGBMのスコアが停滞してきたので,この時点で他の学習モデルの検討を始めました.

アンサンブル学習

LightGBM以外の学習モデルには,TransformerをベースにしたSelf Attentive Knowledge Tracing(SAKT)を選択しました.その理由は,公開されているNotebookの中で,LightGBMなどの勾配ブースティング系とSAKTのアンサンブル学習で大きくスコアを伸ばしているものがちらほら見つかったためです.SAKTは以下の提案論文と公開Notebookを読みながら実装しました.

arxiv.org

SAKTの実装は完了したものの,LightGBMのスコアには到達しませんでした.しかし,アンサンブル学習を試してみると一気にスコアが伸び,そのときの順位で100位くらいで銀メダル圏内に到達しました.アンサンブルの方法は単純な加重平均で,この時点でのスコア的にはAUCでLightGBMが0.778,SAKTが0.772で5:5のアンサンブルをすると一気に0.786まで上がりました.
このようにモデルのアンサンブル学習を行うと,単体のモデルの最高スコアを大きく上回ることが起きうるというのは「Kaggle Ensembling Guide」を読んで知っていたのですが,どうも直感的には低いスコアの方に引きづられそうな気がしていました.ただ,結果を見ると,これがアンサンブル学習の力か…となりました.また,ここでスコアを大きく向上させられた大きな要因として,アンサンブル学習を行ったのが勾配ブースティング系のモデルと,ニューラルネットワーク系のモデルの異なる種類のモデルをアンサンブルさせたことだと思われます.これもアンサンブル学習にはモデルの多様性が重要であるとよく語られることではあるのですが,本などで学んだこともこうやって実際に体験すると,より知識として定着するので実践から学ぶのは良いなと改めて思いました. あとはコンペ終了まで,ひたすら特徴量エンジニアリングを頑張ってLightGBMとSAKTのアンサンブルでスコアを出すというのを繰り返していました.

コンペ終了:Leaderboardの公開

2021/01/08の午前9時にコンペが終了し,最終的なLeaderboardが公開されました, 結果は3406チーム中139位でなんとかソロで銀メダルを獲得できました.

f:id:hirotsuru314:20210108230540p:plain

コンペ全体を通しての感想と今後

初コンペを終えた感想としては,Kaggleはとにかく学びが多くて,とにかく楽しいので,これからもKaggleを続けていこうと思いました.当初は初コンペでメダルが取れたら嬉しいなーくらいの気持ちでしたが,いざ銀メダルを獲得すると,嬉しさよりも,上位陣との歴然の差を目の当たりにした悔しさの方が勝っていました. 銀メダルというと一見あとちょっと頑張れば金メダルだと思えるかもしれませんが,そこには大きな壁があり,コンペ期間中一度も金メダルにかすりそうな瞬間はありませんでした.さらにコンペ終了後に公開された上位の解法にざっと目を通しましたが,どんな頭しているんだと思うくらいすごくてさらに悔しい気持ちがこみ上げてきました.ただ,こういった悔しい経験がどんどん自分を強くしていくので勝てなくても参加するのが大事だと思いました.これからコンペの上位解法をしっかり読んで勉強します.
今後は,自分もコンペで優勝争いをしたいし,いずれは最高位であるGrandmasterにもなりたいです.そのためにもまずは次のコンペでもメダルをとってKaggle Expert(銅メダル以上2個)を目指します.
長くなりましたが最後まで読んでいただきありがとうございました.Kaggleやっていきましょう!

グラフィカルモデルに基づく因果探索手法の調査

最近,因果推論や因果探索に興味を持ち,勉強している.というのも最近,ゆううきさん と一緒に分散システムの異常の原因を即時に診断するための研究を進めている.原因を診断するためのアプローチとして,サーバやコンテナ等から取得できる様々なメトリック(CPU使用率やメモリ使用率など)を(グラフ理論における)ノードとして,因果グラフを構築することを考えている.メトリック同士の単なる「相関」ではなく,結果と原因の関係である「因果」を捉えようとするアプローチである.例えば,システムの障害が発生した場合,相関だけでは,AとBが関連がありそうというところまでしか言えないが,因果を特定できると理想的には,Aの原因はBであるといった議論ができるため,有用だと考えている.
実際に,前述のような因果グラフを構築して障害の原因を特定しようというアプローチは,以下の例に挙げるようにここ数年で増えている印象がある.

今回の記事は,因果探索手法の中でもグラフィカルモデルに基づくアプローチに絞って,まとめていく.以下のサーベイ論文がよくまとまっていたので全面的に参考にした.

Review of Causal Discovery Methods Based on Graphical Models(2019)

グラフィカルモデル

グラフィカルモデルとは,確率モデルをグラフを用いて記述したものである.ここでいうグラフとは, x軸, y軸があるグラフではなく,グラフ理論のことであり,ノード(頂点)とエッジ(枝)の集合で構成されるものである.グラフィカルモデルを用いることで,どの確率変数同士が直接的に関係しているかを視覚的に表現できる.
上記でエッジと書いた部分は,単なる線と矢印の二つの表現があり,それぞれを用いたグラフを無向グラフおよび有向グラフと呼ぶ.変数間の因果を表現する際に利用されるのは主に後者の有向グラフである.有向グラフの代表的な例としてベイジアンネットワークが挙げられる.

因果探索

因果探索とは,データから因果関係(因果グラフ)を推定する手法である.前述のグラフィカルモデルで考えると,考慮すべき確率変数(ノード)のデータは手元にある状態で,因果の構造,すなわちどのようにエッジが繋がっているかを推定するということになる.多くの分野において,データの根底にある因果関係を見つけ出し,それを利用することは非常に重要であるため,因果探索は,ゲノミクス,臨床医学,疫学など幅広い分野に応用されている. 次から,因果探索の具体的なアプローチについてまとめていく.

制約ベースの因果探索

制約ベースのアプローチでは,まず,観測変数にどのような条件付き独立性が成り立つかを推測する.次に推測された条件付き独立性を制約として,それを満たすような因果グラフを探索する.ここでは,代表的な制約ベースのアルゴリズムとして,PCアルゴリズム[1]を説明する.

PCアルゴリズム

PCアルゴリズムの処理の流れを以下の図を用いて説明する.前提として,PCアルゴリズムでは,未観測共通原因(潜在的な交絡因子)がないことを仮定している.

f:id:hirotsuru314:20201007095925p:plain
PCアルゴリズムによる因果グラフの推定

  1. 完全無向グラフを構築する.
  2. 条件付き独立性が成立するエッジをグラフから除外する.
    ここで,条件付き独立性を検索する際は,条件として与える変数の数を0から1ずつ増やしながら,条件付き独立性が成立しうる全ノードの組合せを対象に行われる.
  3. V-structureに基づきエッジの方向を決定する.
    V-structureでは,三つの変数A,B,Cが「A-B-C」でA-C間にエッジがなく,A,Cを条件付き独立とする条件ノード群にBが含まれていないならば,「A→B←C」の関係が成立するとして,無向グラフ中の一部の因果の方向を推定する.
  4. オリエンテーションルールに基づきエッジの方向を決定する.
    オリエンテーションルールは有向非巡回グラフ(DAG)の定義から導かれるグラフの構造に関するルールである.具体的なルールについては「Estimating High-Dimensional Directed Acyclic Graphs with the PC-Algorithm」にまとめてある.

PCアルゴリズムは,一般にノード数が増加すると上記2の条件付き独立性を検索する際のノードの組合せの数が増大し,計算時間がかかるという問題がある.これに対して,条件付き独立性の検定をGPUで並列化して高速に処理するといった提案[2]もあるようだ.また,PCアルゴリズムは未観測共通原因がある場合のFast causal inferece(FCI)アルゴリズム[3]などへ拡張されている.

スコアベースの因果探索

スコアベースのアプローチでは,同じ条件付き独立性を与える因果グラフの集合であるマルコフ同値類ごとに,モデルのよさを測るスコアをつける.ここでは,代表的なスコアベースのアルゴリズムとして,Greedy Equivalence Search(GES)アルゴリズム[4]を説明する.

GESアルゴリズム

GESアルゴリズムでは,PCアルゴリズムのように完全無向グラフから始めるのではなく,空のグラフから始まり,必要なエッジを追加したり,不要なエッジを削除していきながら因果グラフを構築する.アルゴリズムの各ステップでは,グラフに有向エッジを追加することで,ベイズ情報量規準(BIC)のようなスコアや統計的仮説検定のZスコアで計測される適合性が向上するかどうかが決定され,適合性を最も向上させるエッジが追加される.スコアをそれ以上改善できなくなった場合,GESアルゴリズムは,どのエッジを削除するとスコアが最も改善されるかをエッジごとに検証していく.このような手順でスコア(例えばBIC)を用いて因果グラフを推定していく.

構造方程式モデルに基づく因果探索(非ガウス性・非線形性の活用)

先に挙げた制約ベース・スコアベースの因果探索は,データの構造に対してできるだけ仮定をおかずに,データのみから因果グラフを推測するアプローチである.これは,後述する構造方程式モデルにおける関数形と誤差変数の分布に仮定をおかないノンパラメトリックアプローチに分類される.仮定をおかないため様々なデータに適用できるという利点があるが,因果の構造を正確に推定できないケースも多い.そこで,分析者の事前知識を仮定としてモデルに取り入れるアプローチがある.ここでは,仮定を取り入れる上で重要な概念である構造方程式モデルについて述べ,非ガウス性の仮定を利用したLiNGAMという手法について述べる.

構造方程式モデル

構造方程式モデルは,データの生成過程を記述するツールである.ここで, x yという二つの変数について,以下のような単純な因果グラフを考える.

f:id:hirotsuru314:20201007131606p:plain
因果グラフのサンプル

このような因果グラフにおける構造方程式を以下のように表す.

 \displaystyle
x=e_x \\
y=f_y(x, e_y) \\

ここで, x yは観測変数で, e_x e_yは未観測の変数であり,誤差変数という. e_y yの値を決定するために寄与しうる x以外の全ての変数をまとめて表した変数という意味をもつ.因果グラフをこのような方程式で表現した場合,分析者は以下の2点に仮定を導入することができる.

  1. 関数 f_yの仮定(例えば線形か非線形か)
  2. 誤差変数 e_x e_yの分布の仮定(例えば,ガウス分布か非ガウス分布か)

上記の二つに仮定を導入するものをパラメトリックアプローチ,関数には仮定をおく一方,誤差変数の分布には仮定をおかないものをセミパラメトリックアプローチと呼ぶ.ここでは,代表的なセミパラメトリックアプローチであるLiNGAMという手法を紹介する.

LiNGAM

LiNGAMはLinear Non-Gaussian Acyclic Modelの略であり,関数形に線形性を,誤差変数の分布に非ガウス分布を仮定するアプローチである.一見,二つの仮定をおいているようにも見えるが,非ガウス分布という仮定はガウス分布でさえなければどんな(連続)分布でも良いということからほぼ仮定をおいていないという解釈ができ,セミパラメトリックアプローチに分類されるようだ.LiNGAMについては,機械学習プロフェッショナルシリーズの「統計的因果探索」に詳しく解説されている.ここではさらっとアルゴリズムの概要をまとめる.
まずLiNGAMにおける仮定を整理する.

  • 関数形が線形
  • 誤差変数が非ガウス連続分布
  • 有向非巡回グラフ(DAG)
  • 誤差変数が独立性(未観測共通原因がないことを意味する)

p個の観測変数 x_1,x_2,...x_pに対するLiNGAMモデルは構造方程式モデルを用いて以下のように表される.

 \displaystyle
x_i=\sum_{i \neq j}^p b_{ij}x_i+e_i

これを行列・ベクトルで表現すると以下のようになる.

 \displaystyle
\mathbf{x}=\mathbf{B}\mathbf{x}+\mathbf{e}

ここで, p次元ベクトル \mathbf{x}は観測変数ベクトル, p× p行列 \mathbf{B}は係数行列, p次元ベクトル \mathbf{e}は誤差変数ベクトルである.そして, \mathbf{e}の成分 e_iは独立で,それぞれ非ガウス連続分布に従う. 係数行列 \mathbf{B}は,観測変数の分布 p(x)に基づいて一意に推定可能であり, \mathbf{B}の成分のゼロ・非ゼロパターンから因果グラフを描くことができる. 係数行列の推定には,独立成分分析(ICA)によるアプローチと回帰分析と独立性評価によるアプローチがある[5]. どちらのアプローチもcdt15/lingamで実装されている.OSSのコードも含めてLiNGAMはもう少し詳細に追いかけたい.

さいごに

今回は,サーベイ論文を参考にグラフィカルモデルを用いた因果探索の手法を整理した.ここでいくつかアルゴリズムを紹介したが,最近機械学習技術を活用した多くの手法が提案されている.例えば,深層強化学習[6],敵対的生成ネットワーク(GAN)[7]を活用したものなどがあり,まだまだ奥が深い.また,重要なイシューとして,確率変数が時系列データの場合にどのように時系列の情報を消失せずにモデリングするかというのも大変興味深い.今後も因果探索について着実に勉強していって,ブログにまとめていきたい.

参考

[1] An Algorithm for Fast Recovery of Sparse Causal Graphs
[2] cuPC: CUDA-Based Parallel PC Algorithm for Causal Structure Learning on GPU
[3] Causal inference in the presence of latent variables and selection bias
[4] Optimal structure identification with greedy search
[5] 統計的因果探索
[6] Causal Discovery with Reinforcement Learning
[7] SAM: Structural Agnostic Model, Causal Discovery and Penalized Adversarial Learning

経路積分量子モンテカルロ法を用いた量子アニーリングのシミュレータ

0または1の状態しかとれないビットを用いた古典的な計算機において量子力学の原理に基づいた量子アニーリングをシミュレートできるのだろうか.
量子アニーリングの理論について学んだ(前回の記事)のちに,このような疑問が生まれた.調べてみると,どうやら手元のPCでシミュレータを作るには経路積分量子モンテカルロ法という手法がよく用いられるようだった. そのため今回は,経路積分量子モンテカルロ法について学び,Pythonでの既存実装に触れてみることで,量子アニーリングについての理解をより深めていくことを目的とする.

目次

経路積分量子モンテカルロ法

今回,古典的な計算機を用いて量子アニーリングをシミュレートするために,経路積分量子モンテカルロ法という手法を用いる. 量子系をシミュレートするために,まず量子ゆらぎを含む横磁場イジングモデルの分配関数 Zを,対応する古典イジングモデルの分配関数で表す. そのとき,これから示す計算を追っていくと,横磁場イジングモデルの分配関数が1次元高い古典の分配関数に対応することがわかる.この結果を使って量子系の平衡状態における性質を古典系のマルコフ連鎖モンテカルロ法で調べる,というのが量子モンテカルロの基本的な考え方である.

以下,理論的な背景を少し詳細に見ていく.
スピン数が Nの場合,横磁場イジングモデルハミルトニアンは以下のように表される.

 \displaystyle
\begin{eqnarray}
\hat{H} & = & \hat{H_0}+\hat{H_q} \\
\hat{H_0} & = & -\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_i^z\hat{\sigma}_j^z-\sum_{i=1}^Nh_i\hat{\sigma}_i^z \\
\hat{H_q} & = & -\Gamma\sum_{i=1}^N\hat{\sigma}_i^x \\
\end{eqnarray}


ここで, \hat{\sigma_i^z}, \hat{\sigma_i^x}はパウリ行列の z,x成分であり,以下の行列で表される.

 \displaystyle
\hat{\sigma_i^z} = 
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix}
,
\hat{\sigma_i^x} = 
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix}


このハミルトニアンを用いて分配関数は以下のように表せる.

 \displaystyle
\begin{eqnarray}
Z & = & Tr[e^{-\beta\hat{H}}] \\
& = & Tr[e^{-\beta(\hat{H_0}+\hat{H_q})}] \\
& = & \sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid{e^{-\beta(\hat{H_0}+\hat{H_q})}}\mid\boldsymbol{\sigma_1}\rangle \\
\end{eqnarray}


ここで, \beta は温度の逆数(逆温度)である.
Suzuki・Trotter分解により,

 \displaystyle
e^{-\beta(\hat{H_0}+\hat{H_q})} = \lim_{M\rightarrow\infty}\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M


と表せることを用いると,

 \displaystyle
Z = \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M\mid\boldsymbol{\sigma_1}\rangle


ここでM個の積の間に以下の恒等演算子を挟み込むと,

 \displaystyle
\boldsymbol{1}  =  \sum_{\boldsymbol{\sigma_k}}\mid\sigma_k\rangle\langle\sigma_k\mid (k = 1, 2, ..., M)


分配関数は以下のように表される.

 \displaystyle
\begin{eqnarray}
Z & = & \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_k},\boldsymbol{\sigma_k'}}\langle\boldsymbol{\sigma_1}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_1'}\rangle\langle\boldsymbol{\sigma_1'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_2}\rangle \\
& \times & \langle\boldsymbol{\sigma_2}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_2'}\rangle\langle\boldsymbol{\sigma_2'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_3}\rangle \\
& \times & ・・・ \\
& \times & \langle\boldsymbol{\sigma_M}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_M'}\rangle\langle\boldsymbol{\sigma_M'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_1}\rangle
\end{eqnarray}


ここで, \mid\boldsymbol{\sigma_k}\rangle N 個のスピン系の直積空間を表す.

 \displaystyle
\mid\boldsymbol{\sigma_k}\rangle=\mid\boldsymbol{\sigma^z_{1,k}}\rangle\otimes\mid\boldsymbol{\sigma^z_{2,k}}\rangle\otimes・・・\otimes\mid\boldsymbol{\sigma^z_{N,k}}\rangle


分配関数の式から, k = 1, 2, ..., M は元の問題の持つ変数が定義された空間とは別の独立な変数のインデックスであり,トロッタ方向または虚時間方向と呼ばれる.
下図はSuzuki・Trotter分解を行なった際にトロッタ方向に次元が追加されている様子を示した概略図である.

f:id:hirotsuru314:20191210223523p:plain:w500
Suzuki-Trotter分解の概略図

これまでの結果から,分配関数を求めるには,以下の2つの値を求めればよいことがわかる.

 \displaystyle
\begin{eqnarray}
\langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle \tag{1} \\
\langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle \tag{2}
\end{eqnarray}


 (1) \hat{H_0}が対角行列であることから,

 \displaystyle
\begin{eqnarray}
\hat{\sigma}_i^z\mid\sigma_k\rangle & = & \hat{\sigma}_{i, k}\mid\sigma_k\rangle \\
\langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle & = & exp\biggl[\frac{\beta}{M}\biggl(\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}+\sum_{i=1}^Nh_i\hat{\sigma}_{i, k}\biggl)\biggl]\prod_{i=1}^N\delta(\hat{\sigma}_{i, k}, \hat{\sigma}_{i, k}')
\end{eqnarray}


となる.次に (2) は以下の式が成立する.

 \displaystyle
\begin{eqnarray}
\langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle = \biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\gamma}{M}\biggl)\biggl]^{\frac{N}{2}}exp\biggl[\frac{1}{2}\log \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1}\biggl)\biggl] \\
\end{eqnarray}


よって,分配関数 Zは以下のように表すことができる.

 \displaystyle
Z = \lim_{M\rightarrow\infty}\biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\Gamma}{M}\biggl)\biggl]^{\frac{NM}{2}}\sum_{\{\hat{\sigma}_{i, k}=\pm1\}}exp\biggl\{\sum_{k=1}^M\biggl[\sum_{i, j=1}^N\frac{\beta J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\frac{\beta h_{i}}{M}\hat{\sigma}_{i, k} \\
+\frac{1}{2}\ln \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}\hat{\sigma}_{i, k+1}\biggl)\biggl] \\


この結果から横磁場イジングモデルの分配関数は対応する磁場なしイジングモデルに1次元トロッタ方向に次元を追加したもので表されることがわかる.対応するハミルトニアンは以下の式で表される.

 \displaystyle
H_{eff} = -\sum_{i, j=1}^N\sum_{k=1}^M\frac{J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\sum_{k=1}^M\frac{h_{i}}{M}\hat{\sigma}_{i, k}-\frac{1}{2\beta}\ln \coth\biggl(\frac{\beta\Gamma}{M}\biggl)\sum_{i=1}^N\sum_{k=1}^M\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1} \\


シミュレーテッドアニーリングとの対比

これまでの計算から,横磁場イジングモデルの分配関数は対応する古典イジングモデルの分配関数に対して虚時間方向に次元を一つ追加したもので表せることがわかった. このことから,量子アニーリングのシミュレータは,一見単にトロッタ数 Mで並列化したシミュレーテッドアニーリングと同じもののように感じるかもしれないが,上記で計算したハミルトニアンを読み解いてみると,そうでないことがわかる.
ハミルトニアン H_{eff}の第1・2項を見ると,同じ層内( kが同一)での相互作用のみであるが,横磁場の効果を表す第三項では,層間の相互作用を持っている.( k k+1の相互作用)このことから時間発展とともに層同士が互いに独立ではなく,相互に作用し合うため,単にトロッタ数分の並列処理をしているわけではないようだ.
どのような相互作用でそれがどう基底状態探索に効いているのかなどについてはまた別途追っていきたい.

計算の流れ

計算の全体的な流れは以下のようになる.

  1. パラメータの初期化
    ・温度  T
    ・横磁場の強さ  \Gamma

  2. スピンの初期化  M層のそれぞれ N個のスピンをランダムに初期化する.

  3. ランダムに層と層中のスピンを選択し,スピンの向きを反転させる.

  4. スピン反転前後のエネルギーの差分を \Delta E とすると,確率 min(1, e^\frac{\Delta E}{T}) でスピンの向きの反転を受け入れる(受け入れない場合は棄却して変更なし)

  5. 任意の回数3,4のステップを繰り返す(モンテカルロステップ)

  6. 横磁場の強さを  \Gamma小さくする.

  7. 任意の回数5,6のステップを繰り返す(アニーリングステップ)

  8. 各層のエネルギーを求め(横磁場の項は除く)一番エネルギーが低い層のスピン配位を解とする

シミュレータを触ってみる

GitHub上ではいくつかの既存実装が存在したが,今回は以下のものを触ってみることにした.

github.com

必要なライブラリ環境をそろえて,リポジトリのトップディレクトリで以下のように実行すると,ヘルプが見れる.ちなみにPython 2系で実装されている.

$ python solve.py qmc -h
usage: Ising Solvers qmc [-h] [-problem PROBLEM] [-steps STEPS] [-dump DUMP]
                         [-spins SPINS] [-P P] [-T T] [-G0 G0] [-Gf GF]
                         [-e0 E0] [-ef EF]

optional arguments:
  -h, --help        show this help message and exit
  -problem PROBLEM  a file containing problem data in three column format, see
                    README for details
  -steps STEPS      number of monte carlo steps (MCS) used, good values depend
                    on problem size
  -dump DUMP        a file to dump move acceptances to, if desired -- slows
                    down performance
  -spins SPINS      initial spin configuration, if desired
  -P P              number of trotter slices (replicas), 20 to 80 is typical
  -T T              ambient temperature in K, 0.015 is typical
  -G0 G0            initial transverse field strength
  -Gf GF            final tranverse field strength, typically near zero
  -e0 E0            initial coupling pre-factor
  -ef EF            final coupling pre-factor

多くの引数があるが,基本的には指定しなくてもデフォルト値が入るようになっていて,必要であれば実行者が上書きできるようにしてある.
マストで渡さないといけない引数は-problem であり,これは,以下のような形式のファイルのパスを渡す必要がある.(sample_dataディレクトリ配下にサンプルが置いてある)

$ head -5 sample_data/ising12.txt
random 12x12 spin ising model with solution energy of -18,972,276
1 2 -28420
2 3 17885
3 4 -146311
4 5 128784

2行以降は以下の形式で

i j J_ij

すなわち,スピン間の相互作用を示す J_{ij}を設定するためのファイルであることがわかる.今回の最適化では実用的な問題が設定してあるのではなく,単にイジングモデルのスピン間にファイルで設定した相互作用を与え,そのイジングモデル基底状態を探索するのみである. また, -dump 引数を指定しておくと,実行の過程がファイル出力される.-steps でアニーリングステップ数が指定できる.
以下のように実行すると,用いたパラメータ等の出力が返ってくる.

$ python solve.py qmc -problem sample_data/ising12.txt -dump output.txt -steps 1000

Ising Spin Glass
+----------------+-------------------------------------------------------------------+
| Field          | Value                                                             |
+----------------+-------------------------------------------------------------------+
| data file      | sample_data/ising12.txt                                           |
| description    | random 12x12 spin ising model with solution energy of -18,972,276 |
| initial config | 0xe696019ff2add3594c219fcbfd9aef9b0dd4                            |
| current config | 0x313fac59f68f16d3ec6a17d6aca964b2025d                            |
| initial energy | 1627456.0000000007                                                |
| current energy | -15899168.000000006                                               |
+----------------+-------------------------------------------------------------------+

Solver: Path-Integral Quantum Monte Carlo
+-----------+------------+
| Parameter | Value      |
+-----------+------------+
| G0        | 3          |
| Gf        | 1e-06      |
| P         | 40         |
| T         | 0.015      |
| dump      | output.txt |
| e0        | 1e-06      |
| ef        | 4          |
| steps     | 1000       |
+-----------+------------+
To replicate starting conditions, run with arguments:
    qmc -dump output.txt -G0 3 -ef 4 -P 40 -steps 1000 -T 0.015 -problem sample_data/ising12.txt -Gf 1e-06 -e0 1e-06 -spins 0xe696019ff2add3594c219fcbfd9aef9b0dd4

出力ファイルを見てみると以下のような形式で出力されている.

$ head output.txt
Path-Integral Quantum Monte Carlo
{'dump': 'output.txt', 'G0': 3, 'ef': 4, 'P': 40, 'steps': 1000, 'T': 0.015, 'Gf': 1e-06, 'e0': 1e-06}
<class 'ising.SpinGlass'>(data_file="sample_data/ising12.txt", spin_configuration="0xe696019ff2add3594c219fcbfd9aef9b0dd4")
3,16.27456,0xe696019ff2add3594c219fcbfd9aef9b0dd4,-1
3.0,15.64962,0xe696019ff2add3594c219fcbfd9aef9b0cd4,7
3.0,13.09112,0xe696019ff2add3594c219fcbfd9aef9b0ddc,28
3.0,13.02468,0xe69e019ff2add3594c219fcbfd9aef9b0dd4,12
3.0,16.31304,0xe696019ff2add3194c219fcbfd9aef9b0dd4,30
3.0,15.9128,0xe696019ff2add3594c219fcbfd9aefbb0dd4,0
3.0,14.19152,0xe697019ff2add3594c219fcbfd9aef9b0dd4,24

4行目以降の結果は,左から横磁場の強さ,全層の中から最も低いエネルギーの値,スピンの状態を16進数で表したもの,何層目のスピンを反転させたのか.
このライブラリの実装では,1アニーリングステップの中で,それぞれの層の中から1スピンずつ反転させている.すなわちモンテカルロステップ数はトロッタ数と同じである.そのモンテカルロステップごとに1行の出力がある.

ステップを重ねるにつれてイジングモデルが低エネルギー状態に向かっていることを確認するために,横軸にアニーリングステップ,縦軸はエネルギーでグラフを書くと以下のようになる.

f:id:hirotsuru314:20191216173856p:plain

アニーリングステップが進むにつれてエネルギーが低い状態に向かっていることから,基底状態を探索している様子が見て取れる.

まとめ

古典的な計算機で量子アニーリングをシミュレートする代表的な手法として,経路積分量子モンテカルロ法がある.経路積分を用いる過程でSuzuki・Trotter展開によって得られるハミルトニアンの分配関数の近似値から,マルコフ連鎖モンテカルロ法を用いることで,量子アニーリングの挙動をシミュレートすることができる. 動作としては,問題とするイジングモデルをトロッタ数分だけコピーして,それらが互いに相互作用しながら基底状態を探索していく. 今回用いたqmcのような既存実装を用いると,実用的な問題ではないものの手元のPCで簡単に量子アニーリングの挙動をシミュレートすることができる.

参考

量子アニーリングの原理について

量子アニーリングという言葉を初めて耳にした時から,ずっと気になっていた.なぜなら学生時代に材料工学を専攻していた私にとって,「量子」と「アニーリング」という言葉はどちらも聞き馴染みのある言葉だったからである. 最近になって,やっと勉強を始めることができ,これから理論的な背景から応用までしっかり学んでいきたいと思っている. せっかくなので勉強した時のメモを少しずつブログにしていく予定である.

目次

量子アニーリングとは

量子アニーリングとは,量子力学の法則を利用して,ある種の難しい問題を解くための枠組みである. その大きな特徴は,組み合わせ最適化問題に対するメタヒューリスティックな解法として利用可能な点である. 組み合わせ最適化問題とは,離散値をとる多くの変数があるときに,それらの目的関数を最小化ないし最大化する問題である. この組み合わせ最適化問題をイジング模型基底状態探索として表現し,量子力学的なゆらぎを利用して探索する方法が量子アニーリングである.

量子アニーリングは,シミュレーテッドアニーリング量子力学版ともとらえることができる.温度の代わりに量子ゆらぎを導入し,古典的な確率過程の代わりに量子力学的な状態の重ね合わせを利用している.
量子アニーリングは1998年に以下の論文で提案されたものである.

arxiv.org

量子アニーリングをハードウェアの動作原理として実装したものには,カナダのベンチャー企業D-Wave Systemsによって開発されたD-Waveマシンなどがある.

イジング模型

量子アニーリングの原理を理解する出発点となるのがイジング模型である. 量子アニーリングでは,最適化問題の目的関数を統計力学のイジング模型で表現し,イジング模型の基底状態(最低エネルギー状態)を探索することで最適解を得る.
イジング模型のハミルトニアンは以下のように表される.

 \displaystyle
\hat{H_0}=-\sum_{i\verb|<|j}J_{ij}\hat{\sigma_i^z}\hat{\sigma_j^z}-\sum_{i=1}^Nh_i\hat{\sigma_i^z}

 \hat{\sigma_i^z}パウリ行列のz成分である. また, iはイジング模型における格子点の番号, J_{ij}は相互作用の強さ, h_iは局所磁場の強さを表している.

イジング模型における格子点は,量子アニーリングマシンにおけるスピンいわゆる量子ビットに対応する.
ここで,量子ビットがもつ上向きスピン +1および下向きスピン -1という状態を以下のようにベクトルで表現する.

 \displaystyle
\mid\uparrow\rangle=
\begin{pmatrix}
1 \\
0 \\
\end{pmatrix}
,
\mid\downarrow\rangle=
\begin{pmatrix}
0 \\
1 \\
\end{pmatrix}

パウリ行列のz成分を上記の2つのベクトルに作用させると,

 \displaystyle
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix}
\begin{pmatrix}
1 \\
0 \\
\end{pmatrix}
=
\begin{pmatrix}
1 \\
0 \\
\end{pmatrix}
,
\begin{pmatrix}
1 & 0 \\
0 & -1 \\
\end{pmatrix}
\begin{pmatrix}
0 \\
1 \\
\end{pmatrix}
=-
\begin{pmatrix}
0 \\
1 \\
\end{pmatrix}


となるから, \binom{1}{0} \hat{\sigma_i^z}固有値1の固有ベクトル \binom{0}{1}固有値-1の固有ベクトルである. したがって, \hat{\sigma_i^z}を上向きスピンと下向きスピンが同時に存在する量子力学的な重ね合わせの行列表現であると解釈できる.

一般の組み合わせ最適化問題は,多くの場合,前述のようなスピン間の相互作用をもつイジング模型の基底状態探索問題として表現できる. 逆に言うと,基底状態が最適解になるように,解きたい最適化問題に合わせて J_{ij} h_iを設定することができれば,量子アニーリング最適化問題を解くことができる.

量子ゆらぎの導入

量子アニーリングでは,イジング模型の相互作用のパラメータや局所磁場のパラメータのみでなく,量子ゆらぎと呼ばれるものを導入する.  \hat{H_0}基底状態を求めるために,量子力学的なゆらぎを利用して状態探索を行うのが量子アニーリングの基本戦略である.
量子ゆらぎを導入したハミルトニアンを以下の式で表す.

 \displaystyle
\hat{H}(t)=\hat{H_0}+\Gamma(t)\hat{H_1}


 \Gammaは量子ゆらぎの制御変数である.
初期状態 t=0においては, \Gammaを十分大きく選ぶことで右辺第1項を無視でき,第2項の量子ゆらぎの項のみとなる.  \Gammaを時間とともにだんだん小さくしていって,終状態 t=Tにおいては, \Gamma(T)=0となり,最適化したい関数の項のみが残る.
ここで,初期ハミルトニアン \hat{H_1}基底状態が自明に見つかるように \hat{H_1}を選ぶことが重要である. よく用いられる量子ゆらぎとして,横磁場があり,横磁場の効果を示すハミルトニアンを以下のように表す.

 \displaystyle
\hat{H_1}=-\sum_{i=1}^N\hat{\sigma_i^x}

 \hat{\sigma_i^z}はパウリ行列のx成分である.  \hat{\sigma_i^z}固有ベクトルは,

 \displaystyle
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix}
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
=
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
,
\begin{pmatrix}
0 & 1 \\
1 & 0 \\
\end{pmatrix}
\begin{pmatrix}
1 \\
-1 \\
\end{pmatrix}
=-
\begin{pmatrix}
1 \\
-1 \\
\end{pmatrix}


 \binom{1}{1} \binom{1}{-1}の二つである.
これらはいずれも \binom{1}{0} \binom{0}{1}を同じ絶対値を持つ係数(1と1あるいは1と-1)で加えたものであるから, \hat{\sigma_i^z}は上向きと下向きのスピンが等確率で重ね合わされた状態を表している.

 \hat{H_1}基底状態は, 2^N個の状態が重ね合わされて,量子力学の意味で同時に存在する. 当然ながら,元のイジング模型のハミルトニアン \hat{H_0}基底状態はあらかじめわかっていないため, 2^N個の全ての状態の同じ確率の重ね合わせを出発点として探索を始めることは理にかなっている.

 \hat{H_1}の自明な基底状態を初期状態として選び, \hat{H}(t)で記述される系を時間に依存したシュレディンガー方程式に従って時間発展させる.

 \displaystyle
i\frac{\partial}{\partial t}\mid\Psi(t)\rangle=\hat{H}(t)\mid\Psi(t)\rangle


 \mid\Psi(t)\rangleは,各時刻における系の状態である.
このとき, \Gamma(t)の時間変化が十分ゆっくりであれば,後述する量子力学の断熱定理が適用できる.

断熱変化

ハミルトニアン \hat{H}(t)は時間変数 tに依存するが,系が断熱的に(ゆっくりと)変化していると見なせる場合には,各時刻における系の状態 \mid\Psi(t)\rangleは, tを固定パラメータとみなしたときの定常状態のシュレディンガー方程式

 \displaystyle
\hat{H}(t)\mid\Phi_t\rangle=E_0(t)\mid\Phi_t\rangle


基底状態 \mid\Phi_t\rangleに十分に近い.実際の状態 \mid\Psi(t)\rangleと定常基底状態 \mid\Phi_t\rangleいずれも規格化されているとして,これらの内積が1に十分近いことを意味している.

 \displaystyle
|\langle\Phi_t\mid\Psi(t)\rangle|^2=1-\delta^2\quad(\delta\ll1)


これは各瞬間の基底状態をほぼ正確にたどっていくということを示している.
こうして,横磁場項 \Gamma(t)\hat{H_1}により引き起こされる量子ゆらぎを \Gamma(t)を通じて適切に制御することにより,自明な初期状態から出発して,非自明な状態である最適問題の解に到達する.以上が量子アニーリングの基本的な原理である.

応用

最後に少しだけ量子アニーリングの応用について触れる.最初に述べた通り,量子アニーリングの応用先として期待されているのは組み合わせ最適化問題である.(実際には機械学習のサンプリングにも有効であるそうだ)
一般的に組み合わせ最適化問題の例としてよく登場するものは,巡回セールスマン問題やナップサック問題などが挙げられる.一方,「量子アニーリングの基礎」の中で,クラスタリングが例に挙げられていたのが個人的には興味深かった. クラスタリング組み合わせ最適化問題として考えたことはなかったが,ある点がどのクラスに属するかというのは直感的にも組み合わせであるし,異なるクラスに属する点の間の距離を最大化するという問題に読みかえれば組み合わせ最適化問題の目的関数として定式化できる. この辺りはもっと調べたいので,以下の論文などを読んで学んでいきたい.

arxiv.org

参考

GMOペパボ株式会社を退職しました

2019年7月24日が最終出社日でした。ペパボでは、ロリポップやへテムルといったレンタルサーバサービスのインフラエンジニアとして働いており、約1年半在籍していました。
今回は、ペパボでの一年半の振り返りと、転職に至った経緯などについて書いていきたいと思います。

ペパボでの一年半の振り返り

私はペパボカレッジの枠でペパボに採用されました。ペパボカレッジとは第二新卒エンジニア向け研修のことで、中途採用でも現場に配属される前にしっかり研修を受けることができます。私はエンジニア歴も浅いし、そもそもインフラの経験がなかったため、このペパボカレッジの応募はまたとないチャンスだと感じました。実際に研修では約一ヶ月間で広大なインフラの分野の知識を爆速で学んでいきます。具体的な研修内容などは以前のブログに書きました。頑張りたい人を応援してくれる本当にいい仕組みです。

研修が終わって、配属後も周りのエンジニアの方々に多大なるサポートをいただき、未経験の私でもなんとかインフラエンジニアとして仕事をすることができました。ペパボの「いるだけで成長できる環境」とはこのことだったのかと身に染みて感じました。 どれだけ成長できたかをざっと振り返ってみると、ペパボに入った時にはサーバーに入ってもlscdくらいしか叩けずにSSHって何?といった状態だった私が、入社後半年ちょいでSSHのプロキシサーバを書くまでに至ったので、おそらくかなり成長させてもらったんだろうなーと思います。

また、働く環境も私にとっては最高でした。私のチームは完全フレックスだったのですが、名ばかりのフレックスではなく、本当にフレックスでした。私は子どもがまだ小さいため、早めに出社して早めに退社したり、趣味の筋トレをしてから遅めに出社したりと本当に自分のリズムで働かせていただきました。残業時間も月に20時間を超えたことはありませんでした。スーパーホワイトな会社で心からおすすめできる会社です。実際に私の高校時代の友人をペパボに誘って、入社しました。

転職に至った経緯

以上のように、学ぶ環境も整っていて、生活リズムに合わせて働けるとても良い環境ではありましたが、一方で、ずっと心の中にもどかしい感情がありました。

私はもともとペパボの採用面接のときから、「いずれはインフラの領域に機械学習を活用する研究がしたい」と明言していました。配属後も常にそのことを目標に仕事をしてきました。インフラエンジニアになりたての時はインフラ関連の勉強でいっぱいいっぱいでしたが、最近ではコンテナのリソース割り当ての最適化の論文を読んだり、サーバの負荷状態の把握に機械学習を用いる構想について登壇をしたりなど、やりたいことの妄想は膨らんでいました。しかしながら、なかなか妄想を実行に移すことができませんでした。それは、私自身の力不足はもちろんですが、ロリポップのように長くやっているサービスの運用は非常に大変であることも要因の一つとしてあり、結局一年半で一度も機械学習を絡めた仕事をすることができませんでした。

また仕事では、私自身が自分の強みだと思っている「深く思考すること」がなかなかできず、リアクティブに急ぎの対応するような仕事が占める割合が非常に多いのが現状でした。もちろん、自分が優秀であれば、与えられた仕事をこなしながらも、やりたいことができるかもしれませんが、私にはできませんでした。そんな時に消防士からエンジニアに夢を持って転職した自分の感情が蘇り、本当にこのままでは何者にもなれないまま終えてしまいそうにだという焦りが抑えきれなくなり、転職を決意しました。

研究者の道を選んだ

2019年8月1日からさくらインターネット研究所の研究員として働かせていただくこととなりました。
正直、このタイミングでエンジニアから研究者になることに関してはかなり悩みました。それは、私が固定観念で「エンジニアとして優秀じゃないと研究者として活躍できない」と思っていたからです。
この固定観念による不安に対する一つの答えを与えてくれたのは、@matsumotoryさんが書いた下のブログでした。研究をいずれやりたいと思っているなら今すぐ研究に飛び込んだ方がいいという考えに変わるきっかけとなりました。

hb.matsumoto-r.jp

実は松本さんには、もう一年以上前に松本さんがペパボに在籍されていた時から、エンジニアとしての自分の理想と現実が乖離していることや、いずれ研究をやりたいけどまだやっていける自信が持てないことなど色々相談させていただいていました。会議室で二時間近く二人でお話させていただいたこともありました。
その時に松本さんがおっしゃった言葉のいくつかは時間をおいてからやっとその意味が自分の中で消化できて、私の人生に大きく影響を与えた言葉となりました。松本さんはペパボ時代から私の学生時代の研究(専攻は材料工学でした)の業績をしっかり見て下さっていて「過去にこれだけ研究で活躍できていた人間が分野が違うからできないなんてことはありえない」と言ってくれていました。その言葉が最終的に私に研究に飛び込む勇気を与えてくれて、研究への挑戦を心に決めました。

そんな中、現在さくらインターネット研究所の研究員をされている@yuuk1tさんや、@kumagalliumさんとお会いする機会に恵まれ、色々と研究のお話をさせていただく中で、研究に対する熱量や考え方に共感し「私も一緒に研究をやりたい」という気持ちが高まって、さくらインターネット研究所に応募させていただきました。

さいごに

お世話になったペパボの皆さんには「辞めてすみませんでした」ではなく「教えてもらったことを活かしてこれからさらに飛躍します」とポジティブな気持ちでいます!
一年半の間本当にありがとうございました!

遺伝的アルゴリズムの並列化とgoroutineによる実装

先日、「遺伝的アルゴリズムをGoで実装してみた」という記事を書きました。
この内容で2019年7月13日(土)に開催されるGo Conference'19 in Fukuokaに登壇させていただくことになったので、開発中のeagoというパッケージをもっと作り込んで行きたいと思います。
今回は、遺伝的アルゴリズム(以下、GA)の計算処理をgoroutineでサクッと並列化した話です。

概要

GAは複数の個体を用いた多点探索のアルゴリズムであるため、本質的に並列化と親和性が高い手法であると言えます。実用的な観点からも短い計算時間で良好な解を得ることが望まれているため、これまでに数多くの並列化手法が提案されています。

今回はGo言語の強みの一つであるgoroutineを使ってGAの処理の一部を並列化してみます。また、擬似問題に対してベンチマークを取り、実装した並列処理が全体の計算時間短縮に有効であることも検証しました。

代表的な並列化手法

マスタースレーブ方式

GAを並列化するための最も単純な手法は適応度計算を並列化することです。マスタースレーブ方式では下図のように、全体の制御および遺伝的操作を行うマスターと、適応度計算を行う複数のスレーブから構成されます。

f:id:hirotsuru314:20190619141143p:plain
マスタースレーブ方式

この方式は、同時に計算できる処理を並列化しただけの極めて単純かつ一般的な手法であり、GA固有の並列化処理ではありません。
今回はこちらをgoroutineを使って実装してみます。

島モデル

島モデルはGA自体を並列化するというアプローチをとっています。集団を複数の部分集団に分割することでGAを同時に実行し、並列処理を実現しています。
その概要を下図に示します。

f:id:hirotsuru314:20190619142435p:plain
島モデル

それぞれの部分集合を初期化して、評価・選択・交叉・突然変異などをそれぞれで行います。ここで本手法の重要な特徴となるのが、ある一定の条件を元に部分集団から一部の個体を選んで交換を行う点です。このプロセスにより、島モデルは並列化手法であるという側面のほかに、個体の多様性を維持するための手法としての側面も併せもちます。

goroutineによる並列化

今回はマスタースレーブ方式を実装してみます。前述の通り、マスタースレーブ方式では適応度計算を並列化するため、適応度関数(目的関数)が非常に複雑で計算処理に時間がかかる場合に有効な並列化であると言えます。
個体群に属するそれぞれの個体の適応度計算は独立して行えるため、この適応度計算をするためのスレーブとしてgoroutineを立ち上げて、並列で計算を行うようにします。マスターについては、それぞれのスレーブの計算結果を受け取るプロセスがいるというよりは、それぞれのスレーブが以下のように定義したIndividual構造体のFitnessに計算結果を書き込むことで、以降の処理で全ての結果を参照できるようにしています。

type Individuals []Individual

type Individual struct {
    Chromosome Genome
    Fitness    float64
    Evaluated  bool
}

実装にはsync.WaitGroupを用いました。
WaitGroupを用いると、独立した複数のタスクを平行で処理して、それら全ての終了を待ち合わせる処理が非常にシンプルな記述で実装できます。

func (indis Individuals) Evaluate(parallel bool) {
    var wg sync.WaitGroup
    if !parallel {
        for i := range indis {
            indis[i].Evaluate()
        }
    } else {
        for i := range indis {
            wg.Add(1)
            go func(i int) {
                defer wg.Done()
                indis[i].Evaluate()
            }(i)
        }
        wg.Wait()
    }
}

上記の Evaluate() の部分が適応度計算の関数です。
WaitGroupはカウンターのようになっていて、Addを呼び出すと引数に渡された整数だけカウンターを増やし、Doneを呼び出すとカウンターを1つ減らします。そしてWaitを呼び出すとカウンターがゼロになるまでブロックします。

パフォーマンス検証

実行時間の測定には、testingパッケージに備わっているベンチマークの仕組みを用いました。
並列化が有効かどうかを検証するために、実用的ではないですが、適応度関数の中で1 ms sleepするようにして、擬似的に計算処理に時間がかかる適応度関数を作りました。

ベンチマークのための全体コードの以下のようになります。

package main

import (
    "github.com/tsurubee/eago"
    "log"
    "testing"
    "time"
)

type Vector []float64

func (V Vector) Initialization() eago.Genome {
    return Vector(eago.InitFloatVector(2, 32, -32))
}

func (V Vector) Fitness() float64 {
    time.Sleep(1 * time.Millisecond) // Heavy process
    var s float64
    for _, x := range V {
        s += x * x
    }
    return s
}

func (V Vector) Mutation() {
    eago.AddNormalFloat(V, 0.5)
}

func (V Vector) Crossover(X eago.Genome) eago.Genome {
    return Vector(eago.BLXalpha(V, X.(Vector), 0.3))
}

func BenchmarkParallelGA(b *testing.B) {
    for i := 0; i < b.N; i++ {
        runGA(true)
    }
}

func BenchmarkNonParallelGA(b *testing.B) {
    for i := 0; i < b.N; i++ {
        runGA(false)
    }
}

func runGA(parallel bool) {
    var v Vector
    ga := eago.NewGA(eago.GAConfig{
        PopulationSize: 30,
        NGenerations:   20,
        CrossoverRate:  0.8,
        MutationRate:   0.01,
        ParallelEval:   parallel,
    })
    ga.PrintCallBack = func() {} // Do not print messages while running
    if err := ga.Minimize(v); err != nil {
        log.Fatal(err)
    }
}

このファイルを *_test.goという形式で保存して、ファイルが存在するディレクトリで以下のようなコマンドを実行することでベンチマークが取得できます。

$ go test -bench .
goos: darwin
goarch: amd64
pkg: github.com/tsurubee/eago/examples/ga
BenchmarkParallelGA-4             50      37193820 ns/op
BenchmarkNonParallelGA-4           2     821147753 ns/op
PASS
ok      github.com/tsurubee/eago/examples/ga    4.365s

適応度関数内でのsleepの時間を1・5・10 msにして、並列化あり(parallel)並列化なし(non-parallel)における処理時間をグラフにすると以下のようになりました。

f:id:hirotsuru314:20190623114700p:plain

今回のベンチマークの条件では、個体数30、世代数20であるため、適応度計算は600回程度行われ、sleepが10 msのときに処理を並列化しない場合は、適応度計算だけでも6000 msかかることになります。
このことを考えると、適応度計算が重い場合に、今回のようなマスタースレーブ方式の並列化が、全体の処理時間の短縮に有効であることがわかります。

今後

Go Conference'19 in Fukuokaまでに以下の2点について取り組んで行きたいと思います。

  • 他の並列化手法の検証(島モデルなど)

  • 実用的な問題への適用

参考

www.morikita.co.jp

www.oreilly.co.jp