Start QUBO with Qiskit¶

In this section, we are going to calculate some expection values in "quantum way".

1. We will first understand how to construct our state basis such as $|100\rangle$ using different ways in 'Define Basis' section.
2. We will introduce how to calculate Hamiltonian and its expectation value.

Define Basis¶

As we mentioned before, we use $\lvert xyz\rangle$, $x,y,z\in (0,1)$ to define the computational basis in our previous example. We can define these value in the following 3 ways.

1. Define one-qubit state $\lvert 0 \rangle$ and $\lvert 1 \rangle$ and computer their tensor products by using statevector
2. Initialize statevector object from an integer such as 0 or 1.
3. A more concise way.

First, let's import some necessary packages

In [1]:

from qiskit import *

from qiskit import *

1. Define one-qubit to a statevector¶

In [2]:

from qiskit.quantum_info import Statevector

# Defining 0
zero = Statevector([1,0]) # Remember that the vector form for |0\ is [1,0]!
print("zero is", zero)

# Defining 1
one = Statevector([0,1]) # Remember that the vector form for |1\ is [0,1]!
print("one is", one)

from qiskit.quantum_info import Statevector # Defining 0 zero = Statevector([1,0]) # Remember that the vector form for |0\ is [1,0]! print("zero is", zero) # Defining 1 one = Statevector([0,1]) # Remember that the vector form for |1\ is [0,1]! print("one is", one)

zero is Statevector([1.+0.j, 0.+0.j],
            dims=(2,))
one is Statevector([0.+0.j, 1.+0.j],
            dims=(2,))

2. Initialize statevector object from an integer such as 0 or 1.¶

In [3]:

# Defining 0
zero = Statevector.from_int(0, dims =2) # 0 @position 0 in dims = 2 array
print("zero is", zero)

# Defining 1
one = Statevector.from_int(1, dims =2) # 1 @position 1 in dims = 2 array
print("one is", one)

# Defining 0 zero = Statevector.from_int(0, dims =2) # 0 @position 0 in dims = 2 array print("zero is", zero) # Defining 1 one = Statevector.from_int(1, dims =2) # 1 @position 1 in dims = 2 array print("one is", one)

zero is Statevector([1.+0.j, 0.+0.j],
            dims=(2,))
one is Statevector([0.+0.j, 1.+0.j],
            dims=(2,))

Now we can construct states with a higher number of qubits by computing tensor products with the tensor method.

In [4]:

# Computing state by using `tensor` method.
psi = one.tensor(zero.tensor(zero)) # |100\
print ("psi is",psi)

# Computing state by using `tensor` method. psi = one.tensor(zero.tensor(zero)) # |100\ print ("psi is",psi)

psi is Statevector([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j,
             0.+0.j],
            dims=(2, 2, 2))

The result should show the following:

psi is Statevector([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], dims=(2, 2, 2))

The results shows 4 means $\lvert 100\rangle$ in binery representation. Remember, we start our at 0 to 7 for 1 to 8 in binery form.

3. A more concise way¶

In [5]:

# A more concise way
psi = one^zero^zero
print ("psi is",psi)
psi.draw("latex")

# A more concise way psi = one^zero^zero print ("psi is",psi) psi.draw("latex")

psi is Statevector([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j,
             0.+0.j],
            dims=(2, 2, 2))

Out[5]:

$$ |100\rangle$$

The result should show it's statevector

psi is Statevector([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], dims=(2, 2, 2))

and $|100\rangle$.

Also, a faster way of constructing the $|100\rangle$ is using, again, the from_int method:

In [6]:

# A faster way to construct |100\ using from_int method
psi = Statevector.from_int(4, dims = 8)
print("psi is", psi)

# A faster way to construct |100\ using from_int method psi = Statevector.from_int(4, dims = 8) print("psi is", psi)

psi is Statevector([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j,
             0.+0.j],
            dims=(2, 2, 2))

The result should show the following:

psi is Statevector([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], dims=(2, 2, 2))

We specify that we are working with three qubits by setting dims = 8 (because we need 8 amplitudes to define a three-qubit state)

After knowing how to create a basis, we can learn how to create a superposition state! In Qiskit, this is very easy, you just use the following command:

In [7]:

import numpy as np

ghz = 1/np.sqrt(2)*(zero^zero^zero)+1/np.sqrt(2)*(one^one^one)
print('ghz is:', ghz)

import numpy as np ghz = 1/np.sqrt(2)*(zero^zero^zero)+1/np.sqrt(2)*(one^one^one) print('ghz is:', ghz)

ghz is: Statevector([0.70710678+0.j, 0.        +0.j, 0.        +0.j,
             0.        +0.j, 0.        +0.j, 0.        +0.j,
             0.        +0.j, 0.70710678+0.j],
            dims=(2, 2, 2))

You should get:

ghz is: Statevector([0.70710678+0.j, 0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j, 0. +0.j, 0.70710678+0.j], dims=(2, 2, 2))

which created a state $\frac{1}{\sqrt{2}}|000\rangle + \frac{1}{\sqrt{2}}|111\rangle$.

!! in Qiskit, ^, which is a tensor operator, has a lower precedence than $+$ in python, therefore, we need to use parenthesis more carefully!

To compute our expectation values, we need not only the state values but also Hamiltonians. So, let's work on the Hamiltonians in the next section!

Hamiltonians¶

Let's work with hamiltonians in this section:

First, we intorduce Pauli from the qiskit.quantum_info object.

