前回は

https://blueqat.com/yuichiro_minato2/270025ec-8478-45a3-992f-1385c9b3e1e4

今回は4回目で、実際に前回までの量子コンピュータの計算をより効率的にツールなどを使って行います。

テンソルネットワーク対応の新しいツールがどんどん出ています。一番有名なのはcuQuantumですが、今回はquimbを使ってみたいと思います。こちらは最初から量子回路シミュレータが搭載されていて前回みたいにゲートを決めなくても自動でやってくれます。

https://quimb.readthedocs.io/en/latest/tensor-circuit.html

!pip install -q quimb cotengra scikit-optimize

今回はまず、80量子ビットの量子もつれ計算をします。

%config InlineBackend.figure_formats = ['svg']

import random
import quimb as qu
import quimb.tensor as qtn

N = 80
circ = qtn.Circuit(N)

randomly permute the order of qubits

regs = list(range(N))
random.shuffle(regs)

準備は済みました。

次は最初の量子ビットに0を適用し、順番にCXゲートを適用してみます。

最後にサンプルを100個取ります。

# hamadard on one of the qubits
circ.apply_gate('H', regs[0])

chain of cnots to generate GHZ-state

for i in range(N - 1):
circ.apply_gate('CNOT', regs[i], regs[i + 1])

from collections import Counter

sample it 100 times, count results:

Counter(circ.sample(100))

結果は、

Counter({'00000000000000000000000000000000000000000000000000000000000000000000000000000000': 56, '11111111111111111111111111111111111111111111111111111111111111111111111111111111': 44})

きちんと取れました。次に量子回路を確認してみます。

circ.psi.draw(color=['PSI0', 'H', 'CNOT'])

できました。続いてQAOAをやってみます。こちらも例題があります。

https://quimb.readthedocs.io/en/latest/examples/ex_tn_qaoa_energy_bayesopt.html

問題設定は、

54-qubit, p=4(depth 12), random 3-regular graphのmaxcutをQAOAで解きます。

p=4なので結構長いですね。量子ビットも50超えてます。

%config InlineBackend.figure_formats = ['svg']

import quimb as qu
import quimb.tensor as qtn

これまでは計算順序は適当に決めてましたが、これを最適化するツールcotengraを導入します。

import cotengra as ctg

opt = ctg.ReusableHyperOptimizer(
methods=['greedy'],
reconf_opts={},
max_repeats=32,
max_time="rate:1e6",
parallel=True,
# use the following for persistently cached paths
# directory=True,
)

次に問題設定です。グラフ問題としてnetworkxで設定します。

import networkx as nx

reg = 3
n = 54
seed = 666
G = nx.random_regular_graph(reg, n, seed=seed)

terms = {(i, j): 1 for i, j in G.edges}

次にQAOAのパラメータを設定します。ステップ数p=4。変分パラメータgamma/betaを設定します。

p = 4
gammas = qu.randn(p)
betas = qu.randn(p)
circ_ex = qtn.circ_qaoa(terms, p, gammas, betas)

簡単に量子回路も描画できます。

circ_ex.psi.draw(color=['PSI0', 'H', 'RZZ', 'RX'])

計算が重たくなりすぎないか、全てのローカルハミルトニアンで、リハーサルをするのがおすすめのようです。これで見通しをつける。結構時間かかりますね。

import tqdm

ZZ = qu.pauli('Z') & qu.pauli('Z')

local_exp_rehs = [
circ_ex.local_expectation_rehearse(weight * ZZ, edge, optimize=opt)
for edge, weight in tqdm.tqdm(list(terms.items()))
]

100%|██████████| 81/81 [27:16<00:00, 20.21s/it]

27分近くかかりました。

import matplotlib.pyplot as plt

with plt.style.context(qu.NEUTRAL_STYLE):
fig, ax1 = plt.subplots()
ax1.plot([rehs['W'] for rehs in local_exp_rehs], color='green')
ax1.set_ylabel('contraction width, W, [log2]', color='green')
ax1.tick_params(axis='y', labelcolor='green')

ax2 = ax1.twinx()
ax2.plot(\[rehs\['C'\] for rehs in local\_exp\_rehs\], color='orange')
ax2.set\_ylabel('contraction cost, $C$, \[log10\]', color='orange')
ax2.tick\_params(axis='y', labelcolor='orange')

これくらいのコストだといいようです。見極めがよくわからない。次にコストの期待値の計算ですね。jaxバックエンドで計算をします。

def energy(x):
p = len(x) // 2
gammas = x[:p]
betas = x[p:]
circ = qtn.circ_qaoa(terms, p, gammas, betas)

ZZ = qu.pauli('Z') & qu.pauli('Z')

ens = \[
circ.local\_expectation(weight \* ZZ, edge, optimize=opt, backend="jax")
for edge, weight in terms.items()
\]
return sum(ens).real

求める期待値はscikit-learn optimizeで最適化します。ベイズ最適化です。２時間くらいかかりました。

from skopt import Optimizer
from skopt.plots import plot_convergence, plot_objective

eps = 1e-6
bounds = (
[(0.0 + eps, qu.pi / 2 - eps)] * p +
[(-qu.pi / 4 + eps, qu.pi / 4 - eps)] * p
)

bopt = Optimizer(bounds)

for i in tqdm.trange(100):
x = bopt.ask()
res = bopt.tell(x, energy(x))

結果の表示です。

with plt.style.context(qu.NEUTRAL_STYLE):
plot_convergence(res);

エネルギーランドスケープも取れるそうです。かっこいい。

with plt.style.context(qu.NEUTRAL_STYLE):
plot_objective(
res,
cmap='RdYlBu_r',
dimensions=[f'\\\\gamma\_{i}' for i in range(p)] + [f'\\\\beta\_{i}' for i in range(p)],
);

今後の量子コンピュータの計算はこうなっていくと思います。以上です。