tkm2261's blog

研究員(OR屋) → データ分析官 → MLエンジニア → ニートがデータ分析諸々書いてます

理研AIP主催Micheal Jordan先生最適化セミナー聴講メモ

お久しぶりです。

今日は10/20に参加した、理研AIPでのMicheal Jordan先生の講演のメモを書いていきたいと思います。

参加を許可頂いた理研AIPの方や司会の福水先生も良い会をありがとうございました。

Talk by Prof. Michael I. Jordan (University of California, Berkeley, USA) | 革新知能統合研究センター

In this post, I am writing about Prof. Micheal Jordan's speak at Riken AIP in Oct. 20th. If you have some problems about this post, please inform to @tkm2261 on Twitter.

いまはニートなので参加したセミナーの報告をする義務は全く無いのですが、やらないのもケツの座りが良くないので簡単なメモだけ書いていきます。

免責:私の拙いリスニング力に基づいて書かれているので、各自必要な箇所は調べて下さい。検索できるようにワードはなるべく残します。

間違っている箇所は積極的にコメントで教えて下さい

私の備忘のためでもあるので、不確かでも不確かなりに聴講すぐの時点のメモを残したいと思います。

タイトルと要約

webサイトの告知と異なり2つのトークを聞くことが出来ました。

  • Part I: Variational, Hamiltonian and Symplectic Perspectives on Acceleration
  • Part II: How to Escape Saddle Points Efficiently

両方共、連続最適化のお話です。

Part. I 要約

  • Nesterovの加速法を始めとした、加速法はBergman Lagrangianを考える事により微分方程式として統一的に議論できる
  • この時、全ての加速法は同じ経路を異なる時間で進んでいると解釈できる。
  • ただし、この解析は経路が連続(おそらくステップサイズが無限小取れる)時の議論で、これを単に離散化(discretization)した時は精度保証ができず、実験的にも不安定だった。
  • コレに対してSympletic integrationを用いて離散化した場合を解析する、Nesterovの加速法よりも速く、先行研究より安定したアルゴリズムが構築出来た?

最後の方は早くて少し聞き取れていません。

少し前の資料ですが、スライドはこれっぽい

https://www.sigmetrics.org/sigmetrics2017/MI_Jordan_sigmetrics2017.pdf

Part. II 要約

  • PCA、Tensor分解など、局所解に悩む場面は多い
  • これらの学習において、如何に鞍点を抜けるか、如何に効率よく鞍点を抜けるかが重要になってくる
  • GDやSGDでの鞍点に対する研究はいくつかある
    • GDは漸近的に局所解を抜ける?
    • GDが局所解を抜けるのに指数時間かかる場合がある。
    • SGDは鞍点を多項式時間で抜ける?
  • 彼らの研究はPerturbation Gradient Descent(PGD)がSecond Order Stationary Point(SOSP)に収束する事を示した。
  • (SOSPは、鞍点でなく停留点である度合いを2次まで考慮したものと思われ、First Order Stationary Pointへの収束を示した先行研究のよりも、より鞍点を避けて、より停留点と思われる点に収束できると思われる。)

おそらく↓が元ネタなので参照下さい。

[1703.00887] How to Escape Saddle Points Efficiently

講演メモ

初めに

  • 最適化とコンピュータ資源にはトレードオフがある

    • 凸緩和でのエラー
    • 並列実行でのエラー
    • 最適化オラクルの性能
    • community complexity?
    • subsampleing
  • 近年のデータ解析の分野でも重要な問題

  • アルゴリズムを開発するだけでなく、下界を得ることは様々な場面で有用。

