Expected behavior
qml.density_matrix should return the density matrix that represents the simulated quantum
state. For a pure state :math:|\psi\rangle, this is
rho = |psi><psi| = outer(psi, conj(psi))
Consider the following one-qubit Clifford circuit:
It prepares
|psi> = (|0> + i|1>) / sqrt(2)
Therefore, the expected density matrix is
[[0.5+0.j 0.0-0.5j]
[0.0+0.5j 0.5+0.j ]]
This matrix is Hermitian, has trace one, and assigns probability 0.5 to both computational
basis states.
Actual behavior
default.clifford returns:
[[ 0.5+0.j 0.0+0.5j]
[ 0.0+0.5j -0.5+0.j ]]
This is not a valid density matrix:
Hermitian: False
Trace: 0j
P(0): 0.49999997
P(1): -0.49999997
The result gives the |1> basis state a negative probability and has total probability zero.
The same circuit on default.qubit returns the expected valid density matrix.
Additional information
The issue is in DefaultClifford._measure_density_matrix in
pennylane/devices/default_clifford.py. The implementation constructs the outer product without
complex-conjugating the second state vector:
state_vector = math.array(tableau_simulator.state_vector(endian="big"))
return math.reduce_dm(math.einsum("i, j->ij", state_vector, state_vector), wires)This computes :math:|\psi\rangle|\psi\rangle^T, rather than the required
:math:|\psi\rangle\langle\psi|.
The problem is hidden for states whose amplitudes are all real, but it produces incorrect results
for valid Clifford states containing complex amplitudes, such as states prepared using the S
gate.
A possible fix is to conjugate the second state-vector operand before constructing the outer
product:
math.einsum("i, j->ij", state_vector, math.conj(state_vector))Source code
import numpy as np
import pennylane as qml
def get_density_matrix(device_name):
dev = qml.device(device_name, wires=1)
@qml.qnode(dev)
def circuit():
qml.Hadamard(0)
qml.S(0)
return qml.density_matrix(wires=[0])
return np.asarray(circuit())
actual = get_density_matrix("default.clifford")
expected = get_density_matrix("default.qubit")
print("default.clifford:")
print(actual)
print("\ndefault.qubit:")
print(expected)
print("\nMatches default.qubit:", np.allclose(actual, expected))
print("Hermitian:", np.allclose(actual, actual.conj().T))
print("Trace:", np.trace(actual))Tracebacks
default.clifford:
[[ 0.49999997+0.j 0. +0.49999997j]
[ 0. +0.49999997j -0.49999997+0.j ]]
default.qubit:
[[0.5+0.j 0. -0.5j]
[0. +0.5j 0.5+0.j ]]
Matches default.qubit: False
Hermitian: False
Trace: 0j
System information
Version: 0.45.0
Platform info: macOS
Python version: 3.13.13
Existing GitHub issues
Expected behavior
qml.density_matrixshould return the density matrix that represents the simulated quantumstate. For a pure state :math:
|\psi\rangle, this isConsider the following one-qubit Clifford circuit:
It prepares
Therefore, the expected density matrix is
This matrix is Hermitian, has trace one, and assigns probability
0.5to both computationalbasis states.
Actual behavior
default.cliffordreturns:This is not a valid density matrix:
The result gives the
|1>basis state a negative probability and has total probability zero.The same circuit on
default.qubitreturns the expected valid density matrix.Additional information
The issue is in
DefaultClifford._measure_density_matrixinpennylane/devices/default_clifford.py. The implementation constructs the outer product withoutcomplex-conjugating the second state vector:
This computes :math:
|\psi\rangle|\psi\rangle^T, rather than the required:math:
|\psi\rangle\langle\psi|.The problem is hidden for states whose amplitudes are all real, but it produces incorrect results
for valid Clifford states containing complex amplitudes, such as states prepared using the
Sgate.
A possible fix is to conjugate the second state-vector operand before constructing the outer
product:
Source code
import numpy as np import pennylane as qml def get_density_matrix(device_name): dev = qml.device(device_name, wires=1) @qml.qnode(dev) def circuit(): qml.Hadamard(0) qml.S(0) return qml.density_matrix(wires=[0]) return np.asarray(circuit()) actual = get_density_matrix("default.clifford") expected = get_density_matrix("default.qubit") print("default.clifford:") print(actual) print("\ndefault.qubit:") print(expected) print("\nMatches default.qubit:", np.allclose(actual, expected)) print("Hermitian:", np.allclose(actual, actual.conj().T)) print("Trace:", np.trace(actual))Tracebacks
System information
Existing GitHub issues