2022-06-04

Denoising Diffusion Probabilistic Models

はじめに

　2020年の以下の論文の式変形を記載する。 arxiv.org

拡散モデル

　この論文で考える「拡散モデル」はマルコフ連鎖モデル（未来の状態は現在の状態だけで決まるモデル）である（下図参照）。

1番左の状態 $x_T$ はガウスノイズ画像、1番右の状態 $x_0$ は物体が明確に写っている画像である。いま、 $x_0$ の状態にわずかなガウスノイズを載せた画像を作り、それを状態 $x_1$ とする。これを $T$ 回繰り返すと $x_T$ になるとする。 $x_0$ から $x_T$ に至る過程をForward過程（拡散過程）と呼ぶ。一方、 $x_T$ から $x_0$ に至る逆向きの過程も考えることができる。この過程は、ガウスノイズ画像が徐々にきれいな画像に変化していく過程である。この過程をReverse過程（生成過程）と呼ぶ。

　拡散過程の各ステップでは画像にガウスノイズを載せるので、 $q(x_t|x_{t-1})$ を次式で与えることができる。

${\displaystyle \begin{align} q(x_t|x_{t-1})=\mathcal{N}(x_t;\sqrt{1-\beta_t}x_{t-1},\beta_tI) \label{eq0} \end{align} }$

ここで、 $x_t$ はベクトル、 $I$ は単位行列、 $\beta_t\in(0,1)$ である。本論文では $\beta_t$ （ハイパーパラメータ）として $10^{-4}\sim10^{-2}$ 程度の値を与えている。 $\beta_t$ が小さな値のとき生成過程 $p_{\theta}(x_{t-1}|x_t)$ もガウス分布関数で近似することができる。

${\displaystyle \begin{equation} p_{\theta}(x_{t-1}|x_{t})=\mathcal{N}\left(x_{t-1};\mu_{\theta}(x_t,t),\sigma_t^2I\right) \end{equation} }$

ここで、 $\mu_{\theta}(x_t,t)$ は $t$ と $x_t$ から計算される平均ベクトル、 $\sigma_{t}^2$ （ハイパーパラメータ）は分散である。拡散モデルでは、 $\mu_{\theta}(x_t,t)$ をニューラルネットワークを用いて学習する。

変分推論

　学習を行うには最小化する関数が必要である。拡散モデルでは次式を最小化する。

${\displaystyle \begin{align} -\ln{p_{\theta}(x_0)}=-\ln{\int dx_{1:T}p_{\theta}(x_{0:T})} \label{eq1} \end{align} }$

ここで、 $dx_{1:T}$ は $dx_1dx_2\cdots dx_T$ を表し、 $p_{\theta}(x_{0:T})$ は $p_{\theta}(x_{0},x_1,\cdots,x_T)$ を表す。マルコフ連鎖モデルなので

${\displaystyle \begin{equation} p_{\theta}(x_{0:T})=p(x_T)\prod_{t=1}^Tp_{\theta}(x_{t-1}|x_t) \end{equation} }$

と書くことができる。これを\eqref{eq1}に代入すると

${\displaystyle \begin{align*} -\ln{p_{\theta}(x_0)}=-\ln{\int dx_{1:T}p(x_T)\prod_{t=1}^Tp_{\theta}(x_{t-1}|x_t)} \end{align*} }$

を得る。ここで、 $q(x_{1:T}|x_0)$ を用いて上式を書き換える。

${\displaystyle \begin{align*} -\ln{p_{\theta}(x_0)}&=-\ln{\int dx_{1:T}p(x_T)\prod_{t=1}^Tp_{\theta}(x_{t-1}|x_t)}\\ &=-\ln{\int dx_{1:T}q(x_{1:T}|x_0)\dfrac{p(x_T)}{q(x_{1:T}|x_0)}\prod_{t=1}^Tp_{\theta}(x_{t-1}|x_t)}\\ &\leq-\int dx_{1:T}q(x_{1:T}|x_0)\ln{\dfrac{p(x_T)}{q(x_{1:T}|x_0)}\prod_{t=1}^Tp_{\theta}(x_{t-1}|x_t)}\equiv L_1 \end{align*} }$

3行目でイェンセンの不等式を用いた。左辺を最小化するには右辺（ $L_1$ ）を最小化すれば良いので、これ以降の議論では $L_1$ の最小化を考える。 $L_1$ を変形すると

${\displaystyle \begin{align} L_1&=-\int dx_{1:T}q(x_{1:T}|x_0)\ln{\dfrac{p(x_T)}{q(x_{1:T}|x_0)}\prod_{t=1}^Tp_{\theta}(x_{t-1}|x_t)}\notag \\ &=-\int dx_{1:T}q(x_{1:T}|x_0)\ln{p(x_T)}-\int dx_{1:T}q(x_{1:T}|x_0)\ln{\dfrac{\prod_{t=1}^Tp_{\theta}(x_{t-1}|x_t)}{q(x_{1:T}|x_0)}} \label{eq2} \end{align} }$

$q(x_{1:T}|x_0)$ は以下のように書くことができる。

${\displaystyle \begin{align*} q(x_{1:T}|x_0)&=q(x_1,\cdots,x_T|x_0)\\ &=q(x_{T}|x_{T-1})\cdots q(x_{1}|x_{0})\\ &=\prod_{t=1}^{T}q(x_{t}|x_{t-1}) \end{align*} }$

これを式\eqref{eq2}に代入すると

${\displaystyle \begin{align*} L_1&=-\int dx_{1:T}q(x_{1:T}|x_0)\ln{p(x_T)}-\int dx_{1:T}q(x_{1:T}|x_0)\ln{\prod_{t=1}^T\dfrac{p_{\theta}(x_{t-1}|x_t)}{q(x_{t}|x_{t-1})}}\\ &=-\int dx_{1:T}q(x_{1:T}|x_0)\ln{p(x_T)}-\int dx_{1:T}q(x_{1:T}|x_0)\sum_{t=1}^T\ln{\dfrac{p_{\theta}(x_{t-1}|x_t)}{q(x_{t}|x_{t-1})}} \end{align*} }$

を得る。上式の両辺に $q(x_0)$ をかけて $x_0$ で積分すると

${\displaystyle \begin{align*} \int dx_0q(x_0)L_1&=-\int dx_{0:T}q(x_{0:T})\ln{p(x_T)}-\int dx_{0:T}q(x_{0:T})\sum_{t=1}^T\ln{\dfrac{p_{\theta}(x_{t-1}|x_t)}{q(x_{t}|x_{t-1})}}\equiv L_2 \end{align*} }$

となる。いま、期待値

${\displaystyle \begin{align*} E_{q(0:T)}[f]=\int dx_{0:T}q(x_{0:T})\;f \end{align*} }$

を定義すると、上式の右辺は

${\displaystyle \begin{align*} L_2=-E_{q(0:T)}\bigl[\ln{p(x_T)}\bigr]-E_{q(0:T)}\Bigl[\sum_{t=1}^T\ln{\dfrac{p_{\theta}(x_{t-1}|x_t)}{q(x_{t}|x_{t-1})}}\Bigr] \end{align*} }$

と書くことができる。 $L_2$ を少し変形して

${\displaystyle \begin{align} L_2&=-E_{q(0:T)}\bigl[\ln{p(x_T)}\bigr]-E_{q(0:T)}\Bigl[\ln{\dfrac{p_{\theta}(x_{0}|x_1)}{q(x_{1}|x_{0})}}\Bigr]-E_{q(0:T)}\Bigl[\sum_{t=2}^T\ln{\dfrac{p_{\theta}(x_{t-1}|x_t)}{q(x_{t}|x_{t-1})}}\Bigr] \label{eq3} \end{align} }$

を得る。ここで

${\displaystyle \begin{align*} q(x_t|x_{t-1})=q(x_t|x_{t-1},x_0)=\dfrac{q(x_{t-1}|x_t,x_0)q(x_t|x_0)}{q(x_{t-1}|x_0)} \end{align*} }$

が成り立つので（後半の等式はベイズの定理）これを式\eqref{eq3}の右辺第3項に代入すると

${\displaystyle \begin{align} E_{q(0:T)}\Bigl[\sum_{t=2}^T\ln{\dfrac{p_{\theta}(x_{t-1}|x_t)}{q(x_{t}|x_{t-1})}}\Bigr] &=E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } \dfrac{ q(x_{t-1}|x_0) }{ q(x_{t}|x_0) } } \Bigr]\notag\\ &= E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } } \Bigr] + E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ q(x_{t-1}|x_0) }{ q(x_{t}|x_0) } } \Bigr] \label{eq4} \end{align} }$

となる。式\eqref{eq4}の第2項は

${\displaystyle \begin{align*} E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ q(x_{t-1}|x_0) }{ q(x_{t}|x_0) } } \Bigr] &= E_{q(0:T)} \Bigl[ \ln{ \dfrac{ q(x_{1}|x_0) }{ q(x_{2}|x_0) }\dfrac{ q(x_{2}|x_0) }{ q(x_{3}|x_0) }\cdots\dfrac{ q(x_{T-1}|x_0) }{ q(x_{T}|x_0) } } \Bigr]\\ &= E_{q(0:T)} \Bigl[ \ln{ \dfrac{ q(x_{1}|x_0) }{ q(x_{T}|x_0) } } \Bigr] \end{align*} }$

と変形される。従って

${\displaystyle \begin{align*} L_2=-E_{q(0:T)}\bigl[\ln{p(x_T)}\bigr]-E_{q(0:T)}\Bigl[\ln{\dfrac{p_{\theta}(x_{0}|x_1)}{q(x_{1}|x_{0})}}\Bigr]-E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } } \Bigr] -E_{q(0:T)} \Bigl[ \ln{ \dfrac{ q(x_{1}|x_0) }{ q(x_{T}|x_0) } } \Bigr] \end{align*} }$

を得る。右辺第2項と第4項をまとめて

${\displaystyle \begin{align*} L_2= -E_{q(0:T)}\bigl[\ln{p(x_T)}\bigr] -E_{q(0:T)}\Bigl[\ln{\dfrac{p_{\theta}(x_{0}|x_1)}{q(x_{T}|x_{0})}}\Bigr] -E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } } \Bigr] \end{align*} }$

右辺第1項と第2項を変形して

${\displaystyle \begin{align} L_2= -E_{q(0:T)}\bigl[\ln{p_{\theta}(x_0|x_1)}\bigr] -E_{q(0:T)}\Bigl[\ln{\dfrac{p_{\theta}(x_{T})}{q(x_{T}|x_{0})}}\Bigr] -E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } } \Bigr] \label{eq5} \end{align} }$

を得る。式\eqref{eq5}の右辺第2項を変形すると

${\displaystyle \begin{align*} E_{q(0:T)}\Bigl[\ln{\dfrac{p_{\theta}(x_{T})}{q(x_{T}|x_{0})}}\Bigr] &=\int dx_0dx_1\cdots dx_Tq(x_0,x_1,\cdots,x_T)\ln{\dfrac{p_{\theta}(x_{T})}{q(x_{T}|x_{0})}}\\ &=\int dx_0dx_T\Bigl[\int dx_1\cdots dx_{T-1}q(x_0,x_1,\cdots,x_T)\Bigr]\ln{\dfrac{p_{\theta}(x_{T})}{q(x_{T}|x_{0})}}\\ &=\int dx_0dx_Tq(x_0,x_T)\ln{\dfrac{p_{\theta}(x_{T})}{q(x_{T}|x_{0})}}\\ &=\int dx_0q(x_0)\Bigl[\int dx_Tq(x_T|x_0)\ln{\dfrac{p_{\theta}(x_{T})}{q(x_{T}|x_{0})}}\Bigr]\\ &=-\int dx_0q(x_0)D_{KL}\bigl(q(x_T|x_0)\;\|\;p(x_T)\bigr) \end{align*} }$

式\eqref{eq5}の右辺第3項を変形すると