Part I: Variational, Hamiltonian and Symplectic Perspectives on Acceleration

  • 加速勾配法はsublinearのオーダーで収束する

  • Nemirovski先生などが、90年台により優れた収束性を持つアルゴリズムはないかと研究。

  • 一般的にモーメント法として知られている

  • 加速法がどういった最適化オラクルにはいるのだろうか

  • まだまだ加速法はミステリアスで解析されていない。

  • 加速法の解釈にはいくつか研究がある

    • Chebyshev polynomial
    • Linear Coupling
    • etc...
  • しかし依然として謎があり、一体加速法の裏には何があるのだろうか?何故、一度下がって上がった方が収束が速いのだろうか?

  • 加速法を連続(continious time perspective)(恐らくステップサイズが無限小取れる)を考えると、Boyd先生の研究により微分方程式として定義できる。

  • Bregman Lagrangianを導入して、最適化問題をBergman Lagrangianの時間で積分した値の最小化として考える。

  • この時、全ての加速法は同じ経路を異なる時間で移動する手法と考える事ができる。(All accelerated methods are traveling the same curve in space-time at different speeds)

  • p=2のときはNesterovの加速法で、P=3のときはaccelerated cubic-regularized Newton’s methodになる。

  • ただこの解析は、連続の場合であり離散化(discretization)をする必要がある。

  • ただしそのまま離散化すると、収束性の理論保証がなく、経験的に不安定なアルゴリズムとなる。

  • そこでSympletic integrationを用いる。物理の分野では有名

  • (ただし、この手法においてはconverce vs quantitiesのトレードオフがあるらしく、スライド無いでは『Is it possible to find discretiztion whose solutions exactly conserve these same quantities? YES!!』と言っていたが、分かる方解説求む)

  • Bergman Hamiltonian上でSympletic integrationを考えるとNesterovの加速法よりも(実験的に)速いアルゴリズムが開発できた。

  • Gradient法と組み合わせると、より安定したアルゴリズムとなった。

Part II:How to Escape Saddle Points Efficiently

  • 悪い局所解はいつも問題である。

  • 如何に鞍点を抜けるか、如何に効率よく鞍点を抜けるかが重要である

  • GDやSGDでの鞍点に対する研究はいくつかある

    • GDは漸近的に局所解を抜ける?
    • GDが局所解を抜けるのに指数時間かかる場合がある。
    • SGDは鞍点を多項式時間で抜ける?
  • PCA、Tensor分解など、局所解に悩む場面は多い

  • 1980年にNesterovが勾配降下法(GD)がFirst Order Stationary Point(FOSP)に収束することを示した。

  • 本研究はPerturbation Gradient Descent(PGD)がSecond Order Stationary Point(SOSP)に収束する事を示した。

  • (FOSPやSOSPは鞍点でなく停留点である度合いを示していると思われ。Second orderだとより停留点っぽいところに収束できる事が示せると思われる)

  • 幾何学的な解釈をすると、高次元の場合、鞍点付近では取れる値の幅がとても薄くなる。(恐らく勾配がなだらかすぎて鞍点判別&抜け出すのが難しい事を言っていると思う。わかった方解説求む)

  • この薄さが問題を難しくしている。

  • 先行研究ではoverdamped Langevin diffusionを用いて収束に必要な反復回数がO(d/¥epsilon2)だったのに対し、本研究ではoverdamped Langevin diffusionを用いてO(¥sqrt(d)/¥epsilon)を達成した。(dは次元、イプシロンは目標精度)

Githubを整備した

Kaggleの汚いコードばっかなので、Bitbucketのコードを公開するのを躊躇ってたが、

このご時世githubにコードが上がっていないのもリスクなので諸々上げました

ニートになって3ヵ月経ったので、ニートについて書く

皆様、お久しぶりです。

以前のエントリでも触れましたが、絶賛ニートを楽しんでおります。

ニートになって早3ヵ月ですが、率直な感想としては、

世界がこんなに輝いていたとは!

っといったところです。

将来を犠牲にした束の間の天国かもしれませんが、体力ある20代に自分の時間を取れるのは最高です

5年もエンジニアしていたら勉強したいことが貯まるのは当たり前なので30歳前後で1年ぐらいニート期間は超おすすめです。

ポエムはエントリの最後の方にするとして、ニートのなり方とニートの私流楽しみ方を共有してニートの布教に努めたいと思います

手続き編

健康保険

健康保険は大きく分けて3つ戦略があると思います。

  1. 任意継続
  2. 国民健康保険
  3. 誰かの扶養に入る

結局私は国民健康保険に入りました。理由としては、

  • 半年働いてたので扶養には入れず
  • ニートが1年以上を想定したので、来年以降に国保が所得ゼロで安くなる
  • 自治体のトレーニング室無料券が付いていた

任意継続と国保は金額大して変わらないことが多いみたいなので、半年以上ニートするなら国保で良いかなぁと思ってます

あと国保の加入は資格喪失日以降でないと自治体で手続き出来ないので注意。

退職日後になるはやで役場に行きましょう。手続きは一瞬です。

