ベイジアンネットワークのノードをテンソルネットワークで表現する

[前回の記事](https://qiita.com/yoshi_mino/items/fa646bda7a002fa4990c)ではベイジアンネットワークをシンプルな方法で量子回路で表現して、
ネットワークが深くなるほど、また親ノードの数が多くなるほど回路が深くなることを確認しました。

今回は試行的にテンソルネットワークでベイジアンネットワークを表現します。
ここで試行的としたのは、ベイジアンネットワークではデータ規模が小さく、テンソルネットワークを導入するメリットが見えてこないためです。

本記事の執筆に当たり、理論面では文献[1]、実装面では文献 [2][3] を参考にしました。
実装に便利なライブラリが存在しないので、Python + numpy を駆使しなければなりません。かなり複雑なコードになってしまったので、本記事ではコードレベルでの実装の詳細は省きます。

# 表現するデータ
ベイジアンネットワークのノードを1つ取り出した、図のような 4 入力 1 出力のノードについて考えます。
各入出力は 0,1 の 2 通りの状態を持つものとします。
このようなノードを複数用意して結合すれば、CANCER のようなベイジアンネットワークを表現することができるはずです。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/b09c7be9-06c4-6be3-5a97-2c04e3299dfb.png" width="200px"/>


# テンソルネットワークの形式
テンソルネットワークにはいろいろな形式がありますが、ここでは最も基本的な Matrix Product State (MPS) 形式を採用します。上部にある 4 つの矢印 ($x_1,x_2,x_3,x_4$) が入力を、下部にある 1 つの矢印 $f(x)$ が出力を表します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/681d5791-fa01-56c1-ce36-a8305e4fbb98.png" width="400px" />

MPS の入力 $x_1,x_2,x_3,x_4$ と出力 $f(x)$ は状態が 0 になる確率と 1 になる確率を表現する 2 次元ベクトルとします。
例えば、状態 0 になる確率が 0.3 で 1 になる確率が 0.7 のとき、ベクトル $(0.3, 0.7)$ になります。

図では便宜上、ノードを円形、テンソルを矩形、入出力を矢印で表記することで区別しています。

# パラメータ
ここで使用する MPS は以下に示す数の変数値を持ちます。
- 赤線部分のボンド次元が 2 の場合、$2^2 + 2^3 \times 3 = 28$
- 赤線部分のボンド次元が 4 の場合、$2^2 + 2^2 \times 4 + 2\times 4^2 + 4 \times 2^2 = 68$

ここからは、ノードと MPS の入出力が一致するように MPS の変数値（パラメータ）を調整していきます。
すべてのテンソルの足の数を同数にするため、テンソル $A$ にダミーの足を追加してベクトル $(1,0)$ を接続しておきます。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/e1ce2c03-8706-056c-a293-d948a0960a80.png" width="400px" />

# テンソルネットワークの変形
パラメータの調整過程において、次の手順に示すような縮約と特異値分解の操作を繰り返します。
これは Density-Matrix Renormalization Group (DMRG) という手法を元にしたもので、このような操作をテンソルの「繰り込み」と言うようです。

1. 2つのテンソル $A,B$ を縮約して、4 脚テンソル $Y$ にします。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/49d01cea-9df2-3577-7009-e770876e1b83.png" width="350px" /> <br>
2. $Y$ を特異値分解して、元の形式に復元します。
$Y$ を行列 $Y^\prime$ に変形して $Y^\prime = Q \Sigma V$ の形式で特異値分解した後、$Q$ と $V$ を 3 脚に変形します。（図の左上）
その後、$\Sigma$ と $V$ を縮約して改めて $Q \rightarrow A, \Sigma V \rightarrow B$ とします。(図の右上)

3. 次に $B, C$ を縮約して改めて $Y$ を作り、同様の操作を行います。最後に $C,D$ を縮約して同様の操作を行うと特異値が右端に移動します。（図の左下）
![mps5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/b083e5d9-a208-9f87-d277-d7406264e745.png)
4. 再び $C,D$ を縮約して特異値分解し、今度は $Q\Sigma \rightarrow C, V \rightarrow D$ とすると、特異値を左に移動させることができます。（図の右下）

この繰り込み操作によって特異値行列 $\Sigma$ を MPS の任意の位置に移動させることができます。特に図の左下の状態は left-canonical といい、$A,B,C$ を次のように縮約すると単位行列になるので計算が簡単になることがあります。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/afbac11b-5719-6ce6-ed4b-207d89d9ae2e.png" width="400px" />

# パラメータの調整方法
ノードの入出力の値を学習データとして、機械学習で行われる方法と同じように勾配降下法で MPS の値を調整していきます。
例えば、テンソル $B,C$ を縮約して 4 脚のテンソル $Y$ にした状態では、以下の計算式で勾配が求められます。
（$f(x)$ と $y$ がベクトルであることを明示するために一部で添字を追加しています。）
<div>
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/d7ff3824-25d8-7aca-40f5-c3b390343019.png" width="350px"/>
</div>

- MPS の出力：$f(x)_l = \sum\_{ijk} A\_{ij}^{x_1} Y^{x_2 x_3}\_{jk} D^{x_4}\_{kl}$
- ノードの出力（学習データ）：$y_l$
- ノードと MPS の出力誤差（MSE）：$L = \frac{1}{2}(f(x) - y)^2$
- $L$ の勾配：
$$
\frac{\partial L}{\partial Y^{x_2 x_3}_{jk}} = \frac{\partial f(x)}{\partial Y^{x_2 x_3}\_{jk}} \frac{\partial L}{\partial f(x)} = \sum\_{il} (f(x) - y)_l A^{x_1}\_{ij} D^{x_4}\_{kl}
$$

この勾配をもとにテンソル $Y$ の値を更新します。
テンソルの繰り込みによって $Y$ の位置を左右に動かすことで、すべてのテンソルの値を調整していきます。

# 実験
実際に MPS のパラメータ調整を行なった後、ノードと MPS の出力値を確認する実験を行いました。

結果、
- MPS のボンド次元が 2 のとき、出力の二乗誤差の総和が $0.294$
- MPS のボンド次元が 4 のとき、出力の二乗誤差の総和が $1.75 \times 10^{-7}$
- MPS 出力の正規化を行なっていないため、出力の 2 つの値を足し合わせた結果が 1 にならず、値が 1 を超える場合がある。（何度も実験を行うと、負になる場合もあった。）

となりました。詳細は下表の通りです。

| 入力 $(x_1, x_2, x_3, x_4)$ | ノード出力 $(y_1, y_2)$ | MPS 出力（ボンド次元 2）     | MPS 出力（ボンド次元 4）    |
| ------------------------- | ------------------ | ------------------- | ------------------ |
| $(0, 0, 0, 0)$            | $(0.2, 0.8)$       | $(0.2000, 0.8159)$  | $(0.1999, 0.7998)$ |
| $(0, 0, 0, 1)$            | $(0.3, 0.7)$       | $(0.1526, 0.8532)$  | $(0.3001, 0.7002)$ |
| $(0, 0, 1, 0)$            | $(0.9, 0.1)$       | $(0.6799, 0.2541)$  | $(0.8999, 0.1000)$ |
| $(0, 0, 1, 1)$            | $(0.25, 0.75)$     | $(0.4939, 0.5272)$  | $(0.2500, 0.7500)$ |
| $(0, 1, 0, 0)$            | $(0.12, 0.88)$     | $(0.2329, 0.7269)$  | $(0.1200, 0.8799)$ |
| $(0, 1, 0, 1)$            | $(0.27, 0.73)$     | $(0.1762, 0.7993)$  | $(0.2700, 0.7300)$ |
| $(0, 1, 1, 0)$            | $(0.5, 0.5)$       | $(0.4009, 0.6300)$  | $(0.5000, 0.5000)$ |
| $(0, 1, 1, 1)$            | $(0.05, 0.95)$     | $(0.2979, 0.8094)$  | $(0.0500, 0.9502)$ |
| $(1, 0, 0, 0)$            | $(0.2, 0.8)$       | $(0.1590, 0.8249)$  | $(0.2000, 0.7998)$ |
| $(1, 0, 0, 1)$            | $(0.3, 0.7)$       | $(0.1224, 0.8284)$  | $(0.3000, 0.7002)$ |
| $(1, 0, 1, 0)$            | $(0.9, 0.1)$       | $(0.9043, 0.0575)$  | $(0.8999, 0.1000)$ |
| $(1, 0, 1, 1)$            | $(0.75, 0.25)$     | $(0.6486, 0.3026)$  | $(0.7500, 0.2500)$ |
| $(1, 1, 0, 0)$            | $(0.05, 0.95)$     | $(0.1553, 0.8772)$  | $(0.0500, 0.9502)$ |
| $(1, 1, 0, 1)$            | $(0.01, 0.99)$     | $(0.1200, 0.8692)$  | $(0.0100, 0.9902)$ |
| $(1, 1, 1, 0)$            | $(0.9, 0.1)$       | $(1.0437, 0.02479)$ | $(0.8999, 0.1000)$ |
| $(1, 1, 1, 1)$            | $(0.85, 0.15)$     | $(0.7460, 0.2567)$  | $(0.8501, 0.1500)$ |

また、[前回の記事](https://qiita.com/yoshi_mino/items/fa646bda7a002fa4990c)で扱った CANCER ネットワークを MPS で表現する実験も行い、ボンド次元 4 ではほとんど誤差なく表現できることも確認しました。

# 課題
実装を簡略化するため、正規化を行いませんでしたが、本来は行うべきです。
その際、入出力のデータを $(\cos{\frac{\pi}{2}x}, \sin{\frac{\pi}{2}x})$ や $(x, 1-x)$ で符号化し、さらに KL ダイバージェンスなどの誤差関数を用いるといった方法が考えられます。

# おわりに
本記事では、ベイジアンネットワークを MPS で表現し、実際にシンプルなベイジアンネットワークであれば誤差なく MPS で表現できることを確認しました。

パラメータ調整の際に、テンソル繰り込みによる方法を採用しましたが、逆誤差伝播法でも調整が可能です。その場合、PyTorch などの強力なライブラリを使えるため実装が容易になります。
一方で、繰り込みの方法には特異値分解によって同時に次元圧縮ができるというメリットがあります。（問題に依存すると思われますが）勾配の計算がシンプルな点もメリットと言えると思います。

最初に書いたように、ベイジアンネットワークでは MPS を使うメリットは見えてきません。
では、物理モデル以外でどんなデータを扱えばメリットが見出だせるのか？
引き続き模索していきたいと思います。


### 参考文献
[1]: 西野 友年 (著)、テンソルネットワーク入門 (KS情報科学専門書)、講談社

[2]: Zhao-Yu Han, Jun Wang, Heng Fan, Lei Wang, Pan Zhang. Unsupervised Generative Modeling Using Matrix Product States. (https://arxiv.org/abs/1709.01662)

[3]: [Unsupervised Generative Modeling using Matrix Product States](https://github.com/congzlwag/UnsupGenModbyMPS)