Use simulator.simulate() to obtain probability distribution instead of simulator.run() (cirq)

So far, I have my code set up to find the probability distribution using simulator.run(), and I was wondering if I can use the simulator.simulate() function instead to obtain a probability distribution.

The context of my question is quantum walk simulations, but my question doesn't require any knowledge of quantum walks. So, I am looking at quantum walk simulations and want to get the probability distribution of finding the walk in a certain position state (which are represented by the computational basis states of the Hilbert space). The way my code is set up now is as follows: I define the quantum walk in terms of a circuit, where I take a measurement at the end of the circuit (i.e. at the end of each 'walk' I measure which (position) state the walk is found in). I then repeat this process multiple times and find which state is measured each time:

simulator = cirq.Simulator()
result = simulator.run(circuit, repetitions=sample_number)
final = result.histogram(key='x')


Then for 50 repetitions the final looks like Counter({5: 20, 8: 5, 2: 25}) for example, which means that out of 50 runs, I find the walk 20 times in the (position) state 5, 5 times in state 8 etc. I then extract two arrays from this, one with the positions: [5,8,2] and one with the counts: [20,5,25]. I then normalise the latter to obtain probabilities. Then I simply plot these two arrays and find the probability distribution of finding the walk in a certain position.

I was wondering if I can similarly obtain a probability distribution using simulator.simulate() instead of simulator.run(). I still have the same circuit which represents the quantum walk, but now I don't take a measurement, and instead I would like to use

simulator = cirq.Simulator()
result = simulator.simulate(circuit)


Printing the result will give me something like

qubits: (q0, q1, ...) (list of qubits used)
output vector: [0.      a.      c.    ...      0.] (amplitudes of all different states)


It seems like this output vector contains all the information necessary to make a probability distribution: I have the amplitudes of each (position) state, and so squaring them will give me the probability of finding the walk in each state.

The only problem is that I am not sure how to actually 'read' this output vector, and find which state the amplitude corresponds to? For example, is there a way to 'translate' the states in the same way the simulator.run() automatically does (so for 4 qubits, state $$|0010\rangle$$ becomes $$|8\rangle$$ in binary, which the final in the run function seems to automatically do)?

Many thanks!

Yes, the probability distribution should be directly accessible from the full state vector simulation's result, i.e. the state vector as you wrote.

You were also right that the amplitudes in the state vector will correspond to the components' binary representation. Some subtlety there is that the bit ordering is the other way around as you wrote it, i.e. 1000=8 and 0010=2, and that this value also depends on the qubit ordering, which you can control with the qubit_order argument. Also, you can use cirq.dirac_notation to double-check and debug your understanding.

In the following example, only the $$|{6}\rangle = |{110}\rangle$$ and $$|{4}\rangle = |{100}\rangle$$ components have amplitudes, so they will be in the 6th and 4th (in 0 based indexing) place in the vector:

import cirq

a,b,c = cirq.LineQubit.range(3)

circuit = cirq.Circuit(cirq.X(a), cirq.H(b))

print("=== Method one, using cirq.Circuit.final_state_vector ===")
print("state vector:")
state1=circuit.final_state_vector(initial_state=0, qubit_order=[a,b,c])
print(state1)
print("state vector in Dirac notation:")
print(cirq.dirac_notation(state1))

print("=== Method two, using cirq.Simulator explicitly ===")
print("state vector:")
state2 = cirq.Simulator().simulate(circuit, qubit_order=[a,b,c], initial_state=0).state_vector()
print(state2)
print("state vector in Dirac notation:")
cirq.dirac_notation(state2)


=== Method one, using cirq.Circuit.final_state_vector ===
state vector:
[0.        +0.j 0.        +0.j 0.        +0.j 0.        +0.j
0.70710678+0.j 0.        +0.j 0.70710678+0.j 0.        +0.j]
state vector in Dirac notation:
0.71|100⟩ + 0.71|110⟩
=== Method two, using cirq.Simulator explicitly ===
state vector:
[0.        +0.j 0.        +0.j 0.        +0.j 0.        +0.j
0.70710677+0.j 0.        +0.j 0.70710677+0.j 0.        +0.j]
state vector in Dirac notation:
0.71|100⟩ + 0.71|110⟩

• Thanks a lot for your answer! It's a lot clearer now! I just have a follow-up question: say my Hilbert space is a composite Hilbert space consisting of position Hilbert space and an auxiliary Hilbert space, i.e. if I have 3 qubits in total, then 2 qubits indicate the position state, and 1 qubit indicates the auxiliary state. Then how can I find the state vector only for the two qubits indicating the position? In other words, I want to find the amplitudes for my state to be in a certain position, regardless of what the auxiliary component is. How can I do that? Thanks in advance! Aug 19, 2022 at 17:38
• I recommend checking out the methods here: github.com/quantumlib/Cirq/blob/v1.0.0/cirq-core/cirq/linalg/… - if you have a product state, you can use cirq.sub_state_vector, but otherwise, you'll have a mixed state in your "position sub-system", so you'll have to deal with that using cirq.partial_trace_of_state_vector_as_mixture for example. Aug 19, 2022 at 17:55
• Thanks for the link, I had a look at the functions, but they doesn't seem to do what I would like. The cirq.sub_state_vector doesn't work as my states are entangled. And the cirq.partial_trace_of_state_vector_as_mixture returns ((0.5, array([0, 0, 0, 1, 0, 0, 0, 0])), (0.5, array([0, 0, 0, 0, 0, 1, 0, 0]))) for the state vector 0.71|0111⟩ + 0.71|1010⟩, whereas I would like it to return `[0, 0, 0, 0.71, 0, 0.71, 0, 0], is this not possible? Aug 19, 2022 at 18:26
• Well, that's the thing, it looks like that your subsystem is in a mixed state with a probability of 50% of being in one of two pure states. That is not the same in general as a pure state with an equal superposition between the two states, as this is a classical probability instead of quantum at this point. I suggest reading 2.4.3 The reduced density operator from Nielsen & Chuang. Aug 19, 2022 at 19:12