年金

年金は特に選択肢ないですが、健康保険と一緒に役場で手続きします。これも一瞬で終わります。

確定拠出年金に加入していた場合は個人向けに切替える必要がありますので手続きしましょう。

住民税

振込用紙が送られて来るので振り込みましょう。一括納付するとかなり通帳削られるので、メンタルを強く持ちましょう。

失業給付

失業給付の申請だけは役場で出来ないので、近くのハローワークに行きましょう。

早く行ったほうが給付が早いのと、申請しておけば途中で就職しても残りの支給額の6割ぐらいの再就職手当をもらうことが出来ます。

次で諸々紹介します。

ハローワーク(失業給付)編

STEP:1 書類の準備

とりあえずハロワに行きましょう。持っていく物は

  • 会社から送られてきた書類全部
  • 身分証明書
  • マイナンバー ※
  • 顔写真(3x2.5) ※
  • 通帳かキャッシュカード ※

※が付いているの物は最初なくてもなんとかなるので、無くてもとりあえずハロワに行きましょう

自己都合退職の場合は3ヵ月の給付制限期間があるので、支給開始までに書類が揃えばokです。

STEP:2 ハロワに行く

とりあえずハロワに行って受付すれば、よしなに案内してくれます。

書類書いて、求職票書いて、失業給付説明会の日程決めてで終わりです。

STEP:3 失業給付説明会に行く

自治体のどこかの大ホールでやります。

特に大きい情報はありませんが、この説明会参加も就職活動1回としてカウントできるので、そのあたりの説明は聞いておきましょう。

STEP:4 初回認定日に行く

手続きして1ヵ月後に初回認定日があります。時間も指定されて、これに行かないと失業給付もらえないので必ず行きましょう。

STEP:5 就活する(ハロワで求人検索)

3ヵ月の給付制限期間中に最低2回は求職活動しないといけません。

他にも求職活動と認められる行為がありますが、ハロワで求人検索するのがてっとり早いです。

注意が必要なのは、PC終わったあと、ハンコを貰わないと実績にならないので注意しましょう。ラジオ体操と思って下さい。

STEP:6 認定日に行く

給付制限が明けた認定日に行くと、自己都合退職でも晴れて給付がもらえます。

勉強編

Kaggle

分析官や機械学習エンジニアとして働いていたので、無職期間中に資格的な何かが欲しくて、Kaggle頑張りました。

Kaggle Masterになりました - tkm2261's blog

一応masterになれば、どこに出しても通用するっぽいので良かったです。

Master要件のGoldメダルはTop10ぐらいに入らないともらえないので、大変ですが見合うリターンはあると思ってます。

ぶっちゃけKaggleが出来るからといって分析官や機械学習エンジニアが務まるわけれはありません。

ただこの界隈は必要スキルセットが曖昧なので有象無象の外野を説得するのに肩書があったほうが便利です。

分析出来るとか、レポーティングできるとか、開発出来るとかは実際に出来ればよく、政治に時間を使わないで済むに越したことはないです。

競プロ

前職で競プロガチ勢たちがひしめき合ってたので折角なので始めました。

最初は怖いところかと思ってましたが、初心者はdivision 2から始まるので安心して始められます。

Div. 2は独断と偏見ですが、こんなイメージなのでpythonが少しかければ参加できます。

  • 1問目: if文とfor文しってる?
  • 2問目: 全列挙探索できる?
  • 3問目以降: 優先度付きキューとかDP使える?グラフしってる?

とりあえず私はcodeforcesで3問解けると喜んでいます。

幾つかサイトを紹介しています。

  • TopCoder

    • 言わずと知れた競プロの代表サイト
    • ただ今は活発でなくcodeforcesが一番人気?
    • SRMに出るとレートが変動していく。
    • クライアントソフトを入れるまでが一苦労だが、topcoder-greedyを入れるとテストコード付きテンプレが自動で生成されるので参加しやすい。
  • codeforces

    • たぶん今一番活発な競プロサイト?
    • 週1ぐらいでコンテストが開催される。
    • topcoder-greedyみたいなツールはないので、自分でテンプレ作る必要があります。他人のパクりましょう

その他にも、日本のAtCoderやCS AcademyやCodeChefあたりが有名らしく、余裕があればやってみたいです。

勉強

