量子情報理論の基本：シュミット分解 #Python

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

はじめに

量子情報理論の分野においては、環境系と注目系のように複合した全体系の中で、各々の系がどう変化していくとか、どんな性質をもっているとかいう観点で、いろいろな研究がなされてきました。そういった研究を進める中で、有用なツール・知見の蓄積がなされてきたわけですが、今回は、その代表選手の一つである「シュミット分解(Schmidt decomposition)」について勉強してみます。だいたいわかってきたところで、量子計算シミュレータqlazyを使って、そこから導き出される性質を確認しつつ、その面白さを噛み締めてみたいと思います。

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

シュミット分解とは

部分系A,Bからなる全体系ABにおける、ある純粋状態を$\ket{\Phi}^{AB}$とすると、

\ket{\Phi}^{AB} = \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}}^{A} \ket{\psi_{k}}^{B}  \tag{1}

を満たす、負でない実数$\{ \lambda_{k} \}$、および系Aと系Bにおける正規直交基底$\{ \ket{\phi_k}^{A} \}$、$\{ \ket{\psi_{k}}^{B} \}$が必ず存在します。この表現のことを「シュミット分解」と言い、$\{ \lambda_{k}\}$のことを「シュミット係数(Schmidt coefficient)」、$\{ \lambda_{k}\}$の中のゼロでないものの個数を「シュミット・ランク(Schmidt rank)」、正規直交基底$\{\ket{\phi_{k}}^{A}\}$、$\{\ket{\psi_{k}}^{B}\}$のことを、系A、系Bに対する「シュミット基底(Schmidt basis)」と言います。

ここで、$\ket{\Phi}^{AB}$の正規化条件から、

\sum_{k=1}^{R} \lambda_k^{2} = 1  \tag{2}

が成り立ちます。

それでは、任意の$\ket{\Phi}^{AB}$が式(1)のように分解できることを証明してみます。

【証明】

$\ket{\Phi}^{AB}$は、ABの純粋状態なので、系Aの正規直交基底$\{ \ket{i}^A\}$と系Bの正規直交基底$\{ \ket{j}^B\}$を使って、一般に以下のように書けます¹。

\ket{\Phi}^{AB} = \sum_{i=1}^{M} \sum_{j=1}^{N} c_{ij} \ket{i}^{A} \ket{j}^{B}  \tag{3}

ここから出発して、式(1)を導ければ証明完了なのですが、それには特異値分解を使います。

いま、式(1)の$c_{ij}$をM行N列の行列$C$の($i,j$)要素とみなすと、M行M列のユニタリ行列$U$、N行N列のユニタリ行列$V$、M行N列の対角行列$D$を用いて、

C = UDV  \tag{4}

のように特異値分解できます²。つまり、

c_{ij} = \sum_{k=1}^{R} u_{ik} d_{kk} v_{kj}  \tag{5}

です。ここで、$u_{ik}, \space d_{kk}, \space v_{kj}$は各々$U,D,V$の要素です。また、$R$は、$D$の対角成分に並ぶゼロでない成分の個数を表します。

式(5)を式(3)に代入します。

\ket{\Phi}^{AB} = \sum_{i=1}^{M} \sum_{j=1}^{N} \sum_{k=1}^{R} u_{ik} d_{kk} v_{kj} \ket{i}^{A} \ket{j}^B  \tag{6}

ここで、

\begin{align}
\ket{\phi_{k}}^{A} &= \sum_{i=1}^{M} u_{ik} \ket{i}^{A} \\
\ket{\psi_{k}}^{B} &= \sum_{j=1}^{N} v_{kj} \ket{j}^{A} \tag{7}
\end{align}

とおくと、

\ket{\Phi}^{AB} = \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}}^{A} \ket{\psi_{k}}^{B}  \tag{1}

となり、式(1)が導けました。また、式(7)より、$\{ \ket{\phi_{k}}^{A} \}$および$\{\ket{\psi_{k}}^{B}\}$は正規直交基底に対するユニタリ変換なので、両方とも正規直交基底です。（証明終）

