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

はじめに

全体系の量子状態が純粋状態の場合、その部分系に注目すると、一般には混合状態になります。その逆に、混合状態があった場合、それに何らかの外部系を加えることで常に純粋状態にできます。この手続きのことを「純粋化(purification)」と言います。今回はそのお話です。純粋化の具体的な手順とそこから導き出される知見の一つを説明した後、量子計算シミュレータqlazyを使って、純粋化の手続きを実際に計算してみたいと思います。

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

純粋化とは

いま注目している系Aで定義された混合状態$\rho_A$があったとします。これにある別の系Rを加えることで、

\rho_A = Tr_R(\ket{\Phi}^{AR} \bra{\Phi}^{AR})  \tag{1}

をみたす、A+R系の純粋状態$\ket{\Phi}^{AR}$を常に作ることができます。この操作のことを「純粋化」と言います。これは純粋に数学的な手続きであり、何らかの物理的な変化や操作に対応しているわけではありません、ということにご注意ください。これを用いて任意の混合状態（ヒルベルト空間上でのベクトルの確率統計的なアンサンブル）をある純粋状態（ヒルベルト空間上の１つのベクトル）に結びつけることができるようになります¹。実際の物理系の時間変化は、一般には混合状態の時間変化（＝CPTPマップ）として見えるわけですが、理論解析を進める際には、純粋状態とユニタリ変換で理論計算を実施しておいて、最後に部分トレースして現実系に落とし込む方が楽な場合が多いので、量子情報理論における一つの重要テクニックとして「純粋化」を覚えておけば良いのだと思います。

それでは、純粋化の手続きを具体的に見ていきます。まず、$\rho_A$は密度演算子でありエルミートなので、以下のように対角化してスペクトル分解できます（dは密度演算子$\rho_A$のランクです）。

\rho_A = \sum_{k=1}^{d} \lambda_k \ket{\psi_k}^{A} \bra{\psi_k}^{A}  \tag{2}

ここで、Aとは別のヒルベルト空間Rにおける適当な正規直交基底$\{ \ket{\mu_k}^{R} \}$を使って、$\ket{\Psi}^{AR}$を、

\ket{\Psi}^{AR} = \sum_{k=1}^{d} \sqrt{\lambda_k} \ket{\psi_k}^{A} \ket{\mu_k}^{R}  \tag{3}

のように作ります。純粋化の手順は以上です（簡単！）。得られた$\ket{\Phi}^{AR}$が、式(1)を本当に満たしているかを確認します。

式(1)の右辺に式(3)を代入します。

\begin{align}
Tr_R(\ket{\Phi}^{AR} \bra{\Phi}^{AR}) &= Tr_R(\sum_{j=1}^{d} \sum_{k=1}^{d} \sqrt{\lambda_j \lambda_k} \ket{\psi_j}^{A} \ket{\mu_j}^{R} \bra{\psi_k}^{A} \bra{\mu_k}^{R}) \\
&= \sum_{i=1}^{d} \sum_{j=1}^{d} \sum_{k=1}^{d} \sqrt{\lambda_j \lambda_k} \bra{\mu_i}^{R} \ket{\psi_j}^{A} \ket{\mu_j}^{R} \bra{\psi_k}^{A} \bra{\mu_k}^{R} \ket{\mu_i}^{R}　\\
&= \sum_{i=1}^{d} \sum_{j=1}^{d} \sum_{k=1}^{d} \sqrt{\lambda_j \lambda_k} \delta_{ij} \delta_{ik} \ket{\psi_j}^{A} \bra{\psi_k}^{A} \\
&= \sum_{k=1}^{d} \lambda_k \ket{\psi_k}^{A} \bra{\psi_k}^{A} \tag{4}
\end{align}

最後の行は式(2)の右辺なので、これは$\rho_A$に等しいです。というわけで、確認できました。

ここで、要注意ポイントを一つ。ヒルベルト空間Rの正規直交系を適当に選ぶと言いましたが、これは本当にどんな正規直交系でも良いという意味です。上の式展開を改めて見ていただくとわかる通り、この選び方にはユニタリの自由度があります。つまり、$\{ \ket{\mu_k}^{R} \}$に任意のユニタリ変換を施したものを新たな$\{ \ket{\mu_k}^{R} \}$と置いても、上の議論は全く同じようにできます。結局、要するに、R系における正規直交系であれば何でも良いということになります。ということは、与えらた混合状態に対する純粋化の結果は一つには決まりません、ということです。覚えておきましょう。

アンシラの大きさ

