$$
\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}}
$$

はじめに

「量子情報理論の基本」と題して、これまでいろんな話題を取り上げてきましたが、量子状態を表現するために、いつでも登場するのが「密度演算子」というものでした。特に、混合状態になっている量子系を考える際に、使わざるを得ない非常に便利な理論的な道具なので、もうこれなしでは生きていけません。とても便利なので、未知の物理系があったときに、「この量子状態はいかなる密度演算子で表せるものであるか？」ということを思わず問いたくなるわけです（なってください！）。実は、「量子状態トモグラフィ(quantum state tomograpy)」という手法を使うと、それが推定できます。今回の記事では、まず、その手法について説明した後、量子計算シミュレータqlazyを使って、任意の量子状態に対する密度演算子を推定してみたいと思います。

参考にさせていただいたのは、以下の文献です。

量子状態トモグラフィとは

混合状態について（復習）

密度演算子を推定する前に、混合状態がいかにして現れてくるのかについて、少し復習しておきます。どんな場合にこの推定手法が役に立つのかを把握しておくためです。大きく分けると以下に示す通り、２パターンあると思います。

[1] 純粋状態の確率的なアンサンブル
[2] 純粋状態の部分系に注目

[1]は、例えば、複数パターンの純粋状態を発射できる装置をアリスが持っていて、それを使って純粋状態を確率的に選択して次々発射する場合に現れる混合状態です。任意に偏光角を設定できる光子発射装置をイメージすれば良いです。あるいは、何らかの純粋状態を測定したのに、その測定結果を忘れてしまった、という場合もこのパターンの混合状態になると思います¹。

純粋状態の確率的なアンサンブルを$\{ p_i, \ket{\psi_i}\}$と表すと、密度演算子$\rho$は、

\rho = \sum_{i} p_i \ket{\psi_i} \bra{\psi_i}  \tag{1}

と表すことができます。

[2]は、例えば、全体として閉じた系（＝純粋状態）の部分系Sに注目することで現れる混合状態です²。部分系Sは、一般にそれ以外の残りの系（環境系E）と相互作用していると考えられるので、開いた系（＝混合状態）になります。あるいは、量子コンピュータで量子回路を実行しているとき、その一部の量子ビットに注目すると、一般に混合状態になっています³。

全体系が$\ket{\Psi}$だったとき、部分系への注目を数式で書くと、

\rho = Tr_E (\ket{\Psi} \bra{\Psi})  \tag{2}

となります。

密度演算子の推定（１量子状態）

さて、混合状態の具体的な現れ方を軽く復習したところで、では、そのときの密度演算子をどうやって推定するかについて説明します。が、ここで一つの前提を置きます。$\rho$という状態がたった1回だけ現れて、測定する機会が1回しかない場合は、残念ながら推定は不可能です。測定によって得られる数値は1個である一方、密度演算子を記述するための数値はたくさんあるので、不良設定問題となり推定できません。従って、$\rho$という状態を繰り返し発生させることができて、それを都度測定できる、という前提とします。その測定結果を多数集めて、総合することで密度演算子を推定するというのが基本的な考え方です。

では、まず、１量子状態に対する密度演算子について考えてみます。

１量子の場合に限りますが、密度演算子はパウリ演算子の線形和として表現することができました⁴。この表現が推定方法を説明する際のベースになるので、それがどうやって導かれたのか、やや詳しく見ていきます。

パウリ演算子は以下のような行列として定義されました。

\sigma_0 = I = 
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix}, \space
\sigma_1 =  
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}, \space
\sigma_2 =
\begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix}, \space
\sigma_3 =
\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}  \tag{3}

この定義から明らかですが、パウリ演算子はエルミート演算子です。かつ、

\frac{1}{2} Tr(\sigma_i \sigma_j) = \delta_{ij}  \tag{4}

という性質があります。ちなみに、ここに登場する$Tr(AB)$は、行列$A,B$をベクトルと見なしたときの一種の内積であり、「ヒルベルト・シュミット内積」と呼ばれることがあります（後で出てくるので覚えておいてください）。

パウリ演算子のこの性質に注意すると、２次元ヒルベルト空間上で定義される任意のエルミート演算子$A$は、$\{ u_i\} \space (i=0,1,2,3)$を任意の実数として、以下のような線形和で表現できることがわかります。