さて、ここで、出発点とした式(3)とシュミット分解の式(1)を改めて眺めてください。証明の過程で何をやったかというと、要はユニタリ変換で部分系Aと部分系Bの基底を変換しているだけです。どのように変換したかというと、考えるべき基底の数をなるべく少なくなるようにした、ということになります³。と同時に、部分系Aと部分系Bの間のエンタングルメント（もつれ）が最大に見えるように変換したとも言えます。とすると、シュミットランク$R$というのは、ある意味、エンタングルメントの大きさを表す指標と考えることができます。シュミットランクが１というのは、エンタングルメントがない状況、つまり、部分系AとBが積状態（＝テンソル積として分割できる状態）になっているということを表しており、シュミットランクが２以上というのは、両者エンタングルしている状況ということです。

部分系の性質

全体系のある部分系をA、残りの部分系をBとすると、全体系で定義されたどんな純粋状態も、式(1)のようにシュミット分解の表式に書き直すことができる、ということがわかりました。この知見を利用すると、部分系に関して以下の性質が導けます。

密度演算子の固有値

全体系の密度演算子$\rho_{AB}=\ket{\Phi}^{AB} \bra{\Phi}^{AB}$の部分トレースをとると、部分系の密度演算子$\rho_{A}, \rho_{B}$は、

\begin{align}
\rho_{A} &= Tr_{B} (\ket{\Phi}^{AB} \bra{\Phi}^{AB}) = \sum_{k=1}^{R} \lambda_{k}^{2} \ket{\phi_{k}}^{A} \bra{\phi_{k}}^{A}  \\
\rho_{B} &= Tr_{A} (\ket{\Phi}^{AB} \bra{\Phi}^{AB}) = \sum_{k=1}^{R} \lambda_{k}^{2} \ket{\psi_{k}}^{B} \bra{\psi_{k}}^{B}  \tag{8}
\end{align}

となります。したがって、$\rho_{A}$の固有値も$\rho_{B}$の固有値も$\{ \lambda_{k}^2 \}$で、同じということになります。また、固有値（＝ユニタリ変換で対角化したときの対角成分）が同じということは、当然ランクも同じで、さらに、その２乗トレースも等しいということになります。つまり、

Tr(\rho_{A}^{2}) = Tr(\rho_{B}^{2})  \tag{9}

です。

このことは、全体系をどのように分割したとしても成り立つわけで、とても面白い性質だと思います。どっちの部分系で固有値を計算したとしても全く同じになるということなので、この固有値というのは、何かAとBとの間の分割面というか、関係性を特徴づけるような指標に思えてきます。そういえば、この固有値（の平方根）はシュミット係数そのものでした。つまり、その関係性というのは、要はエンタングルメントのことではないでしょうか、ということが、じんわりとわかってきます。

また、ある純粋状態に対するシュミット係数（あるいはシュッミットランク）を知りたいとすると、想定している部分系（どっちの部分系でも良い）で部分トレースをとって、その密度演算子の固有値を求めれば良い、ということもわかってきます（一つの計算方法として）。

ユニタリの自由度

部分系Aで定義される局所的なユニタリ変換を$U_{A}$、部分系Bで定義される局所的なユニタリ変換を$U_{B}$とすると、それによって全体系の純粋状態$\ket{\Phi}^{AB}$は、

\begin{align}
(U_{A} \otimes U_{B}) \ket{\Phi}^{AB} &= (U_{A} \otimes U_{B}) \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}}^{A} \ket{\psi_{k}}^{B} \\
&= \sum_{k=1}^{R} \lambda_{k} (U_{A} \ket{\phi_{k}}^{A}) (U_{B} \ket{\psi_{k}}^{B}) \\
&= \sum_{k=1}^{R} \lambda_{k} \ket{\phi_{k}^{\prime}}^{A} \ket{\psi_{k}^{\prime}}^{B} \tag{10}
\end{align}

のように変換されます。ここで、

\begin{align}
\ket{\phi_{k}^{\prime}}^{A} &= U_{A} \ket{\phi_{k}}^{A} \\
\ket{\psi_{k}^{\prime}}^{B} &= U_{B} \ket{\psi_{k}}^{B} \tag{11}
\end{align}

とおきました。

