量子回路のテンソルネットワークシミュレーション 〜特定の状態の振幅の計算〜

はじめに

最近では、様々な量子コンピュータの実機が各種クラウドサービスから利用可能となっており、少しずつ一般のユーザでも量子コンピュータを利用しやすくなってきています。しかしながら、費用面の負担が大きく、気軽に使えるものではありません。また、現在の量子コンピュータはノイズも多く、計算結果にはエラーが含まれます。そのため、量子アルゴリズムの勉強や研究をする際には依然として古典コンピュータによるシミュレーションもよく用いられています。

本記事では、量子回路のシミュレーション手法として、近年注目を浴びているテンソルネットワークを用いたシミュレーション手法について、細かい理論的な説明は割愛して、簡単な実装例を紹介し、その性能を体験することを目的にします。

テンソルネットワークの活用方法は様々あるのですが、今回は最も簡単な、特定の状態の確率振幅を計算する手法について実装してみます。

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

テンソルネットワーク

テンソルネットワークは、大規模なテンソルを用いる計算の効率化を行う手法の一種で、計算時間やメモリ使用量の削減を目的とした技術です。
詳細な説明は本記事では割愛しますが、大規模なテンソルを小規模なテンソルの積に分解する、計算の順序を工夫する、適切な近似を導入する、などを行うことで、上記の目的を達成します。

量子コンピュータの文脈でテンソルネットワークが広く注目を浴びるようになった一つのきっかけとしては、2021年に発表された論文[3]だと考えています。その内容としては、2019年にGoogleが量子超越を達成したと主張したランダム量子回路サンプリングと呼ばれるタスクを、テンソルネットワークを駆使することで、スパコンで300秒ほどで実行した、というものです。

あらゆる量子回路に対して高速化が保証されるというわけではないのですが、適用先によっては、以前から使われている一般的な状態ベクトルシミュレーション手法に比べて、大幅な高速化・省メモリ化が可能です。

量子回路とテンソルネットワーク

まず、量子回路について簡単におさらいします。

一般に、量子回路は$n$個の量子ビットに対し、様々な種類のゲート操作を適用し、最終的に得られた量子状態に対して観測を繰り返すことで、出力を測定します。
このとき、古典コンピュータ上では量子状態を正規化された$2^n$次元のベクトルとして表現することができます。そして、各ゲート操作は、このベクトルに対して、ゲート操作に対応する行列をかけることに相当します。この計算を行うことで、量子コンピュータのシミュレーションが可能なのですが、$n$が増加するにつれ、状態ベクトルのサイズが指数的に大きくなり、この状態を保有することが難しくなるため、古典コンピュータでの計算は難しいとされています。

ここでテンソルネットワークの出番です。当然のことですが、ベクトルや行列というのもテンソルの一種です。また、実際の量子回路では、順番にゲート操作を作用させていくことになりますが、シミュレーションの上では、計算結果が変わらない範囲で、計算を分割したり、計算順序を工夫することができ、計算が効率的になることがあります。

実装例

ここからは、実際に簡単な量子回路を題材に、テンソルネットワークを組んでみて、Python上で実行してみます。
ここでは、Googleの開発したTensorNetwork[4]というライブラリを利用します。

■分析環境

Python 3.7.6
tensornetwork 0.4.6

ベル状態の回路

最初に、ベル状態を作成する回路を作ってみます。この回路は、$\lvert 00 \rangle$と$\lvert 11 \rangle $の二つの状態が50％ずつで出現します。

まずは$\lvert 00 \rangle $の振幅を計算するテンソルネットワークを作成してみます。

図として表現すると、以下のような形状になります。

今回は初期状態と振幅を求めたい状態が等しいので、少し分かりづらいですが、一番左側の二つのノードが初期状態、一番右側の二つのノードが振幅を計算したい状態に対応しています。すなわち、通常の量子回路の終端部分に、振幅を計算したい状態を接続した形になっています。

ノード間をつなぐエッジの番号は、理解の手助けのため、後ほど紹介するコード内のコメントに対応づけて振ってあります。

少しだけ数式による補足をしておくと、回路の出力を$\lvert \psi \rangle = U \lvert 00 \rangle$、振幅を求めたい状態を$\lvert 00 \rangle $とすると、このテンソルネットワークの縮約を計算することで、$ \langle 00 \lvert \psi \rangle$が得られ、回路の出力における$\lvert 00 \rangle $の振幅を求めることができます。

