量子回路のテンソルネットワークシミュレーション　〜縮約密度行列の計算〜

# はじめに

この記事は「[量子回路のテンソルネットワークシミュレーション　〜特定の状態の振幅の計算〜](https://qiita.com/kazuo_watanabe/items/97a30806f40386cb8e7f)」の続編です。

量子アルゴリズムの中には計算に用いる全ての量子ビットの観測を必要とせず、特定の量子ビットのみを観測し、その出力を得ることで、目的を達成することができるものも存在します。

今回の記事では、テンソルネットワークを活用した量子回路のシミュレーション手法の内、特定の量子ビットに着目し、その出力を測定するような手法について実装します。

本記事の内容は、東京大学「量子ソフトウェア」寄付講座[[1]](#参考文献)内で実施された第3回量子ソフトウェアハンズオン(産学協働ゼミ)[[2]](#参考文献)の演習資料を元に作成しています。


# Quantum Bit String Comparator


今回題材とするのは、Quantum Bit String Comparator[[3]](#参考文献)と呼ばれるアルゴリズムで、二つの$n$ビットの量子ビット列$\lvert a \rangle ,\lvert b \rangle$を入力として受け取り、その大小関係を比較するアルゴリズムとなっています。このとき、$\lvert a \rangle ,\lvert b \rangle$は重ね合わせ状態でも問題ありません。

このアルゴリズムでは、入力ビット長に関わらず、2量子ビットのみを観測すればよく、この値によって、大小関係を判定できます。

具体的には、以下のような出力が得られます。
- $a>b$のとき$\lvert 10 \rangle $
- $b>a$のとき$\lvert 01 \rangle $
- $a=b$のとき$\lvert 00 \rangle $

重ね合わせ状態のときは、その出現確率に合わせて出力が得られます。  
たとえば、$\lvert a \rangle = \lvert 10 \rangle$, $\lvert b \rangle = \frac{1}{\sqrt{2}} (\lvert 00 \rangle + \lvert 11 \rangle)$の場合、得られる出力は$\frac{1}{\sqrt{2}} (\lvert 10 \rangle + \lvert 01 \rangle)$となります。

ただし、このアルゴリズムは量子加速するアルゴリズムではないです。

## 量子回路の概要

回路は二つの入力量子ビット列$\lvert a \rangle,\lvert b \rangle$の各量子ビットの値を比較していくような形状となっています。以下に、入力ビット長が3の場合の例を示します。

![回路の概要](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/8d6abe15-7d83-0e77-500c-6148cec380c3.png)

（図は文献[[3]](#参考文献)中 Fig.1を参考に作成）

図中の$\lvert O_1 O_2 \rangle $が出力に相当します。また、図中の白丸で表現されたTofforiゲートは、制御ビットが$\lvert 00 \rangle $のときに、標的ビットを反転させるものです。

今回はこの回路を元にテンソルネットワークを作成します。

## テンソルネットワークでの表現

今回の場合、観測したい2量子ビットの状態のみを得られれば十分であり、その他の量子ビットの状態については計算する必要がありません。

数式を用いた理論の説明は本記事では割愛しますが、以下のような元の回路を左右反転させてくっつけた形状のテンソルネットワークを構築することで、この2量子ビットに関する縮約密度行列というものを計算することができ、この行列の対角成分を確認することで、求めたい確率を得ることができます。

![テンソルネットワークのイメージ図](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/a14f0f0a-15d2-91cc-ce86-589acc0687c3.png)


# 実装例

それでは上記のテンソルネットワークを元に、Python上での実装を行います。

通常、量子コンピュータ上でトフォリゲートを実装する場合、複数の2量子ビット以下のゲートの組み合わせに分解する必要がありますが、今回のシミュレーションでは、$2^3 \times 2^3$の大きな行列として、そのまま扱うこととします。

まずは、$\lvert a \rangle = \lvert 100 \rangle$、$\lvert b \rangle = \lvert 000 \rangle$の簡単な場合を実装します。

■分析環境
```
Python 3.7.6
tensornetwork 0.4.6
```

最初に、基本的なゲートをnumpyのarray形式で定義しておきます。

```Python
import numpy as np
import tensornetwork as tn

# パラメータ設定
bit_length = 3 # 入力するビット列の長さ
n_qubits = bit_length * 5 - 1 # 回路に必要な量子ビット数

# 各種ゲートの定義
h_gate = np.array([[1, 1], [1, -1]]) / np.sqrt(2)

x_gate = np.array([[0, 1], [1, 0]])

cx_gate = np.array([[1, 0, 0, 0], 
                    [0, 1, 0, 0], 
                    [0, 0, 0, 1], 
                    [0, 0, 1, 0]])

toffori_gate = np.array([[1, 0, 0, 0, 0, 0, 0, 0],
                         [0, 1, 0, 0, 0, 0, 0, 0],
                         [0, 0, 1, 0, 0, 0, 0, 0],
                         [0, 0, 0, 1, 0, 0, 0, 0],
                         [0, 0, 0, 0, 1, 0, 0, 0],
                         [0, 0, 0, 0, 0, 1, 0, 0],
                         [0, 0, 0, 0, 0, 0, 0, 1],
                         [0, 0, 0, 0, 0, 0, 1, 0]])
```

まずは、初期状態に対応するノードを定義します。今回は、実装を簡単にするため、入力するビット列は$\lvert 000 \rangle$からゲート操作を通じて用意するのではなく、あらかじめノードにセットしておきます。

```Python
initial_nodes = []

# ここでは直接、ノードに入力する初期状態を変えている
# ここを変更することで、初期状態を変更可能
for i in range(n_qubits):
    if i == 0: # 添字0がaの最上位ビット、添字1がbの最上位ビットになる 
        initial_nodes.append(tn.Node(np.array([0, 1]))) 
    else:
        initial_nodes.append(tn.Node(np.array([1, 0])))
```

ただし、このやり方では、量子もつれの状態は表現できません。量子もつれの状態を入力したい場合には、ゲートのノードを作る必要があります。


次に、アルゴリズム本体の実装をしていきます。以下では見通しをよくするために、ゲート操作を3区間に分けて考え、一つずつ実装していきます。

以下にQiskitで描画した今回の回路の図を示します。barrierで区切っている範囲ごとに実装していきます。

![量子回路のイメージ図](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/fb39fee2-8b83-a0a0-e0f2-9b949cf3e2dd.png)


```Python
first_layer_gate_nodes = []

for i in range(bit_length):
    first_layer_gate_nodes.append(tn.Node(x_gate)) # トフォリゲート前のX

    first_layer_gate_nodes.append(tn.Node(toffori_gate.reshape(2,2,2,2,2,2))) # 1つ目のトフォリゲート

    first_layer_gate_nodes.append(tn.Node(x_gate)) # トフォリゲート間のXゲート(1ビット目)
    first_layer_gate_nodes.append(tn.Node(x_gate)) # トフォリゲート間のXゲート(2ビット目)

    first_layer_gate_nodes.append(tn.Node(toffori_gate.reshape(2,2,2,2,2,2))) # 2つ目のトフォリゲート

    first_layer_gate_nodes.append(tn.Node(x_gate)) # トフォリゲート後のXゲート
```

```Python
second_layer_gate_nodes = []

for i in range(bit_length - 1):
    second_layer_gate_nodes.append(tn.Node(x_gate)) 
    second_layer_gate_nodes.append(tn.Node(x_gate))
    second_layer_gate_nodes.append(tn.Node(toffori_gate.reshape(2,2,2,2,2,2)))
    second_layer_gate_nodes.append(tn.Node(x_gate))
    second_layer_gate_nodes.append(tn.Node(x_gate))
```

```Python
third_layer_gate_nodes = []

for i in range(bit_length - 1):
    third_layer_gate_nodes.append(tn.Node(toffori_gate.reshape(2,2,2,2,2,2)))
    third_layer_gate_nodes.append(tn.Node(toffori_gate.reshape(2,2,2,2,2,2)))
```

反対方向から同じ形状のテンソルネットワークで挟む必要があるので、同じものをもう1セット作成します。


```Python
import copy
initial_nodes_r = copy.deepcopy(initial_nodes)
first_layer_gate_nodes_r = copy.deepcopy(first_layer_gate_nodes)
second_layer_gate_nodes_r = copy.deepcopy(second_layer_gate_nodes)
third_layer_gate_nodes_r = copy.deepcopy(third_layer_gate_nodes)
```

つづいて、これらをつなぐエッジを定義していきます。
同じ処理を反対側のネットワークでもやりたいので、関数化しています。

```Python
def make_edges(initial_nodes, first_layer_gate_nodes, second_layer_gate_nodes, third_layer_gate_nodes):
    # initial_nodes -> first_layer_nodes
    for i in range(bit_length):
        tn.connect(initial_nodes[i*5][0], first_layer_gate_nodes[i*6 + 1][0])
        tn.connect(initial_nodes[i*5 + 1][0], first_layer_gate_nodes[i*6][0])
        tn.connect(initial_nodes[i*5 + 2][0], first_layer_gate_nodes[i*6 + 1][2])
        tn.connect(initial_nodes[i*5 + 3][0], first_layer_gate_nodes[i*6 + 4][2])

        if i != bit_length - 1:
            tn.connect(initial_nodes[i*5 + 4][0], second_layer_gate_nodes[i*5 + 2][2])

    # in first_layer_nodes
    for i in range(bit_length):
        tn.connect(first_layer_gate_nodes[i*6 + 0][1], first_layer_gate_nodes[i*6 + 1][1])
        tn.connect(first_layer_gate_nodes[i*6 + 1][3], first_layer_gate_nodes[i*6 + 2][0])
        tn.connect(first_layer_gate_nodes[i*6 + 1][4], first_layer_gate_nodes[i*6 + 3][0])

        tn.connect(first_layer_gate_nodes[i*6 + 2][1], first_layer_gate_nodes[i*6 + 4][0])
        tn.connect(first_layer_gate_nodes[i*6 + 3][1], first_layer_gate_nodes[i*6 + 4][1])

        tn.connect(first_layer_gate_nodes[i*6 + 4][3], first_layer_gate_nodes[i*6 + 5][0])

    # first_layer_nodes -> second_layer_nodes
    for i in range(bit_length - 1):
        tn.connect(first_layer_gate_nodes[i*6 + 1][5], second_layer_gate_nodes[i*5 + 0][0])
        tn.connect(first_layer_gate_nodes[i*6 + 4][5], second_layer_gate_nodes[i*5 + 1][0])

    # in second_layer_nodes
    for i in range(bit_length - 1):
        tn.connect(second_layer_gate_nodes[i*5 + 0][1], second_layer_gate_nodes[i*5 + 2][0])
        tn.connect(second_layer_gate_nodes[i*5 + 1][1], second_layer_gate_nodes[i*5 + 2][1])

        tn.connect(second_layer_gate_nodes[i*5 + 2][3], second_layer_gate_nodes[i*5 + 3][0])
        tn.connect(second_layer_gate_nodes[i*5 + 2][4], second_layer_gate_nodes[i*5 + 4][0])

    # second_layer_nodes -> third_layer_nodes
    for i in range(bit_length - 1):
        tn.connect(second_layer_gate_nodes[i*5 + 2][5], third_layer_gate_nodes[i*2 + 0][1])
        tn.connect(second_layer_gate_nodes[i*5 + 3][1], third_layer_gate_nodes[i*2 + 0][2])
        tn.connect(second_layer_gate_nodes[i*5 + 4][1], third_layer_gate_nodes[i*2 + 1][2])

    # first_layer_nodes -> third_layer_nodes
    tn.connect(first_layer_gate_nodes[-2][5], third_layer_gate_nodes[-1][0])
    tn.connect(first_layer_gate_nodes[-5][5], third_layer_gate_nodes[-2][0])

    # in third_layer_nodes
    for i in range(bit_length - 1):
        tn.connect(third_layer_gate_nodes[i*2 + 0][4], third_layer_gate_nodes[i*2 + 1][1])

        if i != 0:
            tn.connect(third_layer_gate_nodes[i*2 + 0][5], third_layer_gate_nodes[i*2 - 2][0])
            tn.connect(third_layer_gate_nodes[i*2 + 1][5], third_layer_gate_nodes[i*2 - 1][0])
```

```Python
make_edges(initial_nodes, first_layer_gate_nodes, second_layer_gate_nodes, third_layer_gate_nodes)
make_edges(initial_nodes_r, first_layer_gate_nodes_r, second_layer_gate_nodes_r, third_layer_gate_nodes_r)
```

最後に、これらの二つのネットワークをつなげます。この際、測定したい部分の足は接続せずに残しておきます。

```Python
for i in range(bit_length):
    tn.connect(first_layer_gate_nodes[i*6 + 5][1], first_layer_gate_nodes_r[i*6 + 5][1])
    tn.connect(first_layer_gate_nodes[i*6 + 4][4], first_layer_gate_nodes_r[i*6 + 4][4])


for i in range(bit_length - 1):
    tn.connect(third_layer_gate_nodes[i*2 + 0][3], third_layer_gate_nodes_r[i*2 + 0][3])
    tn.connect(third_layer_gate_nodes[i*2 + 1][3], third_layer_gate_nodes_r[i*2 + 1][3])
    tn.connect(third_layer_gate_nodes[i*2 + 1][4], third_layer_gate_nodes_r[i*2 + 1][4])
```

以上でテンソルネットワークの構築は完了です。

それでは、出来上がったテンソルネットワークに対して縮約を取ります。今回は接続されていない足が残っているので、この部分をoutput_edge_orderに指定します。

```Python
nodes = (initial_nodes + first_layer_gate_nodes + second_layer_gate_nodes + third_layer_gate_nodes + 
         initial_nodes_r + first_layer_gate_nodes_r + second_layer_gate_nodes_r + third_layer_gate_nodes_r)

output_edges = []

output_edges.append(third_layer_gate_nodes[0][5])
output_edges.append(third_layer_gate_nodes[1][5])

output_edges.append(third_layer_gate_nodes_r[0][5])
output_edges.append(third_layer_gate_nodes_r[1][5])

result = tn.contractors.auto(nodes=nodes, output_edge_order=output_edges)
```

得られた計算結果を$4\times4$の行列に整形することで、縮約密度行列が得られます。

```Python
print(result.tensor.reshape(4,4))

# 左上から00,01,10,11に相当
# array([[0, 0, 0, 0],
#        [0, 0, 0, 0],
#        [0, 0, 1, 0],
#        [0, 0, 0, 0]])
```

動作の確認のため、別の入力でも試してみます。
ここでは、$\lvert a \rangle = \frac{1}{\sqrt{2}}(\lvert 000 \rangle+\lvert 100 \rangle)$、$\lvert b \rangle = \frac{1}{\sqrt{2}}(\lvert 000 \rangle+\lvert 001 \rangle)$として実行した場合の結果を確認してみます。これは、$\lvert a \rangle$の最上位ビットと、$\lvert b \rangle$の最下位ビットにアダマールゲートを適用した場合に相当します。

この場合、初期状態は以下のようになります。

```Python
initial_nodes = []

for i in range(n_qubits):
    if i == 0: # aの最上位ビット
        initial_nodes.append(tn.Node(np.array([1, 1]) / np.sqrt(2))) 
    elif i == n_qubits - 3: # bの最下位ビット
        initial_nodes.append(tn.Node(np.array([1, 1]) / np.sqrt(2))) 
    else:
        initial_nodes.append(tn.Node(np.array([1, 0])))
```

後は同じようにネットワークを定義して、縮約をとることで、以下の結果が得られます。

```Python
print(result.tensor.reshape(4,4))

# array([[0.25, 0.  , 0.  , 0.  ],
#        [0.  , 0.25, 0.  , 0.  ],
#        [0.  , 0.  , 0.5 , 0.  ],
#        [0.  , 0.  , 0.  , 0.  ]])
```

# 量子ビット数と計算時間の変化

作成したテンソルネットワークに対して、入力ビット長を変化させながら`tn.contractors.auto()`の部分の計算時間を確認してみます。量子ビット数は入力ビット長の5倍から1を引いたものになります。実行環境はMacbook Pro（CPU：2GHz クアッドコア Intel Core i5, メモリ：16GB）です。

|入力ビット長|量子ビット数|計算時間|
|----------|----------|-------|
| 10 | 49 | 0.09s |
| 100 | 499 | 1.82s |
| 200 | 999 | 8.76s |
| 300 | 1499 | 17.1s |
| 400 | 1999 | 28.3s |
| 500 | 2499 | 59.3s |
| 600 | 2999 | 89s |
| 700 | 3499 | 122s |
| 800 | 3999 | 191s |

![計算時間の変化](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/30ee8b7b-8caa-b5f1-6660-a2d18d8c4158.png)


量子ビット数が非常に多くても、現実的な時間で計算できている様子が確認できました。

# まとめ

今回は、特定の量子ビットのみを測定する必要のあるアルゴリズムの一例として、Quantum Bit String Comparatorを題材としたテンソルネットワークシミュレーションの実装を行いました。

今回の実験設定では、量子ビット数がかなり多くなっても問題なく計算できることがわかりました。入力する状態には制限を加えているので、量子もつれを含む状態を入力する場合に計算時間がどうなるのかは別途確認してみる必要がありそうです。

以上です。



# 参考文献

1. [東京大学「量子ソフトウェア」寄付講座](http://qsw.phys.s.u-tokyo.ac.jp/)

2. [東京大学「量子ソフトウェア」寄付講座 - 第3回量子ソフトウェアハンズオン(産学協働ゼミ)](http://qsw.phys.s.u-tokyo.ac.jp/joint202303)

3. [David Sena Oliveira, and Rubens Viana Ramos. "Quantum bit string comparator: circuits and applications." Quantum Computers and Computing, 7.1 (2007): 17-26.](https://www.researchgate.net/publication/228574906_Quantum_bit_string_comparator_Circuits_and_applications)