A = u_0 \sigma_0 + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 = \sum_{i=0}^{3} u_i \sigma_i \tag{5}

簡単なので、ちょっと証明してみます。

【証明】

まず、式(3)の右辺がエルミートであるかどうかです。

\begin{align}
A &=
u_0
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix} +
u_1
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix} +
u_2
\begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix} +
u_3
\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}  \\
&= 
\begin{pmatrix}
u_0 & 0 \\
0 & u_0
\end{pmatrix} +
\begin{pmatrix}
0 & u_1 \\
u_1 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & -i u_2 \\
i u_2 & 0
\end{pmatrix} +
\begin{pmatrix}
u_3 & 0 \\
0 & -u_3
\end{pmatrix}  \\
&=
\begin{pmatrix}
u_0 + u_3 & u_1 - iu_2 \\
u_1 + iu_2 & u_0 - u_3
\end{pmatrix}  \tag{6}
\end{align}

となり、明らかに$A^{\dagger}=A$なので、$A$はエルミートです。

次に、逆の証明です。つまり、任意のエルミート演算子が、式(5)のように表現できるかどうかです。

$A$がエルミートであるとすると、任意の実数$a,b,c,d$を使って、

A =
\begin{pmatrix}
a & b-ic \\
b+ic & d
\end{pmatrix}  \tag{7}

のように表せます。ここで、どんな$a,b,c,d$を持ってきたとしても、

\begin{align}
a &= u_0 + u_3 \\
b &= u_1 \\
c &= u_2 \\
d &= u_0 - u_3 \tag{8}
\end{align}

を満たす実数値$u_0,u_1,u_2,u_3$は必ず存在するので、式(7)は、

\begin{align}
A &=
\begin{pmatrix}
u_0 + u_3 & u_1 - iu_2 \\
u_1 + iu_2 & u_0 - u_3
\end{pmatrix} \\
&=
u_0
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix} +
u_1
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix} +
u_2
\begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix} +
u_3
\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}  \\
&= u_0 \sigma_0 + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 \\
&= \sum_{i=0}^{3} u_i \sigma_i \tag{9}
\end{align}

となり、式(5)が証明できました。

（証明終）

ここで、式(5)が意味していることについて、少し別の言い方をしてみます。すなわち、２次元ヒルベルト空間上で定義されたエルミート演算子は、４つのパウリ演算子をヒルベルト・シュミット内積の意味での直交基底としたときに、それによって張られる線形空間上のベクトルと見なすことができます、と言えます。

ん？「演算子（＝行列）」と言っているのに「ベクトル」とはどういうこと？？、という疑問の声が聞こえてきそうですが、ここは頭を柔らかくして、ついてきてください。ざっくり言うと、内積が定義できて直交基底があって、その線形和で表わせるものがあれば、何でも「ベクトル」と言って良いのです、と思うことにしましょう。

で、１量子の密度演算子$\rho$は、まさに、２次元ヒルベルト空間上でのエルミート演算子だったので、たったいま証明した式(5)のように展開できます。

\rho = u_0 \sigma_0 + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 = \sum_{i=0}^{3} u_i \sigma_i \tag{10}

ただし、密度演算子には$Tr(\rho)=1$という制約があり、かつ、$Tr(\sigma_1)=Tr(\sigma_2)=Tr(\sigma_3)=0$なので、$u_0 = Tr(\rho) = 1$となります。つまり、変数が１個減り、

\rho = I + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 = I + \sum_{i=1}^{3} u_i \sigma_i \tag{11}

と書けます。ここで、式(4)のパウリ演算子の性質に注意すると、

u_i = \frac{1}{2} Tr(\sigma_i \rho)  \tag{12}

が成り立つので、式(11)は、結局、

\begin{align}
\rho &= \frac{1}{2} (Tr(\rho) I + Tr(\sigma_1 \rho) \sigma_1 + Tr(\sigma_2 \rho) \sigma_2+ Tr(\sigma_3 \rho) \sigma_3) \\
&= \frac{1}{2} \sum_{i=0}^{3} Tr(\sigma_i \rho) \sigma_i  \tag{13}
\end{align}