純粋化がわかったところで、そこから導出できる一つの知見について説明してみます。以前の記事でCPTPマップを説明した際、注目系を含む大きな閉じた全体系（例えば宇宙全体）で定義された純粋状態に対しユニタリ変換して部分系に注目（部分トレース）することで、現実の物理系の状態変化がCPTPマップとして（あるいはKraus表現として）記述できるということでした。では、現実の物理系があったときに、どのくらいの大きさの補助系（よく「アンシラ(ancilla)系」と言ったりします）を用意しておけば十分と言えるのでしょうか、という疑問が湧いてきます。頭の中で形式的な理論構築だけをやっている分には、あまり気にしなくても良いかもしれませんが、理論に基づき現実の計算をやろうとすると、途端に問題になります。宇宙全体もしくはとてつもなく大きな次元の空間を用意しないといけないとなると、現実的に困ってしまいます。が、実は心配ご無用でして、注目系の次元数を２乗した次元の空間をアンシラ系として用意しておけば十分！ということが、純粋化の知識を使えばわかります。

いま、注目系Aのヒルベルト空間の次元をNとしておきます。初期状態$\rho_{ini}$を純粋状態とすると、周りの環境系との相互作用によって、最終状態$\rho_{fin}$は混合状態になります。先程説明した純粋化の手続きによると、混合状態の密度演算子のランクに等しい次元の環境系Eを用意すれば純粋状態にすることができるのでした。つまり、全体系A+Eの最終状態は、$\rho_{fin}$を対角化する基底$\{ \ket{\psi_k}^{A} \}$と環境系の適当な基底$\{ \ket{\mu_k}^{E} \}$を使って、

\ket{\Phi_{fin}}^{AE} = \sum_{k=1}^{d} \lambda_k \ket{\psi_k}^A \ket{\mu_k}^{E}  \tag{5}

という形の純粋状態にすることができます。一方で、初期状態は純粋状態なので、この基底$\{ \ket{\psi_k}^{A} \}$を使って、

\ket{\Phi_{ini}}^{AE} = \sum_{k=1}^{d} c_k \ket{\psi_k}^A \ket{0}^{E}  \tag{6}

と表現することができます²。N次元のヒルベルト空間上で定義されたN個以下の直交基底同士は、ユニタリ変換で互いに変換することができるので、式(6)から式(5)に変換するユニタリ変換は必ず存在します。そして、式(5)から得られる密度演算子を環境系について部分トレースすると、$\rho_{fin}$になることは簡単に示せます（式変形は省略）。ということで、環境系として必要になる次元はランクd、つまり高々N次元あれば十分ということになります。

ここまでの議論の前提は、初期状態が純粋状態の場合でした。一般には、混合状態の場合も考えないといけません。注目系Aの次元Nにおいて定義された混合状態に対して、これを純粋化するための参照系Rの次元Nが必要になります。さらに環境系Eの次元Nが必要になります（全体系はA+R+Eになる）。つまり、追加で必要な次元はNとNのテンソル積なので、Nの２乗次元です。これは運が悪かった最悪の場合の見積もりです。日頃の行いが良ければ（笑）、もっと少なくて済みます（嘘！）。

純粋化のシミュレーション

それでは、純粋化をシミュレータで再現してみます。初期状態としてランダムな密度演算子を用意し、純粋化の計算を実行した結果、本当にそれが純粋状態になっており、さらに追加したアンシラをトレースアウトすることで、元の混合状態にきちんと戻ることを確認してみます。全体のPythonコードは以下です。

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

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

MIN_DOUBLE = 0.000001

def random_vector(qubit_num):

    dim = 2**qubit_num
    vec_ini = np.array([0.0]*dim)
    vec_ini[0] = 1.0
    mat = unitary_group.rvs(dim)
    vec = np.dot(mat, vec_ini)

    return vec

def random_mixed_state(qubit_num, mixed_num):

    # random probabilities
    r = [random.random() for _ in range(mixed_num)]
    s = sum(r)
    prob_mixed = [r[i]/s for i in range(mixed_num)]

    # random quantum states
    vec_A = [random_vector(qubit_num) for _ in range(mixed_num)]
    qs_A = [QState(vector=vec_A[i]) for i in range(mixed_num)]

    # random density operator
    de_A = DensOp(qstate=qs_A, prob=prob_mixed)

    # free memory
    for i in range(mixed_num):
        qs_A[i].free()
    
    return de_A
    
def computational_basis(qubit_num):

    dim = 2**qubit_num
    basis = np.array([[1 if j==i else 0 for j in range(dim)] for i in range(dim)])

    return basis

def eigen_values_vectors(densop):

    matrix = densop.get_elm()
    eigvals, eigvecs = eigh(matrix, eigvals_only=False)
    eigvecs = np.conjugate(eigvecs.T)

    eigvals_list = [eigvals[i] for i in range(len(eigvals)) if abs(eigvals[i]) > MIN_DOUBLE]
    eigvecs_list = [eigvecs[i] for i in range(len(eigvecs)) if abs(eigvals[i]) > MIN_DOUBLE]
            
    eigvals = np.array(eigvals_list)
    eigvecs = np.array(eigvecs_list)

    return eigvals, eigvecs

