ベイジアンネットワークを量子回路で表現する

本記事では単なる確率分布ではない、構造化されたデータを量子回路で表現する方法を紹介します。

代表的なデータ構造として、ベイジアンネットワークを取り扱います。
ベイジアンネットワークは、因果関係を視覚的に表現するループのない有向グラフです。

ベイジアンネットワークには多くの種類がありますが、
ここでは以下の図のようなシンプルな例で量子表現の方法を示します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/01e9e072-57db-fc98-3527-b7fda29c176f.png" width="35%">

A, B, C の 3 つの確率変数があり、それぞれ 0 または 1 の 2 種類の状態をとるものとします。

A と B が 1 になる確率には前提条件がありません。
一方、C が 1 になる確率は、A, B の状態に依存して決定される条件付き確率になります。
図において、この因果関係は矢印の向きで表現されています。

以下は [論文 Quantum circuit representation of Bayesian networks](https://arxiv.org/abs/2004.14803) に記載の手法を再現したものです。

# 親ノードを持たないノード (A, B) の表現
ノードの状態が 1 になる確率 $P(1)$、0になる確率 $P(0)$ として、回転角を定め、
$R_y$ ゲートを作用させることで表現できます。

回転角は
$$
2 \times \arctan(\sqrt{\frac{P(1)}{P(0)}})
$$
になります。ただし、$P(0)=0$ のときは回転角 $\pi$ とします。

（例）30%で状態が 1 になるとき
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/a06f08eb-e50d-6bf1-c330-6e85d865000d.png" width="50%">

# 親ノードを持つノード (C) の表現
親ノードに相当する量子ビットを制御ビットとする、制御 $R_y$ ゲートを作用させます。
回転角は $\sum_{i,j\in 0,1} P(C=1 | A=i, B=j)$ によって決定されます。
「親ノードの数＝制御ビットの数」になりますので、例では 2 個の制御ビットを持つ $CCR_y$ ゲートを使います。

# ノード (C) の量子回路
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/db35ee10-ecc8-88c6-bd19-cce293beab12.png">

$q_0$ がノード A、$q_1$ がノード B、$q_2$ がノード C に相当します。先にノード A、B を $R_Y(\theta_A), R_Y(\theta_B)$ でエンコードしておく必要があります。

$R_Y(\theta_C|00)$ は $A=0, B=0$ のときに作用する回転ゲートで、回転角は $P(C=1|A=0, B=0)$ によって決定されます。（残りの回転ゲートも同様）
この量子回路で $q_2$ を測定することで $P(C=1)$ を得ることができます。（ただし、厳密には量子振幅推定が必要）

親ノード数が一般に $n$ の場合、制御ビットを n 個持つ制御回転ゲートを $2^n$ 個使用する回路になります。そのため、回路規模はベイジアンネットワークの大きさに対して、指数的に大きくなってしまいます。

# ベイジアンネットワークの具体例

肺ガンに関する因果関係を表したベイジアンネットワーク (CANCER) を使って実装を行います。
このデータは [bnlearn](https://www.bnlearn.com/bnrepository/discrete-small.html) から入手することができます。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/23bef163-611c-c370-57e2-8831a34d6a71.png" width="60%">


## ノードの説明 
| ノード名      | 値                         |
| :-------- | ------------------------- |
| Pollution | 0 = 大気汚染度が低い、1 = 大気汚染度が高い |
| Smoker    | 0 = 非喫煙者、1 = 喫煙者          |
| Cancer    | 0 = ガン患者でない、1 = ガン患者      |
| Xray      | 0 = X線検査異常なし、1 = X線検査異常あり |
| Dyspnoea  | 0 = 呼吸困難症状なし、1 = 呼吸困難症状あり |

## 各状態の発生確率（実務的にはサンプルデータから推定する）

| ノード名      |  0  |  1  |
| :-------- | :-: | :-: |
| Pollution | 0.9 | 0.1 |
| Smoker    | 0.7 | 0.3 |

|          (pollution, Smoker)           | (0, 0) | (0, 1) | (1, 0) | (1, 1) |
| :------------------------------------: | :----: | :----: | :----: | :----: |
| $P($Cancer = 0 \| pollution, Smoker$)$ | 0.999  |  0.97  |  0.98  |  0.95  |
| $P($Cancer = 1 \| pollution, Smoker$)$ | 0.001  |  0.03  |  0.02  |  0.05  |

|           (Cancer)            | (0) | (1)  |
| :---------------------------: | :-: | :--: |
|   $P($Xray = 0 \| Cancer$)$   | 0.8 | 0.1  |
|   $P($Xray = 1 \| Cancer$)$   | 0.2 | 0.9  |
| $P($Dyspnoea = 0 \| Cancer$)$ | 0.7 | 0.35 |
| $P($Dyspnoea = 1 \| Cancer$)$ | 0.3 | 0.65 |

例えば、大気汚染度が低い地域に住む喫煙者が、ガン患者である確率は 0.03、
ガン患者のうち、呼吸困難症状がある確率は 0.65 になっています。

このベイジアンネットワークを表現する量子回路と、１万回測定を行ったヒストグラムは下の図のようになります。

![7-5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/c980f445-5b5c-14a9-10d0-986498008d45.png)

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/5dddfba7-9b78-1bff-65b6-c65401c9d01d.png" width="80%">

ヒストグラムで "00010" は 1443 回観測されているので、約 14.4% だと分かります。
検証のために実際に計算してみると、
$$
P(Pollution=0) \times P(Smoker=1) \times P(Cancer=0 | Pollution=0, Smoker=1)
$$

$$
\times P(Xray=0 | Cancer=0) \times P(Dyspnoea=0 | Cancer=0)
$$

$$
 = 0.9 \times 0.3 \times 0.97 \times 0.8 \times 0.7 \approx 0.1467
$$
で、ほぼ一致していることが分かります。

また、量子ビット $q_4$ だけ測定すれば、呼吸困難の（条件付きでない）確率を知ることができます。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2603676/f25cdbda-077d-958d-3299-4ac3db82bf6b.png" width="80%">

# 実装

qiskit を使って実装を行いました。まず、必要なライブラリを読み込みます。

``` python
import math
import numpy as np

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit_aer import AerSimulator
from qiskit.visualization import plot_histogram

n_shots = 10_000 # シミュレーターでのサンプリング回数
backend_sim = AerSimulator() # シミュレーターの用意
```

続いて、制御$R_y$ゲートを実装します。
制御ビットが1個の場合、2個の場合、3個以上の場合で分けています。
制御ビットが3個以上の場合には、アンシラを使用します。
``` python
def multicontrol_ry(controls, theta) -> tuple[QuantumCircuit, list, int]:
	"""
	制御ビットが3個以上のマルチ制御Ryゲート
	
	Parameters:
	controls (int) : 制御ビットの数
	theta (float) : Ryゲートの回転角
	Returns:
	QuantumCircuit : 量子回路（量子ビット数は制御ビット数の2倍）
	0 から controls-1 --> 制御ビット
	controls から 2*controls-2 --> アンシラ
	2*controls-1 --> ターゲット
	"""
	assert controls >= 3 # 制御ビット数が2のときはアンシラ不要
	circ = QuantumCircuit(controls*2)
	
	circ.ccx(0, 1, controls)
	for i in range(controls - 2):
		circ.ccx(2 + i, controls + i, controls + 1 + i)
	circ.ry(theta/2, controls*2 - 1)
	circ.cx(controls*2 - 2, controls*2 - 1)
	circ.ry(-theta/2, controls*2 - 1)
	circ.cx(controls*2 - 2, controls*2 - 1)
	for i in reversed(range(controls - 2)):
		circ.ccx(2 + i, controls + i, controls + 1 + i)
	circ.ccx(0, 1, controls)
	
	return circ

def ccry(theta) -> QuantumCircuit:
	"""
	CCRy ゲート
	"""
	circ = QuantumCircuit(3)
	circ.ry(theta/2, 2)
	circ.ccx(0, 1, 2)
	circ.ry(-theta/2, 2)
	circ.ccx(0, 1, 2)
	return circ
  
def cry(theta) -> QuantumCircuit:
	"""
	CRy ゲート
	"""
	circ = QuantumCircuit(2)
	circ.ry(theta/2, 1)
	circ.cx(0, 1)
	circ.ry(-theta/2, 1)
	circ.cx(0, 1)
	return circ
```

親ノードを1個持つノードをエンコードする回路を実装します。

``` python
def encode_1(probs) -> QuantumCircuit:
	"""
	親ノード1個のエンコードをする。2量子ビット回路になる。
	
	Parameters:
	probs (list) : 親ノードをA、子ノードをBとして、P(B=1|A=0)、P(B=1|A=1) に対応する確率のリスト
	Returns:
	QuantumCircuit : 量子回路
	"""
	circ = QuantumCircuit(2)
	
	circ.x(0)
	if probs[0] != 1:
		circ.append(cry(2*math.atan(math.sqrt(probs[0]/(1 - probs[0])))), [0, 1])
	else:
		circ.append(cry(np.pi), [0, 1])
	circ.x(0)
	if probs[1] != 1:
		circ.append(cry(2*math.atan(math.sqrt(probs[1]/(1 - probs[1])))), [0, 1])
	else:
		circ.append(cry(np.pi), [0, 1])
	
	return circ
```

親ノードを2個以上持つノードをエンコードする回路を実装します。

``` python
def encode_2state_bayesian(n, probs) -> QuantumCircuit:
	"""
	ベイジアンネットワークにおいて、n個の親ノードを持つノードをエンコードする。
	n+1 個の量子ビットでネットワークを表現し、エンコードの過程で (n+1) ビットの制御Ry回転を実装するために、
	(n-1) 個のアンシラを必要とする。合計で量子ビット数 2n の回路になる。
	
	q_0 から q_(n-1) までの量子ビットは既にエンコードが済んでいる前提で、
	|0> 状態にある量子ビット q_n を回転させる。
	
	Parameters:
	thetas (list) : P(q_n = 1)となる確率が格納された長さ 2^n のリスト（回転角ではない）
	リストの内容は次の通り
	P(q_n = 1 | q_0=q_1=...q_(n-1)=0), P(q_n = 1 | q_0=q_1=...q_(n-2)=0,q_(n-1)=1), ...
	Returns:
	QuantumCircuit : 量子回路 q_0 から q_(n-1) まで親ノード、q_n がエンコード対象となるノード、
	q_(n+1) 以降がアンシラ。n < 3 ではアンシラなし
	"""
	assert n > 1 and len(probs) == 2**n
	if n == 2:
		circ = QuantumCircuit(3)
	else:
		circ = QuantumCircuit(2*n)
	for state in range(2**n):
		# Xゲート
		for i in range(n):
			if state & 2**i == 0:
				circ.x(n - i - 1)
		# 制御Ryゲート
		if probs[state] != 1:
			angle = 2*math.atan(math.sqrt(probs[state]/(1 - probs[state])))
		else:
			angle = np.pi
		if n == 2:
			mc = ccry(angle)
			circ.append(mc, [0, 1, 2])
		else:
			mc = multicontrol_ry(n, angle)
			# append の2番目の引数は、回路 mc の各ビットを回路 circ の何番のビットに接続するかを示す。
			circ.append(mc, list(range(n)) + list(range(n + 1, 2*n)) + [n])
		
		# Xゲートを元に戻す
		for i in range(n):
			if state & 2**i == 0:
				circ.x(n - i - 1)
	
	return circ
```

これら関数を使って、CANCER ベイジアンネットワークをエンコードします。

``` python
q = QuantumRegister(5)
c = ClassicalRegister(5)
circ = QuantumCircuit(q, c)
  
# Pollution
circ.ry(2*math.atan(math.sqrt(0.1/0.9)), 0)
# Smoker
circ.ry(2*math.atan(math.sqrt(0.3/0.7)), 1)
  
# Cancer
probs = [0.001, 0.03, 0.02, 0.05]
sp_circ = encode_2state_bayesian(2, probs)
circ.append(sp_circ, [0, 1, 2])
  
# Xray
Xray = [0.2, 0.9]
circ_Xray = encode_1(Xray)
circ.append(circ_Xray, [2, 3])
  
# Dyspnoea
Dyspnoea = [0.3, 0.65]
circ_Dyspnoea = encode_1(Dyspnoea)
circ.append(circ_Dyspnoea, [2, 4])
  
circ.measure([q[0], q[1], q[2], q[3], q[4]], [0, 1, 2, 3, 4])
```

シミュレーターを使って実行し、ヒストグラムを取得します。
``` Python
result_ideal = backend_sim.run(circ.decompose().decompose(), shots=n_shots).result()
plot_histogram(result_ideal.get_counts(0))
```

# 最後に

本記事ではベイジアンネットワークのような構造化されたデータを量子回路で表現できることを示しました。
しかし、紹介した表現方法は、多数のゲート演算が必要で
量子コンピュータのエラーに対する考慮もないため、実機で動作させることはできません。

実機で意味のある計算を行うためには、異なる手法を探求しなければなりません。