となります。

さて、ようやく密度演算子の推定方法について語るときがやって来ました。

式(13)をよく見てください。$i=1,2,3$に対する$Tr(\sigma_i \rho)$は、オブザーバブル$\sigma_i$を測定したときの期待値でした。$\sigma_1,\sigma_2,\sigma_3$は、各々X軸、Y軸、Z軸方向のスピンのようなものです。従って、飛んできた量子に対して、X軸方向の測定を何度も実行して、その測定値(+1,-1)を記録して期待値を求め、次に、Y軸方向の測定を何度も実行して同様に期待値を求め、さらに、Z軸方向の測定を何度も実行して期待値を求めて、式(13)に代入します。$Tr(\sigma_0 \rho)$の$\sigma_0$は単位行列で、$\rho$のトレースは1なので有無を言わさず、この項は1にしておけば良いです。ということで、式(13)は完全に計算できます。

いま、アリスが物理的に発射した量子状態を何となくイメージしながら説明しましたが、別の例として、量子コンピュータの一つの量子ビットに関する密度演算子を推定したい場合もあると思います。基本的には同様ですが、量子ビットの測定方向を都度切り替えるところだけだけ要注意です。Z方向の測定しかサポートされていないシミュレータやクラウドサービスをお使いの場合は、Z方向の測定の前に、測定方向を変える量子ゲートを噛ませる必要があります。

X方向の測定は、

--H--M--

Y方向の測定は、

--S+--H--M--

Z方向の測定は、

--M--

でいけます。ブロッホ球をイメージしていただければ、単に軸の回転を実行しているだけということがわかると思います。

密度演算子の推定（２量子状態）

２量子状態の密度演算子についても、いままでの議論を拡張することができます。

１量子の場合の密度演算子は、４つのパウリ演算子を基底とする線形空間のベクトルと思って良い、という説明をしました。２量子状態の場合、この線形空間をもう一つ持ってきて直積空間をつくって、その上でベクトルを定義すれば良いです。この直積空間の基底は４つのパウリ演算子のテンソル積で規定されます。つまり、

\{\sigma_i \otimes \sigma_j \} \space (i,j = 0,1,2,3)  \tag{14}

という16個の4X4行列が、この直積空間の基底となります。容易にわかる通り、この基底は、

\frac{1}{2^2} Tr[(\sigma_i \otimes \sigma_j)(\sigma_k \otimes \sigma_l)] = \delta_{ik} \delta_{jl} \tag{15}

という直交条件を満たしますので、先程と同様に、この２量子状態の密度演算子は、

\rho = \sum_{i}^{3} \sum_{j}^{3} u_{ij} \sigma_i \otimes \sigma_j \tag{16}

のようにパウリ演算子を使った基底で展開できます。式(15)に注意すると、この$u_{ij}$は、

u_{ij} = \frac{1}{2^2} Tr[(\sigma_i \otimes \sigma_j) \rho]  \tag{17}

と計算できますので、式(16)は、

\rho = \frac{1}{2^2} \sum_{i=0}^{3} \sum_{j=0}^{3} Tr[(\sigma_i \otimes \sigma_j) \rho] \cdot \sigma_i \otimes \sigma_j  \tag{18}

となります。

密度演算子を推定するには、２量子の測定を各々の方向を切り替えながら実行する必要があります。パウリ演算子の組み合わせは4X4=16通りあるので、その各々の場合について、期待値を求めます。具体的には、まず、パウリ演算子の一つの組み合わせについて、飛んできた２つの量子各々の測定値(+1,-1)を得ます。その積を２量子全体の測定値として記録します。これを何度も繰り返して、期待値を得ます。これで、パウリ演算子の一つの組み合わせについての期待値が得られるので、これをパウリ演算子の全組み合わせ（16パターン）で実行します。16個の期待値が得られますので、それを式(17)に代入することで、密度演算子の推定が完了します。ただし、すべてのパウリ演算子が単位演算子である場合が一つ出てきますので、これは期待値を計算するまでもないです。有無を言わさず１にしておきます。以上です。

密度演算子の推定（n量子状態）

さらに、n量子状態の密度演算子についても、同様の議論でOKです。密度演算子は、

