# Very simple matrix operation (Markov chain) in terms of Qiskit

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]])
v0=np.array([1,0,0])
v1=v0.dot(M)
print(v1)
v3=v0.dot(matrix_power(M, 3))
print(v3)


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.

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: 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]])
v0=np.array([1,0,0])
v1=v0.dot(M)
print(v1)
v2=v0.dot(np.linalg.matrix_power(M, 2))
print(v2)
v3=v0.dot(np.linalg.matrix_power(M, 3))
print(v3)


output:

[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, qreg_q)
markov_step.ccx(qreg_q, qreg_q, qreg_q)
markov_step.cx(qreg_q, qreg_q)

# 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, qreg_q)
markov_step.ccx(qreg_q, qreg_q, qreg_q)
markov_step.cx(qreg_q, qreg_q)

# 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, qreg_q)
markov_step.ccx(qreg_q, qreg_q, qreg_q)
markov_step.cx(qreg_q, qreg_q)

# Swap
markov_step.swap(qreg_q, qreg_q)
markov_step.swap(qreg_q, qreg_q)
markov_step.swap(qreg_q, qreg_q)

# Initialise circuit
circuit = QuantumCircuit(qreg_q,creg_c)

# Initialise state (1,0,0)
circuit.x(0)

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

# 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)

plot_histogram(counts)


And the output: 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.