Go to Qiita Advent Calendar 2024 Top
LoginSignup
47
26

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

量子コンピューターAdvent Calendar 2019

Day 4
@j-i-k-o(Kohji Nishimura)inJij Inc.

IBM-Q の53qubit量子コンピュータのシミュレータを作ってみた

Last updated at Posted at 2019-12-03

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

はじめに

こんにちは、@j-i-k-o と申します。
この記事は量子コンピューター Advent Calendar 2019の四日目の記事です。

普段は量子アニーリングに関する研究だったり、株式会社JijにてOpenJijの開発を行っていたりしますが、気分転換に量子ゲートのほうにも触れてみました。
テンソルネットワークと呼ばれる手法を用いて量子コンピュータのシミュレータを作ってみた話をします。

背景

最近Googleが量子超越を達成したというニュースが話題ですが、IBMもオンラインから利用可能な53qubitの量子コンピュータを発表しました。(参考リンク: IBM、53量子ビットの新型量子コンピュータを発表)

この量子コンピュータの構造は下の図で表され、線でつながっている量子ビット間にゲートを作用させることができます。
(出典: https://www.ibm.com/blogs/research/2019/09/quantum-computation-center)

53Q_topology.jpg

例えば、9番目の量子ビットに関して、隣接している5, 8, 10番目の量子ビット間にSWAPゲートやCNOTゲート等を作用させることができます。

ただ残念ながら今の所私は53qubitの量子コンピュータを使えないので、以下に書くテンソルネットワークを利用した手法でシミュレータを作ります。

テンソルネットワーク

古典コンピュータで量子計算が非常に難しい理由の一つとして、（真面目に計算しようとすると）必要となる量子状態を扱うのにサイズに応じて指数的に増大する要素を持ったベクトルを取り扱わなければならないことが挙げられます。例えば量子ビット数5の場合、波動関数は

\ket{\psi} = \sum_{\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5}\Psi^{{\sigma_1 \sigma_2 \sigma_3 \sigma_4 \sigma_5}}\ket{\sigma_1 \sigma_2 \sigma_3 \sigma_4 \sigma_5}

のように書けますが、全ての0と1の組合せを考慮すると$\Psi$は32個の要素を持つベクトルとなります。5量子ビットであればこの程度ですが量子ビット数が例えば50ぐらいになると$\Psi$の要素数は1125兆8999億684万2624となり、スパコンで取り扱うのも難しくなってくるメモリサイズになります。

しかし、仮に初期状態$\ket{0000\cdots}$から始めるとして、2量子間ゲートを作用させない場合、例えば
download.png
のように全部1量子ビットだけに作用するゲートの場合、異なる量子ビットに相関が発生し得ない (エンタングルしない)ので、波動関数はもっとシンプルに

\ket{\psi} = \sum_{\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5}c^{\sigma_1}\ket{\sigma_1}c^{\sigma_2}\ket{\sigma_2}c^{\sigma_3}\ket{\sigma_3}c^{\sigma_4}\ket{\sigma_4}c^{\sigma_5}\ket{\sigma_5}

のようにお互いに相関がない形で書いてもよいはずです。このときの波動関数に存在するパラメータの数は、各$c^{\sigma}$が2要素のベクトルであるので合計10個のみとなります。量子ビット数が50のときでもパラメータ数は100しかなく、上で書いた1125兆ナントカとはえらい違いです。

実際の計算には当然2量子間ゲート (CNOTやSWAPなど)も存在するためここまで話は上手くいかないですが、量子ビット間の相関も取り入れるために$c^{\sigma}$にもっと自由度を追加して、例えば

\ket{\psi} = \sum_{a_1=1}^{\chi}\sum_{a_2=1}^{\chi}\sum_{a_3=1}^{\chi}\sum_{a_4=1}^{\chi}\sum_{a_5=1}^{\chi}\sum_{\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5}A^{\sigma_1}_{a_1}A^{\sigma_2}_{a_1 a_2}A^{\sigma_3}_{a_2 a_3}A^{\sigma_4}_{a_3 a_4}A^{\sigma_5}_{a_4}\ket{\sigma_1 \sigma_2 \sigma_3 \sigma_4 \sigma_5}

のように$a_1, a_2$のような自由度を追加して$c$の代わりにテンソル$A$を導入し、$\chi$の大きさを調節することにより扱える相関の度合いをコントロールできます。相関が発生していない場合は$\chi=1$に相当しますし、$\chi=2$にすれば例えば2量子ビット間に相関がある

\ket{\psi} = \frac{1}{\sqrt{2}}(\ket{00} + \ket{11})

みたいな量子状態も扱えるようになります。$\chi$を大きくすればするほど波動関数の表現力も増えていきます。$N$量子ビットの波動関数を表すのに必要なパラメータ数が$2N\chi^2$程度で済むので、元の波動関数をおおよそ再現できる程度の$\chi$の値を選んでおけば圧倒的に効率よく古典コンピュータで波動関数を扱うことができるというわけです。

式だと少々見づらいので、この状態を図で表すと、
2019-12-02-222301_526x226_scrot.png
こんな具合になります。各四角がテンソルを表していて、丸が$\sigma_1, \sigma_2, \cdots$のような量子ビット、テンソル間の直線が$a_1, a_2$などの自由度に相当します。

このようにテンソル間の相関を用いて波動関数を書き表したものをテンソルネットワーク と呼びます。1
歴史的にはこの手法で1次元のそこそこ大きいサイズの量子系の基底状態を多項式オーダーで計算できる密度行列繰り込み群 (DMRG) 2の開発から始まり、今ではユニタリ時間発展を扱うTEBD (Time-evolving block decimation)などテンソルネットワークを用いた様々な手法が開発されてきました。

より詳しく知りたい方はこの文献 3がおすすめです。かなり読み応えがありますが、どう実装すればよいかなどが全て記述されています。

シミュレータ

いよいよシミュレータを作ってみましょう。github上のコードはこちらにあります。
また、テンソルを内部で扱う際にITensorというライブラリを用いているので事前にダウンロードが必要です。
わりと短時間で作ったものなのでREADMEにある以上のドキュメントは作っていません...

テンソルネットワーク状態の生成

まず初めにテンソルネットワーク状態を生成するのですが、IBM Qの構造に合わせて相関を導入してあげればよさそうなので、量子コンピュータの構造の図の量子ビットが位置する場所にテンソルを入れていき、一部だけ図にすると
2019-12-02-234302_765x514_scrot.png

このようなネットワークを持った波動関数を作ってあげればよさそうです。 ibmq.hに一つ一つネットワークを導入している様子が書いてあります。45

//generate tensors
generate_link(1,2);
generate_link(2,3);
generate_link(3,4);
generate_link(4,5);

generate_link(1,6);
generate_link(5,7);
generate_link(6,8);
generate_link(7,12);

generate_link(8,9);
generate_link(9,10);
generate_link(10,11);
generate_link(11,12);

generate_link(8,13);
generate_link(12,14);
generate_link(13,15);
generate_link(14,16);
generate_link(15,17);
generate_link(16,19);

generate_link(10,18);

generate_link(17,20);
generate_link(19,21);
generate_link(20,22);
generate_link(21,23);
generate_link(22,24);
generate_link(23,28);

generate_link(18,26);

generate_link(24,25);
generate_link(25,26);
generate_link(26,27);
generate_link(27,28);

generate_link(24,29);
generate_link(28,30);
generate_link(29,31);
generate_link(30,35);

generate_link(31,32);
generate_link(32,33);
generate_link(33,34);
generate_link(34,35);

generate_link(31,36);
generate_link(35,37);
generate_link(36,38);
generate_link(37,39);
generate_link(38,40);
generate_link(39,42);

generate_link(33,41);

generate_link(40,43);
generate_link(42,44);
generate_link(43,45);
generate_link(44,46);
generate_link(45,47);
generate_link(46,51);

generate_link(41,49);

generate_link(47,48);
generate_link(48,49);
generate_link(49,50);
generate_link(50,51);

generate_link(47,52);
generate_link(51,53);

ここで、各量子ビットのつながりを図示すると次のようになります。シミュレータ上のインデックスと上の画像のIBM Qの量子コンピュータのインデックスとは対応してないので注意してください。
2019-12-04-073441_555x629_scrot.png

また、2量子ビットゲートの計算で後で使うため、1番目と2番目の量子ビットに対応するテンソルを掛け算しておき、

\sum_{a=1}^{\chi} A^{\sigma_1}_{ae} A^{\sigma_2}_{ab} = \Psi^{\sigma_1 \sigma_2}_{eb}

で置き換えておきます。これで、波動関数の状態は
2019-12-02-234822_759x514_scrot.png
こんな具合になっています。式で書くとかなりややこしいので、図のありがたみがよく分かります。

量子ゲートの掛け算

今回はテンソルネットワーク状態の時間発展で使われる、参考文献 3の7.3.1節にある**tDMRG (time-dependent DMRG)**アルゴリズムを使います。また、1量子ビットに作用するゲート (Xゲートとか)も2量子ビットに作用するゲートを使って$X \otimes I$のようにかけるため、このアルゴリズムでは2量子ビットに作用するゲートだけ扱います。

基本的なアルゴリズムの流れとしては、量子ゲートを作用させたい量子ビットまで$\Psi$を"移動"させて、量子ゲートを掛け算する操作を繰り返すことになります。この"移動"について、以下、$\sigma_2, \sigma_3$に2量子ビット間ゲート$U^{(\sigma_2' \sigma_3'),(\sigma_2 \sigma_3)}$を作用させる例を用いて説明します。

まず、$\Psi$を特異値分解することにより、

\Psi^{\sigma_1 \sigma_2}_{eb} = \sum_{a}U^{\sigma_1}_{ea}S_{a a}(V^{\dagger})^{\sigma_2}_{ab}

という記述を得ます。ここで、$S$は対角行列で成分は全て非負となります。$U$や$V^{\dagger}$はユニタリ行列です。

次に、特異値分解の性質上、$S_{aa}$は$a=0$に近い成分が一番値が大きく、そこから$S_{00} \geq S_{11} \geq \cdots$と0に向かって減少していくため、$\Psi$の中の和について、ある程度大きい$a$については打ち切ってしまってもよさそうです。そこで、打ち切り誤差となる

\epsilon = 1 - \frac{\sum_{a=1}^{\chi}S_{aa}}{\sum_{a}S_{aa}}

が一定量以下となる$\chi$を選び、

\Psi^{\sigma_1 \sigma_2}_{eb} \simeq \sum_{a=1}^{\chi}U^{\sigma_1}_{ea}S_{a a}(V^{\dagger})^{\sigma_2}_{ab}

と近似します。上で述べた「元の波動関数をおおよそ再現できる程度の$\chi$の値を選んでおけば圧倒的に効率よく古典コンピュータで波動関数を扱うことができる」のトリックはここにあります。

そして、一番目の量子ビットに対応するテンソル$A^{\sigma_1}$を

A^{\sigma_1}_{ea} = U^{\sigma_1}_{ea}

で定義しなおし、$\Psi$を

\Psi^{\sigma_2 \sigma_3}_{ac} = \sum_{b}S_{a a}(V^{\dagger})^{\sigma_2}_{ab}A^{\sigma_3}_{bc}

で定義し直せば、$\Psi$が$\sigma_2, \sigma_3$にまたがる演算子となったため、"移動"に成功したと言えます。
このときの波動関数を図に表すと、
2019-12-03-004937_762x506_scrot.png
となっています。ここまで特異値分解の誤差を除き波動関数の情報は何も変わっていないです。
後はユニタリ演算を作用させて、

\Psi^{\sigma_2' \sigma_3'}_{ac} = \sum_{\sigma_2,\sigma_3}U^{(\sigma_2' \sigma_3'),(\sigma_2 \sigma_3)} \Psi^{\sigma_2 \sigma_3}_{ac}

とすればゲート演算ができます。この処理をひたすら繰り返すことで量子ゲートをどんどん作用させることができます。

ibmq.hにあるshift_toという関数に上の処理が記述されています。隣接する量子ビットに$\Psi$を"移動"させる処理です。
また、applyという関数で量子ゲートを作用させています。

結果の取得

この方法で生成される波動関数は直接テンソルを見ても何がどうなっているのかサッパリ分からないと思うので、$\braket{\psi_1}{\psi_2}$のような値 (overlap)を計算する必要があります。6 この量を用いることで2つの波動関数の重なり具合を調べることができます。

IBM Qの量子コンピュータの構造で説明しようとすると複雑になりすぎるので、最初に紹介した
2019-12-02-222301_526x226_scrot.png
このテンソルネットワークで説明します。

\ket{\psi_1} = \sum_{a_1=1}^{\chi}\sum_{a_2=1}^{\chi}\sum_{a_3=1}^{\chi}\sum_{a_4=1}^{\chi}\sum_{a_5=1}^{\chi}\sum_{\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5}A^{\sigma_1}_{a_1}A^{\sigma_2}_{a_1 a_2}A^{\sigma_3}_{a_2 a_3}A^{\sigma_4}_{a_3 a_4}A^{\sigma_5}_{a_4}\ket{\sigma_1 \sigma_2 \sigma_3 \sigma_4 \sigma_5} \\
\ket{\psi_2} = \sum_{b_1=1}^{\chi}\sum_{b_2=1}^{\chi}\sum_{b_3=1}^{\chi}\sum_{b_4=1}^{\chi}\sum_{b_5=1}^{\chi}\sum_{\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5}B^{\sigma_1}_{b_1}B^{\sigma_2}_{b_1 b_2}B^{\sigma_3}_{b_2 b_3}B^{\sigma_4}_{b_3 b_4}B^{\sigma_5}_{b_4}\ket{\sigma_1 \sigma_2 \sigma_3 \sigma_4 \sigma_5} \\

これら2つの波動関数でoverlapとなる$\braket{\psi_1}{\psi_2}$を計算すると、

\braket{\psi_1}{\psi_2} = \sum_{a_4,b_4,\sigma_5}\left(\sum_{a_3,b_3,\sigma_4}\left(\sum_{a_2,b_2,\sigma_3}\left(\sum_{a_1,b_1,\sigma_2}\left(\sum_{\sigma_1}(A^{\dagger})^{\sigma_1}_{a_1}B^{\sigma_1}_{b_1}\right)(A^{\dagger})^{\sigma_2}_{a_1 a_2}B^{\sigma_2}_{b_1 b_2}\right)(A^{\dagger})^{\sigma_3}_{a_2 a_3}B^{\sigma_3}_{b_2 b_3}\right)(A^{\dagger})^{\sigma_4}_{a_3 a_4}B^{\sigma_4}_{b_3 b_4}\right)(A^{\dagger})^{\sigma_5}_{a_4}B^{\sigma_5}_{b_4}

と書けることが分かります。
内側の和から順番に取っていけば和を実行する際に扱わなければならない行列のサイズは高々$\chi^2 \times \chi^2$ぐらいで済むので、実用的な速度で計算できます。順序を気にせずに適当に和をとろうとすると扱わなければならない行列のサイズがどんどん増大してしまうため、注意が必要です。ibmq.hのoverlap関数に処理が書いてあります。7

実際に使ってみる

まずは小手調べに、3量子ビットのGHZ状態を作ってみます。回路は
circuit.png
こんな具合です。

\ket{\psi} = \frac{1}{\sqrt{2}}(\ket{000}-\ket{111})

の状態が生成されるはずなので、$\ket{000}$と$\ket{111}$それぞれとのoverlapを計算してみます。(ついでに検算のために$\braket{\psi}{\psi} = 1$となるかどうかも計算しています。)
以下、実行コードです。q[0],q[1],q[2]をシミュレータ上ではそれぞれ7,11,12番目の量子ビットとして扱っています。4

#include "itensor/all.h"
#include "itensor/mps/siteset.h"
#include <stdio.h>
#include <fstream>
#include <cmath>
#include "ibmq.h"

using namespace itensor;
using namespace ibmq;

int main(int argc, char const* argv[]){

    std::array<std::pair<double, double>, IBMQPeps::NUM_BITS> init_param;

    //start with |0> (53qubit)
    for(auto& elem : init_param){
        elem = std::make_pair(1.0, 0.0);
    }

    //回路の初期化 (ここでテンソルネットワーク状態を作る)
    IBMQPeps circuit(init_param);


    //Below is the demonstration of genearating GHZ state

    //the default cursor is located on qubit number 1 and 2
    //shift the cursor to qubit number (7,12)
    //$\Psi$を7,12番目の量子ビットまで"移動"させている。(打ち切り誤差は1e-5ぐらいにしています)
    circuit.shift_to(3, {"Cutoff", 1E-5});
    circuit.shift_to(4, {"Cutoff", 1E-5});
    circuit.shift_to(5, {"Cutoff", 1E-5});
    circuit.shift_to(7, {"Cutoff", 1E-5});
    circuit.shift_to(12, {"Cutoff", 1E-5});

    //$\Psi$を移動させながら量子ゲートを掛けていく
    //apply Hadamard and X to gate (7,12)
    circuit.apply(H(circuit.site(7))*X(circuit.site(12)));
    //shift the cursor to qubit number (11,12)
    circuit.shift_to(11, {"Cutoff", 1E-5});
    //apply Hadamard to gate 11 (Idは恒等演算子なので、12番目の量子ビットにはなにもしません。)
    circuit.apply(H(circuit.site(11))*Id(circuit.site(12)));
    //apply CNOT to gate (11, 12)
    circuit.apply(CNOT(circuit.site(11), circuit.site(12)));
    //shift the cursor to qubit number (7,12)
    circuit.shift_to(7, {"Cutoff", 1E-5});
    //apply CNOT to gate (7, 12)
    circuit.apply(CNOT(circuit.site(7), circuit.site(12)));
    //apply Hadamard to gate (7,12)
    circuit.apply(H(circuit.site(7))*H(circuit.site(12)));
    //shift the cursor to qubit number (11,12)
    circuit.shift_to(11, {"Cutoff", 1E-5});
    //apply Hadamard to gate 11
    circuit.apply(H(circuit.site(11))*Id(circuit.site(12)));
    //the result should be bell state (1/sqrt(2))(|000> + |111>)
    
    //shift the cursor to qubit number (7,12)
    circuit.shift_to(7, {"Cutoff", 1E-5});
    circuit.shift_to(5, {"Cutoff", 1E-5});
    circuit.shift_to(4, {"Cutoff", 1E-5});
    circuit.shift_to(3, {"Cutoff", 1E-5});
    circuit.shift_to(2, {"Cutoff", 1E-5});
    circuit.shift_to(1, {"Cutoff", 1E-5});


    //to show GHZ state is generated, calc the overlap between |0...000....0> and |0...111....0> where 000 and 111 are located on the qubit (7,11,12).
    
    //ここから、|0...000....0>状態と|0...111....0>も生成してみます。
    
    //|0...000....0>
    IBMQPeps circuit000(init_param, circuit.site());

    //|0...111....0>
    IBMQPeps circuit111(init_param, circuit.site());
    //shift the cursor to qubit number (7,12)
    circuit111.shift_to(3, {"Cutoff", 1E-5});
    circuit111.shift_to(4, {"Cutoff", 1E-5});
    circuit111.shift_to(5, {"Cutoff", 1E-5});
    circuit111.shift_to(7, {"Cutoff", 1E-5});
    circuit111.shift_to(12, {"Cutoff", 1E-5});
    //flip the qubit number (7,12)
    circuit111.apply(X(circuit.site(7))*X(circuit.site(12)));
    //shift the cursor to qubit number (11,12)
    circuit111.shift_to(11, {"Cutoff", 1E-5});
    //flip the qubit number 11
    circuit111.apply(X(circuit.site(11))*Id(circuit.site(12)));
    
    myvector<ITensor> op(IBMQPeps::NUM_BITS);
    for(auto i : range1(IBMQPeps::NUM_BITS)){
        op[i] = Id(circuit.site(i));
    }
    
    //overlapの計算
    Print(overlap(circuit, op, circuit000));
    Print(overlap(circuit, op, circuit111));
    //自分自身とのoverlap (1にならなきゃおかしい)
    Print(overlap(circuit, op, circuit));
    return 0;
}

実行した結果、

$ ./mps_mylib
overlap(circuit, op, circuit000) = (-0.707107,0)
overlap(circuit, op, circuit111) = (0.707107,0)
overlap(circuit, op, circuit) = (1,0)

となるので正しそうです。8なんか全体の位相が-1倍されている気がしますが、物理量には特に影響ないので大丈夫でしょう。

3量子ビットでは面白くないので、10数量子ビットぐらい使ってひたすらSWAPゲートを噛まして量子状態を転送させまくってみましょう。

SWAPゲートはCNOT演算子3つを
circuit (1).png
このように組めばいいのでa番目とb番目の量子ビットをswapさせる関数は

void swap(IBMQPeps& circuit, size_t a, size_t b){
    circuit.apply(CNOT(circuit.site(a), circuit.site(b)));
    circuit.apply(CNOT(circuit.site(b), circuit.site(a)));
    circuit.apply(CNOT(circuit.site(a), circuit.site(b)));
}

となります。12番目の量子ビットにアダマール演算子Hをかけ、SWAPゲートで十数量子ビットぐらい経由しながら量子状態を転送させ続け、7番目の量子ビットまで動かした回路circuitと、直接7番目の量子ビットにアダマール演算子Hをかけたcircuit111とoverlapを取ってみます。

#include "itensor/all.h"
#include "itensor/mps/siteset.h"
#include <stdio.h>
#include <fstream>
#include <cmath>
#include "ibmq.h"

using namespace itensor;
using namespace ibmq;

void swap(IBMQPeps& circuit, size_t a, size_t b){
    circuit.apply(CNOT(circuit.site(a), circuit.site(b)));
    circuit.apply(CNOT(circuit.site(b), circuit.site(a)));
    circuit.apply(CNOT(circuit.site(a), circuit.site(b)));
}

int main(int argc, char const* argv[]){

    std::array<std::pair<double, double>, IBMQPeps::NUM_BITS> init_param;

    //start with |0> (53qubit)
    for(auto& elem : init_param){
        elem = std::make_pair(1.0, 0.0);
    }

    IBMQPeps circuit(init_param);


    //Below is the demonstration of SWAP gates

    //the default cursor is located on qubit number 1 and 2
    //shift the cursor to qubit number (7,12)
    circuit.shift_to(3, {"Cutoff", 1E-5});
    circuit.shift_to(4, {"Cutoff", 1E-5});
    circuit.shift_to(5, {"Cutoff", 1E-5});
    circuit.shift_to(7, {"Cutoff", 1E-5});
    circuit.shift_to(12, {"Cutoff", 1E-5});

    //apply H to gate (12)
    circuit.apply(Id(circuit.site(7))*H(circuit.site(12)));
    //shift the cursor to qubit number (11,12)
    circuit.shift_to(11, {"Cutoff", 1E-5});
    swap(circuit, 11, 12);
    circuit.shift_to(10, {"Cutoff", 1E-5});
    swap(circuit, 10, 11);
    circuit.shift_to(9, {"Cutoff", 1E-5});
    swap(circuit, 10, 9);
    circuit.shift_to(8, {"Cutoff", 1E-5});
    swap(circuit, 8, 9);
    circuit.shift_to(6, {"Cutoff", 1E-5});
    swap(circuit, 6, 8);
    circuit.shift_to(1, {"Cutoff", 1E-5});
    swap(circuit, 1, 6);
    circuit.shift_to(2, {"Cutoff", 1E-5});
    swap(circuit, 1, 2);
    circuit.shift_to(3, {"Cutoff", 1E-5});
    swap(circuit, 2, 3);
    circuit.shift_to(4, {"Cutoff", 1E-5});
    swap(circuit, 3, 4);
    circuit.shift_to(5, {"Cutoff", 1E-5});
    swap(circuit, 4, 5);
    circuit.shift_to(7, {"Cutoff", 1E-5});
    swap(circuit, 5, 7);
    //この時点で元々12番目にあった量子ビットは7番目の量子ビットへと転送されています。

    
    //circuit.shift_to(5, {"Cutoff", 1E-5});
    circuit.shift_to(4, {"Cutoff", 1E-5});
    circuit.shift_to(3, {"Cutoff", 1E-5});
    circuit.shift_to(2, {"Cutoff", 1E-5});
    circuit.shift_to(1, {"Cutoff", 1E-5});

    //直接7番目にHを作用させる
    IBMQPeps circuit111(init_param, circuit.site());
    //shift the cursor to qubit number (7,12)
    circuit111.shift_to(3, {"Cutoff", 1E-5});
    circuit111.shift_to(4, {"Cutoff", 1E-5});
    circuit111.shift_to(5, {"Cutoff", 1E-5});
    circuit111.shift_to(7, {"Cutoff", 1E-5});
    circuit111.shift_to(12, {"Cutoff", 1E-5});
    //apply H to the qubit number (7)
    circuit111.apply(H(circuit.site(7))*Id(circuit.site(12)));
    myvector<ITensor> op(IBMQPeps::NUM_BITS);
    for(auto i : range1(IBMQPeps::NUM_BITS)){
        op[i] = Id(circuit.site(i));
    }

    Print(overlap(circuit, op, circuit111));

	return 0;
}

実行してみると

$ ./mps_mylib 
overlap(circuit, op, circuit111) = (1,0)

と結果が一致するのでまあ問題なさそうです。
一見するとSWAPゲートを作用させただけなので対して面白くは見えないですが、テンソルネットワークを用いた手法だと100量子ビットを経由しようが1000量子ビットを経由しようが何の問題も起こりません。秒もかからずに答えが出てきます。(愚直に波動関数を計算する従来の手法だとこの計算でも$2^{1000} \times 2^{1000}$ぐらいのサイズの行列の演算が必要になります。無理です。)

まとめ

本記事ではテンソルネットワークを用いた手法でIBMの量子コンピュータのシミュレータを作りました。

一見するとテンソルネットワークを用いれば量子コンピュータなんてなくても古典コンピュータで全てシミュレートできるような気もしてきますが別にそんなこともなく、各量子ビット間の相関が大きい (強くエンタングルメントしている)場合には打ち切り誤差$\epsilon$を小さくするために$\chi$をかなり大きな値にしなければならず、結果的に古典コンピュータでも計算が困難になる場合は多々あります。9

参考までに、IBMの量子コンピュータの外周全て(大体40量子ビットぐらい)にランダムな2量子ビットユニタリー演算を作用させるシミュレーションを何度も行わせた結果、その回数 (横軸)に応じた計算時間 (縦軸)は大体こんな感じになりました。
2019-12-04-021419_588x357_scrot.png
横軸はおおよそ量子回路のdepthに対応する量なので、この結果によると大体depthが20ぐらいならシミュレータでなんとかなるみたいです。(かなり大雑把に調べた結果なので大して厳密な結果ではないです。)

しかしながらそれでも愚直に$2^N$の要素を持つベクトルを扱って行列演算するよりは遥かに効率がいいため、そこそこ大きなサイズの量子コンピュータのシミュレータとして実用的であると思います。

今回の記事ではまだ大した量子計算をさせていないので、これからいろいろ遊んでみようと思います。

追記 (2021/11/03)

Googleの量子超越に関する回路に関しても、テンソルネットワークベースのアルゴリズムで304秒で解けたというニュースが流れてきました。10
基本的な手法は上記の手法であると思いますので、今後詳細を読んでみようと思います。

  1. この例で扱っている波動関数は1次元的に量子ビット間に相関を導入しており、テンソルネットワークの中でも特に**行列積状態 (Matrix Product State)**と呼ばれる名前がついています。 ↩

  2. Steven R. White, "Density matrix formulation for quantum renormalization groups", Phys. Rev. Lett. 69, 2863 (1992) ↩

  3. https://arxiv.org/abs/1008.3477v2 ↩ ↩2

  4. シミュレーション内で用いている量子ビットのインデックスはIBM Qの量子コンピュータの図中のインデックスとは対応してないので注意してください ↩ ↩2

  5. generate_linkで後からリンクを作成しているため、実をいうとIBMだけでなくちゃっかり任意の量子コンピュータのトポロジーに対応してたりします。 ↩

  6. 後にこれ書いてから、overlapよりも、$n$番目の量子状態を測定した結果が$q_n$となる確率$P(q_n)$を求めるほうがシミュレータとして自然なのではと気づきましたが、これは一般に密度行列を用いて$P(q_n) = \mathrm{Tr}(\ket{q_n}\bra{q_n}\otimes \hat{\rho})$とかけるはずなのでどちらにせよ演算を工夫すれば下に書いたように簡単に計算できます。 ↩

  7. 私の実装したoverlap関数はより一般に$\bra{\psi_1}\hat{M}\ket{\psi_2}$ のような量を計算できるように作ってあります。$\hat{M}$を恒等演算子にしておけば$\braket{\psi_1}{\psi_2}$になりますし、$n$番目の量子ビットだけ$\ket{q_n}\bra{q_n}$、残りを恒等演算子にすれば$\bra{\psi}\hat{M}\ket{\psi}$が上の注釈で書いた$P(q_n)$になります。 ↩

  8. (a,b)は複素数で$a+bi$を表しています。 ↩

  9. おそらくGoogleの量子超越で用いられれた回路もこうなっているのでしょう。まだ実験とかしていないので何とも言えないですが... ↩

  10. https://arxiv.org/abs/2110.14502 ↩

47
26
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
47
26

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?