\rho = \frac{1}{2^n} \sum_{\mu_1=0}^{3} \sum_{\mu_2=0}^{3} \cdots \sum_{\mu_{n}=0}^{3} Tr[(\sigma_{\mu_1} \otimes \sigma_{\mu_2} \otimes \cdots \otimes \sigma_{\mu_{n}})\rho] \cdot \sigma_{\mu_1} \otimes \sigma_{\mu_2} \otimes \cdots \otimes \sigma_{\mu_{n}} \tag{19}

のように表すことができます。ちょっと複雑に見えますが、よく見れば、じんわりわかってくると思います。

密度演算子の推定方法は、

Tr[(\sigma_{\mu_1} \otimes \sigma_{\mu_2} \otimes \cdots \otimes \sigma_{\mu_{n}})\rho]   \tag{20}

という期待値を求めれば良いです。n量子なのでパウリ演算子の組み合わせは4のn乗です。その各々の組み合わせについて、各量子ビットに割り当てられたパウリ演算子に対応した方向の測定を実行して、n量子の測定値（各量子の測定値の積）を得て、期待値を求めれば良いです。

以上のように推定する手法が「量子状態トモグラフィ」です。

推定シミュレーション

実装

それでは、任意の量子状態アンサンブルを作成し、量子状態トモグラフィによって密度演算子を推定してみたいと思います⁵。そんなまだるっこしいことをやらなくても、量子計算シミュレータqlazyには、量子状態アンサンブルから行列計算によって密度演算子をダイレクトに取得できる機能が搭載されています。が、今回は、量子状態トモグラフィの勉強なので、一旦忘れて、量子状態アンサンブルの測定のみによって、密度演算子を推定してみます。その推定の精度を測るための参照値（真の値）として、qlazyで計算した密度演算子を使います。精度は忠実度（フィデリティ）で評価してみます（忠実度もqlazyの機能として搭載されています）。

それでは、全体のPythonコードを示します。

【2021.9.5追記】qlazy最新版でのソースコードはここに置いてあります。

import random
import math
import numpy as np
from scipy.stats import unitary_group
from qlazypy import QState,DensOp

def random_qstate_ensemble(qubit_num, mimxed_num):

    qstate = []
    dim = 2**qubit_num
    for _ in range(mixed_num):
        vec_ini = np.zeros(dim)
        vec_ini[0] = 1.0
        mat = unitary_group.rvs(dim)
        vec = np.dot(mat, vec_ini)
        qstate.append(QState(vector=vec))
        
    prob = [random.random() for _ in range(mixed_num)]
    total = sum(prob)
    prob = [p/total for p in prob]

    return (qstate,prob)

def get_pauli_index(index, total):
    # ex) index = 9 = 1001 -> [10,01] -> [2,1] --reverse-> [1,2] = pauli_index
    #     '1' means X for 0th-qubit, '2' means Y for 1st-qubit

    pauli_index = [0]*int(math.log2(total)/2)
    count = 0
    while index > 0:
        pauli_index[count] = index%4
        index = index // 4
        count += 1
    pauli_index.reverse()

    return pauli_index

def make_pauli_product(index, total, pauli_mat):

    pauli_index = get_pauli_index(index, total)
    pauli_prod_dim = 2**len(pauli_index) 
    pauli_prod = np.array([1.0])
    for pid in pauli_index:
        pauli_prod = np.kron(pauli_prod, pauli_mat[pid])

    return pauli_prod
    
def make_densop(expect, qubit_num, pauli_mat):

    dim = 2**qubit_num
    measure_num = len(expect)
    matrix = np.zeros((dim,dim))
    for i in range(measure_num):
        pauli_prod = make_pauli_product(i,measure_num,pauli_mat)
        matrix = matrix + expect[i] * pauli_prod
    matrix = matrix / (2.0**qubit_num)
    
    return DensOp(matrix=matrix)

def calc_expect(prb):

    N = len(prb)
    val = np.zeros(N)
    for index in range(N):
        bin_list = [int(s) for s in list(format(index, 'b'))]
        val[index] = (-1)**sum(bin_list)  # odd -> -1, even -> +1

    return np.dot(prb, val)

