この辺りは技術発展も速く、ドキュメントもないと思うので、最新版で使い方を見ていきます。テンソルネットワークでは、目的意識をもって計算をするのが大事でした。今回は状態ベクトルからサンプリングまでを実装してみたいと思います。Google Cirqは量子コンピュータSDKでもちょっと使い方が特殊なところもあるので、blueqatSDKなどで簡単に計算できる方法をサポートすべきですが、手間の問題で良しなにやってみます。
現在Google ColabはPythonのバージョンが3.8なので、簡単にcuQuantumが入るようです。こちらのブログ記事では状態ベクトルシミュレータの使い方を解説しています。
https://blueqat.com/yuichiro_minato2/09eb8212-dfd6-4330-9032-e9c1cd4667d8
状態ベクトルシミュレータは状態ベクトルと呼ばれる、量子状態をベクトルの形で取り出せるので、そこからさまざまな計算に派生できます。一方今回利用するテンソルネットワークは計算順序が異なり、最初に行うべき計算に向けて量子回路を構築します。同じ量子回路でもサンプルを取るのと、期待値を計算するのと、状態ベクトルを取る（量子ビット数が小さいとき限定）では最終的に計算する回路のサイズが異なります。
最初にインストールしてcirq経由でやってみましょう。cuQuantumのgithubリポジトリにサンプルがあるので、そのあたりを良しなにまとめてみました。
!pip install cuquantum-python cirq
cuquantum + cirq だけでしたらcuquantum-pythonとcirqだけ入れれば十分です。次にツールを読み込みますが、cuquantumからの読み込みはテンソルネットワークで利用する縮約と量子回路をEinsumに変換するものです。
import cirq
import cupy as cp
import numpy as np
from cuquantum import contract, CircuitToEinsum
Google cirqは恵まれていて、Einsum形式に変換してもらえてますね。Qiskitも最近対応しています。手順は単純です。今回は状態ベクトルを得たいので、量子回路を作ってそれを状態ベクトルタイプに変換を咬ませます。
#量子ビット数を決めます
n_qubits = 3
#回路を初期化
circuit = cirq.Circuit()
#量子ビットの並びを決める必要があり、今回は一列に配置
qubits = cirq.LineQubit.range(n_qubits)
#量子ゲートを適用
circuit.append(cirq.H(qubits[0]))
circuit.append(cirq.H(qubits[1]))
circuit.append(cirq.H(qubits[2]))
#回路を確認
print(circuit)
#量子回路をアインシュタインの縮約記法によりテンソルネットワークに変換
converter = CircuitToEinsum(circuit, dtype='complex128', backend=cp)
回路を出力しましたのでアダマールゲートが並んでいます。
0: ───H───
1: ───H───
2: ───H───
テンソルネットワークに変換されていますので、あとは計算するだけです。
expression, operands = converter.state_vector()
sv = contract(expression, *operands)
sv
見てみると2**3の状態ベクトル形式で取り出せていますね。
array([[[0.35355339+0.j, 0.35355339+0.j],
[0.35355339+0.j, 0.35355339+0.j]],
\[\[0.35355339+0.j, 0.35355339+0.j\],
\[0.35355339+0.j, 0.35355339+0.j\]\]\])
興味としてexpressionの中身を見てみます。
expression
出力は、
a,b,c,da,eb,fc->def
これを見る限りは、abcが量子ビットの初期状態の|0>の腕を表していて、daが0番目のHゲート、ebが1番目のHゲート、fcが2番目のHゲートとなっていて、最終的にテンソルからdefの腕三本が出ています。
最終的にはsvを出力すると、状態ベクトルが出力されます。次に単一振幅を求めてみます。これ単体ではあまり利用する機会はないかもですが、その他の重要な計算に派生するので一応見ておきます。なんと便利機能で求めたいビット列と量子回路を組み合わせてテンソルネットワーク回路を用意してくれました。
#振幅を求めたいビット列を指定
bitstring = '010'
#便利機能。。。
expression, operands = converter.amplitude(bitstring)
print(f'einsum expression:\n{expression}\n')
出力結果を見てみましょう。
einsum expression:
a,b,c,da,eb,fc,d,e,f->
これをみると、さっきの腕がdefだったところにd,e,fがそれそれ010のビット列がつながってスカラー量が求まるのが確認できます。実際に計算ですが、
#振幅と出現確率を求めてますね
amplitude = contract(expression, *operands)
probability = abs(amplitude) ** 2
print(f'for bitstring {bitstring}, amplitude: {amplitude}, probability: {probability}\n')
単一振幅の値が求まり、その二乗で出現確率を求めています。
for bitstring 010, amplitude: (0.353553390593274+0j), probability: 0.12500000000000017
次にもチョイ役立ちそうなハミルトニアンの期待値を求めてみます。こちらはおなじみですね。NVIDIAさんが準備してくれています。pauli演算子の文字列で準備します。基本的にハミルトニアンの単項が出るので、それを求めます。
pauli_string = 'ZII'
expression, operands = converter.expectation(pauli_string, lightcone=True)
expec = contract(expression, *operands)
print(f'expectation value for {pauli_string}: {expec}')
出力結果は、
expectation value for ZII: (-3.322727519965453e-18+0j)
数値が出ました。Hゲートを選んでしまったのでハミルトニアンZの計算結果はちょっと微妙です。試しにXにしてみます。
pauli_string = 'XII'
expression, operands = converter.expectation(pauli_string, lightcone=True)
expec = contract(expression, *operands)
print(f'expectation value for {pauli_string}: {expec}')
expectation value for XII: (1.0000000000000004+0j)
まぁ、あっていそうです。ローカルハミルトニアンの期待値の計算は比較的簡単です。最後にサンプリングと行きたいところですが、TNではおそらくサンプリングが一番面倒です。こちらのサンプルでは密度行列から求めています。今回も密度行列を求めるTN操作が用意されていたので、それを使って周辺確率（？）を求めます。この辺りはあんまり詳しくないので今後見てみたいと思います。
qubits = sorted(circuit.all_qubits())
where = (qubits[0], qubits[1], qubits[2])
fixed = {}
expression, operands = converter.reduced_density_matrix(where, fixed=fixed)
rdm = contract(expression, *operands)
n_sites = len(where)
sh = 2 ** n_sites
prob = abs(rdm.reshape(sh, sh).diagonal())**2
marginal = prob.reshape((2,)*n_sites) / prob.sum() # normalization
# print the probabilities
print(marginal)
求まりました。量子ビット数が増えるとどうなってしまうのでしょうか。。。
[[[0.125 0.125]
[0.125 0.125]]
[[0.125 0.125]
[0.125 0.125]]]
とりあえずサンプルを取ってみます。
import itertools
import matplotlib.pyplot as plt
def create_samples(marginal, nsample):
# this function generates samples on CPU, since we are sampling over Python strings
nqubits = marginal.ndim
bitstrings = [''.join(bitstring) for bitstring in itertools.product('01', repeat=nqubits)]
ps = cp.asnumpy(marginal.flatten())
return np.random.choice(bitstrings, size=nsample, p=ps)
nsample = 10000
samples = create_samples(marginal, nsample)
keys, counts = np.unique(samples, return_counts=True)
print(f'\n{nsample} samples on {where} assuming:')
# plot the result; note the keys are already sorted by numpy.unique()
plt.bar(keys, counts)
plt.show()
なんか求まりましたね。
ちょっとまとめなおしただけですが、意外と今後この辺りは需要がありそうなので、勉強会とブログ記事を増やしてみます。