def purification(de_A):

    # representation of the input state with orthogonal basis of system A
    coef,basis_A = eigen_values_vectors(de_A)
    rank = len(basis_A)

    # make computational basis of system R
    qnum_R = int(math.log2(rank))
    if 2**qnum_R != rank:
        qnum_R += 1
    basis_R = computational_basis(qnum_R)

    # orthogonal basis of system A+R
    basis_AR = [np.kron(basis_A[i],basis_R[i]) for i in range(rank)]

    # representation of the purified state with orthogonal basis of system A+R
    vec_AR = np.array([0]*len(basis_AR[0]))
    for i in range(rank):
        vec_AR = vec_AR + (math.sqrt(coef[i]) * basis_AR[i])

    qs_AR = QState(vector=vec_AR)
    de_AR = DensOp(qstate=[qs_AR], prob=[1.0])

    # free memory
    qs_AR.free()

    return de_AR

def difference(de_1, de_2):

    mat_1 = de_1.get_elm()
    mat_2 = de_2.get_elm()
    vec_1 = np.ndarray.flatten(mat_1)
    vec_2 = np.ndarray.flatten(mat_2)
    diff = np.linalg.norm(vec_1-vec_2)

    return diff

if __name__ == '__main__':

    qnum_A = 1
    mixed_num = 3
    
    print("=== system A (mixed state) ===")
    de_A = random_mixed_state(qnum_A,mixed_num)
    print("* density operator:")
    print(de_A)
    print("* square trace:", de_A.sqtrace())
    
    print("=== system A+R (pure state) ===")
    de_AR = purification(de_A)
    print("* square trace:", de_AR.sqtrace())

    print("=== system A (partial system of A+R) ===")
    de_AR_pat = de_AR.partial(id=list(range(qnum_A)))
    print("* density operator:")
    print(de_AR_pat)
    print("* square trace:", de_AR_pat.sqtrace())

    print("=== result ===")
    if difference(de_A, de_AR_pat) < MIN_DOUBLE:
        print("* OK!")
    else:
        print("* NG!")

    de_A.free()
    de_AR.free()
    de_AR_pat.free()

何をやっているか、順に説明します。

qnum_A = 1
mixed_num = 3

まず、ランダムな量子状態を生成するためのパラメータ設定です。系Aの量子ビット数qnum_Aと何個の純粋状態を混合させるかを表すmixed_numを設定します。

print("=== system A (mixed state) ===")
de_A = random_mixed_state(qnum_A,mixed_num)
print("* density operator:")
print(de_A)
print("* square trace:", de_A.sqtrace())

関数random_mixed_stateでランダムな密度演算子を作成します。関数定義の中身を見てください。

# random probabilities
r = [random.random() for _ in range(mixed_num)]
s = sum(r)
prob_mixed = [r[i]/s for i in range(mixed_num)]

# random quantum states
vec_A = [random_vector(qubit_num) for _ in range(mixed_num)]
qs_A = [QState(vector=vec_A[i]) for i in range(mixed_num)]

# random density operator
de_A = DensOp(qstate=qs_A, prob=prob_mixed)

まず、混合割合（確率）を示すprob_mixedというリストを作成しまず。次に、考えている系Aにおけるランダムな純粋状態のリストを作成します。これは前回の記事で説明した関数random_vectorを使います。中でscipyのランダムユニタリ行列を生成する関数を召喚しています（説明省略）。最後に、得られた純粋状態のリストと混合割合のリストを使って、密度演算子de_Aを作成します。

main部の説明に戻ります。このde_Aの要素をprintで表示して、さらに２乗トレースを表示します。次の処理は、

print("=== system A+R (pure state) ===")
de_AR = purification(de_A)
print("* square trace:", de_AR.sqtrace())

です。関数purificationで、混合状態de_Aを純粋化し、その２乗トレースを表示します。ここで正常に純粋化がなされていれば１と表示されるはずです。関数の中身を見てみます。ここが、本日の本題となる計算部分になります。

# representation of the input state with orthogonal basis of system A
coef,basis_A = eigen_values_vectors(de_A)
rank = len(basis_A)

関数eigen_values_vectorsで、de_Aをスペクトル分解します。coefがその係数（固有値）を表す１次元配列で、basis_Aがその基底を表す２次元配列（固有ベクトルの配列）になります。eigen_values_vectorsの中でscipyの固有値問題を解く関数を召喚しています。詳細は関数定義を見てください（説明省略）。固有ベクトルの数がランクになるので、それを変数rankに格納します。

# make computational basis of system R
qnum_R = int(math.log2(rank))
if 2**qnum_R != rank:
    qnum_R += 1
