simulate_moment_steps return invalid density matrices

[CirqUser199]

When I try to simulate the density matrix through a circuit using simulate_moment_steps, I have some density matrices that have a trace different than 1 for some steps. For instance, when I execute the simple code below, the last line returns 1/2. If I did no mistake, this means that the simulator returns a density matrix that has trace less than 1. From what I can see, this could be caused by the insert strategy (maybe it interferes with the moment by moment simulator). Did you already hear about such an issue?

import cirq
import numpy as np
from cirq.circuits import InsertStrategy
q0, q1, q2, q3, q4, q5, q6 = cirq.LineQubit.range(7)

example_circuit=cirq.Circuit()
example_circuit.append([cirq.rx(3*np.pi/7)(q0),cirq.rx(3*np.pi/5)(q1),cirq.rx(np.pi/4)(q2)])
example_circuit.append([cirq.ry(1*np.pi/5)(q0),cirq.ry(3*np.pi/2)(q1),cirq.ry(2*np.pi/5)(q2)])
example_circuit.append([cirq.H(q3),cirq.H(q4)], strategy=InsertStrategy.NEW_THEN_INLINE)
example_circuit.append([cirq.CNOT(q3,q5),cirq.CNOT(q4,q6)], strategy=InsertStrategy.NEW_THEN_INLINE)

sim = cirq.DensityMatrixSimulator()
example_densities=[]

for i, step in enumerate(sim.simulate_moment_steps(example_circuit)):
    example_densities.append(step.density_matrix())
    
np.trace(example_densities[2])

[cduck]

It appears that the simulator is re-using and modifying the matrix returned by .density_matrix() during later steps. This was probably intended behavior to improve speed but can cause strange usability issues like you found. If others have done this too I would consider it a bug.

Try this instead. You get the correct result if you use the matrix before the next simulation step.

for i, step in enumerate(sim.simulate_moment_steps(example_circuit)):
    print(i, np.trace(step.density_matrix()))

0 (1+0j)
1 (1-1.1175871e-08j)
2 (1-1.1175869e-08j)
3 (1-1.117587e-08j)

If you do need to store all the density matrices (I would not recommend this for circuits much larger than your example), just make copies of the matrix:

    example_densities.append(np.copy(step.density_matrix()))

[CirqUser199]

Thanks a lot for your quick answer!

[dabacon]

The density matrix is now a copy by default. This should fix your original code by default. If you need the more optimized version you can request it by setting the copy boolean when accessing the density matrix.