Classical shadow -古典情報から量子状態を復元する-

ここでは、最近提案された　classical shadow　と呼ばれる量子アルゴリズムを紹介します。
名前がカッコイイですね！

参考文献は以下。

What is classical shadow?

量子ビット数$N$の量子状態は、$(2^{N})^{2}=4^{N}$個の成分を持つ密度行列$\rho$により記述できます。
https://qiita.com/SamN/items/ecbae603041317511969
(状態ベクトル表示は密度行列表示の特別な場合とみなせます。)
密度行列の成分は、各量子ビットに対する物理量${X,Y,Z,I}$の期待値からなります。
4通りの物理量の$N$量子ビットに対する組み合わせなので$4^{N}$通りあるわけです。
期待値である以上は複数回の測定が必要になります。
量子状態は一般には測定することで破壊されてしまうので、まったく同じ量子状態$\rho$を何度も生成して何度も（壊しながら）測ることが必要です。
量子状態を知ることは、膨大な測定回数が必要そうですね。
さらにまずいことに、密度行列の成分というのは量子ビット数$N$に対して指数的に多くあるのでした($4^{N}個$)。
行列の成分が指数的に多くあり、かつ1つの成分を知るのですら複数回の量子状態生成→測定（破壊）が必要となると、ぞっとしますね。

本題のclassical shadow は、この”膨大な回数の測定”をなんとか回避して、
限られた数の測定で量子状態を（ある精度でよいから）推定できないか？　というコンセプトです。
実は結果的にはこのアルゴリズムはうまくいっていないのですが、、アイデアが面白いので紹介しましょう。

classical shadow　のアプローチ

ある量子回路の出力の量子状態、即ち密度行列$\rho$を推定したいとします。
classical shadow　では、密度行列の成分を1つずつ測定するのを諦めます。
その代わりに、密度行列の成分をランダムに1つ選択してかつ1度だけ測定(shots=1)します。
密度行列の成分をランダムに1つ選択することは、各量子ビットに対して測定する物理量$O_{rand.} \in {X,Y,Z,I}$の割り当てを決めることと同じです。
${X,Y,Z,I}$測定は、$Z$測定の前に対応するゲート操作$U_{rand.}$を挟むこととも同じです。
(量子コンピュータ実機では$Z$測定しかできないので、このような置き換えをよく使います)

得られた古典的な測定結果$b \in [0,1]^{N}$は、量子論の性質から、その量子状態の射影になっています。
これが"classical shadow"の由来です。
得られた射影$b$を、古典コンピュータ上で”逆射影”します。
もちろん射影自体が不可逆な操作なので、この”逆射影”は量子状態を再現してはくれないのですが、元の量子状態（の密度行列)に対する部分的な情報は含まれていると期待します。
”逆射影”されたものは密度行列と等しい次元を持ち、「密度行列のスナップショット $\hat{\rho}$」と呼ばれます。
この「密度行列のスナップショット $\hat{\rho}$」を得る操作を、最初の「密度行列の成分をランダムに1つ選択」をしなおして、繰り返します。
すると密度行列のスナップショットの集まり{$\hat{\rho_{1}}$,$\hat{\rho_{2}}$,...}ができます。
これをclassical shadowと呼びます。