現在は、強化学習系を学んでいて、DeepMindが公開したStarCraft2のAIをなんとか動かせないかやっています。

語るほどわかってきたら何処かに書きます

英語

英語が苦手すぎて東工大を選んだ私としては鬼門ですが、このご時世英語なしには無理になってきたので頑張っています。

TOEFL100点を目標に頑張っています、今にして思うとIELTSの方が良かったかもしれません。

TOEICは日本でしか通用しないので止めました。ただTOEFL一回230ドルは痛い。。。

実はニートになる前の2年前からTOEFLには取り組んでいて、現在スコアが92点までは来ています。大学時代のTOEICは640点でした。

最初はListeningが呪文にしか聞こえず、体質的に英語無理なんじゃないかと真剣に悩んだ時もありましたが、2年経つと8割は分かるようになったので継続は偉大です。

Writingはフィリピンで一ヶ月間にエッセイを1日3本書いてから、書くのが苦になくなりました。ただTOEFLのWritingスコア安定しないので精進してます

SpeakingもLimitedの評価を貰うセクションがあり精進してます。Speakingに関しては、何より瞬間英作文が効きました。

https://www.amazon.co.jp/dp/4860641345

これをやると、言いたいことを日本語で思いついてから1秒で英語になるので会話を成立させることが出来ます。そのうち無意識に出てくるようになります

とりあえず英語が出来るようになりたい人は、これを毎日やるのが良いと思ってます。

健康編

運動の習慣をつける

ニートは下手すると1日1歩も外に出ないので、最低限通勤で使用していた体力すら低下しかねないので、意識的に運動を入れたほうがよいです。

私の場合はこんな感じで、どれかを毎日やるようにしてます

  • ランニング
    • 3kmぐらいの短いコース
    • ダルい日は最低限ランニングぐらいはしようと思える短いコース設計
  • 自転車
  • ジム
    • 後述

ジムの使い方を覚える

ニートになって良かった事の上位としてジムの使い方を覚えた事があります。

ランニングマシンやバイクもあって運動もできるので、外出嫌いでもとりあえずジムに行ってしまえば自然と運動するので良いです。

私のようなヒョロガリは最初ジムに行くのが怖かったので、いくつかTipsを書いておきます

  • ベンチプレスをやろう
    • ジムに来てベンチプレスやるとテンションが上がるのでオススメです。ジム来た感が出ます。
    • YouTubeでフォームを調べてから行きましょう。
    • 初回はジムの人に使い方を聞きましょう
    • 棒だけの20kgからスタートして徐々に上げて行きましょう
    • すぐ終わる初心者なんて誰も気にしないので堂々行きましょう
  • ダンベル中心のメニューを組む
    • ベンチプレスとか器具系は空いていないことがあるので、ダンベルをメインにする。ヒョロガリには十分な負荷なはず
    • ダンベルは重さを変えれば空いていない事はないので、同じメニューが安定して出来て良い
    • スクワットは脂肪燃焼効果が高くて良いみたい
  • 筋トレ -> 有酸素運動(ラン・バイク)
    • これで1時間から2時間ぐらいは潰せる
    • 最後に風呂入っていくと一仕事した感が出る
  • 最初は月額でなく毎回払い
    • 最初はそんなに行かないので月額はもったいないかも
    • 自治体のトレーニング室が安くておすすめ

メンタル編

ニートなので何かやられるワケでは無いですが、ニート故に生活リズムを崩して遅起きしてボーっとして一日終わって自己嫌悪みたいなのが多発しました。

解決策としては、体がダルイと思ったら運動する!を実践するのが大事

運動や筋トレは思考をポジティブにしてくれます。

この習慣が身につくだけでもニートやるかいがあるかと思います

3ヵ月終わってみて

ここまでは非常に日々を楽しんでいます。

そもそも退職した理由が、会社の仕事と自分のやりたい事の方向が合わなかった事なので、退職してやりたい事が出来てれば楽しくわけないかなぁと思ってます。

無職期間を挟む事は、かなりリスキーなので万人に勧められるワケではないし、私も辞める前は相当悩みました。

ただ新卒から転職を挟んで、自分のスキルや感覚が大きく間違って無いことがわかってきて、なあなあで仕事するよりまとまった時間で勉強したほうが将来のキャリアの伸びが良さそうだと考えて実行してしまいました。