def estimate_densop(prob,qstate,shots):

    pauli_mat = [np.eye(2),                   # = I
                 np.array([[0,1],[1,0]]),     # = X
                 np.array([[0,-1j],[1j,0]]),  # = Y
                 np.array([[1,0],[0,-1]])]    # = Z
    
    mixed_num = len(prob)
    qubit_num = qstate[0].qubit_num
    measure_num = 4**qubit_num

    for i in range(mixed_num):

        expect = {}
        for index in range(measure_num):

            pauli_index = get_pauli_index(index,measure_num)

            if index == 0:
                expect[index] = 1.0
                continue
                    
            qs = qstate[i].clone()

            for qid in range(len(pauli_index)):

                if pauli_index[qid] == 0:
                    continue
                elif pauli_index[qid] == 1:
                    qs.h(qid)
                elif pauli_index[qid] == 2:
                    qs.s_dg(qid).h(qid)
                else:
                    pass

            frq = qs.m(shots=shots).frq
            prb = np.array([f/shots for f in frq])
            expect[index] = calc_expect(prb)

            qs.free()

        de_tmp = make_densop(expect, qubit_num, pauli_mat)
        
        if i == 0:
            de_tmp.mul(prob[i])
            densop_est = de_tmp.clone()
        else:
            de_tmp.mul(prob[i])
            densop_est.add(de_tmp)

    de_tmp.free()

    return densop_est
    
if __name__ == '__main__':

    # settings
    qubit_num = 2
    mixed_num = 4
    shots = 100

    # quantum state ensemble (original)
    qstate, prob = random_qstate_ensemble(qubit_num, mixed_num)

    # density operator (original)
    densop_ori = DensOp(qstate=qstate, prob=prob)

    # density operator estimation only from quantum state ensemble
    # (quantum state tomography)
    densop_est = estimate_densop(prob,qstate,shots)

    print("** density operator (original)")
    densop_ori.show()
    print("** density operator (estimated)")
    densop_est.show()
    print("** fidelity =",densop_ori.fidelity(densop_est))

    for q in qstate:
        q.free()
    densop_ori.free()
    densop_est.free()

何をやっているか、説明します。main処理部を見てください。

# settings
qubit_num = 2
mixed_num = 4
shots = 1000

これはシミュレーションのための設定値です。qubit_numは前提とする量子の数、mixed_numは混合する純粋状態の数、shotsは測定の回数を表します。

# quantum state ensemble (original)
qstate, prob = random_qstate_ensemble(qubit_num, mixed_num)

ここで、ランダムに量子状態アンサンブルを作成しています。qstateは量子状態４つからなるリストです。probは４つの純粋状態をどんな割合で混ぜるかを決める確率のリストです。純粋状態も確率もどちらもランダムに決めています。詳細は関数定義を見てください。

# density operator (original)
densop_ori = DensOp(qstate=qstate, prob=prob)

これは、qlazyを使って、今考えている密度演算子の真の値を計算している部分です。最後に精度評価するための参照値として使います。

# density operator estimation only from quantum state ensemble
# (quantum state tomography)
densop_est = estimate_densop(prob,qstate,shots)

ここが、今回の肝になる密度演算子推定の部分になりますので、少し詳しく説明します。関数estimate_densopの中身を見てください。

def estimate_densop(prob,qstate,shots):

    pauli_mat = [np.eye(2),                   # = I
                 np.array([[0,1],[1,0]]),     # = X
                 np.array([[0,-1j],[1j,0]]),  # = Y
                 np.array([[1,0],[0,-1]])]    # = Z
    
    mixed_num = len(prob)
    qubit_num = qstate[0].qubit_num
    measure_num = 4**qubit_num
    ...

最初のpauli_matはパウリ行列を設定しているところです。mixed_num,qubit_numは引数から改めて求めています。qubit_numからパウリ演算子の組み合わせ数が4**qubit_numで求まります(これを変数measure_numに格納しています)。

次のforループは、

for i in range(mixed_num):

    expect = {}
    for index in range(measure_num):

        pauli_index = get_pauli_index(index,measure_num)
	...

のように２重になっています。外側のforループは、混合されている各々の純粋状態ごとに密度演算子を推定するためのものです。後で、確率の重みをつけて各々を足し合わせて、最終的な密度演算子とします⁵。