basis_R = computational_basis(qnum_R)

参照系Rとして必要な量子ビットは上で得たランクのlog2にしておけば良いので、それを計算しています。ただし、ランクが２のべきになっていない場合、切り捨てで整数化された量子ビット数ではランクの次元を表せないため、１つ増やしています。この量子ビット数qnum_Rを使って、関数computational_basisで計算基底を表す２次元配列（基底ベクトルの配列）作成します（関数定義参照、説明省略）。先程説明した通り、これは何でも良いので、一番簡単な計算基底としました。

# orthogonal basis of system A+R
basis_AR = [np.kron(basis_A[i],basis_R[i]) for i in range(rank)]

注目系Aの基底と参照系Rの基底が得られたため、式(5)の形を作るための準備として、A+Rの基底を作成します。両者のテンソル積³をとれば良いので、それを計算します。numpyにクロネッカー積を計算する関数kronがあるので、これを使います（同じことなので）。

# representation of the purified state with orthogonal basis of system A+R
vec_AR = np.array([0]*len(basis_AR[0]))
for i in range(rank):
    vec_AR = vec_AR + (math.sqrt(coef[i]) * basis_AR[i])

式(5)の純粋状態を表すベクトルvec_ARを計算します。

qs_AR = QState(vector=vec_AR)
de_AR = DensOp(qstate=[qs_AR], prob=[1.0])

vec_ARに基づき、量子状態(QState)のインスタンスqs_ARを生成し、qs_ARに基づき、密度演算子(DensOp)のインスタンスを生成します。これで純粋化ができました。

main部に戻ります。

print("=== system A (partial system of A+R) ===")
de_AR_pat = de_AR.partial(id=list(range(qnum_A)))
print("* density operator:")
print(de_AR_pat)
print("* square trace:", de_AR_pat.sqtrace())

純粋状態から部分系Aの密度演算子を取り出します。DensOpクラスのpartialメソッドに量子ビット番号のリストを引数として与えればできます（中で部分トレースをやっています）。結果を変数de_AR_patに格納します。de_AR_patの要素をprintで表示し、さらに２乗トレースも表示します。

print("=== result ===")
if difference(de_A, de_AR_pat) < MIN_DOUBLE:
    print("* OK!")
else:
    print("* NG!")

最後に、初期密度演算子de_Aと純粋化後に部分トレースをやった結果の密度演算子de_AR_patが等しいことを確認します。関数differenceでは、両者の要素間の差の絶対値の総和を計算しています（関数定義参照、説明省略）⁴。両者一致していれば「OK!」と表示されます。

それでは、実行結果を示します。

=== system A (mixed state) ===
* density operator:
[[ 0.4835702+0.j         -0.1852386+0.20585944j]
 [-0.1852386-0.20585944j  0.5164298+0.j        ]]
* square trace: 0.65392277
=== system A+R (pure state) ===
* square trace: 1.0
=== system A (partial system of A+R) ===
* density operator:
[[ 0.4835702+0.j         -0.1852386+0.20585944j]
 [-0.1852386-0.20585944j  0.5164298+0.j        ]]
* square trace: 0.65392277
=== result ===
* OK!

というわけで、混合状態を適当に用意し純粋化の計算を実行した結果、２乗トレースが１になり、確かに純粋状態になることが示せました。さらに、追加した参照系をトレースアウトすることで、元の混合状態に戻ることもわかりました。量子ビット数qnum_Aと何個の純粋状態を混合させるかを表すmixed_numをいろいろ変化させて何度も実行してみましたが、すべて同様の結果でした。

特に驚くべきこともなく、予想通りの結果ですが、とりあえず自分の理解確認とシミュレータのデバッグも兼ねて（汗）やってみました。

おわりに

量子情報理論の基本の基として取り上げていない話題がまだあります。「状態間の距離・忠実度」とか「エントロピー」あたりは、最低限おさえておきたい話題なので、今後順に勉強していく予定です。これらをベースに、さらに「量子誤り訂正」「量子暗号」といった大きなお山の麓あたりをうろうろしながら、基本事項のアレコレをつまみ食いしてみたいと、今のところ思っています。引き続き、マイペース＆気まぐれに記事を書いていくので、予定変更するかもしれませんが、その際はご容赦を。

以上

この意味で、Rのことを「参照系」と言ったりするようです。 ↩
環境系は何でも良いのですが、$\ket{0}^{E}$をもってきました。 ↩
ここは間違えないように。２つのベクトルの要素を繋げれば良いわけではありませんので、念のため。 ↩
密度演算子同士が一致しているかどうかを見たかっただけなので簡易的に（勝手に）距離をこのように定義しました。密度演算子同士を比較するための指標として「トレース距離」とか「忠実度」という由緒正しいものがあったりしますが、これについては別途取り上げます。 ↩