それに次に就職するときに苦労するようなら、変なイキリがとれてそれはそれで良い気がしています。

ニートになるのはかなり決断が必要ですが、

エンジニアは30歳前後で1年ぐらいニートしたほうが良いと思ってます。

Kaggler Slack作りました

Twitterでは拡散しましたが、検索で来る人もいそうなので宣伝

https://kaggler-ja.herokuapp.com/

オープンチャンネルなのでメール入れれば誰でも入れます

誰もチェックとかしてないので気軽にどうぞー

Kaggle Masterになりました

記念に残しておく

f:id:tkm2261:20170822102936p:plain

https://www.kaggle.com/tkm2261

細かすぎて伝わらないLightGBM活用法 (callback関数)

皆様tkm2261です。この頃連投が続いてますが、

最近まで参加していたInstacart Market Basket Analysis | Kaggleで色々やったので残しておこうと思います。

このcallback関数は便利ですが、Kaggleなどでヘビーに使う人以外ここまでしないと思うので活躍するかは微妙なところです。。。

LightGBMのtrain関数を読み解く

xgboostもそうですが、lightgbmにもtrain()という関数がありLightGBMユーザはこれを使って学習を実行します。

scikit-learn APIも内部ではこの関数を呼んでいるので同じです。

この引数にcallbacksというのがあり、殆どのユーザは使っていないと思います。(私も今回初)

このcallbacksの活用法が今回のトピックになります。

train関数の実装を見てみると, 192行目ぐらいでBoostingを回しています。

LightGBM/engine.py at master · Microsoft/LightGBM · GitHub

# lightgbm/engine.py: 192
    for i in range_(init_iteration, init_iteration + num_boost_round):
        for cb in callbacks_before_iter:
            cb(callback.CallbackEnv(model=booster,
                                    params=params,
                                    iteration=i,
                                    begin_iteration=init_iteration,
                                    end_iteration=init_iteration + num_boost_round,
                                    evaluation_result_list=None))

        booster.update(fobj=fobj)

        evaluation_result_list = []
        # check evaluation result.
        if valid_sets is not None:
            if is_valid_contain_train:
                evaluation_result_list.extend(booster.eval_train(feval))
            evaluation_result_list.extend(booster.eval_valid(feval))
        try:
            for cb in callbacks_after_iter:
                cb(callback.CallbackEnv(model=booster,
                                        params=params,
                                        iteration=i,
                                        begin_iteration=init_iteration,
                                        end_iteration=init_iteration + num_boost_round,
                                        evaluation_result_list=evaluation_result_list))
        except callback.EarlyStopException as earlyStopException:
            booster.best_iteration = earlyStopException.best_iteration + 1
            evaluation_result_list = earlyStopException.best_score
            break

意外に思われるかもしれませんが、boostingのfor文はpython側で回ってます。

early stoppingもこっちで例外として書いてあるのでcallback関数さえかければループ毎にかなり動的に書くことが出来ます。

例えばearly stoppingの条件をより複雑にしたり、loggingをちゃんと仕込んで普通は標準出力に吐かれて保存が難しい学習の様子をファイルに残すことが出来ます。

train関数に渡したcallback関数はcallbacks_after_iterに渡るのでこちらをこれから見ていきます。

callback関数の仕様

custom objectiveやcustom metricと異なりcallbackの仕様は英語でもドキュメントがありません。

ただ実装は素直なので見ていきます。

文字で説明してもアレなので、先に私のコンペの実装を見せます。

def get_pred_metric(pred, dtrain):
    """予測値無理やり取るmetric
    """
    return 'pred', pred, True