内側のforループは、パウリ演算子の組み回せ数分、期待値を計算するためのものです。各々の組み合わせに対応した期待値を格納するための辞書expectを最初に空辞書として初期化しておきます。先程のmeasure_numの数分、繰り返すのですが、その組み合わせ番号が変数indexに入ります。これは単なる整数値なので、各量子ビットの測定方向を制御するためには扱いにくいです。ということで、関数get_pauli_indexで、パウリ演算子番号のリストに変換します。整数値を４進数に直して、各桁をリストにしてリターンします。詳細は関数定義を見てください。

次に、そのforループの中身です。

if index == 0:
    expect[index] = 1.0
    continue
        
qs = qstate[i].clone()

for qid in range(len(pauli_index)):

    if pauli_index[qid] == 0:
        continue
    elif pauli_index[qid] == 1:
        qs.h(qid)
    elif pauli_index[qid] == 2:
        qs.s_dg(qid).h(qid)
    else:
        pass

frq = qs.m(shots=shots).frq
prb = np.array([f/shots for f in frq])
expect[index] = calc_expect(prb)

qs.free()

index=0の場合は、有無を言わさず期待値を１にします。そうでない場合は、

qs = qstate[i].clone()

で、i番目の純粋状態をコピーしてテンポラリとしてqsに格納します⁶。

次に、またforループが出てきました。

for qid in range(len(pauli_index)):
...

これは、パウリ演算子番号に相当した測定方向で各量子ビットを順に測定するため、測定前にかける量子ゲートのためのものです。中身は上記、見ての通りです。そのループが終わったら、

frq = qs.m(shots=shots).frq
prb = np.array([f/shots for f in frq])
expect[index] = calc_expect(prb)

qs.free()

としています。これで、実際の測定をやっています。shotは1000だったので1000回測定します⁷。結果、頻度リストが得られるので、確率にするため総和が１になるように正規化します。それを関数calc_expectに渡して、期待値を計算します。いま２量子状態を想定しているので、確率は４つになります。その各々は、特定のパウリ演算子の組を表しています。特定のパウリ演算子の組に対して、２量子の測定値は各々の積として計算できます（+1または-1になります）。４つの確率と４つの測定値が決まれば、期待値は計算できます。関数calc_expectは、その処理をやっています。詳細は関数定義を見てください。qsを使い終わったら、qs.free()でメモリを解放します⁸。

de_tmp = make_densop(expect, qubit_num, pauli_mat)

if i == 0:
    de_tmp.mul(prob[i])
    densop_est = de_tmp.clone()
else:
    de_tmp.mul(prob[i])
    densop_est.add(de_tmp)

４つの期待値がすべて求まったら、関数make_densopで、一旦、この純粋状態に対する密度演算子を計算します（make_densopの詳細は関数定義を見てください）。それを混合数（４つとしました）の分だけ、確率の重みをつけて足し合わせます。if文はその処理です。外側のループを抜けると、最終的なdensop_estが求まったことになるので、これをリターンします。

main処理部に戻ってきて、

print("** density operator (original)")
densop_ori.show()
print("** density operator (estimated)")
densop_est.show()
print("** fidelity =",densop_ori.fidelity(densop_est))

これで、本物の密度演算子(densop_ori)の要素と、推定した密度演算子(densop_est)の要素を順に表示し、最後に両者の忠実度を計算して、表示しています。

結果

シミュレーションの結果を示します。

上のプログラムは２量子状態の推定のためのものでしたが、まず最初は、１量子状態の推定です。パラメータは、

qubit_num = 1
mixed_num = 2
shots = 100

としました。実行したところ、

** density operator (original)
elm[0][0] = +0.3270-0.0000*i : 0.1069 |++
elm[0][1] = +0.1508-0.2138*i : 0.0684 |++
elm[1][0] = +0.1508+0.2138*i : 0.0684 |++
elm[1][1] = +0.6730+0.0000*i : 0.4530 |++++++
** density operator (estimated)
elm[0][0] = +0.3224+0.0000*i : 0.1040 |++
elm[0][1] = +0.1584-0.1887*i : 0.0607 |++
elm[1][0] = +0.1584+0.1887*i : 0.0607 |++
elm[1][1] = +0.6776+0.0000*i : 0.4591 |++++++
** fidelity = 0.9996155234243419