これより、局所的なユニタリ変換で、シュミット分解の性質（シュミット係数やシュミットランク）は不変であるということがわかります。つまり、部分系AとBの間に存在するエンタングルメントの性質は、局所ユニタリ変換で不変に保たれるということです。こういう不変性は、いろいろ理論計算をする上で、何かと便利に使えるものなので、覚えておくと良いと思います。

シミュレータで確認

それでは、上で述べてきた知見の一部をシミュレータで確認してみたいと思います。全体系の任意の純粋状態に対する部分系の密度演算子の固有値が、全体系をどういう形に２分割したとしても、同じということだったので、本当にそうなるかを確認します。

全体のPythonコードは以下です。

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

import random
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_qstate(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)
    qs = QState(vector=vec)

    return qs

def random_qubit_id(qubit_num):

    id = list(range(qubit_num))
    random.shuffle(id)
    qubit_num_A = random.randint(1,qubit_num-1)
    qubit_num_B = qubit_num - qubit_num_A
    id_A = id[:qubit_num_A]
    id_B = id[qubit_num_A:]

    return id_A,id_B

def eigen_values(densop):

    matrix = densop.get_elm()
    eigvals = eigh(matrix, eigvals_only=True)
    eigvals_out = [eigvals[i] for i in range(len(eigvals)) if eigvals[i] > MIN_DOUBLE]

    return eigvals_out

if __name__ == '__main__':

    # whole quantum state
    qubit_num = 5
    qs = random_qstate(qubit_num)
    de = DensOp(qstate=[qs], prob=[1.0])  # pure state

    # partial density operators (system A and B)
    id_A, id_B = random_qubit_id(qubit_num)
    de_A = de.patrace(id=id_B)
    de_B = de.patrace(id=id_A)

    # eigen-values of density operators (system A and B)
    eval_A = eigen_values(de_A)
    eval_B = eigen_values(de_B)

    print("== system A ==")
    print("- qubit id =", id_A)
    print("- square trace = ", de_A.sqtrace())
    print("- eigen values =", eval_A)
    print("- rank =", len(eval_A))

    print("== system B ==")
    print("- qubit id =", id_B)
    print("- square trace = ", de_B.sqtrace())
    print("- eigen values =", eval_B)
    print("- rank =", len(eval_B))

    qs.free()
    de.free()
    de_A.free()
    de_B.free()

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

# whole quantum state
qubit_num = 5
qs = random_qstate(qubit_num)
de = DensOp(qstate=[qs], prob=[1.0])  # pure state

全体系の量子ビット数を5にして、量子状態をランダムに設定します。関数random_qstateの中で実行しているので、その中身を見てみます。

def random_qstate(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)
    qs = QState(vector=vec)

    return qs

指定された量子ビット数に対応した状態の次元数をdimとして、numpyの配列（ベクトル）を適当に決めます。第0成分のみ1.0、それ以外を0.0としました。ここで、ランダムな量子状態にするために、scipyの関数を召喚します。

from scipy.stats import unitary_group
...

mat = unitary_group.rvs(dim)

です。これは、ランダムなユニタリ行列を生成するための関数です。ここにドキュメントがあります⁴。生成されたユニタリ行列をvec_iniに適用して、できたランダムベクトルに基づき、量子状態QStateクラスのインスタンスを作り、リターンします。

もとのmain部に戻ります。できあがったランダム量子状態を用いて、密度演算子DensOpクラスのインスタンスを生成します。これで、全体系に対する純粋状態の密度演算子ができあがりました。

# partial density operators (system A and B)
id_A, id_B = random_qubit_id(qubit_num)
de_A = de.patrace(id=id_B)
de_B = de.patrace(id=id_A)

全体系の量子ビット（5量子ビット）を部分系Aと部分系Bにランダムに分割して、各々の量子ビット番号リストを変数id_A,id_Bに格納しています。関数random_qubit_idで実行しています。説明は省きます。関数定義を見てください。この量子ビット番号リストを使って、全体系の密度演算子の部分トレースをとります。結果を各々、変数de_A,de_Bに格納します。

# eigen-values of density operators (system A and B)
eval_A = eigen_values(de_A)
eval_B = eigen_values(de_B)