def callback(env):
    """ligthgbm callback関数(instacartコンペ)
   
    :param lightgbm.callback.CallbackEnv env: 学習中のデータ
    """

    # 10回毎に実行
    if (env.iteration + 1) % 10 != 0:
        return

    clf = env.model                 # 学習中モデル
    trn_env = clf.train_set         # 学習データ (ラベルだけが入ってる※後述)
    val_env = clf.valid_sets[0]     # 検証データ (ラベルだけが入ってる※後述)

    # 無理やり予測値をとる
    preds = [ele[2] for ele in clf.eval_train(get_pred_metric) if ele[1] == 'pred'][0]
    # ラベル取得Cython用に型を指定
    labels = trn_env.get_label().astype(np.int)

    # ユーザ毎のしきい値で予測を0-1にしてデータ毎に正解なら1, 不正解なら0を返す関数
    # (list_idxはgroupでグローバルでアクセス)
    res = f1_group_idx(labels, preds, list_idx).astype(np.bool)
   
    # 間違ってるデータのウェイトを上げるて正解のデータのウェイトを下げる
    weight = trn_env.get_weight()
    if weight is None:
        weight = np.ones(preds.shape[0])
    weight[res] *= 0.8
    weight[~res] *= 1.25

    trn_env.set_weight(weight)

    # 検証データのmetricを計算
    preds = [ele[2] for ele in clf.eval_valid(get_pred_metric) if ele[1] == 'pred'][0]
    labels = val_env.get_label().astype(np.int)
    t = time.time()
    res = f1_group(labels, preds, list_idx)
    sc = np.mean(res)

    logger.info('cal [{}] {} {}'.format(env.iteration + 1, sc, time.time() - t))

受け取る引数はひとつ (lightgbm.callback.CallbackEnv)

ただの名前付きタプルです。lightgbmユーザなら名前で大体中身が分かるはずです。

唯一よく知らないのがmodel(Booster)かと思うので、この後見ていきます。

begin_iterationとend_iterationがあるのはinit_model指定して学習を途中からやる場合があるからです。

CallbackEnv = collections.namedtuple(
    "LightGBMCallbackEnv",
    ["model",
     "params",
     "iteration",
     "begin_iteration",
     "end_iteration",
     "evaluation_result_list"])

LightGBM/callback.py at master · Microsoft/LightGBM · GitHub

データのアクセスはAPI経由

callback関数最大の難関は、lightgbmの学習データをpythonから直接アクセス出来ないことです。

なぜかというと、学習はCでやるためCのオブジェクトでデータを持ってるためです。

メモリ的にPython側でも持つわけには行かないので仕方ありません。

もっと良いやり方があるかもしれませんが、学習データにはcustom metricでアクセスするのが一番楽そうでした。

なのでget_pred_metricみたいな一切計算しないで予測値を返すだけのcustom metricを実装して無理やり取得します。

preds = [ele[2] for ele in clf.eval_train(get_pred_metric) if ele[1] == 'pred'][0]
preds = [ele[2] for ele in clf.eval_valid(get_pred_metric) if ele[1] == 'pred'][0]

Booster.eval_trainとBooster.eval_testで学習データと検証データのそれぞれの予測値がとれます。

データは複数渡せるので、最後に先頭のだけ取ります。

ラベルは普通にget_label()でとれます。

その他の事はBoosterの実装を見てください。

LightGBM/basic.py at master · Microsoft/LightGBM · GitHub

あとは好きに実装

私の場合は、metricの計算が遅いので10回反復毎にしたり、反復毎にデータのウェイトを変えたり、ログに書いたりといった活用をしました。

ただ反復毎にデータのウェイトを変えるのはC側の変更も必要なので注意下さい。

自前early stoppingのやり方

↑のboosting反復の実装みれば分かる通り、lightgbm.callback.EarlyStopExceptionを上げるだけです。

callback関数内で好きに実装してraiseしましょう。

class EarlyStopException(Exception):
    """Exception of early stopping.
    Parameters
    ----------
    best_iteration : int
        The best iteration stopped.
    """
    def __init__(self, best_iteration, best_score):
        super(EarlyStopException, self).__init__()
        self.best_iteration = best_iteration
        self.best_score = best_score

LightGBM/callback.py at master · Microsoft/LightGBM · GitHub

ただ未検証なので、ちゃんと巻き戻るかとか検証したら教えて下さい。。。

一応書いたけど。。。

ここまで学習を制御すること無いので普通使わないですが、忘れそうなので備忘で記事にしました。

道具としてのCython

皆様tkm2261です。今日は道具としてのCythonと題して、

使うことに特化してCythonの解説をしたいと思います。

きっかけはKaggle

最近まで参加していたInstacart Market Basket Analysis | Kaggle

どうしても高速化したい処理があり実装しました。

Numbaも良いのですが、速くなるときとならないときがあり、JITよりもコンパイルが速いのは自明なのでCythonにしました。

Cythonを使うとき