${\displaystyle \begin{align*} E_{q(0:T)} \Bigl[ \sum_{t=2}^T \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } } \Bigr] &= \sum_{t=2}^T \int dx_0\cdots dx_T q(x_0,\cdots,x_T) \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } }\\ &= \sum_{t=2}^T \int dx_0dx_{t-1}dx_t q(x_0,x_{t-1},x_t) \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } }\\ &= \sum_{t=2}^T \int dx_0dx_{t}q(x_0,x_t)\int dx_{t-1} q(x_{t-1}|x_0,x_t) \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } }\\ &= -\sum_{t=2}^T \int dx_0dx_{t}q(x_0,x_t) D_{KL}\bigl(q(x_{t-1}|x_t,x_0)\;\|\;p_{\theta}(x_{t-1}|x_t) \bigr) \end{align*} }$

を得る。いま

${\displaystyle \begin{align*} E_{q(0)}[f]&\equiv \int dx_0 q(x_0)\;f\\ E_{q(0,t)}[f]&\equiv \int dx_0dx_t q(x_0,x_t)\;f \end{align*} }$

を定義すると式\eqref{eq5}は次式となる。

${\displaystyle \begin{align} L_2= -E_{q(0:T)}\bigl[\ln{p_{\theta}(x_0|x_1)}\bigr] &+E_{q(0)}\Bigl[D_{KL}\bigl(q(x_T|x_0)\;\|\;p(x_T)\bigr)\Bigr]\notag\\ &+\sum_{t=2}^TE_{q(0,t)}\Bigl[D_{KL}\bigl(q(x_{t-1}|x_t,x_0)\;\|\;p_{\theta}(x_{t-1}|x_t)\bigr)\Bigr] \label{eq7} \end{align} }$

この式が論文中の式(5)である。

$q(x_t|x_0)$ の計算

　式\eqref{eq0}で $\alpha_t=1-\beta_t$ とおく。

${\displaystyle \begin{align*} q(x_t|x_{t-1})=\mathcal{N}\bigl(x_t;\sqrt{\alpha_t}x_{t-1},(1-\alpha_t)I\bigr) \end{align*} }$

いま

${\displaystyle \begin{align*} q(x_1|x_0)&=\mathcal{N}\bigl(x_1;\sqrt{\alpha_1}x_{0},(1-\alpha_1)I\bigr)\\ q(x_2|x_1)&=\mathcal{N}\bigl(x_2;\sqrt{\alpha_2}x_1,(1-\alpha_2)I\bigr) \end{align*} }$

を考え、次式を計算する。

${\displaystyle \begin{align} q(x_2|x_0)&=\int dx_1q(x_2|x_1)q(x_1|x_0)\notag\\ &=\int dx_1\mathcal{N}\bigl(x_2;\sqrt{\alpha_2}x_1,(1-\alpha_2)I\bigr)\mathcal{N}\bigl(x_1;\sqrt{\alpha_1}x_{0},(1-\alpha_1)I\bigr)\notag\\ &\propto\int dx_1 \exp{ \Bigl\{ -\dfrac{1}{2} \Bigl[ \dfrac{\left(x_1-\sqrt{\alpha_1}x_0 \right)^T\left(x_1-\sqrt{\alpha_1}x_0\right)}{1-\alpha_1} \Bigr. \Bigr.\\ \Bigl. \Bigl. +\dfrac{\left(x_2-\sqrt{\alpha_2}x_1\right)^T\left(x_2-\sqrt{\alpha_2}x_1\right)}{1-\alpha_2} \Bigr] \Bigr\} } \label{eq6} \end{align} }$

3行目は $x_1$ に依存する項だけを残した。指数関数の肩の $[\cdots$ ]に注目して

${\displaystyle \begin{align*} [\cdots]=a\;x_1^Tx_1-b^Tx_1-x_1^Tb+\dfrac{\alpha_1}{1-\alpha_1}x_0^Tx_0+\dfrac{1}{1-\alpha_2}x_2^Tx_2 \end{align*} }$

ここで

${\displaystyle \begin{align*} a&\equiv\dfrac{1}{1-\alpha_1}+\dfrac{\alpha_2}{1-\alpha_2}\\ b&\equiv\dfrac{\sqrt{\alpha_1}}{1-\alpha_1}x_0+\dfrac{\sqrt{\alpha_2}}{1-\alpha_2}x_2 \end{align*} }$

と置いた。 $a$ はスカラー、 $b$ はベクトルである。さらに計算すると

${\displaystyle \begin{align*} [\cdots]=a\left(x_1-\dfrac{b}{a}\right)^T\left(x_1-\dfrac{b}{a}\right)-\dfrac{b^Tb}{a} +\dfrac{\alpha_1}{1-\alpha_1}x_0^Tx_0 +\dfrac{1}{1-\alpha_2}x_2^Tx_2 \end{align*} }$

を得る。上式を式\eqref{eq6}に代入し、 $x_1$ について積分すると

${\displaystyle \begin{align*} q(x_2|x_0)&\propto \exp{\Bigl\{-\dfrac{1}{2}\Bigl[ -\dfrac{b^Tb}{a} +\dfrac{\alpha_1}{1-\alpha_1}x_0^Tx_0 +\dfrac{1}{1-\alpha_2}x_2^Tx_2 \Bigr]}\Bigr\}\\ &\propto \exp{\Bigl\{-\dfrac{1}{2}\Bigl[ \dfrac{\left(x_2-\sqrt{\alpha_1\alpha_2}x_0\right)^T\left(x_2-\sqrt{\alpha_1\alpha_2}x_0\right)}{1-\alpha_1\alpha_2} \Bigr]}\Bigr\}\\ \end{align*} }$

となる。2行目への変形では $x_2$ に依存する項だけを残した。以上より

${\displaystyle \begin{align*} q(x_2|x_0)=\mathcal{N}\bigl(x_2;\sqrt{\bar{\alpha}_2}\;x_{0},(1-\bar{\alpha}_2)I\bigr) \end{align*} }$

が成り立つことがわかる。ここで

${\displaystyle \begin{align*} \bar{\alpha}_2=\alpha_1\alpha_2 \end{align*} }$

とした。同じ計算を繰り返していけば

${\displaystyle \begin{align} q(x_t|x_0)&=\mathcal{N}\bigl(x_t;\sqrt{\bar{\alpha}_t}\;x_{0},(1-\bar{\alpha}_t)I\bigr)\label{eq9}\\ \bar{\alpha}_t&=\prod_{s=1}^t\alpha_s\notag \end{align} }$

を得る。これが論文の式(4)である。また、 $t$ を大きくしていくと $\bar{\alpha}_t$ は小さくなるので $q(x_t|x_0)$ は標準正規分布に近づいていくことが分かる。

$q(x_{t-1}|x_t,x_0)$ の計算

　ベイズの定理より

${\displaystyle \begin{align*} q(x_{t-1}|x_t,x_0)=\dfrac{q(x_{t-1}|x_0)q(x_t|x_{t-1},x_0)}{q(x_t|x_0)} \end{align*} }$

ここに\eqref{eq0}と\eqref{eq9}を代入する。

${\displaystyle \begin{align*} q(x_{t-1}|x_t,x_0)&=\dfrac{\mathcal{N}\bigl(x_{t-1};\sqrt{\bar{\alpha}_{t-1}}\;x_{0},(1-\bar{\alpha}_{t-1})I\bigr)\mathcal{N}(x_t;\sqrt{\alpha_t}x_{t-1},\beta_tI) }{\mathcal{N}\bigl(x_t;\sqrt{\bar{\alpha}_t}\;x_{0},(1-\bar{\alpha}_t)I\bigr)}\\ &\propto \exp{\Bigl\{-\dfrac{1}{2}\Bigl[ \dfrac{\left(x_{t-1}-\sqrt{\bar{\alpha}_{t-1}}\;x_0\right)^T \left(x_{t-1}-\sqrt{\bar{\alpha}_{t-1}}\;x_0 \right)}{1-\bar{\alpha}_{t-1}} \\ \;\;\;\;\;\;+\dfrac{\left(x_{t}-\sqrt{\alpha_{t}}\;x_{t-1}\right)^T \left(x_{t}-\sqrt{\alpha_{t}}\;x_{t-1} \right)}{\beta_t} \Bigr]\Bigr\}} \end{align*} }$

2行目では $x_{t-1}$ に依存する項だけを残した。指数関数の肩の $[\cdots$ ]に注目して

${\displaystyle \begin{align*} [\cdots]=a\;x_{t-1}^Tx_{t-1}-b^Tx_{t-1}-x_{t-1}^Tb+\cdots \end{align*} }$

ここで

${\displaystyle \begin{align*} a&\equiv\dfrac{1}{1-\bar{\alpha}_{t-1}}+\dfrac{\alpha_t}{\beta_t}\\ b&\equiv\dfrac{\sqrt{\bar{\alpha}_{t-1}}}{1-\bar{\alpha}_{t-1}}x_0+\dfrac{\sqrt{\alpha_t}}{\beta_t}x_{t} \end{align*} }$

と置いた。 $a$ はスカラー、 $b$ はベクトルである。さらに変形すると

${\displaystyle \begin{align*} [\cdots]=a\left(x_{t-1}-\dfrac{b}{a}\right)^T\left(x_{t-1}-\dfrac{b}{a}\right)+\cdots \end{align*} }$

ここに上の $a$ と $b$ を代入すると

