ゆるふわクオンツの日常

WIndowsでGPU版pytorchを使えるようにするCUDA関連のメモ

f:id:dw_dw_dt:20210605001920j:plain

GPU(RTX 3060 ti)を積んだwindowsにて

pytorchを使えるようにセットアップしたので

その作業メモとなります。

事前に必要なもの

Visual Studio(正確にはVS build tools)が必要となります。

手順

基本的には、以下で述べる3つのソフトウェアをインストールすればOKです。

ただし、ライブラリのバージョンなどを踏まえて

あえてバージョンを落としてインストールしたりするのが面倒ではあります。

step.1 NVIDIAドライバーの入手

www.nvidia.co.jp

上記のページから、GPUに対応するインストーラーを入手してください。

今回の3060tiの場合は以下のようになりました。

f:id:dw_dw_dt:20210604235331p:plain

このNVIDIAドライバーはあまりバージョン対応など考えず

とりあえずデフォルトの最新版を入手で大丈夫です。

※日ごろからGeforce Experienceアプリで最新ドライバにしている人はこのプロセスは不要かもしれません。

step.2 CUDA toolkitの入手

これはまず、自分が使いたいライブラリの対応状況を確認する必要があります。

今回はpytorchですので、CUDA 11.1 が必要となりそうです。

f:id:dw_dw_dt:20210604235652p:plain

で、先ほどのNVIDIAドライバーがCUDA 11.1に対応していることを

念のために以下リンク先で確認。

NVIDIAドライバーは後方互換性があるみたいなので、

最新のやつなら古いCUDAにも対応しているみたいです)

docs.nvidia.com

確認が取れたら、以下のアーカイブから欲しいバージョンのCUDAを取得。

(私のバージョンは11.1.1でした)

developer.nvidia.com

step.3 cuDNNの入手

以下のリンク先からcuDNNを入手します。

進めていき、step.2で取得したCUDAのバージョンと対応するcuDNNを選べばOKです。

developer.nvidia.com

step.4 インストール

これまでで入手したものを順にインストールします。

まずはNVIDIAドライバー。

次にCUDA toolkit。ちなみにここでVisual Studioが必要になります。

Visual Studioのバージョンが違う的なメッセージが出たりしますが、

それは無視して進めて問題は出なかったです。

で、最後にcuDNNのzipファイルを展開します。

展開後は、zipの中身のcudaフォルダを適当なCドライブ直下などに配置し

f:id:dw_dw_dt:20210605003256p:plain

ユーザー環境変数のPathにC:/cuDNN/cuda/bin通して完了です。

動作確認

以上でGPUが使えるようになったと思いますのでpytorchで動作確認です。

pytorch.org

あとは公式サイトに従って、ライブラリをインストールします。

conda install pytorch torchvision torchaudio cudatoolkit=11.1 -c pytorch -c conda-forge

で、インストールが完了したら、

import torch
print(torch.cuda.is_available())

を実行し、エラーがなければOKです。

今回の関連図書

東芝SBMの疑似量子アニーリング?でポートフォリオ最適化メモ

f:id:dw_dw_dt:20210417013820j:plain

※今回の記事は専門じゃないことを記載しているので、間違いがあるかもですが、その際はコメントいただけると幸いです。

さて、今更感がありますが、東芝のシミュレーテッド分岐マシン(以下SBM)なるソルバーを使ってみました。

www.global.toshiba

このSBMというのは量子コンピュータ(イジング方式)が得意な計算を

普通のコンピュータでできるようにしてみたアルゴリズムのことだそうです。

イジング方式は組み合わせ最適化問題が得意(超早い)ということで

実用性はいったん置いておいて、とりあえずポートフォリオ最適化を

整数計画にアレンジした際のやり方をメモしました。

ノリ的にはalpacaさんがarXivにあげている以下の論文と似ているかと思います。

https://arxiv.org/pdf/2009.08412.pdf

なお、SBMは以下のAWS marketplacesにて販売されています(FPGAではなくGPU版)。

aws.amazon.com

