So I'm trying to replicate the results from Huang's paper, following Pennylane's tutorial. But instead of using Pennylane (in particular their recently implemented ClassicalShadow class), I'm trying to implement it in Qutip for individual purposes.
The generator function I'm implementing is, in essence, a line-by-line copy of the one in Pennylane's tutorial but in Qutip. The problem is, however, whenever I try to reconstruct any particular state as a benchmark for the code, the distance metrics are not improving at all.
This is the code
import numpy as np
import qutip as qt
from qutip.qip.operations.gates import x_gate, y_gate, z_gate
from qutip.measurement import measure, measure_observable
def gen_cshadow(state, csh_size, num_qubits):
'''
Qobj, int, int, list --> array
This function takes as input a density matrix, prepared beforehand, and generates a minimal representation of it
called a classical shadow, which is a numpy array containing the information of the measurement result and the list index of
the applied single-qubit unitaries. The shape of the output array is determined by both csh_size,
the cardinality of the classical shadow, and num_qubits, the number of qubits in the multi-qubit system. This code is
an implementation of the protocol proposed in Nat. Phys. 16, 1050–1057 (2020) (10.1038/s41567-020-0932-7)
state: Density matrix, a qutip Qobj of shape (2**num_qubits, 2**num_qubits) and dimensions [[2,..., 2], [2,...,2]]
where the lists are of length num_qubits.
csh_size: Integer, defines the cardinality of the classical shadow.
num_qubits: Integer, defines the number of qubits of the system.
'''
# Initializing single-qubit Clifford gates (Random Pauli Basis)
unitaries = [x_gate, y_gate, z_gate]
# Uniform random sampling of the Pauli measurements
uts_ids = np.random.randint(0, len(unitaries), size = (csh_size, num_qubits), dtype = int)
# Initializing outcome array of shape csh_size x num_qubits
outcomes = np.zeros((csh_size, num_qubits))
for snapshot in range(csh_size):
# Extract the Pauli basis measurement
observables = [unitaries[uts_ids[snapshot, qb]](N = num_qubits, target = qb) for qb in range(num_qubits)]
outcomes[snapshot] = [measure_observable(state, observables[i])[0] for i in range(num_qubits)]
return outcomes, uts_ids
def snapshot_state(out_list, obs_list):
'''
DOCSTRING - TO DO
'''
num_qubits = len(out_list)
zero_state = np.array([[1, 0], [0, 0]])
one_state = np.array([[0, 0], [0, 1]])
phaze_z = np.array([[1, 0], [0, -1j]], dtype = complex)
hadamard = (1/np.sqrt(2))*np.array([[1, 1], [1, -1]])
identity = np.eye(2)
rotations = [hadamard, hadamard @ phaze_z, identity]
snapshot = [1]
for outcome, basis in zip(out_list, obs_list):
state = zero_state if outcome == 1 else one_state
U = rotations[basis]
local_state = 3 * (U.conj().T @ state @ U) - identity
snapshot = np.kron(snapshot, local_state)
return snapshot
def reconstruction(cshadow):
'''
DOCSTRING - TO DO
'''
cs_size, num_qubits = cshadow[0].shape
out_lists, obs_lists = cshadow
# State reconstruction
shadow_state = np.zeros((2 ** num_qubits, 2 ** num_qubits), dtype = complex)
for i in range(cs_size):
shadow_state += snapshot_state(out_lists[i], obs_lists[i])
return shadow_state / cs_size
if name == 'main':
bell_state = qt.bell_state()
bell_state = qt.ket2dm(bell_state)
shadow = gen_cshadow(bell_state, 1000, 2)
sh_state = reconstruction(shadow)
print(sh_state)
pass
Now let me explain what I think the code should be doing. Just like the tutorial, the first function is a generator of measurement outcomes for the individual qubits and the IDs of the respective observables. The idea of the for loop within that function is to provide a uniformly sampled measurement basis for each of the qubits per snapshot, thus the observables list is just creating the appropriate $2^n \times 2^n$ matrices ($n$ is the number of qubits) needed for each single-qubit operation.
The following line should be simulating a single-shot measurement (projective according to source) of the state $\rho$ in the measurement basis of the respective unitary, which I think is the same thing they do in the tutorial. The rest of the code is simply the reconstruction of the state using numpy.
I have not been able to pinpoint what is wrong with the code, the only thing I can think of is that the first function is not generating what I think. Any help would be highly appreciated.
