I am trying to approach Markov chains as a use case for Quantum Computing. For this I took the simple introductory case from Wikipedia (Markow-Kette):

  • $v_0=\begin{pmatrix} 1 & 0 & 0 \end{pmatrix}$
  • M=$\begin{pmatrix} 0 & 0.2 & 0.8\\ 0.5 & 0 & 0.5\\ 0.1 & 0.9 & 0\\ \end{pmatrix}$
  • $v_1=v_0M$
  • $v_3=v_0M^3$

As a starter/"appetizer" I implemented this situation with the following 3-liner using Numpy and got very straight forwardly the correct result v3=[0.37 0.126 0.504] which is in line with the result from the Wikipedia Site:

M = np.array([[0, 0.2, 0.8], [0.5, 0, 0.5], [0.1, 0.9, 0]])
v3=v0.dot(matrix_power(M, 3))

Now I'm stuck porting the whole thing to Qiskit, where I am really would be satisfied with a simulator-based solution:

qiskit.IBMQ.save_account('your_token', overwrite=True)
n_wires = 1
n_qubits = 1
provider = qiskit.IBMQ.get_provider('ibm-q')
backend = Aer.get_backend('qasm_simulator')

Asking Google led me to Matrix product state simulation method, but it seems to be not obvious how to apply this to my simple problem. A simple nudge in the right direction would be really appreciated.


1 Answer 1


I made an implementation, I'm not sure whether it has advantages or how it generalizes. I hope it might steer you (or us) in the right direction.

The approach is as follows:

You have a matrix $M$, then:

  1. Create a circuit with $2N=6$ qubits
    • The first 3 qubits represent 'being in state $i$' ($i$ is a state in $N$) (i.e. $|001>$ is state $0$, $|010>$ is state $1$, $|100>$ state $2$)
    • The last 3 qubits represent 'going to state $j$' ($j$ is a state in $N$)
  2. Controlled go from state $i$ to state $j$ with probability $M_{i,j}$ *
  3. Make sure to not go to state $j'$ when you're going to state $j$
  4. Make sure the naming of states checks out (i.e. $|001>$ state $0$, $|010>$ is state $1$, $|100>$ is state $2$)
  5. Repeat steps 2-4 for all $i$
  6. Swap the last 3 qubits with the first 3 qubits, (in words, these were the states that you're 'going to' and now they are the state 'you are in').
  7. Reset the last 3 qubits.

The circuit now looks like this:

markov step

With this circuit, you can do the same steps as you proposed before, so in your case you have:

import numpy as np
M = np.array([[0, 0.2, 0.8], [0.5, 0, 0.5], [0.1, 0.9, 0]])
v2=v0.dot(np.linalg.matrix_power(M, 2))
v3=v0.dot(np.linalg.matrix_power(M, 3))


[0.  0.2 0.8]
[0.18 0.72 0.1 ]
[0.37  0.126 0.504]

In Qiskit, this now is the following **:

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, transpile
from numpy import pi

# Inilialise registers
qreg_q = QuantumRegister(6, 'q')
creg_c = ClassicalRegister(3, 'c')

# Create Markov Step as a circuit
markov_step = QuantumCircuit(qreg_q)

# Create the Markov Step
# From state 0 to state 1 and 2
markov_step.cu3(2*np.arccos(np.sqrt(M[0,1])), pi/2, pi/2, qreg_q[0], qreg_q[4])
markov_step.ccx(qreg_q[4], qreg_q[0], qreg_q[5])
markov_step.cx(qreg_q[0], qreg_q[4])

# From state 1 to state 0 and 2
markov_step.cu3(2*np.arccos(np.sqrt(M[1,2])), pi/2, pi/2, qreg_q[1], qreg_q[5])
markov_step.ccx(qreg_q[5], qreg_q[1], qreg_q[3])
markov_step.cx(qreg_q[1], qreg_q[5])

# From state 2 to state 0 and 1
markov_step.cu3(2*np.arccos(np.sqrt(M[2,0])), pi/2, pi/2, qreg_q[2], qreg_q[3])
markov_step.ccx(qreg_q[3], qreg_q[2], qreg_q[4])
markov_step.cx(qreg_q[2], qreg_q[3])

# Swap
markov_step.swap(qreg_q[0], qreg_q[3])
markov_step.swap(qreg_q[1], qreg_q[4])
markov_step.swap(qreg_q[2], qreg_q[5])

# Initialise circuit
circuit = QuantumCircuit(qreg_q,creg_c)

# Initialise state (1,0,0)

# Do the markov step n times
n = 3
for _ in range(n):
    for ins in markov_step:
        circuit.append(ins[0], ins[1], ins[2])

# Measure outcome
circuit.measure(qreg_q[:3], creg_c)

And you can run it by

from qiskit.visualization import plot_histogram
backend = provider.get_backend('ibmq_qasm_simulator')
job = backend.run(circuit)
result = job.result()
counts = result.get_counts(circuit)


And the output:

output counts

Which is approximately equal to the exact answers [0.37 0.126 0.504].

This approach definitely isn't perfect and I'm quite sure optimizations can be made (e.g. not using 1 qubit per state, but using the full $2^N$ possible states) and I'm not sure how to go to larger state spaces. But it's the first step!

Some notes:

* note: Step 2 is not trivial. I implemented it as a controlled $X$-rotation. An $X$-rotation is (according to the qiskit-textbook) given by

\begin{equation} R_x(\theta) = \begin{pmatrix} \cos(\theta /2) & -i \sin(\theta /2) \\ -i \sin(\theta /2) & \cos(\theta /2) \end{pmatrix}, \end{equation}

and it brings the $|0>$ state to $\cos(\theta /2) |0> - i \sin(\theta /2) |1>$. Now, for example, we want the target qubit to be in state $|0>$ with probability 0.2 (when starting the target is in state 0). The probability of finding the target in state $|0>$ is $|\cos(\theta /2)|^2$ and this must be equal to $0.2$. Thereby, it can be found that $\theta = 2 \arccos{\sqrt{0.2}}$.

** note: Implementing the circuit is a bit annoying because not all gates are available in all backends. Sometimes, you have to 'translate' the gates. A controlled $X$-rotation (CRX) is also a controlled U3 gate, where $CR_x(\theta) = CU3(\theta, \pi/2,\pi/2)$. Some backend also doesn't allow for the CU3 gate, but they do allow the multi-controlled gate .mcu3. In that case, just put your control qubit in a list by putting square brackets around it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.