量子アニーリング

余裕で専門外なのでざっくりとした話しかできないのですが

そもそも量子コンピュータというのは2種類に分類できるもので

  1. ゲート方式

  2. イジング方式

という2種類があるみたいです。

で、ゲート方式というのは汎用的なマシンのことで

イジング方式は組み合わせ最適化に特化した量子コンピュータのようです。

IBM-Qやgoogleマシンがゲート方式で、D-WAVEがイジング方式みたいです。

で、そのイジング方式は、量子アニーリングという手法で

量子焼きなまし法 - Wikipedia

組み合わせ最適化を解くのですが、それを普通のPCで解けるようにしたのが

今回触ってみる東芝様のSBMというものになります。

ちなみに量子コンピュータそのものについては、

東大が公開している以下のmaterialが役に立ちそうです。

utokyo-icepp.github.io

問題設定

あまり筋が良いとは思えませんが、

いわゆるポートフォリオ最適化に対してSBMを適用してみたいと思います。

通常のポートフォリオ最適化は

 w \in \mathbb{R}を最適ウェイトとして  max_{w} \mu^{T} w - \lambda w^{T} \Sigma w

みたいなものを解くわけですが、

SBMで解くためには、 以下の形(イジング問題)に変換する必要があります。

 x を成分が 0または 1からなるベクトルとして  min_{x} x^{T}Z x

この ZSBMインスタンスに投げることで、答えのresponseが返ってくるイメージです。

問題の変形

さて、上で見た通り通常のポートフォリオ最適化から

SBMで受け取れるように変形する必要があります。

その際、ウェイト wを求める問題から数量 xを求める問題に変わってしまいますが

それは特に問題ないかなと思います。

2進展開

たとえば、12という数字は、以下のように2進展開できます。

 12 = 0 * 2^{0} + 0 * 2^{1} + 1 * 2^{2} + 1 * 2^{3}

このテクニックを使うことで、

12という数字は (0,0,1,1)という01ベクトルと対応付けることができそうです。

すると、例えば期待リターンが3の証券を12個保有した際の

期待リターンは以下のようにあらわすことができます。

 3 * (1, 2, 4, 8) (0, 0, 1, 1)^{T} = 36

さらに、複数証券のときは、 d = (1,2,4,8)として

期待リターンが3の証券を12個、3.5の証券を10個保有した際は

 (3, 3.5)C(0,0,1,1,0,1,0,1)^{T}=71

ただしCは対角成分に dが並んでいる2*8行列

のようにあらわすことができそうです。

このように考えることで、

 \mu^{T} C x

という形で、証券をxだけ保有した際の期待リターンが求まる。

ただし、これではまだ二次形式になっていないのでSBMには投げれない。

そこで、01ベクトルはその成分を2乗しても変わらないことを利用して

 \mu^{T} C x = x^{T} diag(\mu) \otimes diag(d) x

とすることで、はれて二次形式となる。

そして、リスクについては

 x^{T}C^{T} \Sigma C x

とすれば、 xだけ保有した際のリスクになるので

 x^{T} (C^{T} \Sigma C - diag(\mu) \otimes diag(d)) x

を最小化することで、いわゆるポートフォリオ最適化と類似の問題に帰着する。

コード

最小化するべき二次形式の真ん中の行列をまずは計算する。

サンプルとして以下のようなものを想定する。

import numpy as np


Z = np.zeros((2,2))
Z[0][0] = 10
Z[0][1] = 3
Z[1][0] = 4
Z[1][1] = 7

次にこれを所定の形式で吐き出す。

from scipy.io import mmwrite
from scipy.sparse import csr_matrix


mmwrite('Z.mtx', csr_matrix(Z), field='real', symmetry='symmetric')

最後にこれをSBMインスタンスに投げつけて、

レスポンスを読み取る

from urllib import request