次に、これを実装していきます。

事前準備として、各ゲート操作に対応する行列を、numpyのarrayで用意しておきます。

import numpy as np
import tensornetwork as tn

h_gate = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
cx_gate = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])

次に、図中のノードに相当するものを一つずつTensorNetworkのNodeクラスを用いて定義していき、リストに格納します。

# 初期状態に相当するノードを格納するリストを作成
initial_state_nodes = list()

for i in range(2):
    initial_state_nodes.append(tn.Node(np.array([1, 0]))) # 0状態 = [1, 0]を入力


# 各ゲートに相当するノードを格納するリストを作成
gate_nodes = list()

gate_nodes.append(tn.Node(h_gate)) # アダマールゲート
gate_nodes.append(tn.Node(cx_gate.reshape(2,2,2,2))) # CNOTゲート、足の数に合わせてreshapeする


# 最終状態に相当するノードを格納するリストを作成
final_state_nodes = list()

for i in range(2):
    final_state_nodes.append(tn.Node(np.array([1, 0]))) # 0状態に対応するベクトル[1, 0]を入力

続いて、各ノードをつなぐエッジを定義します。足の番号の順番は、下図のように、左上を0として、左上、左下、右上、右下の順番に振られていきます。

エッジをつなぐ操作は、tn.connect(a, b)という形で記述できます。入力は、つなぎたいノードの足になります。

簡易的な記述方法として、a ^ bという書き方もできます。

# 0 : 1量子ビット目の初期状態から、Hゲートの入力
tn.connect(initial_state_nodes[0][0], gate_nodes[0][0])

# 1 : Hゲートの出力から、CXゲートの入力の1量子ビット目
tn.connect(gate_nodes[0][1], gate_nodes[1][0])

# 2 : 2量子ビット目の初期状態から、CXゲートの入力の2量子ビット目
tn.connect(initial_state_nodes[1][0], gate_nodes[1][1])

# 3 : CXゲートの出力の1量子ビット目から、1量子ビット目の最終状態
tn.connect(gate_nodes[1][2], final_state_nodes[0][0])

# 4 : CXゲートの出力の2量子ビット目から、2量子ビット目の最終状態
tn.connect(gate_nodes[1][3], final_state_nodes[1][0])

最後に、作成されたネットワークに対して縮約を計算します。テンソルネットワークの計算においては、縮約を取る順番が計算速度に大きく影響するのですが、ここではライブラリの機能に任せて自動で決めてもらいます。

以下のコードのautoの部分を書き換えることで、縮約を取る順番を決めるアルゴリズムを自分で指定することもできます。

# すべてのノードをまとめたリストを作成
all_nodes = initial_state_nodes + gate_nodes + final_state_nodes

# ノードのリストを引数に、contractorsクラスから適当な関数を呼ぶ
result = tn.contractors.auto(nodes=all_nodes)

得られた結果から、以下のようにして振幅を取り出せます。

print(result.tensor)
# 0.7071067811865475

print(result.tensor ** 2)
# 0.4999999999999999

別の状態の確認

上の結果を見ると、うまくいっているように見えます。念の為、他の状態の振幅も正しく表示されることを確認します。これは、最終状態の部分のコードを修正することで、計算できます。

一度定義したノードには、すでにエッジが貼られている状態になっているので、ノードとエッジの定義は改めてすべて実行する必要がありますが、すべてを書くと冗長になるため以下では変更点と結果だけ記載します。

$\lvert 11 \rangle$の場合

# 最終状態に相当するノードを格納するリストを作成
final_state_nodes = list()

for i in range(2):
    final_state_nodes.append(tn.Node(np.array([0, 1]))) # 1状態 = [0, 1]を入力

print(result.tensor)
# 0.7071067811865475

print(result.tensor ** 2)
# 0.4999999999999999

$\lvert 01 \rangle$の場合

# 最終状態に相当するノードを格納するリストを作成
final_state_nodes = list()

final_state_nodes.append(tn.Node(np.array([1, 0]))) # 0状態 = [1, 0]を入力
final_state_nodes.append(tn.Node(np.array([0, 1]))) # 1状態 = [0, 1]を入力

print(result.tensor)
# 0.0

print(result.tensor ** 2)
# 0.0

より大規模な回路

ベル状態を作成する回路での実験を通して、うまく計算できていそうなことがわかりました。それではさらに大きな回路として、$n$量子ビットのもつれ状態を作成する回路を作ってみます。