となりました。忠実度が0.9996...なので、結構な精度で推定できたと思います。

次に、２量子状態の推定です。パラメータは、

qubit_num = 2
mixed_num = 4
shots = 1000

としてみました。密度演算子の要素は長くなるので、忠実度だけ示すと、

** fidelity = 0.9814099144563044

となり、少し精度は落ちました。が、とりあえず、量子状態トモグラフィのお勉強という目的は達成できたと思うので、このあたりで良しとしておきます⁹。

おわりに

今回、状態を推定するための「量子状態トモグラフィ」を取り上げましたが、「状態」ではなく「プロセス」の方を推定する「量子プロセス・トモグラフィ」という手法もあります。物理系の変化を表すクラウス演算子を推定する手法で、基本的には、今回の「状態」トモグラフィを繰り返し使う手法のようです。が、未勉強なので、またの機会とさせてください。

さて、次回は、量子通信関係の話題を考えています。古典情報ではなく量子状態そのものを量子通信チャンネルに乗せて通信する際の手法である「データ圧縮」を取り上げてみたいと思っています（予定は未定ですが）。

以上

ちなみに、このような測定のことを「非選択的な測定」と言います。量子情報の教科書を読むと、「忘れてしまう」とか「測定結果をなくしてしまう」という言い方によく出会います。そんな間抜けな実験家はいないと思いますが、言葉の綾と思いましょう（自分の学生時代、時々、非選択的な測定をやってしまっていましたが、汗）。逆に、測定結果をしっかり覚えておく測定のことを「選択的な測定」と言います。その場合、系の状態は純粋状態になります。 ↩
量子状態の面白いところですね。状態が確定した純粋状態の全体系を分割すると、「もつれ」がほどけて混合状態になります。別の言い方をすると、エントロピーがゼロの系を単に分割するだけで、なぜか部分系にエントロピーが発生するということでもあります！ ↩
量子計算は、全体としては純粋状態に対するユニタリ変換なので、全体の量子状態を見ると常に純粋状態を保ったまま推移していきますが、その一部の量子ビットに注目すると、混合状態になっています。ただし、注目している部分量子ビットが、他の量子ビットと独立＝もつれていない＝積状態になっている場合は純粋状態です。 ↩
「量子チャネル」を取り上げた以前の記事で、密度演算子のブロッホ表示について説明しましたが、そのことです。 ↩
今回のシミュレーションは、量子状態アンサンブルから現れる混合状態の推定としました。が、純粋状態の部分系に注目することで現れる混合状態の推定も、プログラムを少し（？）変えることでできると思います。 ↩ ↩²
シミュレータなので気軽にコピーしちゃっていますが、実際の実験または本物の量子コンピュータでやる場合、この部分は量子発生装置から何度も同じものを発生させたり、あるいは同じ量子回路を何度も通してその状態を得るという操作が必要です。 ↩
ここでも、シミュレータならではの処理で実装をサボっています。リアルな実験ではqsは本当に何度も発生させて測定しないといけないです。 ↩
同じ量子状態を何度もループして使う場合、この処理は大事です。qlazyは中でC言語のライブラリを使っていて、メモリ確保もCの世界でやっています。そのため、量子状態と密度演算子については明示的にメモリ解放すべし、という仕様にしています。Pythonの__del__に任せて、明示的な解放をしない仕様を試したこともあったのですが、Python側でどんなタイミングでメモリ解放されるのかよくわからず、怪しい動きをする場合があったので、こんな仕様にしています。 ↩
量子情報工学に次のような記載がありました。「１量子ビットの場合はこのように比較的簡単だが、多量子ビットの状態を決めるためには量子ビット数の２乗のオーダーで増加する密度演算子の成分を決めなければならない。そのため、必要な測定の数が増加する。さらに、状態がたまたま純粋状態であると、密度行列のランクが１になる（１個の対角要素を残して他の要素がすべて０になる）ため、不要なパラメータも推定してしまう状況になる」とのことで、なかなか難しいようです。実際にシミュレーションをしてみると、こういう難しさも実感できます。 ↩