固有値問題を解いて固有値のリストを求めます。関数eigen_valuesの中で、scipyの関数を使って実行しています。eighはエルミート行列に対する固有値を求める関数です。

from scipy.linalg import eigh
 ...

def eigen_values(densop):

    matrix = densop.get_elm()
    eigvals = eigh(matrix, eigvals_only=True)
    eigvals_out = [eigvals[i] for i in range(len(eigvals)) if eigvals[i] > MIN_DOUBLE]

    return eigvals_out

eigen_values関数の最初の行で、密度演算子DensOpクラスのget_elm()メソッドを使って行列表現（numpyの２次元配列）を取得しています。それを使って、固有値問題を解いています。

メイン部に戻ります。

print("== system A ==")
print("- qubit id =", id_A)
print("- square trace = ", de_A.sqtrace())
print("- eigen values =", eval_A)
print("- rank =", len(eval_A))
 ...

部分系Aの量子ビット番号、２乗トレース、固有値、ランクを表示します。部分系Bも同様に表示します。

さて、実行結果です。以下のようになりました。

== system A ==
- qubit id = [4, 1, 2]
- square trace =  0.38437171
- eigen values = [0.0822181470819045, 0.1464507832125754, 0.2143642754307335, 0.5569667940374802]
- rank = 4
== system B ==
- qubit id = [0, 3]
- square trace =  0.38437171
- eigen values = [0.08221814299667664, 0.14645078689065036, 0.2143642733317681, 0.5569667867809049]
- rank = 4

部分系Aに3量子ビット、部分系Bに2量子ビットが各々割り当てられました。ランクは4で固有値を見ていただくとわかる通り、きちんと一致しました！⁵また、２乗トレースは１以下なので、言わずもがなですが混合状態になっています。

もう一回やってみます。

== system A ==
- qubit id = [2, 3, 0, 4]
- square trace =  0.53145101
- eigen values = [0.3745986253398182, 0.6254013753462571]
- rank = 2
== system B ==
- qubit id = [1]
- square trace =  0.53145101
- eigen values = [0.37459862198601285, 0.6254013780139871]
- rank = 2

というわけで、今度は、4量子ビットと1量子ビットに分割されました。ランクは2となり、この場合も固有値は一致しました！全体系の量子ビット数も変えつつ、何度も実行してみましたが、すべての場合で固有値は一致しました！

これは、それほど自明なことではないので、面白がっていただけると幸いです。純粋状態の密度演算子があって、それを２つの部分系に分けたとすると、必ず成り立つ性質です。

おわりに

今回は、純粋状態を部分系にわけるお話でした。その部分系は一般に混合状態になります。では、逆に混合状態に何らかの部分系を追加することで、純粋状態にできるでしょうか。答えはイエスで、それを実行する手続きのことを「純粋化」と言います。今回の記事で一緒に説明しようと思ったのですが、長くなりました。というわけで、次回は「純粋化」の予定です。

以上

これは大丈夫ですよね。任意の状態は適当な基底を使った重ね合わせで表現できます、と言っているだけです。 ↩
特異値分解は、通常$C = UDV^{\dagger}$と書かれたりします。添字を使って式展開するときに、ダガーがあると複素共役とか添字の順番を逆にしないと、、とかちょっと余計な注意をしないといけないので、$V^{\dagger} \rightarrow V$と置き換えました。 ↩
データサンプルの分布特徴を利用して、大量の説明変数をなるべく減らすという統計分析の世界でよく使われる「主成分分析」に似ています。 ↩
中で何をやっているか理解できていませんが（汗）、Haar測度に基づく由緒正しいランダムユニタリ行列を計算しているようです。今回は単に適当な量子状態を作りたかっただけなので、ここまでちゃんとしたランダム性は必要ないですが、ランダムな量子状態をどうやったら作れるんだろうとちょっと調べていたら、ひっかかってきたので試しに使ってみました。採用している手法は、以下の論文に記載されているようです。参考まで。How to generate random matrices from the classical compact groups ↩
小数点以下７桁以降は微妙に違っていたりしますが、計算誤差なので無視します ↩