classical shadow、すなわち「密度行列のスナップショットの集まり」は、それぞれは推定したい密度行列の部分的な情報(特に、選択した密度行列の成分周辺の情報）しか持ちませんが、これを寄せ集めて平均化してやれば、なんとなく張り合わされてもともとの密度行列が浮かび上がってきそうに思いませんか？
$\rho \approx mean[{\hat{\rho_{1}},\hat{\rho_{2}},...}]$

これがclassical shadow　の基本的なアプローチです。

しかしclassical shadow が厳密な方法に比べて効率的である理由はいまいち直感的に分かりません。
実際のところ、classical shadow における”張り合わせ”は、「スナップショットの集まりに含まれる要素のうち、たまたま欲しい条件に近いものを抜き出して使い、他は捨てているだけと同じである」ということも指摘されています。
classical shadow が実際に測定回数を減らすことができるケースは極めて限定的です。
https://pennylane.ai/qml/demos/tutorial_diffable_shadows.html

とはいえ、古典情報から量子状態の部分的な情報のアンサンブルを作って、平均化することで量子状態を復元するというアイデア自体は面白いです。

実装(Qiskit ver.)

classical shadow の実装は数式を見ていただくのが一番です。
https://pennylane.ai/qml/demos/tutorial_classical_shadows.html

ここでは、上記のpennylane SDK による実装をQiskit SDKに置き換えてみたものを掲載します。
乱択　→　射影測定　→　逆発展　→　密度行列スナップショットを取得　→　最後に平均
という流れを感じてください。
以下では2量子ビットのBell状態の推定をしています。
"shadow size"がスナップショットの個数です。

#必要なモジュールのインポート
from qiskit import IBMQ, QuantumCircuit
from qiskit import execute, Aer
from qiskit.providers.aer import QasmSimulator
import numpy as np

def classical_shadow(shadow_size=100):
    random_paulis = np.random.randint(0,3,(shadow_size,2)) # (shadow_size,num_of_qubits). 0~2 mean pauli X,Y and Z measurements.
    random_unitaries = np.zeros((shadow_size,2,2,2),dtype="complex") # (shadow_size,num_of_qubits). 0~2 mean pauli X,Y and Z measurements.

    gate_h = 1/np.sqrt(2)*np.array([[1,1],[1,-1]])
    gate_p = np.array([[1,0],[0,-1j]])
    gate_i = np.array([[1,0],[0,1]])
    for i_shadow in range(shadow_size):
        i_paulis = random_paulis[i_shadow]
        for i_qubit in range(2):
            if i_paulis[i_qubit]==0:
                random_unitaries[i_shadow][i_qubit] = gate_h # X basis
            elif i_paulis[i_qubit]==1:
                random_unitaries[i_shadow][i_qubit] = gate_h@gate_p # Y basis
            elif i_paulis[i_qubit]==2:
                random_unitaries[i_shadow][i_qubit] = gate_i # Z basis
            else:
                print("error")

    rho_all = np.zeros((4,4),dtype="complex")
    for i_shadow in range(shadow_size):
        sample_unitary = random_unitaries[i_shadow]
        # bell state preparation
        circuit = QuantumCircuit(2) 
        circuit.h(0)
        circuit.cx(0,1) #Controlled-NOT. The controlled and target bit are 0th and 1st one, respectively.
        for i in range(2):
            circuit.unitary(sample_unitary[i],[i]) # Rotate to X or Y or Z basis measurement
        circuit.measure_all()
        circuit.draw('mpl')
        # execute on qasm simulator for a single shot
        results = execute(circuit, backend=Aer.get_backend('qasm_simulator'), shots=1).result()
        # from sample to classical vector
        binary = list(results.data()["counts"].keys())[0] # measure bits
        if binary=='0x0':
            classical_vec0 = [1,0] #0
            classical_vec1 = [1,0] #0
        elif binary=='0x1':
            classical_vec0 = [0,1] #1
            classical_vec1 = [1,0] #0
        elif binary=='0x2':
            classical_vec0 = [1,0] #0
            classical_vec1 = [0,1] #1
        elif binary=='0x3':
            classical_vec0 = [0,1] #1
            classical_vec1 = [0,1] #1
        classical_vec0 = np.array(classical_vec0).reshape([2,1]) # reconstructed projected state vector of 0-th qubit
        classical_vec1 = np.array(classical_vec1).reshape([2,1])

        classical_dm0 = classical_vec0*classical_vec0.conj().transpose() # reconstructed projected density matrix of 0-th qubit
        classical_dm1 = classical_vec1*classical_vec1.conj().transpose()

        rho0 = 3*sample_unitary[0].conj().transpose()@classical_dm0@sample_unitary[0]-gate_i # back propagated sample_unitary[1] matrix of 0-th qubit
        rho1 = 3*sample_unitary[1].conj().transpose()@classical_dm1@sample_unitary[1]-gate_i
        rho = np.kron(rho1,rho0) 
        rho_all += rho
    rho_all = rho_all/shadow_size
    return rho_all

推定したい密度行列が本当の密度行列とどれぐらい”近い”のかを評価するために、
行列間の距離関数を定義します。

ground_truth = np.array([[0.5,0,0,0.5],[0,0,0,0],[0,0,0,0],[0.5,0,0,0.5]])
def matrix_distance(ground_truth,rho):
    diff = ground_truth-rho
    distance = np.sqrt(np.trace(diff.conj().transpose()@diff))
    return np.abs(distance)

では、密度行列スナップショットの個数を変えながら、classical shadowによる推定を行います。個数に対して推定エラー（行列距離）がどうなるか見ます。

dist = []
shadow_size_list = [1,10,20,30,50,100,200,400,1000]
for _shadow_size in shadow_size_list:
    _rho = classical_shadow(_shadow_size)
    _dist = matrix_distance(ground_truth,_rho)
    dist.append(_dist)

plt.plot(shadow_size_list,dist,'.-')
plt.xlabel("shadow_size",fontsize=16)
plt.ylabel("estimation error",fontsize=16)

スナップショットの個数を十分とれば、ある程度の精度で密度行列が復元できています！

実際に、密度行列の成分（は複素数なので、その絶対値）を二次元プロットして見せます。

classical shadow　で推定されたものは、

import matplotlib.pyplot as plt
plt.imshow(np.abs(_rho),cmap="gray")

本当の密度行列は、

plt.imshow(np.abs(ground_truth),cmap="gray")

ということで、確かに一致してきていることがわかります。
classical shadowにより、量子系の部分的な情報（古典情報）だけから量子状態が推定できたのです。
なかなか面白いアルゴリズムでしたね。