res = request.urlopen(
    request.Request(url='***', data=open('Z.mtx', 'rb'), headers={'Content-Type': 'application/octet-stream'}
)
result = json.load(res.read())

となります。

とはいえ

まぁポートフォリオ最適化だと

整数計画に精緻化するメリットが希薄なのと、

変数の数が2進展開で爆増してしまって

結局スピード面の優位性がそがれるので、

あまり良い事例ではないように感じられます。

東芝社の論文内で述べられている為替の三角裁定や

ダルマキャピタルとの共同研究の株 stat arbのほうが筋が良さそうですね。

ではまた!

今回の関連図書

そこそこわかりやすく記載されていた気がします。

windowsでpymc3を使おうとしたときのエラー対応

f:id:dw_dw_dt:20210321222125j:plain

Pythonでのベイズを勉強しようと思い、

ちょうどpymc3対応だった『Pythonによるベイズ統計学入門』

のコードを動かそうとした際のエラー対応メモとなります。

なお、コードは

github.com

をご参照ください。

エラー1: module 'arviz' has no attribute 'geweke'

f:id:dw_dw_dt:20210321221425p:plain

pymc3のインストール後に import pymc3 でコケるエラー

これは

discourse.pymc.io

に書かれているもので、pymc3 インストール時のarvizのバージョンがダメなようです。

conda install -c conda-forge arviz=0.11.0 すれば問題は解消しました。

エラー2

f:id:dw_dw_dt:20210321221305p:plain

コード実行時に上記のメッセージとともに動かなくなる現象。

これは

qiita.com

に書かれているように、

実行時の処理を main() のような関数で包めばちゃんと動くようになりました。

今回の関連図書

pymc3でのコード例がある本はこの本くらいかなぁと思い買いました。

R & rstanだともっといっぱい本があるので、pymc3に強いこだわりがなければ

そちらでやるのもありだとは思います。

秘書問題で学ぶ100%婚活必勝講座(大嘘)

f:id:dw_dw_dt:20210130153934j:plain

今回は、「モーテル問題」「秘書問題」「浜辺の美女問題」

と言われる問題を確認してみたいと思います。

秘書問題

秘書を1人雇いたいとする。  N人が応募してきている。
応募者には順位が付けられ、複数の応募者が同じ順位になることはない。
無作為な順序で1人ずつ面接を行う。
毎回の面接後、その応募者を採用するか否かを即座に決定する。
その応募者を採用するか否かは、それまでの面接に基づいて決定する。
不採用にした応募者を後から採用することはできない。
このような状況で、最良の応募者を選択することが問題の目的である。

秘書問題 - Wikipedia

モデリング

まず、並んでいる N人を 1,2,...,Nの数字の順列とみなします。

すると、1からNまでの数字を並べた順列は N!パターン存在します。

そのN!個の順列によって構成される集合を\Omega、その要素を \omegaとする。

 \omega=(\omega_{1},\omega_{2},...,\omega_{N})

また、t人目の応募者と会った時点で得られている情報を\mathscr{F}_{t}とする。

例えば、\mathscr{F}_{1}=\{Ω,\emptyset\}\mathscr{F}_{2}=\{Ω,\emptyset, \{ \omega |  \omega_{2} >  \omega_{1} \} , \{ \omega |  \omega_{1} >  \omega_{2} \} \} となります。

さて、候補者を選んだ際の効用の設計を以下で行います。

順位が下から k番目の候補者を選んだ際の効用をf(k)とする。
今回は、もっとも良い候補者を選んだ場合のみ効用が発生する場合、つまりf(N)=1k!=Nではf(k)=0となるようなf(k)を想定する。

ここで、  X_{m}を、 m人目の候補者を選んだ際に、その人が全体で何番目かを表す確率変数とします。

つまり、 X_{n}n番目の成分への射影であり、 X_{n}(\omega)=\omega_{n}となります。

例えば、常に5番目の候補者を採用するという戦略は X_{5}であり、

その際の効用は f(X_{5}(\omega))という確率変数であらわされることになる。

しかし、常に5番目の候補者を選択するという戦略のは直感的にもよくなさそうな気がします。

そこで、 X_{m}の添え字の mが機動的に変化する設定を考えたくなる。

そういった場合に役に立つのが最適停止時刻の考え方です。

確率空間 (\Omega, \mathscr{F}, P)上にフィルトレーション\mathscr{F}_{t}と可積分な確率過程 X_{t}が与えられているとする。
このとき、\mathscr{G}を停止時刻全体(つまり \sigma \in \mathscr{G}として\{ \omega |  \sigma(\omega) \le  N \} \in \mathscr{F}_{N})とすると、
 \max_{\sigma \in \mathscr{G}}E[X_{\sigma}]を満たす停止時刻 \sigmaが最適停止時刻である。

この考え方を先の X_{m}に導入します。

つまり、ある時点までに分かっている情報に基づいて採用かどうかを決めることは

停止時刻として表現できるので、

今回の問題は、停止時刻を T(\omega)として
 \max_{T(\omega) \in \mathscr{G}} E[f(X_{T(\omega)}(\omega))]を満たすT(\omega)を求める

という問題に帰着します。

最適停止問題の解き方

本題に入る前に、最適停止問題 \max_{\sigma \in \mathscr{G}}E[ X_{\sigma} ]の一般的な解き方を確認する。

天下り的ですが、 \{Z_{n}\}^{N}_{1}

  •  Z_{N}=X_{N}

  •  Z_{N}, Z_{N-1},...と求まったとき、 Z_{n-1}=E[Z_{n}|\mathscr{F}_{n-1}] \vee X_{n-1}と定める

この時、以下の性質が成り立つ。

  1.  \{Z_{n}\}は優マルチンゲール

  2.  X_{n} \le Z_{n}

  3.  Z_{n}は 1., 2.を満たす中で最小

つまり、この \{Z_{n}\} \{X_{n}\}より値が大きい優マルチンゲールのうちで

最小のものということなのです(スネル包)。

で、この Z_{n}を用いて作られる以下の停止時刻 \sigma_{0}(\omega)が最適停止時刻になっている

最適停止問題 \max_{\sigma \in \mathscr{G}}E[ X_{\sigma} ]における解は \sigma_{0}(\omega)=min(n| X_{n}(\omega)=Z_{n}(\omega))

具体例

コイントスを10回やります。最後に出た目の期待値を大きくするにはいつサイコロを振るのをやめればよいか?

 X_{n} n回目に出た目とします。

 Z_{n}を求めていくと、

 Z_{10}=X_{10}

 Z_{9}=E[Z_{10}|\mathscr{F}_{9}] \vee X_{9}=max(3.5, X_{9})

 Z_{8}=E[Z_{9}|\mathscr{F}_{8}] \vee X_{8}=max(4.25, X_{8})

 Z_{7}=E[Z_{8}|\mathscr{F}_{7}] \vee X_{7}=max(4.6..., X_{7})

 Z_{6}=E[Z_{7}|\mathscr{F}_{6}] \vee X_{6}=max(4.9..., X_{6})

 Z_{5}=E[Z_{6}|\mathscr{F}_{5}] \vee X_{5}=max(5.1..., X_{5})

...

というふうに求まってっていくことがわかります。

で、肝心の停止時刻はというと X_{n}の値を見ながら判定していくことになります。

あとコイントスが2回残っていれば、最後の1回の(10回目の)期待値が3.5なので、

次の(9回目の)コイントスで4以上、つまり Z_{9}=X_{9}なら \sigma_{0}=9となりそこでストップ、3以下なら Z_{9}>X_{9}となり \sigma_{0}=10となる

同様にして、残りコイントスが3~5回なら、次に出た目が4以下なら続行、5以上ならストップ

残りコイントスが6回以上であれば、次に出る目が6以外は続行

というのが最適戦略となる

また、 Z_{n}が、そのまま続けた際の期待値と X_{n}との比較になっていたことがわかります。

 \sigma_{0}の最適性の証明

ここでは \sigma_{0}が最適となっていること、つまり E[X_{\sigma}] \le E[X_{\sigma_{0}}]を確認します。

まず \{Z_{n\wedge\sigma_{0}}\} \mathscr{F}_{n}マルチンゲールであることを示す。

 Z_{n\wedge\sigma_{0}}

 =Z_{1} + \sum_{k=2}^{n}1_{\{k \le \sigma_{0}\}}(Z_{k}-Z_{k-1})

 =Z_{1} + \sum_{k=2}^{n}1_{\{k \le \sigma_{0}\}}(Z_{k}-E[Z_{k}|  \mathscr{F}_{k-1} ])

よって E[Z_{n\wedge\sigma_{0}} -Z_{(n-1)\wedge\sigma_{0}}|  \mathscr{F}_{n-1} ]=0より結論を得る。

 \mathscr{F}_{n}マルチンゲールであることから

 E[Z_{1}] = E[Z_{1\wedge\sigma_{0}}] = E[Z_{N\wedge\sigma_{0}}] = E[Z_{\sigma_{0}} ]

とわかる。

ここで X_{n} \le Z_{n} Z_{n}の優マルチンゲール性より

 E[X_{\sigma}] \le E[Z_{\sigma} ] \le E[Z_{1}] = E[Z_{\sigma_{0}} ]

この時、 \sigma_{0}において X_{\sigma_{0}} = Z_{\sigma_{0}}なので

 E[X_{\sigma}] \le E[Z_{\sigma} ] \le E[Z_{1}] = E[Z_{\sigma_{0}} ]= E[X_{\sigma_{0}} ]

したがって、 E[X_{\sigma}] \le E[X_{\sigma_{0}}]とわかった

今回の最適停止問題の解き方の概要

まずはじめに、 \{f(X_{n}(\omega))\} \mathscr{F}_{n}適合とは限らないので

 Y_{n}=E[ f(X_{n}(\omega) |   \mathscr{F}_{n}]とおく。

 Y_{n} \mathscr{F}_{n}適合なので、 X_{n}の代わりに Y_{n}を考えればよい。

すると、先ほどの議論から T_{0}=min\{n| Y_{n}=Z_{n}\}となる T_{0}を求めればよいとわかる。

 Y_{n}=E[ f(X_{n})| \mathscr{F}_{n}] = \sum_{k=1}^{N}f(k)E[1_{\{X_{n}=k\}}| \mathscr{F}_{n}]

であるので、 A_{n,i}=\{\omega|\omega_{n}が\omega_{1},...,\omega_{n}のうち、下から数えてi番目\}とおく。

ここで A_{n,i} \in  \mathscr{F}_{n}であり、 A_{n,i} \cap A_{n,j} \neq  \emptysetかつ \Omega = \cup_{i=1}^{n}A_{n,i}が成立する。

補題1
 \alpha_{k;n,i}を、 i \le k \le N-n+i \le Nのとき {}_{k-1}C_{i-1}  {}_{N-k}C_{n-i} /{}_{N}C_{n}とし、
それ以外の時は0とする。この時以下が成立する。
 E[1_{\{X_{n}=k\}} |  \mathscr{F}_{n}] = \sum_{i=1}^{n}\alpha_{k;n,i}1_{A_{n,i}}

補題1より、 a_{n,i}= \sum_{k=1}^{N}f(k)\alpha_{k;n,i}とおくと以下が成立する。

 Y_{n}
 =\sum_{k=1}^{N}f(k)E[1_{\{X_{n}=k\}}| \mathscr{F}_{n}]
 =\sum_{i=1}^{n}\sum_{k=1}^{N}f(k)\alpha_{k;n,i}1_{A_{n,i}}
 =\sum_{i=1}^{n}a_{n,i}1_{A_{n,i}}

これで Y_{n}についてかなりわかりやすい形になったので、次に Z_{n}を考える。

命題1
 \{c_{n}\}_{n=0}^{N} c_{N}=0とし、帰納的に c_{n-1}=\frac{1}{n}\sum_{i=1}^{n}(a_{n,i} \vee c_{n} )で定める。
この時以下が成立する。
 Z_{n}=\sum_{i=1}^{n}(a_{n,i} \vee c_{n})1_{A_{n,i}}

これで、 T_{0}=min\{ n| Y_{n} = Z_{n} \} Y_{n},Z_{n}の形が特定できた。

少し書き下すことで、 T_{0}=min\{ n| a_{n,i} \ge c_{n} \}となることがわかる。

実はここまでの文脈では、 f(k)についての情報は使っていない。

そこで、最良選択の f(k)について、さらに内容を確認していく。

最良選択の場合、 a_{n,i}=\alpha_{N;n,i}となる。

そして \alpha_{N;n,i} i=nのときに  {}_{N-1}C_{n-1}  {}_{N-N}C_{n-n} /{}_{N}C_{n} = \frac{n}{N}となり、

それ以外の時は0になる。

これを用いて c_{n}を計算すると以下のようになる。

補題2
 \frac{1}{n+1}+...+\frac{1}{N-1} \le 1のとき、 c_{n}=\frac{n}{N}(\frac{1}{n}+\frac{1}{n+1}+...+\frac{1}{N-1})が成り立つ。
 \frac{1}{n}+...+\frac{1}{N-1} > 1のとき、 c_{n} > \frac{n}{N}となり c_{1}=...=c_{n}が成り立つ。

 k_{N}  \frac{1}{k_{N}+1}+...+\frac{1}{N-1} \le 1 かつ  \frac{1}{k_{N}}+...+\frac{1}{N-1} > 1とすると

 n \le k_{N}では c_{n} > \frac{k_{N}}{N} > a_{n,i}となり、停止条件にヒットすることはない。

そして n k_{n}より大きいと、 c_{n}が減少し、

たとえば、 c_{k_{N}+1}=\frac{k_{N}+1}{N}(\frac{1}{k_{N}+1}+...+\frac{1}{N-1}) \le \frac{k_{N}+1}{N} となる。

一方で a_{n,i}は、それまでで最良の時だけ \frac{n}{N} (\ge c_{n}) をとるので、

この問題での最適な戦略は、
 k_{N}人と会ったあとは、今までで一番いい候補者に遭遇したらその人でdoneすればいい

ということですね!

例えば10人の候補者がいる場合は、初めの3人(k_{N})を見送った後、

それまでで一番いい人が現れたらその人に決めればよいということですね!

補題1の証明

よくあるインディケーター関数と確率に関する E[1_{A}] = P(A)の条件付きverで、

 E[1_{\{X_{n}=k\}}| \mathscr{F}_{n}] = \sum_{\sigma \in  S_{n}}P(\{X_{n}=k\}|A_{\sigma})1_{A_{\sigma}}

ここで、 S_{n}というのは n次対称群( \{1,2,...,n\}の順列全体からなる群)とする。

すると、 P(\{X_{n}=k\}|A_{\sigma})というのは、

「いままでの n人の順位付けが与えられたうえで、 n番目にあった人はその中で i番目の順位でした。

ではその n番目にあった人は全体の N人のうちk番目である確率はいくつでしょう?」

という問題の答えになる。

すると、 i \le k \le N-(n-i)が満たされなければ P(\{X_{n}=k\}|A_{\sigma}) = 0となることがわかる。

上記の条件が満たされている場合について考える。

全体の並び方はN!個となる。

このうち、n番目に会った人が、初めのn人のうちでi番目となる並び方は {}_{N}C_{n}(n-1)!(N-n)!となる。

次に、 n番目に会った人が、初めのn人のうちでi番目かつ全体でk番目となる並び方は {}_{k-1}C_{i-1}  {}_{N-k}C_{n-i} (n-1)!(N-n)!となる。

したがって、 P(\{X_{n}=k\}|A_{\sigma}) = \alpha_{k;n,i}

命題1の証明

帰納法で示す。

n=Nのときは定義より Z_{N}=Y_{N}=\sum_{i=1}^{N}a_{N,i}1_{A_{N,i}}となり成立。

 nの時成立すると仮定すると

 Z_{n-1}

 = E[ Z_{n}|  \mathscr{F}_{n-1}] \vee Y_{n-1}

 = \sum_{i=1}^{n}(a_{n,i }\vee c_{n}) E[1_{A_{n,i}} |  \mathscr{F}_{n-1}] \vee Y_{n-1}

 = \frac{1}{n}\sum_{i=1}^{n}(a_{n,i }\vee c_{n}) \vee Y_{n-1}

 = c_{n-1} \vee Y_{n-1}

 = c_{n-1} \vee \sum_{i=1}^{n-1}a_{n-1,i}1_{A_{n-1,i}}

 = \sum_{i=1}^{n-1}(c_{n-1} \vee a_{n-1,i}) 1_{A_{n-1,i}}

よって帰納法より示された。

補題2の証明

帰納法で示す。

定義から c_{N}=0,  c_{N-1}=\frac{1}{N}(0+0+...+\frac{N}{N})=\frac{1}{N},  c_{N-2}=\frac{1}{N-1}(\frac{N-2}{N}+\frac{N-1}{N})=\frac{N-2}{N}(\frac{1}{N-2}+\frac{1}{N-1})

よって、 n=N-2の時には結論は成り立つ。

 c_{n+1}で成り立つと仮定すると、

 c_{n}

 = \frac{1}{N+1}\sum_{i=1}^{n+1}(a_{n+1,i} \vee c_{n+1})

 = \frac{1}{N+1}(c_{n+1}*n+(\frac{n+1}{N} \vee c_{n+1}))

 = \frac{1}{n+1}(\frac{n+1}{N}(\frac{1}{n+1}+...+ \frac{1}{N-1})*n + \frac{n+1}{N})

 = \frac{n}{N}(\frac{1}{n}+\frac{1}{n+1}+...+\frac{1}{N-1})

よって成立。

まとめ

この問題って、結論はよく知られている割にその途中式を見かけないので写経してみました。

現実問題では、Nが少なすぎたり0だったりで

最適戦略もくそもなかったりすることがよくありますが、

winner takes allよろしく、十分大きな Nを獲得できている人は

今回のような戦略をとってみるとよいかもしれませんね。

ではまた。

今回の関連図書

今回はこちらの本の内容となります。確率論の本ですが、今回の確率制御っぽい内容にも触れられています。

正則化でノイズの影響を減らすしくみ

f:id:dw_dw_dt:20201119203355j:plain

ノイズの影響が大きいデータを扱う際(以下のmediumなど)は、

正則化などで対応することが多いと思いますが、

具体的にどういうロジックでノイズの影響を緩和しているのかを、

Ridgeを題材に確認したいと思います。

medium.com

Ridge回帰

確率ベクトル (\tilde{X}, \tilde{Y})の標本 (x_{1},y_{1}),...,(x_{n},y_{n})が得られたとします。

いわゆる線形な回帰モデル y_{i} = x_{i} \beta + \epsilon_{i} \epsilonN(0,\sigma^{2}))を想定し、

真のモデルがこの回帰モデルで表現できる場合を考えます。

最小二乗法による\betaの推定では、

x_{i}を縦に並べた Xy_{i}を縦に並べたYを用いて、

損失関数を || Y-X\beta ||_{2} と想定し、

これを最小化する \beta_{OLS}=(X'X)^{-1}X'Yを求めることになります。

一方、

Ridge回帰はというと、

寝室関数を||Y-X\beta||_{2}+\lambda ||\beta||_{2} と見て、

これを最小化する \beta_{Ridge}=(X'X+\lambda I)^{-1}X'Yとなります。

さて、

得られた \beta_{OLS}, \beta_{Ridge}とノイズの関係性を確認するため、

 Y = X\beta + \epsilonを代入します。

すると、 \beta_{OLS}=(X'X)^{-1}X'(X\beta + \epsilon)となるわけですが、

この式から、ノイズ\epsilonがあと1だけ大きかったら

 \beta_{OLS}がどれくらい違っていたか等が分かります。

ただ、これだと少し見づらいので

特異値分解 X = UΛV'=\sum{\gamma_{i}u_{i}v_{i}'}を用いることで

 \beta_{OLS} = \beta + \sum\frac{v_{i} u_{i}'}{\gamma_{i}}\epsilon

 \beta_{Ridge}=(\sum\frac{\gamma_{i}^{2}}{\gamma_{i}^{2}+\lambda}v_{i}v_{i}')\beta+\sum\frac{\gamma_{i}v_{i} u_{i}'}{\gamma_{i}^{2}+\lambda} \epsilon

と変形します。

すると、

 \beta_{OLS}は、期待値的には\betaに一致することがわかります。

ただし、\epsilonの分散が大きかったり、Xにゼロに近い特異値などがあった場合は、

推定結果があまり信頼できないものになっている可能性があります。

一方で \beta_{Ridge}は、期待値が \betaじゃなくなっていますが(不偏推定じゃない)

その代わりに \epsilonの分散や特異値による影響を \lambdaで緩和していることがわかります。

ベイズ的に解釈すると

さて、

Ridgeによってノイズの影響が大きい時でも抗えることがわかりましたが

これは実はベイズ的に考えると自然に辿り着くという話があります。

はい、それでは自然に辿り着く様を確認していきましょう。

 y_{i} = x_{i} \beta + \epsilon_{i}において、 \betaの事前分布をN(0,\alpha I)とします。

( \sigma^{2}の事前分布は逆ガンマとすることが多いですが、今回は無視します)

すると、尤度は p(Y|\beta)=N(X\beta, \sigma^{2}I)となります。

ここで、天下り式ですが \betaは共役事前分布になっているので、

事後分布が正規分布になっていることがわかります。

そこで、事後分布を p(\beta|Y)=N(\bar{\beta}, \Gamma^{-1})とし、この\bar{\beta}, \Gamma^{-1}を求めていきます。

 p(\beta|Y) \sim p(Y|\beta)p(\beta)において

両辺の正規分布の指数の中を比較すると

 p(\beta|Y)の指数関数の中身は

 -\frac{1}{2}(\beta-\bar{\beta})'\Gamma(\beta-\bar{\beta})=-\frac{1}{2}\beta' \Gamma \beta+\beta' \Gamma \bar{\beta} + constとなり

 p(Y|\beta)p(\beta)の指数関数の中身は

 -\frac{1}{2}(\frac{1}{\alpha}\beta'I\beta+\frac{1}{\sigma^{2}}(Y-X\beta)'I(Y-X\beta))=-\frac{1}{2}\beta'(\frac{1}{\alpha}I+\frac{1}{\sigma^{2}}X'X)\beta+\frac{1}{\alpha}\beta'X'Y+const

であることから係数比較して

 \Gamma=\frac{1}{\alpha}I+\frac{1}{\sigma^{2}}X'X

 \bar{\beta}=(\frac{\sigma^{2}}{\alpha}I+X'X)^{-1}X'Y

となります。

ここで、\betaの事後分布の期待値である\bar{\beta}

実は\beta_{Ridge}と同じ形をしていることがわかります。

解釈としては

ベイズは事前分布を得られたデータ(尤度)によって調整していく手法ですが、

いわゆる普通のOLSなどはデータのみによって推定を行います。

その違いは、事前分布の考慮によってデータのバリアンスの影響を緩和するかどうかということでしょうか。

最尤法は \ln{p(Data|\beta)}を最大にする\betaの推定であり、

MAP推定では \ln{p(Data|\beta)}+ \ln{p(\beta)}を最大にする\betaの推定ということです。

(ここでは最尤法との比較のためにMAP推定を持ち出しましたが、

\betaの事後分布による期待値とするのが一般的ですかね)

具体的にpの関数に正規分布を考えると、

それぞれ上で見たOLSやRidgeの損失関数と大体一致することがわかりますね。

ちなみに、分散の事前分布をいじってスパースベイズにするとLassoみたいな

スパースな解が得られるようになるみたいです。

ではまた〜

今回の関連書籍

今回の話が載ってました。