${\displaystyle \begin{align*} [\cdots]&=\dfrac{1-\bar{\alpha}_t}{(1-\bar{\alpha}_{t-1})\beta_t} \left(x_{t-1}-\left(\dfrac{\beta_t\sqrt{\bar{\alpha}_{t-1}}}{1-\bar{\alpha}_t}x_0+\dfrac{\sqrt{\alpha_t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_t\right)\right)^T\\ &\;\;\;\;\;\;\;\;\;\;\;\;\times \left(x_{t-1}-\left(\dfrac{\beta_t\sqrt{\bar{\alpha}_{t-1}}}{1-\bar{\alpha}_t}x_0+\dfrac{\sqrt{\alpha_t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_t\right)\right)\\ &=\dfrac{\left(x_{t-1}-\tilde{\mu}_t(x_t,x_0)\right)^T\left(x_{t-1}-\tilde{\mu}_t(x_t,x_0)\right)}{\tilde{\beta}_t} \end{align*} }$

となる。ここで

${\displaystyle \begin{align} \tilde{\mu}_t(x_t,x_0)&\equiv \dfrac{\beta_t\sqrt{\bar{\alpha}_{t-1}}}{1-\bar{\alpha}_t}x_0+\dfrac{\sqrt{\alpha_t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_t\label{eq11}\\ \tilde{\beta}_t&\equiv \dfrac{(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_t}\beta_t\notag \end{align} }$

と置いた。上式を指数関数の肩に戻すと

${\displaystyle \begin{align*} q(x_{t-1}|x_t,x_0) &\propto \exp{\Bigl\{-\dfrac{1}{2}\Bigl[ \dfrac{\left(x_{t-1}-\tilde{\mu}_t(x_t,x_0)\right)^T\left(x_{t-1}-\tilde{\mu}_t(x_t,x_0)\right)}{\tilde{\beta}_t} \Bigr]\Bigr\}} \end{align*} }$

となる。従って

${\displaystyle \begin{align*} q(x_{t-1}|x_t,x_0) =\mathcal{N}\left(x_{t-1};\tilde{\mu}_t(x_t,x_0),\tilde{\beta}_tI\right) \end{align*} }$

を得る。これらが論文中の式(6)と式(7)に相当する。

最小化すべき関数の計算

　最小化すべき関数\eqref{eq7}は以下であった。

${\displaystyle \begin{align*} L_2= -E_{q(0:T)}\bigl[\ln{p_{\theta}(x_0|x_1)}\bigr] &+E_0\Bigl[D_{KL}\bigl(q(x_T|x_0)\;\|\;p(x_T)\bigr)\Bigr]\notag\\ &+\sum_{t=2}^TE_{0,t}\Bigl[D_{KL}\bigl(q(x_{t-1}|x_t,x_0)\;\|\;p_{\theta}(x_{t-1}|x_t)\bigr)\Bigr] \end{align*} }$

右辺第2項は学習から決めるパラメータ $\theta$ に依存しないので無視できる。

${\displaystyle \begin{align} L_3\equiv -E_{q(0:T)}\bigl[\ln{p_{\theta}(x_0|x_1)}\bigr] +\sum_{t=2}^TE_{0,t}\Bigl[D_{KL}\bigl(q(x_{t-1}|x_t,x_0)\;\|\;p_{\theta}(x_{t-1}|x_t)\bigr)\Bigr] \label{eq16} \end{align} }$

最初に $L_3$ の第2項を考える。期待値の中のKullback Leibler divergenceを取り出して

${\displaystyle \begin{align*} D_{KL}\bigl(q(x_{t-1}|x_t,x_0)\;\|\;p_{\theta}(x_{t-1}|x_t)\bigr)= -\int dx_{t-1} q(x_{t-1}|x_0,x_t) \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } } \end{align*} }$

ここに

${\displaystyle \begin{align*} q(x_{t-1}|x_t,x_0)&=\mathcal{N}(x_{t-1};\tilde{\mu}_t(x_t,x_0),\tilde{\beta}_tI)\\ p_{\theta}(x_{t-1}|x_t)&=\mathcal{N}(x_{t-1};\mu_{\theta}(x_t,x_0),\sigma_t^2I) \end{align*} }$

を代入すると

${\displaystyle \begin{align} &-\int dx_{t-1} q(x_{t-1}|x_0,x_t) \ln{ \dfrac{ p_{\theta}(x_{t-1}|x_t) }{ q(x_{t-1}|x_{t},x_0) } }\notag\\ &=\int dx_{t-1} q(x_{t-1}|x_0,x_t)\Bigl\{\dfrac{1}{2}\Bigl[ \dfrac{\left(x_{t-1}-\mu_{\theta}\right)^T\left(x_{t-1}-\mu_{\theta}\right)}{\sigma_t^2} -\dfrac{\left(x_{t-1}-\tilde{\mu}_{t}\right)^T\left(x_{t-1}-\tilde{\mu}_{t}\right)}{\tilde{\beta}_t} \Bigr] \Bigr\}+C \label{eq10} \end{align} }$

となる。 $C$ は $x_{t-1}$ に依存しない項である。 $\tilde{\beta}_t=\sigma_t^2$ の場合を考えて $\{\cdots\}$ の中を計算すると

${\displaystyle \begin{align*} \{\cdots\}&= \dfrac{1}{2\sigma_t^2}\Bigl[ \left(x_{t-1}-\mu_{\theta}\right)^T\left(x_{t-1}-\mu_{\theta}\right) -\left(x_{t-1}-\tilde{\mu}_{t}\right)^T\left(x_{t-1}-\tilde{\mu}_{t}\right) \Bigr]\\ &= \dfrac{1}{2\sigma_t^2}\Bigl[ -\left(x_{t-1}-\tilde{\mu}_{t}\right)^T\left(\mu_{\theta}-\tilde{\mu}_{t}\right) -\left(\mu_{\theta}-\tilde{\mu}_{t}\right)^T\left(x_{t-1}-\tilde{\mu}_{t}\right) +\left(\tilde{\mu}_t-\mu_{\theta}\right)^T\left(\tilde{\mu}_t-\mu_{\theta}\right) \Bigr] \end{align*} }$

上式を\eqref{eq10}に戻して $x_{t-1}$ について積分すると、上式の第1項と第2項はガウス関数と掛け算されて奇関数となるので積分の寄与はゼロになる。以上より

${\displaystyle \begin{align*} D_{KL}\bigl(q(x_{t-1}|x_t,x_0)\;\|\;p_{\theta}(x_{t-1}|x_t)\bigr)= \int dx_{t-1} q(x_{t-1}|x_0,x_t) \Bigl[\dfrac{1}{2\sigma_t^2}\left(\tilde{\mu}_t-\mu_{\theta}\right)^T\left(\tilde{\mu}_t-\mu_{\theta}\right)\Bigr] \end{align*} }$

を得る。ただし定数項は落とした。 $L_3$ の第2項に戻すと

${\displaystyle \begin{align*} \sum_{t=2}^T&E_{q(0,t)}\Bigl[ \int dx_{t-1} q(x_{t-1}|x_0,x_t) \Bigl[\dfrac{1}{2\sigma_t^2}\left(\tilde{\mu}_t-\mu_{\theta}\right)^T\left(\tilde{\mu}_t-\mu_{\theta}\right)\Bigr]\Bigr]\\ &= \sum_{t=2}^T \int dx_0dx_tq(x_0,x_t)\int dx_{t-1} q(x_{t-1}|x_0,x_t) \Bigl[\dfrac{1}{2\sigma_t^2}\left(\tilde{\mu}_t-\mu_{\theta}\right)^T\left(\tilde{\mu}_t-\mu_{\theta}\right)\Bigr]\\ &= \sum_{t=2}^T \int dx_{0:t}q(x_{0:t}) \Bigl[\dfrac{1}{2\sigma_t^2}\left(\tilde{\mu}_t-\mu_{\theta}\right)^T\left(\tilde{\mu}_t-\mu_{\theta}\right)\Bigr]\\ &= \sum_{t=2}^T E_{q(0:t)}\Bigl[ \dfrac{1}{2\sigma_t^2}\left(\tilde{\mu}_t-\mu_{\theta}\right)^T\left(\tilde{\mu}_t-\mu_{\theta}\right) \Bigr]\equiv L_4 \end{align*} }$

を得る。これが論文中の式(8)である。先に求めた式

${\displaystyle \begin{align*} q(x_t|x_0)=\mathcal{N}\bigl(x_t;\sqrt{\bar{\alpha}_t}\;x_{0},(1-\bar{\alpha}_t)I\bigr) \end{align*} }$

に再パラメータ化トリックを適用して

${\displaystyle \begin{align} x_t&=\sqrt{\bar{\alpha}_t}\;x_0+\sqrt{1-\bar{\alpha}_t}\;\epsilon\label{eq14}\\ \epsilon&\sim\mathcal{N}(\epsilon;0,I)\notag \end{align} }$

これを $x_0$ について解くと

${\displaystyle \begin{align} x_0=\frac{1}{\sqrt{\bar{\alpha}}_t}\left(x_t-\sqrt{1-\bar{\alpha}_t}\;\epsilon\right) \label{eq12} \end{align} }$

となる。ここで、上で求めた式\eqref{eq11}を考える。

${\displaystyle \begin{align*} \tilde{\mu}_t= \dfrac{\beta_t\sqrt{\bar{\alpha}_{t-1}}}{1-\bar{\alpha}_t}x_0+\dfrac{\sqrt{\alpha_t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_t \end{align*} }$

右辺の $x_0$ に式\eqref{eq12}を代入すると

${\displaystyle \begin{align*} \tilde{\mu}_t&= \dfrac{\beta_t\sqrt{\bar{\alpha}_{t-1}}}{1-\bar{\alpha}_t}x_0+\dfrac{\sqrt{\alpha_t}(1-\bar{\alpha}_{t-1})}{1-\bar{\alpha}_{t}}x_t\\ &= \frac{1}{\sqrt{\alpha_t}}\left( x_t(x_0,\epsilon)-\frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\epsilon \right) \end{align*} }$

となる。従って $L_4$ は

${\displaystyle \begin{align} L_4= \sum_{t=2}^T E_{q(0:t)} \left[ \dfrac{1}{2\sigma_t^2} \left\| \frac{1}{\sqrt{\alpha_t}}\left( x_t-\frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\epsilon \right) -\mu_{\theta} \right\|^2 \right]\label{eq13} \end{align} }$

となる。これが論文中の式(10)である。上式が最小になるのは

${\displaystyle \begin{align*} \mu_{\theta} =\frac{1}{\sqrt{\alpha_t}}\left( x_t-\frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\epsilon \right) \end{align*} }$

のときである。 $\mu_{\theta}$ は学習で決めるパラメータであるから上式右辺には学習で決める項がなければならない。そこで、新たな量 $\epsilon_{\theta}$ を導入する。

${\displaystyle \begin{align*} \mu_{\theta} =\frac{1}{\sqrt{\alpha_t}}\left( x_t-\frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\epsilon_{\theta}(x_t,t) \right) \end{align*} }$

$\epsilon_{\theta}$ は $x_t$ と $t$ から計算される関数である。上式を式\eqref{eq13}に代入すると

${\displaystyle \begin{align} L_4 = \sum_{t=2}^T E_{q(0:t)}\left[ \dfrac{\beta_t^2}{2\sigma_t^2\alpha_t(1-\bar{\alpha}_t)} \left\| \epsilon-\epsilon_{\theta}(x_t,t) \right\|^2 \right] \label{eq15} \end{align} }$

を得る。ここで、 $x_t$ は式\eqref{eq14}で与えられていることに注意する（次式）。

${\displaystyle \begin{align*} x_t&=\sqrt{\bar{\alpha}_t}\;x_0+\sqrt{1-\bar{\alpha}_t}\;\epsilon\\ \epsilon&\sim\mathcal{N}(\epsilon;0,I)\notag \end{align*} }$

式\eqref{eq15}は、標準正規分布で生成したノイズ $\epsilon$ を $x_t$ と $t$ を用いて予測する目的関数になっている。式\eqref{eq15}を実際に計算する際は、 $x_0$ と $\epsilon$ でサンプリングして期待値を取る。

${\displaystyle \begin{align*} L_4 = \sum_{t=2}^T E_{x_0,\epsilon}\left[ \dfrac{\beta_t^2}{2\sigma_t^2\alpha_t(1-\bar{\alpha}_t)} \left\| \epsilon-\epsilon_{\theta}(x_t,t) \right\|^2 \right] \end{align*} }$

さらに、ノルムの係数も落とし、ステップ $t$ についてもサンプリングを取る。

${\displaystyle \begin{align*} L_4 &= E_{t,x_0,\epsilon}\left[ \left\| \epsilon-\epsilon_{\theta}(x_t,t) \right\|^2 \right]\\ &= E_{t,x_0,\epsilon}\left[ \left\| \epsilon-\epsilon_{\theta}(\sqrt{\bar{\alpha}_t}\;x_0+\sqrt{1-\bar{\alpha}_t}\;\epsilon,t) \right\|^2 \right] \end{align*} }$

これが論文中の式(14)である。最後に式\eqref{eq16}の第1項を考える。

${\displaystyle \begin{align*} -E_{q(0:T)}\bigl[\ln{p_{\theta}(x_0|x_1)}\bigr] &= -\int dx_{0:T}q(x_{0:T})\ln{p_{\theta}(x_0|x_1)}\\ &= -\int dx_1q(x_1)\int dx_{0}q(x_{0}|x_1)\ln{p_{\theta}(x_0|x_1)} \end{align*} }$

$p_{\theta}(x_0|x_1)$ は、生成過程の一番最後のステップである。 $x_0$ を各画素から構成される $D$ 次元のベクトルと考える。

${\displaystyle \begin{align*} &-\int dx_{0}q(x_{0}|x_1)\ln{p_{\theta}(x_0|x_1)}\\ &= -\int dx_0^1\cdots dx_0^D\;q(x_0^1,\cdots,x_0^D|x_1) \ln{p_{\theta}(x_0^1,\cdots,x_0^D|x_1)}\\ &= -\int dx_0^1\cdots dx_0^D\;q(x_0^1,\cdots,x_0^D|x_1) \sum_{i=1}^D\ln{p_{\theta}(x_0^i|x_1)}\\ &= -\sum_{i=1}^D\int dx_0^i\;q(x_0^i|x_1) \ln{p_{\theta}(x_0^i|x_1)}\\ &= -\sum_{i=1}^D\int dx_0^i\;q(x_0^i|x_1) \ln{\mathcal{N}(x_0^i;\mu_{\theta}^i(x_1,1),\sigma_1^2})\\ &= -\sum_{i=1}^D\int_{\mu_{\theta}^i-\delta}^{\mu_{\theta}^i+\delta} dx_0^i \ln{\mathcal{N}(x_0^i;\mu_{\theta}^i(x_1,1),\sigma_1^2}) \label{eq17} \end{align*} }$

最終行では積分範囲を平均値 $\mu_{\theta}^i$ の近傍に制限した（下図参照）。これが論文中の式(13)である（恐らく）。

まとめ

　拡散モデルの2020年の論文Denoising Diffusion Probabilistic Modelsの式変形を追ってみた。恐らく誤りや勘違いがあると思うので指摘いただければ幸いです。

2021-11-21

QSのAcademic Reputationについて

大学

はじめに

　世界大学ランキングのひとつに「QS World University Rankings」がある。 www.topuniversities.com このランキングは6個の指標（後述）の総合点を用いて作られている。今回は6個の中の1つであるAcademic Reputation（学者からの評価）だけを用いて国内の大学を順位付けする。

Academic Reputation

　QSに使われる6個の指標の内訳は以下の通り。

f:id:seiya_kumada:20211121141515p:plain — 表1. 6個の指標

今回はこれらの指標のうち、「Academic Reputation」に注目する。世界中の大学の教員に自分の専門分野でトップと目される研究をしている大学名を上げてもらい、その結果をスコア化したものが「Academic Reputation」である。100点満点のスコアである。

国内の大学の順位

　Academic Reputationだけで並べ替えたときの国内大学の順位は以下の通り。

f:id:seiya_kumada:20211121210226p:plain — 表2. 国内大学の点数

f:id:seiya_kumada:20211121210333p:plain — 図1. 国内大学の点数

上のグラフを見て分かることは以下の通り。

東大は常に満点である。

京大は常に満点に近いが満点を取れない。

国内のトップ大学は3つに分けることができる。東大・京大グループ（A）、阪大・東工大・東北大グループ（B)、その他の旧帝と早慶グループ（C)。

Aグループは常に100点あるいは100点に近い。

BグループとCグループは2019年以降あまり点数が変わらない。

2012年の頃に比べると阪大の点数が随分と落ちている。

慶応大より早大のほうが常に点数は高い。

大学入試の偏差値は阪大より東工大のほうが高いが、QSのAcademic Reputationでは阪大のほうが常に高い。

2022年のアジア地域の順位

　Academic Reputationだけで並べ替えたときのアジア地域の順位は以下の通り。

f:id:seiya_kumada:20211121150751p:plain — 表3. アジア地域の順位（2022年）

2022年の満点大学

　2022年のAcademic Reputationで満点をとった大学は以下の通り。

f:id:seiya_kumada:20211121151204p:plain — 表4. 満点大学

まとめ

　どの世界大学ランキングを見ても、日本の大学の地位が低下している。しかし、QSのAcademic Reputationを見る限り学術的評価は高い。旧帝や早慶は、研究大学として学術的評価を伸ばすことを目標にしてほしい。Academic Reputation以外の指標は気にする必要はないと思う。

2021-09-18

cpplinqサンプル集

C++

はじめに

　ずいぶん昔にGoogleのbloggerに作ったcpplinqのサンプル集が、bloggerの仕様が変わって表示されなくなりました。こちらに転載します。

サンプル集

#include "cpplinq_samples.hpp"
#include <cpplinq.hpp>
#include <iostream>
#include <array>
#include <boost/range/algorithm.hpp>
#include <random>
#include <boost/format.hpp>

void sample_where()
{
    std::cout << "* cpplinq::where\n";
    auto vs = std::array<int, 5>{};
    // make an array of [1,2,3,4,5]
    std::iota(std::begin(vs), std::end(vs), 1);

    std::cout << " even: ";
    cpplinq::from(vs)
        >> cpplinq::where([](const auto& v){ return v % 2 == 0; })
        >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
    
    std::cout << " odd: ";
    cpplinq::from(vs)
        >> cpplinq::where([](const auto& v){ return v % 2 == 1; })
        >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_reference_wrapper()
{
    std::cout << "* std::reference_wrapper/cpplinq::ref\n";
    auto ls = std::list<int>(10);
    // make a list of [-4,-3,...]
    std::iota(std::begin(ls), std::end(ls), -4);
    
    std::cout << " 0: ";
    boost::for_each(ls, [](const auto& l){ std::cout << l << " "; });
    std::cout << std::endl;
    
    // Two expressions have the same meaning.
    auto vs = std::vector<std::reference_wrapper<int>>{
        std::begin(ls), std::end(ls)};
//    auto vs = cpplinq::from(ls) >> cpplinq::ref()
//        >> cpplinq::to_vector();
    
    std::shuffle(std::begin(vs), std::end(vs), std::mt19937{std::random_device{}()});
    
    std::cout << " 1: ";
    boost::for_each(vs, [](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
    
    boost::for_each(ls, [](auto& l){ l *= 2; });
    
    std::cout << " 2: ";
    cpplinq::from(ls)
        >> cpplinq::for_each([](const auto& l){ std::cout << l << " "; });
    std::cout << std::endl;
    
    std::cout << " 3: ";
    boost::for_each(vs, [](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_select()
{
    std::cout << "* cpplinq::select\n";
    auto ps = std::vector<std::pair<int, std::string>>{
        std::make_pair(1, "sabaton"),
        std::make_pair(2, "metallica"),
        std::make_pair(3, "slayer"),
    };
    
    std::cout << " 0: ";
    auto vs = cpplinq::from(ps)
        >> cpplinq::select([](const auto& p){ return p.second; })
        >> cpplinq::to_vector();
    boost::for_each(vs, [](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}
    
void sample_select_many()
{
    std::cout << "* cpplinq::select_many\n";
    auto vs = std::array<int, 3>{1, 2, 3};
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::select_many([](const auto& v){
            auto dst = std::vector<int>{};
            dst.reserve(3);
            dst.emplace_back(v);
            dst.emplace_back(10 * v);
            dst.emplace_back(100 * v);
            return cpplinq::from(dst);})
        >> cpplinq::for_each([](const auto& u){ std::cout << u << " "; });
    std::cout << std::endl;
}
    
void sample_select_many_flatten()
{
    std::cout << "* cpplinq::select_many(flatten)\n";
    auto vss = std::array<std::array<int, 3>, 3>{
        {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}
    };
    std::cout << " 0: ";
    cpplinq::from(vss)
        >> cpplinq::select_many([](const auto& vs){
            return cpplinq::from(vs);})
        >> cpplinq::for_each([](const auto& u){ std::cout << u << " "; });
    std::cout << std::endl;
}

void sample_take()
{
    std::cout << "* cpplinq::take\n";
    auto vs = std::array<int, 100>{};
    std::iota(std::begin(vs), std::end(vs), 1);
    
    std::cout << " 0: ";
    cpplinq::from(vs) >> cpplinq::take(5)
        >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_orderby_descending()
{
    std::cout << "* cpplinq::orderby_descending\n";
    auto vs = std::array<std::pair<int, std::string>, 3>{{
        {std::make_pair(3, "tokyo" )},
        {std::make_pair(1, "osaka" )},
        {std::make_pair(2, "nagoya")},
    }};
    
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::orderby_descending([](const auto& v){ return v.first; })
        >> cpplinq::for_each([](const auto& u){
            std::cout << boost::format("(%1%, %2%) ") %u.first %u.second; });
    std::cout << std::endl;
}

void sample_orderby_ascending()
{
    std::cout << "* cpplinq::orderby_ascending\n";
    auto vs = std::array<std::pair<int, std::string>, 3>{{
        {std::make_pair(3, "tokyo" )},
        {std::make_pair(1, "osaka" )},
        {std::make_pair(2, "nagoya")},
    }};
    
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::orderby_ascending([](const auto& v){ return v.first; })
        >> cpplinq::for_each([](const auto& u){
            std::cout << boost::format("(%1%, %2%) ") %u.first %u.second; });
    std::cout << std::endl;
}

void sample_orderby_true()
{
    std::cout << "* cpplinq::orderby(true)\n";
    auto vs = std::array<std::pair<int, std::string>, 3>{{
        {std::make_pair(3, "tokyo" )},
        {std::make_pair(1, "osaka" )},
        {std::make_pair(2, "nagoya")},
    }};
    // The result is equal to that of orderby_ascending.
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::orderby([](const auto& v){ return v.first; }, true)
        >> cpplinq::for_each([](const auto& u){
            std::cout << boost::format("(%1%, %2%) ") %u.first %u.second; });
    std::cout << std::endl;
}

void sample_orderby_false()
{
    std::cout << "* cpplinq::orderby(false)\n";
    auto vs = std::array<std::pair<int, std::string>, 3>{{
        {std::make_pair(3, "tokyo" )},
        {std::make_pair(1, "osaka" )},
        {std::make_pair(2, "nagoya")},
    }};
    // The result is equal to that of orderby_descending.
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::orderby([](const auto& v){ return v.first; }, false)
        >> cpplinq::for_each([](const auto& u){
            std::cout << boost::format("(%1%, %2%) ") %u.first %u.second; });
    std::cout << std::endl;
}

void sample_thenby()
{
    std::cout << "* cpplinq::thenby\n";
    auto vs = std::array<std::pair<int, std::string>, 4>{{
        {std::make_pair(3, "tokyo" )},
        {std::make_pair(1, "osaka" )},
        {std::make_pair(2, "nagoya")},
        {std::make_pair(3, "tanaka")},
    }};
    // The result is equal to that of orderby_descending.
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::orderby([](const auto& v){ return v.first; },  false)
        >> cpplinq::for_each([](const auto& u){ std::cout << boost::format("(%1%, %2%) ") %u.first %u.second; });
    std::cout << std::endl;

    std::cout << " 1: ";
    cpplinq::from(vs)
        >> cpplinq::orderby([](const auto& v){ return v.first; },  false)
        >> cpplinq::thenby( [](const auto& v){ return v.second; }, true)
        >> cpplinq::for_each([](const auto& u){ std::cout << boost::format("(%1%, %2%) ") %u.first %u.second; });
    std::cout << std::endl;
}

void sample_take_while()
{
    std::cout << "* cpplinq::take_while\n";
    auto vs = std::array<int, 100>{};
    std::iota(std::begin(vs), std::end(vs), 1);
    
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::take_while([](const auto& v){ return v < 5; })
        >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_skip_while()
{
    std::cout << "* cpplinq::skip_while\n";
    auto vs = std::array<int, 10>{};
    std::iota(std::begin(vs), std::end(vs), 1);
    
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::skip_while([](const auto& v){ return v < 5; })
        >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_join()
{
    std::cout << "* cpplinq::join\n";

    struct vocal
    {
        int         id_;
        std::string name_;
    };
    
    struct band
    {
        int         id_;
        std::string name_;
    };
    
    auto bands = std::array<band, 4>
    {{
        {1, "Nightwish"},
        {2, "Amaranthe"},
        {3, "Epica"},
        {4, "Arch Enemy"},
        
    }};
    
    auto vocals = std::array<vocal, 3>
    {{
        {3, "Simone"},
        {1, "Floor"},
        {4, "Alissa"},
    }};
    
    std::cout << " 0: ";
    cpplinq::from(bands)
        >> cpplinq::join(
                cpplinq::from(vocals),
                [](const auto& band)  {return band.id_;},
                [](const auto& vocal) {return vocal.id_;},
                [](const auto& band, const auto& vocal) { return std::make_pair(band, vocal); })
        >> cpplinq::for_each([](const auto& p){
            std::cout << boost::format("(%1%, %2%) ") %p.first.name_ %p.second.name_; });
    std::cout << std::endl;

    std::cout << " 1: ";
    cpplinq::from(vocals)
        >> cpplinq::join(
                cpplinq::from(bands),
                [](const auto& vocal) {return vocal.id_;},
                [](const auto& band)  {return band.id_;},
                [](const auto& vocal, const auto& band) { return std::make_pair(vocal, band); })
        >> cpplinq::for_each([](const auto& p){
            std::cout << boost::format("(%1%, %2%) ") %p.first.name_ %p.second.name_; });
    std::cout << std::endl;

}

void sample_concat()
{
    std::cout << "* concat\n";
    auto vs = cpplinq::range(0, 5);
    auto us = cpplinq::range(5, 5);
    std::cout << " 0: ";
    vs >> cpplinq::concat(us) >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_reverse()
{
    std::cout << "* reverse\n";
    auto vs = cpplinq::range(0, 5);
    std::cout << " 0: ";
    vs >> cpplinq::reverse() >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;
}

void sample_distinct()
{
    std::cout << "* distinct\n";
    auto vs = std::array<int, 9>{5, 4, 3, 2, 1, 2, 3, 4, 5};
    std::cout << " 0: ";
    cpplinq::from(vs) >> cpplinq::distinct() >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;
}
    
void sample_union_with()
{
    std::cout << "* union_with\n";
    auto vs = std::array<int, 3>{5, 4, 3};
    auto us = std::array<int, 3>{3, 5, 1};
    
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::union_with(cpplinq::from(us))
        >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;

    std::cout << " 1: ";
    cpplinq::from(us)
        >> cpplinq::union_with(cpplinq::from(vs))
        >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;
}

void sample_intersect_with()
{
    std::cout << "* intersect_with\n";
    auto vs = std::array<int, 3>{5, 4, 3};
    auto us = std::array<int, 3>{3, 5, 1};
    
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::intersect_with(cpplinq::from(us))
        >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;

    std::cout << " 1: ";
    cpplinq::from(us)
        >> cpplinq::intersect_with(cpplinq::from(vs))
        >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;
}

void sample_except()
{
    std::cout << "* except\n";
    auto vs = std::array<int, 3>{5, 4, 3};
    auto us = std::array<int, 3>{3, 5, 1};
    
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::except(cpplinq::from(us))
        >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;

    std::cout << " 1: ";
    cpplinq::from(us)
        >> cpplinq::except(cpplinq::from(vs))
        >> cpplinq::for_each([](const auto& p){ std::cout << p << " "; });
    std::cout << std::endl;
}

void sample_to_vector()
{
    std::cout << "* to_vector\n";
    std::vector<int> vs = cpplinq::range(0, 3) >> cpplinq::to_vector();
    std::cout << " 0: ";
    boost::for_each(vs, [](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}
    
void sample_to_list()
{
    std::cout << "* to_list\n";
    std::list<int> vs = cpplinq::range(0, 3) >> cpplinq::to_list();
    std::cout << " 0: ";
    boost::for_each(vs, [](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_to_map()
{
    std::cout << "* to_map\n";
    auto vs = std::vector<std::pair<int, std::string>>{
        std::make_pair(1, "tarja"),
        std::make_pair(2, "annette"),
        std::make_pair(3, "floor"),
    };
    
    std::map<int, std::pair<int, std::string>> ms = cpplinq::from(vs)
        >> cpplinq::to_map([](const auto& v){ return v.first; });
    std::cout << " 0: ";
    boost::for_each(ms, [](const auto& p){
        std::cout << "(" << p.first << " " << p.second.second << ") ";
    });
    std::cout << std::endl;
}

void sample_to_lookup()
{
    std::cout << "* to_lookup\n";
    auto vs = std::vector<std::pair<int, std::string>>{
        std::make_pair(1, "fear factory"),
        std::make_pair(2, "soilwork"),
        std::make_pair(3, "arch enemy"),
        std::make_pair(3, "scar symmetry"),
    };
    
    auto ms = cpplinq::from(vs)
        >> cpplinq::to_lookup([](const auto& v){ return v.first; });
    std::cout << " 0: ";
    ms[3]
        >> cpplinq::select([](const auto& m){ return m.second; })
        >> cpplinq::for_each([](const auto& n){ std::cout << n << " "; });
    std::cout << std::endl;
}

void sample_sequence_equal()
{
    std::cout << "* sequence_equal\n";
    auto vs0 = std::vector<int>{1, 2, 3};
    auto vs1 = std::vector<int>{1, 2, 3};
    auto is_same = cpplinq::from(vs0) >> cpplinq::sequence_equal(cpplinq::from(vs1));
    std::cout << " 0: ";
    std::cout << std::boolalpha << is_same << std::endl;
    
    auto vs2 = std::vector<int>{3, 2, 1};
    is_same = cpplinq::from(vs0) >> cpplinq::sequence_equal(cpplinq::from(vs2));
    std::cout << " 1: ";
    std::cout << is_same << std::endl;
}
    
void sample_first()
{
    std::cout << "* first\n";
    auto vs = std::vector<int>{1, 2, 3};
    auto f = cpplinq::from(vs) >> cpplinq::first();
    std::cout << " 0: ";
    std::cout << f << std::endl;
    auto g = cpplinq::from(vs)
        >> cpplinq::first([](const auto& v){ return v % 2 == 0; });
    std::cout << " 1: ";
    std::cout << g << std::endl;
    try {
        auto h = cpplinq::from(vs)
            >> cpplinq::first([](const auto& v){ return v > 9; });
        std::cout << " 2: ";
        std::cout << h << std::endl;
    } catch (const std::exception& error) {
        std::cout << " 3: ";
        std::cout << error.what() << std::endl;
    }
}

void sample_first_or_default()
{
    std::cout << "* first_or_default\n";
    auto vs = std::vector<int>{1, 2, 3};
    auto f = cpplinq::from(vs) >> cpplinq::first_or_default();
    std::cout << " 0: ";
    std::cout << f << std::endl;
    auto g = cpplinq::from(vs)
        >> cpplinq::first_or_default([](const auto& v){ return v % 2 == 0; });
    std::cout << " 1: ";
    std::cout << g << std::endl;
    auto h = cpplinq::from(vs)
        >> cpplinq::first_or_default([](const auto& v){ return v > 9; });
    std::cout << " 2: ";
    std::cout << h << std::endl;
}

void sample_last_or_default()
{
    std::cout << "* last_or_default\n";
    auto vs = std::vector<int>{1, 2, 3, 4, 5};
    auto f = cpplinq::from(vs) >> cpplinq::last_or_default();
    std::cout << " 0: ";
    std::cout << f << std::endl;
    auto g = cpplinq::from(vs) >> cpplinq::last_or_default([](const auto& v){ return v % 2 == 0; });
    std::cout << " 1: ";
    std::cout << g << std::endl;
    auto h = cpplinq::from(vs) >> cpplinq::last_or_default([](const auto& v){ return v > 9; });
    std::cout << " 2: ";
    std::cout << h << std::endl;
}

void sample_element_at_or_default()
{
    std::cout << "* element_at_or_default\n";
    auto vs = std::vector<int>{1, 2, 3, 4, 5};
    auto f = cpplinq::from(vs) >> cpplinq::element_at_or_default(0);
    std::cout << " 0: ";
    std::cout << f << std::endl;
    auto g = cpplinq::from(vs) >> cpplinq::element_at_or_default(3);
    std::cout << " 1: ";
    std::cout << g << std::endl;
    auto h = cpplinq::from(vs) >> cpplinq::element_at_or_default(100);
    std::cout << " 2: ";
    std::cout << h << std::endl;
}
    
void sample_range()
{
    std::cout << "* range\n";
    auto start = 10;
    auto count = 3;
    std::cout << " 0: ";
    cpplinq::range(start, count)
        >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_repeat()
{
    std::cout << "* repeat\n";
    auto element = "tokyo";
    auto count = 3;
    std::cout << " 0: ";
    cpplinq::repeat(element, count)
        >> cpplinq::for_each([](const auto& v){ std::cout << v << " "; });
    std::cout << std::endl;
}

void sample_empty()
{
    std::cout << "* empty\n";
    auto is = cpplinq::empty<int>() >> cpplinq::to_vector();
    std::cout << " 0: ";
    std::cout << std::boolalpha << is.empty() << std::endl;
    
    auto ss = cpplinq::empty<std::string>() >> cpplinq::to_list();
    std::cout << " 1: ";
    std::cout << std::boolalpha << is.empty() << std::endl;
}

void sample_singleton()
{
    std::cout << "* singleton\n";
    std::cout << " 0: ";
    cpplinq::singleton(1)
        >> cpplinq::for_each([](const auto& v){ std::cout << v << std::endl; });
}

void sample_generate()
{
    std::cout << "* generate\n";
    auto mt = std::mt19937{std::random_device{}()};
    auto rnd = std::uniform_int_distribution<>{0, 10};
    std::cout << " 0: ";
    cpplinq::generate([&rnd, &mt](){
        auto value = rnd(mt);
        return (value != 0) ? cpplinq::to_opt(value) : cpplinq::to_opt<int>(); })
        >> cpplinq::for_each([](const auto& x){ std::cout << x << " "; });
    std::cout << std::endl;
}

void sample_any()
{
    std::cout << "* any\n";
    auto vs = std::array<int, 3>{1, 4, 8};
    std::cout << " 0: ";
    bool flag = cpplinq::from(vs)
        >> cpplinq::any([](const auto& v){ return v % 2 == 1; });
    std::cout << std::boolalpha << flag << std::endl;
}

void sample_all()
{
    std::cout << "* all\n";
    auto vs = std::array<int, 3>{2, 4, 8};
    std::cout << " 0: ";
    bool flag = cpplinq::from(vs)
        >> cpplinq::all([](const auto& v){ return v % 2 == 0; });
    std::cout << std::boolalpha << flag << std::endl;
    std::cout << " 1: ";
    flag = cpplinq::from(vs)
        >> cpplinq::all([](const auto& v){ return v % 2 == 1; });
    std::cout << std::boolalpha << flag << std::endl;
}

void sample_contains()
{
    std::cout << "* contains\n";
    auto vs = std::array<int, 3>{2, 4, 8};
    std::cout << " 0: ";
    bool flag = cpplinq::from(vs) >> cpplinq::contains(2);
    std::cout << std::boolalpha << flag << std::endl;

    struct info
    {
        int id_;
        std::string name_;
    };
    
    auto infos = std::vector<info>{
        info{3, "tokyo"},
        info{2, "osaka"},
        info{5, "nagoay"},
    };
    
    flag = cpplinq::from(infos)
        >> cpplinq::contains(infos[0], [](const auto& x, const auto& y){
            return x.id_ == y.id_;
        });
    std::cout << " 1: ";
    std::cout << flag << std::endl;
}

void sample_count()
{
    std::cout << "* count\n";
    auto vs = std::array<int, 5>{1, 2, 3, 4, 8};
    auto c = cpplinq::from(vs) >> cpplinq::count();
    std::cout << " 0: ";
    std::cout << c << std::endl;
    c = cpplinq::from(vs) >> cpplinq::count([](const auto& v){ return v % 2 == 0; });
    std::cout << " 1: ";
    std::cout << c << std::endl;
}

void sample_sum()
{
    std::cout << "* sum\n";
    auto vs = std::array<int, 5>{1, 2, 3, 4, 8};
    auto s = cpplinq::from(vs) >> cpplinq::sum();
    std::cout << " 0: ";
    std::cout << s << std::endl;
    s = cpplinq::from(vs) >> cpplinq::sum([](const auto& v){ return v * v; });
    std::cout << " 1: ";
    std::cout << s << std::endl;
}
    
void sample_min()
{
    std::cout << "* min\n";
    auto vs = std::array<int, 5>{1, 2, 3, 4, 8};
    auto m = cpplinq::from(vs) >> cpplinq::min();
    std::cout << " 0: ";
    std::cout << m << std::endl;
    m = cpplinq::from(vs) >> cpplinq::min([](const auto& v){ return 2 * v; });
    std::cout << " 1: ";
    std::cout << m << std::endl;
}

void sample_max()
{
    std::cout << "* max\n";
    auto vs = std::array<int, 5>{1, 2, 3, 4, 8};
    auto m = cpplinq::from(vs) >> cpplinq::max();
    std::cout << " 0: ";
    std::cout << m << std::endl;
    m = cpplinq::from(vs) >> cpplinq::max([](const auto& v){ return 2 * v; });
    std::cout << " 1: ";
    std::cout << m << std::endl;
}

void sample_avg()
{
    std::cout << "* avg\n";
    auto vs = std::array<double, 5>{1, 2, 3, 4, 8};
    auto m = cpplinq::from(vs) >> cpplinq::avg();
    std::cout << " 0: ";
    std::cout << m << std::endl;
    m = cpplinq::from(vs) >> cpplinq::avg([](const auto& v){ return 2 * v; });
    std::cout << " 1: ";
    std::cout << m << std::endl;
}

void sample_pairwise()
{
    std::cout << "* pairwise\n";
    auto vs = std::array<double, 5>{1, 2, 3, 4, 8};
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::pairwise()
        >> cpplinq::for_each(
            [](const auto& p){
                std::cout << boost::format("(%1%, %2%) ") %p.first %p.second;
            });
    std::cout << std::endl;
}

void sample_zip_with()
{
    std::cout << "* zip_with\n";
    auto vs = std::array<std::string, 3>{"nightwish", "amaranthe", "kamelot"};
    auto us = std::array<std::string, 5>{"finland", "sweden", "america", "japan", "canada"};
    std::cout << " 0: ";
    cpplinq::from(vs)
        >> cpplinq::zip_with(cpplinq::from(us))
        >> cpplinq::for_each(
            [](const auto& p){
                std::cout << boost::format("(%1%, %2%) ") %p.first %p.second;
            });
    std::cout << std::endl;
}

2021-09-18

boost/range/adaptorサンプル集

C++

はじめに

　ずいぶん昔にGoogleのbloggerに作ったboost/range/adaptorのサンプル集が、bloggerの仕様が変わって表示されなくなりました。こちらに転載します。

サンプル集

#include <boost/range/irange.hpp>
#include <boost/range/adaptor/adjacent_filtered.hpp>
#include <boost/range/adaptor/copied.hpp>
#include <boost/range/adaptor/filtered.hpp>
#include <boost/range/adaptor/indexed.hpp>
#include <boost/range/adaptor/indirected.hpp>
#include <boost/range/adaptor/map.hpp>
#include <boost/range/adaptor/replaced.hpp>
#include <boost/range/adaptor/replaced_if.hpp>
#include <boost/range/adaptor/reversed.hpp>
#include <boost/range/adaptor/sliced.hpp>
#include <boost/range/adaptor/strided.hpp>
#include <boost/range/adaptor/tokenized.hpp>
#include <boost/range/adaptor/type_erased.hpp>
#include <boost/range/adaptor/transformed.hpp>
#include <boost/range/adaptor/uniqued.hpp>
#include <boost/range/algorithm/copy.hpp>
#include <boost/range/combine.hpp>
#include <boost/range/algorithm/for_each.hpp>
#include <iostream>
#include <boost/regex.h>

void test_adjacent_filtered()
{
    const auto range = {1, 1, 1, 2, 2, 3, 4, 5, 5};
    const auto result = range | boost::adaptors::adjacent_filtered([](auto a, auto b){ return a != b; });
    boost::copy(result, std::ostream_iterator<int>(std::cout, ", "));
    std::cout << std::endl;
}

void test_copied()
{
    const auto range = std::vector<int>{10, 20, 30, 40, 50};
    const auto result = range | boost::adaptors::copied(1, 3);
    boost::copy(result, std::ostream_iterator<int>(std::cout, ", "));
    std::cout << std::endl;
}

void test_filtered()
{
    const auto range = boost::irange(0, 10);
    const auto result = range | boost::adaptors::filtered([](auto v){ return v % 2 == 0; });
    boost::copy(result, std::ostream_iterator<int>(std::cout, ", "));
    std::cout << std::endl;
}

void test_indexed()
{
    const auto range = boost::irange(100, 110);
    for ( const auto& v : range | boost::adaptors::indexed(3) ) {
        std::cout << v.index() << " --- " << v.value() << std::endl;
    }
}

void test_indirected()
{
    const auto range = std::vector<std::shared_ptr<int>> {
        std::make_shared<int>(10),
        std::make_shared<int>(20),
        std::make_shared<int>(30),
    };
    boost::copy(range | boost::adaptors::indirected, std::ostream_iterator<int>(std::cout, ", "));
    std::cout << std::endl;
    
}

void test_map()
{
    const auto range = std::map<int, std::string> {
        {1, "tokyo"},
        {2, "osaka"},
        {3, "nagoya"},
    };
    
    boost::copy(range | boost::adaptors::map_keys, std::ostream_iterator<int>(std::cout, ", "));
    std::cout << std::endl;
    
    boost::copy(range | boost::adaptors::map_values, std::ostream_iterator<std::string>(std::cout, ", "));
    std::cout << std::endl;
}

void test_replaced()
{
    const auto range = std::string{"tokyo"};
    boost::copy(range | boost::adaptors::replaced('o', 'O'), std::ostream_iterator<char>(std::cout, ""));
    std::cout << std::endl;
}

void test_replaced_if()
{
    const auto range = std::string{"tokyo"};
    boost::copy(range | boost::adaptors::replaced_if([](auto c){ return c == 'o'; }, 'O'), std::ostream_iterator<char>(std::cout, ""));
    std::cout << std::endl;
}

void test_reversed()
{
    const auto range = std::vector<int>{1, 2, 3, 4, 5};
    boost::copy(range | boost::adaptors::reversed, std::ostream_iterator<int>(std::cout, " "));
    std::cout << std::endl;
}

void test_sliced()
{
    const auto range = std::string{"tokyo"};
    boost::copy(range | boost::adaptors::sliced(1, 3), std::ostream_iterator<char>(std::cout, ""));
    std::cout << std::endl;
}

void test_strided()
{
    const auto range = std::string{"!o!y!"};
    boost::copy(range | boost::adaptors::strided(2), std::ostream_iterator<char>(std::cout, ""));
    std::cout << std::endl;
}

void test_tokenized()
{
    auto range = std::string{"[100,200];[300,400];[500,600]"};
    for ( const auto& m : range | boost::adaptors::tokenized(boost::regex(R"(\[\d+,\d+\])")) ) {
        std::cout << m << std::endl;
    }
}

void test_transformed()
{
    auto range = std::vector<int> {1, 2, 3, 4, 5};
    boost::copy(range | boost::adaptors::transformed([](auto i){ return 2 * i; }), std::ostream_iterator<int>(std::cout, " "));
    std::cout << std::endl;
}

void test_uniqued()
{
    auto range = std::vector<int> {1, 1, 3, 4, 4};
    boost::copy(range | boost::adaptors::uniqued, std::ostream_iterator<int>(std::cout, " "));
    std::cout << std::endl;
}

void test_combine()
{
    std::cout << "> combine\n";
    auto v1 = std::vector<int>{1, 2, 3};
    auto v2 = std::vector<std::string>{"tokyo", "osaka", "nagoya"};
    boost::for_each(boost::combine(v1, v2), [](const auto& p){
        std::cout << boost::get<0>(p) << " " << boost::get<1>(p) << std::endl;}
    );
}

int main(int, char* argv[])
{
    test_adjacent_filtered();
    test_copied();
    test_filtered();
    test_indexed();
    test_indirected();
    test_map();
    test_replaced();
    test_replaced_if();
    test_reversed();
    test_sliced();
    test_strided();
    test_tokenized();
    test_transformed();
    test_uniqued();
    test_combine();
    return 0;
}

出力は以下の通り。

> adjacent_filtered
1, 2, 3, 4, 5, 
> copied
20, 30, 
> filtered
0, 2, 4, 6, 8, 
> indexed
3 --- 100
4 --- 101
5 --- 102
6 --- 103
7 --- 104
8 --- 105
9 --- 106
10 --- 107
11 --- 108
12 --- 109
> indirected
10, 20, 30, 
> map
1, 2, 3, 
tokyo, osaka, nagoya, 
> replaced
tOkyO
> replaced_if
tOkyO
> reversed
5 4 3 2 1 
> sliced
ok
> strided
!!!
> tokenized
[100,200]
[300,400]
[500,600]
> transformed
2 4 6 8 10 
> uniqued
1 3 4 
> combine
1 tokyo
2 osaka
3 nagoya

2021-08-28

2つの自然数が互いに素である確率

数学

はじめに

　Youtubeで面白い数学の問題が紹介されていたのでまとめておく。

問題

　問題文は以下の通り。

任意の2つの自然数が互いに素である確率を求めよ。

解答

　ある自然数を考えたとき、それが $n$ の倍数である確率は $1/n$ である。いま、2つの自然数 $A, B$ を取り出したとき、それらが同時に $n$ の倍数にならない確率 $P_n$ は

${\displaystyle \begin{equation} P_n=1-\dfrac{1}{n^2} \end{equation} }$

となる。これを用いると、 $A, B$ が互いに素になる確率 $P$ は

${\displaystyle \begin{equation} P=P_2 P_3 P_5 P_7 \cdots=\prod_{p:{\rm prime}}P_p\tag{1} \end{equation} }$

と書くことができる。ここで、 $p$ は素数である。式(1)の計算をさらに進める。

${\displaystyle \begin{eqnarray} P &=& \prod_{p:{\rm prime}}P_p\\ &=& \dfrac{1}{\prod_{p:{\rm prime}}\left(\dfrac{1}{P_p}\right)}\tag{2}\\ \end{eqnarray} }$

を得る。この分母を考える。

${\displaystyle \begin{eqnarray} \prod_{p:{\rm prime}}\left(\dfrac{1}{P_p}\right) &=& \prod_{p:{\rm prime}}\left(\dfrac{1}{1-\dfrac{1}{p^2}}\right)\\ &=& \left(\dfrac{1}{1-\dfrac{1}{2^2}}\right)\left(\dfrac{1}{1-\dfrac{1}{3^2}}\right)\left(\dfrac{1}{1-\dfrac{1}{5^2}}\right)\cdots\tag{3}\\ &=& 1+\dfrac{1}{2^2}+\dfrac{1}{3^2}+\dfrac{1}{4^2}+\dfrac{1}{5^2}+\cdots\tag{4}\\ &=& \dfrac{\pi^2}{6}\tag{5} \end{eqnarray} }$

式(3)から式(4)への変形にゼータ関数のオイラー積表示を用いた。また、式(4)から式(5)への変形はバーゼル問題と呼ばれる計算である。この結果を式(2)へ代入して

${\displaystyle \begin{equation} P=\dfrac{6}{\pi^2} \end{equation} }$

を得る。

まとめ

　問題文はシンプルであるが、その解答には

バーゼル問題
ゼータ関数のオイラー積

が現れるとても興味深い問題である。

参考動画

２つの自然数が互いに素である確率

2021-08-15

量子線形回帰

量子コンピュータ

はじめに

　線形回帰を量子コンピュータで解くアルゴリズムを説明する。

線形回帰

　 ${\bf x} \in \mathbb{R}^{M}$ , ${\bf x}=\{x_j\}_{j=1}^{M}$ と $y\in\mathbb{R}$ の組 $({\bf x},y)$ が $N$ 個与えられたとき

${\displaystyle \begin{equation} y={\bf w}^T{\bf x} \end{equation} }$

を満たすベクトル ${\bf w}$ を求めることが目的である（ $x_1=1$ としてあり、バイアス項は考慮されているとする）。最小化すべき損失関数は

${\displaystyle \begin{eqnarray} L({\bf w}) &=& \dfrac{1}{2}\sum_{i=1}^N\left( {\bf w}^T{\bf x}_i-y_i \right)^2\\ &=& \dfrac{1}{2}\sum_{i=1}^N\left( \sum_{j=1}^{M}w_jx_{ij} -y_i \right)^2 \end{eqnarray} }$

である。 $w_k$ で偏微分した値を0と置くと

${\displaystyle \begin{eqnarray} \dfrac{\partial L({\bf w})}{\partial w_k} &=& \sum_i\left( \sum_jw_jx_{ij}-y_i \right)\sum_l\delta_{lk}x_{il}\\ &=& \sum_i\left( \sum_jw_jx_{ij}-y_i \right)x_{ik}\\ &=&0 \end{eqnarray} }$

を得る。いま行列 $X=\{x_{ij}\}$ を定義すると上式より

${\displaystyle \begin{eqnarray} && \sum_i\sum_jw_jX_{ij}X_{ik}=\sum_iy_iX_{ik}\\ &\iff& \sum_i\sum_jX^T_{ki}X_{ij}w_j=\sum_iX^T_{ki}y_i\\ &\iff& \sum_j\left(X^TX\right)_{kj}w_j=\sum_iX^T_{ki}y_i\\ &\iff& \left(X^TX{\bf w}\right)_k=\left(X^T{\bf y}\right)_k\\ \end{eqnarray} }$

を得る。任意の $k$ について上式が成り立つから

${\displaystyle \begin{equation} X^TX{\bf w}=X^T{\bf y} \end{equation} }$

となり、これより

${\displaystyle \begin{equation} {\bf w}=\left(X^TX\right)^{-1}X^T{\bf y}\tag{1} \end{equation} }$

を得る。 ${\bf w}\in\mathbb{R}^M,{\bf y}\in\mathbb{R}^N$ である。ところで、任意の行列は特異値分解できるから（特異値分解の詳細はこちらの説明を見てほしい）

${\displaystyle \begin{equation} X=U\Sigma V^T\tag{2} \end{equation} }$

と書くことができる。ここで、 $X\in\mathbb{R}^{N\times M},U\in\mathbb{R}^{N\times N},\Sigma\in\mathbb{R}^{N\times M},V\in\mathbb{R}^{M\times M}$ である。 $U$ と $V$ は直交行列であり

${\displaystyle \begin{eqnarray} U^TU=UU^T&=&I\\ V^TV=VV^T&=&I \end{eqnarray} }$

を満たす。また、 $r\leq M$ として

${\displaystyle \begin{equation} \Sigma={\rm diag}\left(\sigma_1,\cdots,\sigma_r,0,\cdots,0\right) \end{equation} }$

と書くことができる（ $\sigma_i> 0,i=1,\cdots,r$ ）。式(2)を使うと

${\displaystyle \begin{equation} \left(X^TX\right)^{-1}X^T=V{\rm diag}\left( \dfrac{1}{\sigma_1},\cdots,\dfrac{1}{\sigma_r},0,\cdots,0 \right)U^{T}\tag{3} \end{equation} }$

を示すことができる。 $V\in\mathbb{R}^{M\times M},U\in\mathbb{R}^{N\times N}$ であり、(3)式右辺の $V,U$ に挟まれた行列は $M$ 行 $N$ 列なので、上の行列は全体として $M$ 行 $N$ 列の行列になる。式(3)を式(1)に代入すれば ${\bf w}$ の解析解を得ることができる。

　未知のデータ $\tilde{{\bf x}}$ が与えられたときの出力値 $\tilde{y}$ は以下のように計算される。

${\displaystyle \begin{eqnarray} \tilde{y}&=&{\bf w}^T\tilde{{\bf x}}\\ &=& \sum_{j=1}^{M}\biggl( V{\rm diag}\left( \dfrac{1}{\sigma_1},\cdots,\dfrac{1}{\sigma_r},0,\cdots,0 \right)U^{T} {\bf y} \biggr)_j\; \tilde{x}_j\\ &=& \sum_{k=1}^{r}\sum_{i=1}^{N}\sum_{j=1}^{M}(V)_{jk}\dfrac{1}{\sigma_k}(U^T)_{ki}y_i\tilde{x}_j\\ &=& \sum_{k=1}^{r}\dfrac{1}{\sigma_k}\left(\sum_{j=1}^{M}v_{jk}\tilde{x}_j\right) \left( \sum_{i=1}^{N}u_{ik}y_i \right)\tag{4} \end{eqnarray} }$

ここで、 $U,V$ の成分を $u_{ij},v_{ij}$ とした。また、式(2)より

${\displaystyle \begin{eqnarray} &&X^TX=V\;{\rm diag}\left(\sigma_1^2,\cdots,\sigma_r^2,0,\cdots,0\right)V^T\\ &\iff&X^TXV=V\;{\rm diag}\left(\sigma_1^2,\cdots,\sigma_r^2,0,\cdots,0\right) \end{eqnarray} }$

成分をとると

${\displaystyle \begin{eqnarray} &&\left(X^TXV\right)_{ij}=\left(V\;{\rm diag}\left(\sigma_1^2,\cdots,\sigma_r^2,0,\cdots,0\right)\right)_{ij}\\ &\iff& \sum_k\left(X^TX\right)_{ik}V_{kj}=\sum_{m=1}^rV_{im}\;\sigma_m^2\;\delta_{mj}\\ &\iff& \sum_{k}\left(X^TX\right)_{ik}V_{kj}=\sigma_j^2V_{ij}\\ &\iff& \big(\left(X^TX\right)V_{:j}\bigr)_i=\sigma_j^2\left(V_{:j}\right)_i \end{eqnarray} }$

ここで、行列 $V$ の第j列の列ベクトル $V_{:j}$ を ${\bf v}_j$ と書くと

${\displaystyle \begin{equation} X^TX{\bf v}_j=\sigma_j^2{\bf v}_j \tag{5} \end{equation} }$

を得る。すなわち、 $X^TX$ の固有値は $\sigma^2_j$ であり、それに対応する固有ベクトルは ${\bf v}_j$ である。

線形回帰の量子アルゴリズム

　ここでは、前節の定式化を量子コンピュータ上で行う手順を示す。まず最初に、入力データ ${\bf x}_i(i=1,\cdots,N)$ を振幅エンコーディングした状態を用意する。

${\displaystyle \begin{eqnarray} |\psi_X\rangle &=& \sum_{i=1}^N|{\bf x}_i\rangle_A\otimes|i\rangle_B\\ &=& \sum_{i=1}^N\sum_{j=1}^Mx_{ij}\;|j\rangle_A\otimes|i\rangle_B\\ \end{eqnarray} }$

ここで、添字の $A,B$ は量子ビットを収めるレジスタであり、独立した2つの状態空間であることを示す。式(2)より

${\displaystyle \begin{eqnarray} |\psi_X\rangle &=& \sum_{i=1}^N\sum_{j=1}^M\sum_{k=1}^r u_{ik}\sigma_kv_{jk} \;|j\rangle_A\otimes|i\rangle_B\\ &=& \sum_{k=1}^r\sigma_k \sum_{j=1}^Mv_{jk}|j\rangle_A\otimes \sum_{i=1}^M u_{ik}|i\rangle_B\\ &=& \sum_{k=1}^r\sigma_k |{\bf v}_k\rangle_A\otimes|{\bf u}_k\rangle_B\tag{6} \end{eqnarray} }$

を得る。ここで、 ${\bf v}_k,{\bf u}_k$ はそれぞれ行列 $V,U$ の第 $k$ 列の列ベクトルである。 $|\psi_X\rangle$ の密度行列は

${\displaystyle \begin{eqnarray} \rho_{AB}&=&|\psi_X\rangle\langle\psi_X|\\ &=&\sum_{i=1}^N\sum_{j=1}^N|{\bf x}_i\rangle_A\langle{\bf x}_j|_A\otimes|i\rangle_B\langle j|_B \end{eqnarray} }$

となり、Bレジスタの空間でトレースを取ると以下のようになる。

${\displaystyle \begin{eqnarray} {\rm Tr}_B\left(\rho_{AB}\right) &=& \sum_{i=1}^N\sum_{j=1}^N|{\bf x}_i\rangle_A\langle{\bf x}_j|_A\otimes \sum_{m=1}^N\langle m|_B |i\rangle_B\langle j|_B|m\rangle_B\\ &=& \sum_{i=1}^N\sum_{j=1}^N|{\bf x}_i\rangle_A\langle{\bf x}_j|_A\otimes \sum_{m=1}^N\delta_{mi}\delta_{mj}\\ &=& \sum_{i=1}^N\sum_{j=1}^N|{\bf x}_i\rangle_A\langle{\bf x}_j|_A\otimes \delta_{ij}\\ &=& \sum_{i=1}^N|{\bf x}_i\rangle_A\langle{\bf x}_i|_A\\ &=& \sum_{i=1}^N|{\bf x}_i\rangle\langle{\bf x}_i|\\ &=& \sum_{i=1}^N\sum_{j=1}^M\sum_{k=1}^Mx_{ij}x_{ik}|j\rangle\langle k|\\ &=& \sum_{i=1}^N\sum_{j=1}^M\sum_{k=1}^M(X^T)_{ji}(X)_{ik}|j\rangle\langle k|\\ &=& \sum_{j=1}^M\sum_{k=1}^M(X^TX)_{jk}|j\rangle\langle k|\\ &\equiv&\rho_A \end{eqnarray} }$

すなわち、行列 $X^TX$ を埋め込んだ密度行列 $\rho_A$ を得ることできる。式(5)より次式が成り立つことも分かる。

${\displaystyle \begin{equation} \rho_A|{\bf v}_k\rangle=\sigma_k^2|{\bf v}_k\rangle \end{equation} }$

　さて、式(6)に、複数ビットからなる3つ目のレジスタ $C$ を加える。

${\displaystyle \begin{equation} |\psi_X\rangle =\sum_{k=1}^r\sigma_k |{\bf v}_k\rangle_A\otimes|{\bf u}_k\rangle_B\otimes|0\rangle_C\tag{7} \end{equation} }$

上の状態の $A$ レジスタにユニタリー行列 $e^{i2\pi\rho_A}$ を施し、 $\rho_A$ の固有値 $\sigma_k^2$ を $C$ レジスタに入れる（量子位相推定）。

${\displaystyle \begin{equation} \sum_{k=1}^r\sigma_k |{\bf v}_k\rangle_A\otimes|{\bf u}_k\rangle_B\otimes|\sigma_k^2\rangle_C \end{equation} }$

4つ目のレジスタ $D$ を加え $C$ レジスタの値を入力とし算術演算を行う（算術演算や次の回転ゲートの詳細な手順はここを参照のこと）。

さらに、1ビットからなる $E$ レジスタを加え、 $D$ レジスタの各ビットを制御ビットとし $E$ レジスタに回転ゲートを施す。

${\displaystyle \begin{equation} \sum_{k=1}^r\sigma_k |{\bf v}_k\rangle_A\otimes|{\bf u}_k\rangle_B\otimes|\sigma_k^2\rangle_C\otimes|\dfrac{1}{\pi}\cos^{-1}{\dfrac{1}{\sigma_k^2}}\rangle_D\otimes \left( \dfrac{1}{\sigma_k^2}|0\rangle_E+\sqrt{1-\left(\dfrac{1}{\sigma_k^2}\right)^2}|1\rangle_E \right) \end{equation} }$

$E$ レジスタを測定し0が観測されたとき、上の状態は次の状態に収縮する。

${\displaystyle \begin{equation} c\sum_{k=1}^r\dfrac{1}{\sigma_k} |{\bf v}_k\rangle_A\otimes|{\bf u}_k\rangle_B\otimes|\sigma_k^2\rangle_C\otimes|\dfrac{1}{\pi}\cos^{-1}{\dfrac{1}{\sigma_k^2}}\rangle_D\otimes |0\rangle_E \end{equation} }$

$c$ は規格化定数

${\displaystyle \begin{equation} c=\dfrac{1}{\sum_k\sigma_k^{-2}} \end{equation} }$

である。各種演算の逆演算を施し $C,D$ レジスタを初期化する。

${\displaystyle \begin{equation} c\sum_{k=1}^r\dfrac{1}{\sigma_k} |{\bf v}_k\rangle_A\otimes|{\bf u}_k\rangle_B\otimes|0\rangle_C\otimes|0\rangle_D\otimes |0\rangle_E \end{equation} }$

$A,B$ レジスタ部分を取り出し、これを $|\psi_1\rangle$ と置く。

${\displaystyle \begin{equation} |\psi_1\rangle=c\sum_{k=1}^r\dfrac{1}{\sigma_k} |{\bf v}_k\rangle_A\otimes|{\bf u}_k\rangle_B \end{equation} }$

ここまでの計算で式(3)を振幅に埋め込んだ状態が得られたことになる。さて、未知データ $\tilde{{\bf x}}$ が与えられたときの出力値 $\tilde{y}$ を計算するため次の状態を用意する。

${\displaystyle \begin{eqnarray} |\psi_2\rangle&=&|\tilde{{\bf x}}\rangle_A\otimes|{\bf y}\rangle_B\\ &=& \sum_{j=1}^M\tilde{x}_j|j\rangle_A\otimes\sum_{i=1}^Ny_i|i\rangle_B \end{eqnarray} }$

内積 $\langle\psi_1|\psi_2\rangle$ を計算すると

${\displaystyle \begin{eqnarray} \langle\psi_1|\psi_2\rangle&=& \left( c\sum_{k=1}^r\dfrac{1}{\sigma_k}\langle{\bf v}_k|_A\otimes\langle{\bf u}_k|_B\right) \left( \sum_{j=1}^M\tilde{x}_j|j\rangle_A\otimes\sum_{i=1}^Ny_i|i\rangle_B \right)\\ &=& c\sum_{k=1}^r\dfrac{1}{\sigma_k} \sum_{j=1}^M \tilde{x}_j \langle{\bf v}_k|_A|j\rangle_A \otimes \sum_{i=1}^Ny_i \langle{\bf u}_k|_B|i\rangle_B\\ &=& c\sum_{k=1}^r\dfrac{1}{\sigma_k} \sum_{j=1}^M v_{jk}\tilde{x}_j \sum_{i=1}^Nu_{ik}y_i \end{eqnarray} }$

を得る。これは式(4)で得た出力値 $\tilde{y}$ を $c$ 倍したものである。 $c$ は求めることができる値なので、内積 $\langle\psi_1|\psi_2\rangle$ を計算できれば、未知データに対する出力を得ることができる。今回は、内積計算の量子アルゴリズムには言及しない。

まとめ

　線形回帰の量子アルゴリズムについて説明した。上の量子位相推定を行う際に $e^{i2\pi\rho_A}$ を $|{\bf v}_k\rangle$ に施す必要がある。このとき、ここで説明したDensity Matrix Exponentiationを利用する。

参考文献

2021-08-07

Density Matrix Exponentiation

量子コンピュータ

はじめに

　Density Matrix Exponentiationは、密度演算子 $\rho$ をハミルトニアンとする時間発展演算子 $\exp{\left(-i\rho T\right)}$ を量子状態に作用させる際に使われるテクニックである。

微小時間への分解

　 $N$ を十分大きな数として、 $T=N\delta$ と置くと

${\displaystyle \begin{equation} \exp{\left(-i\rho T\right)}=\bigl(\exp{\left(-i\rho\delta\right)}\bigr)^N \end{equation} }$

と書ける。つまり、微小時間 $\delta$ の時間発展を $N$ 回繰り返すことを考える。時間発展させる量子状態を $|\psi\rangle$ とすると、 $\delta$ 後の状態は

${\displaystyle \begin{equation} |\psi^{\prime}\rangle=\exp{\left(-i\rho\delta\right)}|\psi\rangle \end{equation} }$

となる。この状態に対応する密度演算子は

${\displaystyle \begin{equation} |\psi^{\prime}\rangle\langle\psi^{\prime}|=\exp{\left(-i\rho\delta\right)}|\psi\rangle\langle\psi|\exp{\left(i\rho\delta\right)} \end{equation} }$

である。 $\sigma^{\prime}=|\psi^{\prime}\rangle\langle\psi^{\prime}|$ 、 $\sigma=|\psi\rangle\langle\psi|$ と置くと

${\displaystyle \begin{equation} \sigma^{\prime}=\exp{\left(-i\rho\delta\right)}\;\sigma\;\exp{\left(i\rho\delta\right)} \end{equation} }$

を得る。 $\exp{\left(\pm i\rho\delta\right)}$ を級数展開して上式に代入すると

${\displaystyle \begin{equation} \sigma^{\prime}=\sigma+i[\sigma, \rho]\delta+\mathcal{O}(\delta^2) \end{equation} }$

となる。 $\delta$ は微小量であるから、これの1次までとり近似する。

${\displaystyle \begin{equation} \sigma^{\prime}\sim\sigma+i[\sigma, \rho]\delta\tag{1} \end{equation} }$

この近似式をユニタリー演算子だけを用いて作る手順が、Density Matrix Exponentiationである。

Density Matrix Exponentiation

　最初に、同じ数の量子ビットからなら2つのレジスタA,Bを考える。Aレジスタに密度演算子 $\rho$ を、Bレジスタに密度演算子 $\sigma$ を振幅エンコーディングする。

${\displaystyle \begin{equation} \rho_A\otimes\sigma_B\tag{2} \end{equation} }$

次に、AレジスタとBレジスタの内容を交換するスワップ演算子

${\displaystyle \begin{equation} S_{AB}=\sum_{i,j}|i\rangle_A\langle j|_A\otimes|j\rangle_B\langle i|_B\tag{3} \end{equation} }$

を考える。これが2つのレジスタ間の状態を交換することは以下のように確認することができる。いま

${\displaystyle \begin{eqnarray} |u\rangle_A&=&\sum_i u_i|i\rangle_A\\ |v\rangle_B&=&\sum_i v_i|i\rangle_B \end{eqnarray} }$

とすると

${\displaystyle \begin{eqnarray} S_{AB}|u\rangle_A\otimes|v\rangle_B&=&\sum_{i,j}|i\rangle_A\langle j|_A\otimes|j\rangle_B\langle i|_B \sum_m u_m|m\rangle_A\otimes\sum_n v_n|n\rangle_B\\ &=& \sum_{i,j}\sum_{m,n}u_m v_n|i\rangle_A \langle j|_A|m\rangle_A\otimes |j\rangle_B \langle i|_B|n\rangle_B\\ &=& \sum_{i,j}\sum_{m,n}u_m v_n|i\rangle_A\;\delta_{j,m}\otimes |j\rangle_B\;\delta_{i,n}\\ &=& \sum_{i,j}u_j v_i|i\rangle_A\otimes |j\rangle_B\\ &=& \sum_{i}v_i|i\rangle_A\otimes \sum_j u_j|j\rangle_B\\ &=& |v\rangle_A\otimes|u\rangle_B \end{eqnarray} }$

となり、確かにAレジスタとBレジスタの内容が交換されたことが分かる。同じ計算を繰り返せば

${\displaystyle \begin{equation} S_{AB}^2=I_{AB}\tag{4} \end{equation} }$

を示すことができる（元の状態に戻るということ）。

　さて、式(2)と(3)を用いた次式を考える。

${\displaystyle \begin{equation} F= e^{-iS_{AB}\;\delta}\left(\rho_A\otimes\sigma_B\right)e^{iS_{AB}\;\delta} \end{equation} }$

式(4)から

${\displaystyle \begin{equation} e^{\pm iS_{AB}\;\delta}=I_{AB}\cos{\delta}\pm iS_{AB}\sin{\delta} \end{equation} }$

が成り立つので、これらを $F$ に代入すると

${\displaystyle \begin{equation} F=\left(\rho_A\otimes\sigma_B\right)\cos^{2}{\delta} +\bigl(S_{AB}\left(\rho_A\otimes\sigma_B\right)S_{AB}\bigr)\sin^{2}{\delta}+i \bigl( \left(\rho_A\otimes\sigma_B\right)S_{AB}-S_{AB}\left(\rho_A\otimes\sigma_B\right) \bigr)\sin{\delta}\cos{\delta} \end{equation} }$

を得る。いま

${\displaystyle \begin{equation} S_{AB}\left(\rho_A\otimes\sigma_B\right)S_{AB}=\sigma_A\otimes\rho_B \end{equation} }$

が成り立つので

${\displaystyle \begin{equation} F=\left(\rho_A\otimes\sigma_B\right)\cos^{2}{\delta} + \left(\sigma_A\otimes\rho_B\right)\sin^{2}{\delta} +i \bigl( \left(\rho_A\otimes\sigma_B\right)S_{AB}-S_{AB}\left(\rho_A\otimes\sigma_B\right) \bigr)\sin{\delta}\cos{\delta} \end{equation} }$

を得る。レジスタAの空間だけでトレースを取る。

${\displaystyle \begin{eqnarray} {\rm Tr}_A\left(F\right) &=& {\rm Tr}_A \left( \rho_A\otimes\sigma_B \right) \cos^{2}{\delta} + {\rm Tr}_A \left( \sigma_A\otimes\rho_B \right) \sin^{2}{\delta}\\ &&\hspace{1cm}+i \left[ {\rm Tr}_A \bigl( \left(\rho_A\otimes\sigma_B\right)S_{AB} \bigr) - {\rm Tr}_A \bigl( S_{AB}\left(\rho_A\otimes\sigma_B\right) \bigr) \right] \sin{\delta}\cos{\delta} \end{eqnarray} }$

いま

${\displaystyle \begin{eqnarray} {\rm Tr}_A \left( \rho_A\otimes\sigma_B \right)&=&\sigma_B\\ {\rm Tr}_A \left( \sigma_A\otimes\rho_B \right)&=&\rho_B\\ {\rm Tr}_A \bigl( \left(\rho_A\otimes\sigma_B\right)S_{AB} \bigr) &=&(\sigma\rho)_B\\ {\rm Tr}_A \bigl( S_{AB}\left(\rho_A\otimes\sigma_B\right) \bigr) &=& (\rho\sigma)_B \end{eqnarray} }$

が成り立つ。右辺はすべてレジスタBの空間での量である。従って

${\displaystyle \begin{equation} {\rm Tr}_A\left(F\right)=\sigma\cos^{2}{\delta}+\rho\sin^{2}{\delta}+i[\sigma,\rho]\cos{\delta}\sin{\delta} \end{equation} }$

を得る（添字 $B$ は省いた）。 $\delta$ を十分小さな量としているので、 $\delta$ の1次までの量で近似すると

${\displaystyle \begin{equation} {\rm Tr}_A\left(F\right)\sim\sigma+i[\sigma,\rho]\delta \end{equation} }$

となる。これは求めたかった式(1)である。

Tech Blog

勉強したことをまとめます。

Denoising Diffusion Probabilistic Models

はじめに

拡散モデル

変分推論

$q(x_t|x_0)$ の計算

$q(x_{t-1}|x_t,x_0)$ の計算

最小化すべき関数の計算

まとめ

QSのAcademic Reputationについて

はじめに

Academic Reputation

国内の大学の順位

2022年のアジア地域の順位

2022年の満点大学

まとめ

cpplinqサンプル集

はじめに

サンプル集

boost/range/adaptorサンプル集

はじめに

サンプル集

2つの自然数が互いに素である確率

はじめに

問題

解答

まとめ

参考動画

量子線形回帰

はじめに

線形回帰

線形回帰の量子アルゴリズム

まとめ

参考文献

Density Matrix Exponentiation

はじめに

微小時間への分解

Density Matrix Exponentiation

参考文献