Cythonはコードの可読性を下げてしまうので、下記のような事を満たすケースか考えた方が良いです。

  • 処理が切り出せる
  • numpy配列の処理(C言語で書けって言われたら書けるレベル)
  • numpyの関数では表現できない。(配列をfor文でアクセスしないといけない)
  • 呼び出し回数が数千回以上

今回のKaggleではこれが当てはまったので紹介します。

FaronさんのF1最適化

このコンペではユーザ毎にアイテムを推薦してそのF1値の平均が指標でした。

つまり、予測値に対して閾値を設けて0-1にする必要があります。

Faronさんは「Ye, N., Chai, K., Lee, W., and Chieu, H. Optimizing F-measures: A Tale of Two Approaches. In ICML, 2012.」のアルゴリズムを実装してカーネルとして公開してくれました。

F1-Score Expectation Maximization in O(n²) | Kaggle

これが絶大な威力を発揮して、コンペはこれの公開後から別の様相を呈しました。

DP (Dynamic Programming)を含んだ実装はCythonの出番

このFaronさんの実装を見てもらえばすぐわかりますが、

いつくかDPテーブルを確保して、for文を回しています。これは超Cython向きです。

使い方 その1 『Cython実装』

実装自体はこの記事最後に付けたので参照下さい。ここではPythonとの相違点をメインに話します。

ファイルは.pyx

とくに考えずに、.pyではなく.pyxとしましょう。

cimport

CythonやNumpyのCライブラリに直接アクセスするためにcimportします

cimport cython
cimport numpy as np

ライブラリやヘッダにパスを通す必要がありますが、それは後ほど

型宣言

cdefを付けて使う前に宣言してください。

cdef int i, j, k

一番上でi,j,kをintで宣言するのは競プロっぽいですね。

実際、競プロ+CythonをやればPythonコードをかなり速く出来ると思います。

配列の場合は次元数と型を宣言してください。

cdef np.ndarray[double, ndim = 2] DP_C
cdef np.ndarray[long, ndim= 1] idx

np.float64で渡すと、Cython側ではlong型になるので注意しましょう。メモリにうるさい場合は呼び出しの方でも小さい整数型を使って下さい。

Cythonでは配列の次元、型指定が速度のキーとなるので必ず指定しましょう。

これやらないと使う意味ないです。

オプション指定

Pythonでは配列の長さを超えてアクセスするとエラーになったり、-1でアクセスすると一番後ろが取れたり便利な機能があります。

裏を返せば、そこのチェックに時間を使っているので、これを切ると高速になります。

@cython.boundscheck(False)
@cython.wraparound(False)

ただしC言語同様にSegmentation Faultを吐くようになります。

普通のPythonも書ける

今回はユーザ分並列したかったのでmultiprocessingも呼んでいます。

付録の感じで普通のPythonコードもpyxファイル内で呼ぶことが出来ます。

ただし、高速化は余りされないので注意しましょう。

※ CythonでGIL外してマルチスレッド処理も可能ですが、超面倒なので極力避けましょう

このf1_groupはlightgbmの評価関数として使える形になっており、これを使ってコンペではチューニングしていますた。

使い方 その2『コンパイル

ここではsetup.pyを使うコンパイルを紹介します。

pyximport使うのもありますが、呼び出し元でcythonを意識したくないのでsetup.pyで.soを作ります。

コンパイルするpyxファイルはutils.pyxとします。

setup.pyの書き方

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy

setup(
    cmdclass={'build_ext': build_ext},
    ext_modules=[Extension("utils", ["utils.pyx"])],
    include_dirs=[numpy.get_include()],
    extra_compile_args=['-fopenmp'],
    extra_link_args=['-fopenmp'],
)

numpy.get_include()で動的に持ってくるのがミソです。

コンパイルの実行

下記を実行してすると、

python setup.py build_ext --inplace

というファイルが作成されます。これを実行したいフォルダに持っていきます。

使い方 その3『呼び出し』

作成されたsoファイルを実行したいフォルダに持ってきて普通にimportできます

from utils import f1_opt

soファイルのドット前の名前でimport出来ます。

速度検証

実際に1万人に対して並列せずに実行して速度を見てみます。

計算時間
Python実装 241.1秒
Cython実装 21.5秒

約12倍の高速化になっています。

使いすぎには注意ですが、Cythonは非常に強力なツールです。

付録: 実装 utils.pyx

スネークケースとかごっちゃになってますがコピペして持ってきたのですみません。。。