In [8]:

from qiskit.quantum_info import Pauli

# Create Z0Z1I matrix mentioned before!
Z0Z1 = Pauli("ZZI")

# Print some results:
print("Z0Z1 is:", Z0Z1)
print("And its matrix is:")
print(Z0Z1.to_matrix())

from qiskit.quantum_info import Pauli # Create Z0Z1I matrix mentioned before! Z0Z1 = Pauli("ZZI") # Print some results: print("Z0Z1 is:", Z0Z1) print("And its matrix is:") print(Z0Z1.to_matrix())

Z0Z1 is: ZZI
And its matrix is:
[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j -1.+0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  1.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  1.+0.j]]

You should get a matrix representation of Z0Z1.

Also, if you find the matrix hard to read, it's also welcome to use the sparse method like below:

In [9]:

# We can also using the sparse method to demonstrate our results
print("The sparse representation of Z0Z1 is")
print(Z0Z1.to_matrix(sparse=True))

# We can also using the sparse method to demonstrate our results print("The sparse representation of Z0Z1 is") print(Z0Z1.to_matrix(sparse=True))

The sparse representation of Z0Z1 is
<Compressed Sparse Row sparse matrix of dtype 'complex128'
	with 8 stored elements and shape (8, 8)>
  Coords	Values
  (0, 0)	(1+0j)
  (1, 1)	(1+0j)
  (2, 2)	(-1+0j)
  (3, 3)	(-1+0j)
  (4, 4)	(-1+0j)
  (5, 5)	(-1+0j)
  (6, 6)	(1+0j)
  (7, 7)	(1+0j)

The main drawback for the Pauli object is that you cannot add them or multiply them by scaler. To compute something like $Z_{0}Z_{1}+Z_{1}Z_{2}$, we need to convert Pauli object to PauliOp first.

In [10]:

from qiskit.quantum_info import SparsePauliOp

H_cut = SparsePauliOp(Pauli("ZZI")) + SparsePauliOp(Pauli("ZIZ"))
print("H_cut is:", H_cut)
print("The sparse representation of H_cut is", H_cut.to_matrix(sparse=True))

from qiskit.quantum_info import SparsePauliOp H_cut = SparsePauliOp(Pauli("ZZI")) + SparsePauliOp(Pauli("ZIZ")) print("H_cut is:", H_cut) print("The sparse representation of H_cut is", H_cut.to_matrix(sparse=True))

H_cut is: SparsePauliOp(['ZZI', 'ZIZ'],
              coeffs=[1.+0.j, 1.+0.j])
The sparse representation of H_cut is <Compressed Sparse Row sparse matrix of dtype 'complex128'
	with 8 stored elements and shape (8, 8)>
  Coords	Values
  (0, 0)	(2+0j)
  (1, 1)	0j
  (2, 2)	0j
  (3, 3)	(-2+0j)
  (4, 4)	(-2+0j)
  (5, 5)	0j
  (6, 6)	0j
  (7, 7)	(2+0j)

Also we can try more complicated method using parsePauliOp(pauli_strings, coef) by specifying Paulis and it's coefficient of the following equation:

$$ 0.5Z_{0}Z_{1}+2Z_{0}Z_{2}-Z_{1}Z_{2}+Z_{1}-5Z_{2} $$

In [11]:

from qiskit.quantum_info import SparsePauliOp

pauli_strings = ["ZZI", "ZIZ", "IZZ", "IZI", "IIZ"]
coef = [-0.5, 2.0, -1.0, 1.0, -1.5]

H_ising = SparsePauliOp(pauli_strings, coef)
print("H_ising is:", H_ising)
print("The sparse representation of H_ising is", H_ising.to_matrix(sparse=True))

from qiskit.quantum_info import SparsePauliOp pauli_strings = ["ZZI", "ZIZ", "IZZ", "IZI", "IIZ"] coef = [-0.5, 2.0, -1.0, 1.0, -1.5] H_ising = SparsePauliOp(pauli_strings, coef) print("H_ising is:", H_ising) print("The sparse representation of H_ising is", H_ising.to_matrix(sparse=True))

H_ising is: SparsePauliOp(['ZZI', 'ZIZ', 'IZZ', 'IZI', 'IIZ'],
              coeffs=[-0.5+0.j,  2. +0.j, -1. +0.j,  1. +0.j, -1.5+0.j])
The sparse representation of H_ising is <Compressed Sparse Row sparse matrix of dtype 'complex128'
	with 8 stored elements and shape (8, 8)>
  Coords	Values
  (0, 0)	0j
  (1, 1)	(1+0j)
  (2, 2)	(1+0j)
  (3, 3)	(-2+0j)
  (4, 4)	(-3+0j)
  (5, 5)	(6+0j)
  (6, 6)	(-4+0j)
  (7, 7)	(1+0j)

Now, we can calculate the expectation value of our H_cut by using command: psi.expectation_value(H_cut).

The expectation value should be (-2+0j).

In [12]:

print("The expectation value is", psi.expectation_value(H_cut))

print("The expectation value is", psi.expectation_value(H_cut))

The expectation value is (-2+0j)

Since the $Z_{0}Z_{1} + Z_{0}Z_{2}$ is the Hamiltonian for the Max-Cut problem for the three-node. The Hamiltonian is represented by $|100\rangle$ (vertex 0 on one set, and 1 and 2 in the other) cuts the two edges of the graph and is an optimal solution.

---