Computing specific observable for a Clifford circuit

Given a clifford circuit, how can I simulate it efficiently to get the 2-RDM (reduced density matrix) $$D^2=\langle a_i^\dagger a_j^\dagger a_k a_l\rangle$$ in the presence of deplorizing error? The simulation will take a lot of time, which seems conter-intuitive for simulation of clifford circuit. It seems that the cost is mainly caused by CircuitSampler. Is there a way to avoid the using of it while can get the noisy expectation result for the 2-RDM or other weighted pauli string-like observables, such as a Hamiltonian, in the Clifford case?

from qiskit import *
import numpy as np
import itertools

#Circuit imports
import matplotlib
# import cirq
import qiskit.quantum_info as qi
from qiskit.quantum_info import random_clifford
matplotlib.use('Agg')
from qiskit_nature.operators.second_quantization import FermionicOp
from qiskit.opflow.converters import CircuitSampler
from qiskit_nature.mappers.second_quantization import JordanWignerMapper
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit.utils import QuantumInstance
from qiskit.opflow import Z, X, I, StateFn, CircuitStateFn
from qiskit.opflow.expectations import PauliExpectation
import qiskit.providers.aer.noise as noise
from qiskit.providers.aer import AerSimulator
np.set_printoptions(threshold=np.inf)

prob_1 = 1e-3  # 1-qubit gate
prob_2 = 1e-2  # 2-qubit gate

# Depolarizing quantum errors
error_1 = noise.depolarizing_error(prob_1, 1)
error_2 = noise.depolarizing_error(prob_2, 2)

# Add errors to noise model
noise_model = noise.NoiseModel()
noise_model.add_all_qubit_quantum_error(error_1, ['s', 'sdg', 'x', 'y', 'z', 'i', 'h'])
noise_model.add_all_qubit_quantum_error(error_2, ['cx', 'cz', 'swap'])
num_qubits = 10
two_rdm = np.zeros((num_qubits,) * 4)
mapper = JordanWignerMapper()
converter = QubitConverter(mapper=mapper, two_qubit_reduction=False)
backend =  AerSimulator(method='stabilizer',
noise_model=noise_model)
quantum_instance = QuantumInstance(backend=backend,
shots=1000,
noise_model=noise_model)
circuit_sampler = CircuitSampler(quantum_instance)
cliff = random_clifford(num_qubits)
cliff = cliff.to_circuit()

for i, j, k, l in itertools.product(range(num_qubits), repeat=4):
s = "+_{i} +_{j} -_{k} -_{l}".format(i=str(i),j=str(j),k=str(k),l=str(l))
fermi_term = FermionicOp(s, register_length=num_qubits)
qubit_term = converter.convert(fermi_term, num_particles=(1,1))

# Evaluate the 2-RDM term w.r.t. the given circuit
temp = ~StateFn(qubit_term) @ CircuitStateFn(primitive=cliff, coeff=1.)
temp = circuit_sampler.convert(
PauliExpectation().convert(temp)
).eval()
two_rdm[i,j,k,l] = temp

• While I do not know a specific answer, I am afraid quite often Clifford circuits with noise are no longer easy to simulate (since noise is effectively a non-Clifford gate.) This paper could be useful: arxiv.org/abs/2001.08373 Commented Dec 29, 2021 at 23:16