cimport cython
import numpy as np
cimport numpy as np
from sklearn.metrics import f1_score


@cython.boundscheck(False)
@cython.wraparound(False)
def f1_opt(np.ndarray[long, ndim=1] label, np.ndarray[double, ndim=1] preds):
    """ ユーザ毎の期待F1値最大化&その時の実際のF1値取得関数

    :param label np.ndarray: 0-1の正解ラベル 
    :param preds np.ndarray: 予測確率ラベル 
    """

    cdef int i, j, k, k1
    cdef double f1, score, f1None, pNone
    cdef long n = preds.shape[0]

    pNone = (1 - preds).prod()

    cdef np.ndarray[long, ndim= 1] idx = np.argsort(preds)[::-1]
    label = label[idx]
    preds = preds[idx]

    cdef np.ndarray[double, ndim = 2] DP_C = np.zeros((n + 2, n + 1), dtype=np.float)

    DP_C[0, 0] = 1.0
    for j in range(1, n):
        DP_C[0, j] = (1.0 - preds[j - 1]) * DP_C[0, j - 1]
    for i in range(1, n + 1):
        DP_C[i, i] = DP_C[i - 1, i - 1] * preds[i - 1]
        for j in range(i + 1, n + 1):
            DP_C[i, j] = preds[j - 1] * DP_C[i - 1, j - 1] + (1.0 - preds[j - 1]) * DP_C[i, j - 1]

    cdef np.ndarray[double, ndim = 1] DP_S = np.zeros((2 * n + 1,))
    cdef np.ndarray[double, ndim = 1] DP_SNone = np.zeros((2 * n + 1,))
    for i in range(1, 2 * n + 1):
        DP_S[i] = 1. / (1. * i)
        DP_SNone[i] = 1. / (1. * i + 1)

    score = -1
    cdef np.ndarray[double, ndim= 1] expectations = np.zeros(n + 1)
    cdef np.ndarray[double, ndim= 1] expectationsNone = np.zeros(n + 1)

    for k in range(n + 1)[::-1]:
        f1 = 0
        f1None = 0
        for k1 in range(n + 1):
            f1 += 2 * k1 * DP_C[k1][k] * DP_S[k + k1]
            f1None += 2 * k1 * DP_C[k1][k] * DP_SNone[k + k1]
        for i in range(1, 2 * k - 1):
            DP_S[i] = (1 - preds[k - 1]) * DP_S[i] + preds[k - 1] * DP_S[i + 1]
            DP_SNone[i] = (1 - preds[k - 1]) * DP_SNone[i] + preds[k - 1] * DP_SNone[i + 1]
        expectations[k] = f1
        expectationsNone[k] = f1None + 2 * pNone / (2 + k)

    if expectations.max() > expectationsNone.max():
        i = np.argsort(expectations)[n] - 1
        tp = label[:i + 1].sum()
        if tp > 0:
            precision = tp / (i + 1)
            recall = tp / label.sum()
            f1 = (2 * precision * recall) / (precision + recall)
        else:
            f1 = 0
    else:
        i = np.argsort(expectationsNone)[n] - 1
        tp = label[:i + 1].sum() if label.sum() != 0 else 1
        if tp > 0:
            precision = tp / (i + 2)
            recall = tp / max(label.sum(), 1)
            f1 = (2 * precision * recall) / (precision + recall)
        else:
            f1 = 0

    return f1

from multiprocessing import Pool


@cython.boundscheck(False)
@cython.wraparound(False)
def f1_group(np.ndarray[long, ndim=1] label, np.ndarray[double, ndim=1] preds, np.ndarray[long, ndim=1] group):
    """ ユーザ平均F1値取得関数

    :param np.ndarray label: 0-1の正解ラベル 
    :param np.ndarray preds: 予測確率ラベル
    :param  np.ndarray group: ユーザの切れ目(lightgbmのgroup定義に準拠)
    """
    cdef int i, start, end, j, s
    cdef double score = 0.
    cdef long m = group.shape[0]
    cdef long n = preds.shape[0]
    start = 0

    p = Pool()
    list_p = []
    for i in range(m):
        end = start + group[i]
        list_p.append(p.apply_async(f1_opt, (label[start:end], preds[start:end],)))
        start = end
    scores = [a.get() for a in list_p]
    p.close()
    p.join()
    return np.mean(scores)