I'm trying to reproduce the basic method of classical shadow, which is based on the tutorial of pennylane. However, I've met some realization problems here when I finish reading the tutorial of pennylane, and trying to finish the method myself, just to check if my understanding is correct, because I'm not sure if my understanding about the inverse of the map $M$ mentioned in the paper is correct.
I will describe the method shortly first, and then showing my code with qiskit, pennylane, and matlab. But all of them failed to have the same effect as the tutorial of pennylane does(increasing the number of measurements, the distance between the state I want to reconstruct and the original state should be more and more close).
The idea of classical shadow(or the process of the algorithm) is kind of simple while the math behind it might be complicated. The process states that for any density matrix $\rho$, we act some unitary matrix $U$ which are chosen randomly from a specific set of the unitary matrix $\mathcal{U}$ on it, i.e., $U\rho U^\dagger$. Then we do one-shot measurement based on the computational basis on $U\rho U^\dagger$. Then the state will collapse into some state $|\hat{b}\rangle$. And then we do the rest of the work in classical data analysis style. First we undo the unitary matrix, i.e., $U^\dagger |\hat{b}\rangle\langle \hat{b}| U$ . Then we do the inverse of the map $\hat{\rho}\equiv M^{-1}(U^\dagger |\hat{b}\rangle\langle \hat{b}| U)$ which can be defined as $M(\rho)\equiv E(U^\dagger |\hat{b}\rangle\langle \hat{b}| U)$, where the $E$ stands for expectation over both unitary matrix and measurement result $|\hat{b}\rangle$. And for a specific choice of unitary set(seems Clifford group, not very clear here), we have a form of the inverse of $M$ states as $$ \hat{\rho} = \bigotimes_{j=1}^n(3U^{\dagger}_j|\hat{b}_j\rangle\langle\hat{b}_j|U_j-\mathbb{I})\tag{1} $$ Then I will introduce my code first. We specify the unitary group into Hadamard gate, phase gate, and identity. Then we do the computational basis measurement and rebuild the classical shadow $\hat{\rho}$ with the help of eq.(1), and then calculate the expectation value of $\hat{\rho}$ by directly divide measurement times.
Following the codes(pennylane, qiskit, Matlab), aiming at construct the classical shadow of bell state:
from networkx.algorithms.centrality import harmonic
from networkx.exception import HasACycle
from networkx.readwrite.sparse6 import write_sparse6
from numpy import dtype
from numpy.random.mtrand import rand
import pennylane as qml
from pennylane import wires
import pennylane.numpy as np
import matplotlib.pyplot as plt
import time
def distance(rho):
return np.sqrt(np.trace(rho.conjugate().transpose() @ rho))
def my_quantum_function(x, y):
unitary = [qml.Hadamard, qml.S, qml.Identity]
qml.Hadamard(wires=0)
qml.CNOT(wires=[0,1])
unitary[y[0]](wires=0)
unitary[y[1]](wires=1)
# all measure in computational basis, i.e., mean value of pauliz, one-shot case
return [qml.expval(qml.PauliZ(0)), qml.expval(qml.PauliZ(1))]
# one-shot case shots = 1 to simulate the measure in computational basis requirement
dev = qml.device('default.qubit', wires=2, shots=1)
circuit = qml.QNode(my_quantum_function, dev)
# generate random number seed for easy replicate the experiment
np.random.seed(666)
# init
phase_z = np.array([[1, 0], [0, 1j]], dtype=complex)
hadamard = qml.Hadamard(0).matrix
identity = qml.Identity(0).matrix
unitary = [hadamard, phase_z, identity]
snapshot = 1000
state0 = np.array([[1,0],[0,0]])
state1 = np.array([[0,0],[0,1]])
record_rho = np.zeros([4,4])
for i in range(snapshot):
randnum = np.random.randint(0,3,size=2)
[res0, res1] = circuit(0,randnum)
# print(circuit.draw())
if res0 == 1:
rho1 = 3*(unitary[randnum[0]].conj().T @ state0 @ unitary[randnum[0]]) - identity
else:
rho1 = 3*(unitary[randnum[0]].conj().T @ state1 @ unitary[randnum[0]]) - identity
if res0 == 1:
rho2 = 3*(unitary[randnum[1]].conj().T @ state0 @ unitary[randnum[1]]) - identity
else:
rho2 = 3*(unitary[randnum[1]].conj().T @ state1 @ unitary[randnum[1]]) - identity
record_rho = record_rho + np.kron(rho1,rho2)
record_rho = record_rho/snapshot
bell_state = np.array([[0.5, 0, 0, 0.5], [0, 0, 0, 0], [0, 0, 0, 0], [0.5, 0, 0, 0.5]])
print(record_rho)
print(distance(record_rho - bell_state))
from math import exp
from qiskit import *
from qiskit import Aer
import numpy as np
import matplotlib.pyplot as plt
from random import randrange
np.random.seed(222)
def one_shot(operator):
sim = Aer.get_backend('aer_simulator')
qc = QuantumCircuit(2)
unitary = [qc.h, qc.sdg, qc.id]
qc.h(0)
qc.cx(0,1)
unitary[operator[0]](0)
unitary[operator[1]](1)
qc.measure_all()
qobj = assemble(qc,shots=1)
result = sim.run(qobj).result().get_counts()
return result
def distance(rho):
'''
calculate distance of two density matrix
'''
return np.sqrt(np.trace(rho.conjugate().transpose().dot(rho)))
hadamard = 1/np.sqrt(2)*np.array([[1,1],[1,-1]])
s_gate = np.array([[1,0],[0,-1j]],dtype=complex)
id = np.identity(2)
unitary = [hadamard,np.dot(hadamard,s_gate),id]
snapshot_num = 1000
state0 = np.array([[1,0],[0,0]])
state1 = np.array([[0,0],[0,1]])
record_rho = np.zeros([4,4])
for i in range(snapshot_num):
randnum = np.random.randint(0,3,size=2)
result = one_shot(randnum)
if result.get('00') == 1:
rho = np.kron(3*np.dot(unitary[randnum[0]].conj().T,state0).dot(unitary[randnum[0]] - id),3*np.dot(unitary[randnum[1]].conj().T,state0).dot(unitary[randnum[1]]) - id)
elif result.get('01') == 1:
rho = np.kron(3*np.dot(unitary[randnum[0]].conj().T,state0).dot(unitary[randnum[0]] - id),3*np.dot(unitary[randnum[1]].conj().T,state1).dot(unitary[randnum[1]]) - id)
elif result.get('10') == 1:
rho = np.kron(3*np.dot(unitary[randnum[0]].conj().T,state1).dot(unitary[randnum[0]] - id),3*np.dot(unitary[randnum[1]].conj().T,state0).dot(unitary[randnum[1]]) - id)
else:
rho = np.kron(3*np.dot(unitary[randnum[0]].conj().T,state1).dot(unitary[randnum[0]] - id),3*np.dot(unitary[randnum[1]].conj().T,state1).dot(unitary[randnum[1]]) - id)
record_rho = record_rho + rho
record_rho = record_rho/snapshot_num
bell_state = np.array([[0.5, 0, 0, 0.5], [0, 0, 0, 0], [0, 0, 0, 0], [0.5, 0, 0, 0.5]])
print(record_rho)
print(distance(record_rho - bell_state))
% Pauli matrix as random unitary
pauli = eye(2);
pauli(:,:,2) = [1 0;0 -1j]; pauli(:,:,3) = 1/sqrt(2)*pauli(:,:,2)'*[1 1;1 -1];
% two computational basis |0) and |1)
state = eye(4);
state0 = [1;0];
state1 = [0;1];
psi = 1/sqrt(2)*[1;0;0;1];
rho = psi*psi';
record_rho = zeros(4);
n = 6000;
for i = 1:n
randnum = randi([1 3],[1 2]);
rhot = kron(pauli(:,:,randnum(1)),pauli(:,:,randnum(2)))*rho*kron(pauli(:,:,randnum(1)),pauli(:,:,randnum(2)))';
prob1 = state(:,1)'*rhot*state(:,1);
prob2 = state(:,2)'*rhot*state(:,2);
prob3 = state(:,3)'*rhot*state(:,3);
prob4 = state(:,4)'*rhot*state(:,4);
% Utilizing if to simulate the quantum measurement
% The inverse process is using the formula of eq(S44) in supplemental
% material of the original paper
if rand < prob1
rhot = 3*pauli(:,:,randnum(1))'*(state0*state0')*pauli(:,:,randnum(1)) - eye(2);
rhot = kron(rhot,3*pauli(:,:,randnum(2))'*(state0*state0')*pauli(:,:,randnum(2)) - eye(2));
elseif rand < prob1 + prob2
rhot = 3*pauli(:,:,randnum(1))'*(state0*state0')*pauli(:,:,randnum(1)) - eye(2);
rhot = kron(rhot,3*pauli(:,:,randnum(2))'*(state1*state1')*pauli(:,:,randnum(2)) - eye(2));
elseif rand < prob1 + prob2 + prob3
rhot = 3*pauli(:,:,randnum(1))'*(state1*state1')*pauli(:,:,randnum(1)) - eye(2);
rhot = kron(rhot,3*pauli(:,:,randnum(2))'*(state0*state0')*pauli(:,:,randnum(2)) - eye(2));
else
rhot = 3*pauli(:,:,randnum(1))'*(state1*state1')*pauli(:,:,randnum(1)) - eye(2);
rhot = kron(rhot,3*pauli(:,:,randnum(2))'*(state1*state1')*pauli(:,:,randnum(2)) - eye(2));
end
record_rho = record_rho + rhot;
end
record_rho = record_rho/n;
sqrt(trace((rho - record_rho)'*(rho - record_rho)))
Fidelity(rho,record_rho)