回路の全体像としては以下のような形になります。

それでは、これをテンソルネットワークとして、実装していきます。先ほどの例と同様に、$\lvert 00....0 \rangle $が出現する確率が50％になることを確認してみましょう。

まずはノードの定義を行います。

# 量子回路のビット数
n_qubits = 100

# ノードの定義
# 初期状態に相当するノードを格納するリストを作成
initial_state_nodes = list()

for i in range(n_qubits):
    initial_state_nodes.append(tn.Node(np.array([1, 0]))) # 0状態 = [1, 0]を入力

# 各ゲートに相当するノードを格納するリストを作成
gate_nodes = list()

gate_nodes.append(tn.Node(h_gate))

for i in range(n_qubits - 1):
    gate_nodes.append(tn.Node(cx_gate.reshape(2,2,2,2)))

# 添字0：Hゲート
# 添字1〜n_qubits-1：CNOTゲート


# 最終状態に相当するノードを格納するリストを作成
final_state_nodes = list()

for i in range(n_qubits):
    final_state_nodes.append(tn.Node(np.array([1, 0]))) # 0状態 = [1, 0]を入力

次に、これらをつなぐエッジを定義します。

# 1量子ビット目の初期状態から、Hゲートの入力
tn.connect(initial_state_nodes[0][0], gate_nodes[0][0])

# Hゲートの出力から、1個目のCXゲートの入力の1量子ビット目
tn.connect(gate_nodes[0][1], gate_nodes[1][0])

# i個目のCXゲートの出力の1量子ビット目から、i+1個目のCXゲートの入力の1量子ビット目
for i in range(1, n_qubits - 1):
    tn.connect(gate_nodes[i][2], gate_nodes[i+1][0])
    
# i個目の量子ビットの初期状態から、i-1個目のCXゲートの入力の2量子ビット目
for i in range(1, n_qubits):
    tn.connect(initial_state_nodes[i][0], gate_nodes[i][1])
    
# 最後のCXゲートの出力の1量子ビット目から、最終状態の1量子ビット目
tn.connect(gate_nodes[-1][2], final_state_nodes[0][0])

# 各CXゲートの出力の2量子ビット目から、各最終状態
for i in range(1, n_qubits):
    tn.connect(gate_nodes[i][3], final_state_nodes[i][0])

最後に、完成したネットワークで縮約を取ると、計算結果が得られます。

all_nodes = initial_state_nodes + gate_nodes + final_state_nodes

result = tn.contractors.auto(nodes=all_nodes)

print(result.tensor)
# 0.7071067811865475

print(result.tensor**2)
# 0.4999999999999999

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

作成したテンソルネットワークに対して、量子ビット数を変化させながら計算時間を観察してみます。ここでは、私の手元のMacbook Pro（CPU：2GHz クアッドコア Intel Core i5, メモリ：16GB）で計算した結果を記載しています。

ネットワークの定義部分はほとんど時間はかかりませんので、縮約計算部分の時間を以下に整理しています。

量子ビット数	計算時間
10	5.96ms
100	52.4ms
1000	1.26s
2000	7.69s
3000	14.4s
4000	24.2s
5000	35.8s

計算時間の伸びが緩やかで、非常に大きな回路でも現実的な時間でシミュレーションできている様子が確認できると思います。

今回は、非常に簡単な回路で実験をおこなっており、対象とする量子回路によって、計算時間の伸び方には違いが生じると思われますので、この点には注意が必要です。

まとめ

今回は簡単な例として、単独の状態の振幅を確認するテンソルネットワークシミュレーションについて、実装例を紹介しました。

ただし、今回のような、特定の状態の振幅のみを求めたいというケースは、実用上はあまりないかもしれません。

この他にも、テンソルネットワークを用いて、サンプリングを行ったり、量子ビットの一部のみを測定するなど、様々な使い方が可能です。

次回は、量子ビットの一部に着目し、その部分の密度行列を求めることで、量子ビットの一部のみを測定するような内容についてまとめようと思います。

参考文献

1. 東京大学「量子ソフトウェア」寄付講座

2. 東京大学「量子ソフトウェア」寄付講座 - 第2回量子ソフトウェア産学協働ゼミ

3. Yong Liu, et al. "Closing the "quantum supremacy" gap: achieving real-time simulation of a random quantum circuit using a new sunway supercomputer." Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. 2021.

4. GitHub - google/TensorNetwork: A library for easy and efficient manipulation of tensor networks.