Building Variational Quantum Eigensolver from scratch

Let us assume that we are given a $2^N\times 2^N$$2^N\times2^N$ dimension matrix (like the one given below).

Today, we will be finding eigenvalues of such matrices using VQE-like circuits for such matrices. While doing so, we will cover the following topics:

1. Hamiltonian in Pauli Basis
2. Calculations of Expectation Value
3. Designing Ansatze
4. Running Simulations
5. Running Noisy Simulations
6. Extension to Excited Energy States
7. Interpreting the Results

 

Necessary Imports

xxxxxxxxxx
 
 
 
 
 
 
 
import numpy as np
 
 
import scipy as sp
 
import itertools
 
import functools as ft
 
%matplotlib notebook
 
import matplotlib.pyplot as plt
 
from mpl_toolkits import mplot3d
 
import re 
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
#!pip install qiskit
 
 
from qiskit import *
 
from qiskit.providers.aer.noise import NoiseModel
 
from qiskit.providers.aer.noise import QuantumError, ReadoutError
 
from qiskit.providers.aer.noise import phase_amplitude_damping_error
 
from qiskit.providers.aer.noise import depolarizing_error
 
from qiskit.providers.aer.noise import thermal_relaxation_error
 
from qiskit.ignis.mitigation.measurement import complete_meas_cal, CompleteMeasFitter
 
 
 
 
 
 

 

Hamiltonian in Pauli Basis

The normalized Pauli matrices $\{{\sigma }_0,{\sigma }_x,{\sigma }_y,{\sigma }_z\}/\sqrt{2}$$\{\sigma_0,\sigma_x,\sigma_y,\sigma_z\}/{\sqrt{2}}$ form an orthogonal basis of ${\mathcal{M}}_2$$\mathcal{M}_2$, this vector space can be endowed with a scalar product called the Hilbert-Schmidt inner product: $⟨A,B⟩=Tr(A^{†}B)$$⟨A,B⟩=Tr(A^{\dagger}B)$. Since the Pauli matrices anticommute, their product is traceless, and since they are Hermitian this implies that they are orthogonal with respect to that scalar product. Hence this property can be used for decomposing a Hamiltonian as:

$H=\sum _{i_1,\dots ,i_n=\{0,x,y,z\}}{h}_{i_1,\dots ,i_n}\cdot \frac{1}{2^n}{\sigma }_{i_1}\otimes \dots \otimes {\sigma }_{i_n}$$H = \sum_{i_{1},\ldots,i_{n} = \{0,x,y,z\}}h_{i_{1},\ldots,i_{n}}\cdot \frac{1}{2^{n}}\sigma_{i_1}\otimes\ldots\otimes\sigma_{i_n}$

$h_{i_1,\dots ,i_n}=\frac{1}{2^n}\text{Tr}(({\sigma }_{i_1}\otimes \dots \otimes {\sigma }_{i_n}{)}^{†}\cdot H)=\frac{1}{2^n}\text{Tr}(({\sigma }_{i_1}\otimes \dots \otimes {\sigma }_{i_n})\cdot H)$$h_{i_{1},\ldots,i_{n}} = \frac{1}{2^{n}}\text{Tr}\big((\sigma_{i_1}\otimes\ldots\otimes\sigma_{i_n})^\dagger \cdot H\big) = \frac{1}{2^{n}}\text{Tr}\big((\sigma_{i_1}\otimes\ldots\otimes\sigma_{i_n}) \cdot H\big)$

 

Decomposing Hamiltonian into Pauli Terms

xxxxxxxxxx
 
 
 
 
 
 
 
def decompose_ham_to_pauli(H):
 
 
    """Decomposes a Hermitian matrix into a linear combination of Pauli operators.
 

 
    Args:
 
        H (array[complex]): a Hermitian matrix of dimension 2**n x 2**n.
 

 
    Returns:
 
        tuple[list[float], list[string], list [ndarray]]: a list of coefficients, 
 
        a list of corresponding string representation of tensor products of Pauli 
 
        observables that decompose the Hamiltonian, and
 
        a list of their matrix representation as numpy arrays
 
    """
 

 
    n = int(np.log2(len(H)))
 
    N = 2 ** n
 

 
    # Sanity Checks
 
    if H.shape != (N, N):
 
        raise ValueError(
 
            "The Hamiltonian should have shape (2**n, 2**n), for any qubit number n>=1"
 
        )
 

 
    if not np.allclose(H, H.conj().T):
 
        raise ValueError("The Hamiltonian is not Hermitian")
 

 
    sI = np.eye(2, 2, dtype=complex)
 
    sX = np.array([[0, 1], [1, 0]], dtype=complex)
 
    sZ = np.array([[1, 0], [0,-1]], dtype=complex)
 
    sY = complex(0,-1)*np.matmul(sZ,sX)
 
    paulis = [sI, sX, sY, sZ]
 
    paulis_label = ['I', 'X', 'Y', 'Z']
 
    obs = []
 
    coeffs = []
 
    matrix = []
 
    
 
    for term in itertools.product(paulis, repeat=n):
 
        matrices = [pauli for pauli in term]
 
        coeff = np.trace(ft.reduce(np.kron, matrices) @ H) / N 
 
        coeff = np.real_if_close(coeff).item()
 
        
 
        # Hilbert-Schmidt-Product
 
        if not np.allclose(coeff, 0): 
 
            coeffs.append(coeff)
 
            obs.append(''.join([paulis_label[[i for i, x in enumerate(paulis) 
 
            if np.all(x == t)][0]]+str(idx) for idx, t in enumerate(reversed(term))]))
 
            matrix.append(ft.reduce(np.kron, matrices))
 

 
    return obs, coeffs , matrix
 

 
 
 
 
 
 

 

Composing Hamiltonian into Pauli Terms

xxxxxxxxxx
 
 
 
 
 
 
 
def compose_ham_from_pauli(terms, coeffs):
 
 
    """Composes a Hermitian matrix from a linear combination of Pauli operators.
 

 
    Args:
 
        tuple[list[float], list[string]]: a list of coefficients, 
 
        a list of corresponding string representation of tensor products of 
 
        Pauli observables that decompose the Hamiltonian.
 

 
    Returns:
 
        H (array[complex]): a Hermitian matrix of dimension 2**n x 2**n.
 
    """
 

 
    pauli_qbs = [re.findall(r'[A-Za-z]|-?\d+\.\d+|\d+', x) for x in terms]
 
    qubits = max([max(list(map(int, x[1:][::2]))) for x in pauli_qbs])
 

 
    N = int(2**np.ceil(np.log2(qubits+1)))
 

 
    sI = np.eye(2, 2, dtype=complex)
 
    sX = np.array([[0, 1], [1, 0]], dtype=complex)
 
    sZ = np.array([[1, 0], [0,-1]], dtype=complex)
 
    sY = complex(0,-1)*np.matmul(sZ,sX)
 
    paulis = [sI, sX, sY, sZ]
 
    paulis_label = ['I', 'X', 'Y', 'Z']
 
    paulis_term = {'I':sI, 'X':sX, 'Y':sY, 'Z':sZ}
 
    hamil = np.zeros((2**N, 2**N), dtype=complex) 
 

 
    for coeff, pauli_qb in zip(coeffs, pauli_qbs):
 
      term_str = ['I'] * N
 
      for term, index in zip(pauli_qb[0:][::2], pauli_qb[1:][::2]):
 
        term_str[int(index)] = term
 
      matrices = [paulis_term[pauli] for pauli in term_str]
 
      term_matrix = np.asarray(ft.reduce(np.kron, matrices[::-1]))
 
      hamil += coeff*term_matrix
 
    
 
    # Sanity Check
 
    if not np.allclose(hamil, hamil.conj().T):
 
        raise ValueError("The Hamiltonian formed is not Hermitian")
 

 
    return hamil
 
 
 
 
 
 

 

Test Hamiltonian

 

xxxxxxxxxx
 
 
 
 
 
 
 
H0 = np.array([
 
 
                [ 0.7056    ,  0.         , -0.       ,  0.        ],
 
                [ 0.        ,  -1.1246    , 0.182     ,  0.        ],
 
                [-0.        , 0.182       ,  0.4318   , -0.        ],
 
                [ 0.        ,  0.         , -0.       ,  0.888     ]
 
               ])
 
H0
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
array([[ 0.7056,  0.    , -0.    ,  0.    ],
 
 
       [ 0.    , -1.1246,  0.182 ,  0.    ],
 
       [-0.    ,  0.182 ,  0.4318, -0.    ],
 
       [ 0.    ,  0.    , -0.    ,  0.888 ]])
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
a, b , c = decompose_ham_to_pauli(H0)
 
 
a, b
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
(['I0I1', 'Z0I1', 'X0X1', 'Y0Y1', 'I0Z1', 'Z0Z1'],
 
 
 [0.2252, 0.3435, 0.091, 0.091, -0.43470000000000003, 0.5716])
 
 
 
 
 
 

 

$H_0=0.2252\times I_0{I}_1+0.3435\times Z_0-0.4347\times Z_1+0.91\times X_0{X}_1+0.91\times Y_0{Y}_1+0.5716\times Z_0{Z}_1$$H_{0} = 0.2252\times I_{0}I_{1} + 0.3435\times Z_{0} - 0.4347\times Z_{1} + 0.91\times X_{0}X_{1} + 0.91\times Y_{0}Y_{1} + 0.5716 \times Z_{0}Z_{1}$

 

xxxxxxxxxx
 
 
 
 
 
 
compose_ham_from_pauli(a, b)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
array([[ 0.7056+0.j,  0.    +0.j,  0.    +0.j,  0.    +0.j],
 
 
       [ 0.    +0.j, -1.1246+0.j,  0.182 +0.j,  0.    +0.j],
 
       [ 0.    +0.j,  0.182 +0.j,  0.4318+0.j,  0.    +0.j],
 
       [ 0.    +0.j,  0.    +0.j,  0.    +0.j,  0.888 +0.j]])
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
assert(np.isclose(compose_ham_from_pauli(a, b).real, H0).all())
 
 
 
 
 
 

 

Given Hamiltonian

xxxxxxxxxx
 
 
 
 
 
 
 
H1 = np.array([[1, 0, 0, 0], 
 
 
             [0, 0, -1, 0], 
 
             [0, -1, 0, 0], 
 
             [0, 0, 0, 1]])
 
H1
 
 
 
 
 
 

array([[ 1,  0,  0,  0],
 
 
       [ 0,  0, -1,  0],
 
       [ 0, -1,  0,  0],
 
       [ 0,  0,  0,  1]])
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 
a, b
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
(['I0I1', 'X0X1', 'Y0Y1', 'Z0Z1'], [0.5, -0.5, -0.5, 0.5])
 
 
 
 
 
 

 

$ H_{1} = 0.5\times (I_{0}I_{1} - X_{0}X_{1} - Y_{0}Y_{1} + Z_{0}Z_{1})$

xxxxxxxxxx
 
 
 
 
 
 
compose_ham_from_pauli(a, b)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
array([[ 1.+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, -1.+0.j,  0.+0.j,  0.+0.j],
 
       [ 0.+0.j,  0.+0.j,  0.+0.j,  1.+0.j]])
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
assert(np.isclose(compose_ham_from_pauli(a, b).real, H1).all())
 
 
 
 
 
 

 

Calculating the Expectation values

For an $n$$n$-qubit quantum state, measuring one qubit corresponds to projecting the quantum state onto one of two half-spaces determined by the unique eigenvalues of our measurement operator.

By convention, performing a computational basis measurement is equivalent to measuring in measuring Pauli Z, which gives us two eigenvectors $|0⟩$$|0\rangle$ and $|1⟩$$|1\rangle$, with corresponding eigenvalues $\pm 1$$\pm 1$. Therefore, doing a computation measurement of a qubit and obtaining $0$$0$ means the state of our qubit is in the $+1$$+1$ eigenstate of the $Z$$Z$ operator. In general, we can then use the definition of the tensor product to perform multi-qubit measurements as will be explained below.

For some qubit $|\psi ⟩$$|\psi\rangle$ we have $|\psi ⟩=\alpha |0⟩+\beta |1⟩$$|\psi\rangle = \alpha |0\rangle + \beta|1\rangle$, we can calculate expectation value of a single-qubit operator $M=\sum _{m}m|m⟩⟨m|$$M = \sum_m m|m\rangle\langle m|$, as $⟨M⟩=⟨\psi |\sum _{m}m|m⟩⟨m|\psi ⟩=\sum _{m}m⟨\psi |m^{†}m|\psi ⟩=\sum _{m}mp(m)=(+1)p(0)+(-1)p(1)=p(0)-p(1)$$\langle M \rangle = \langle \psi | \sum_{m} m|m\rangle\langle m| \psi\rangle = \sum_m m \langle \psi |m^{\dagger} m| \psi\rangle = \sum_m m\space p(m) = (+1)p(0) + (-1)p(1) = p(0)-p(1)$. Similarly, for a two-qubit operator $M$$M$ we will get $⟨M⟩=p(00)+p(11)-p(01)-p(10)$$\langle M \rangle = p(00) + p(11) - p(01) - p(10)$.

 

Change of Basis

Therefore, to measure a qubit, we can use any $2\times 2$$2\times2$ operator that is a unitary transformation of $Z$$Z$. That is, we could also use a measurement operator $M=U^{†}ZU$$M=U^{\dagger}ZU$, where $U$$U$ is an unitary operator and $M$$M$ gives two unique outcomes of a measurement in its $\pm 1$$\pm 1$ eigenvectors. Thereofre, to measure in a basis other than the computational (or $Z$$Z$), we need to find the corresponding unitary equivalence. Similar to the one-qubit case, all two-qubit Pauli-measurements can be written as $(U_1^{†}\otimes U_2^{†})[Z\otimes Z](U_1\otimes U_2)$$(U_{1}^{\dagger}\otimes U_{2}^{\dagger}) [Z\otimes Z] (U_{1}\otimes U_{2})$ for $2\times 2$$2\times 2$ unitary matrices $U_1$$U_1$ and $U_2$$U_2$.

For $X$$X$, and $Y$$Y$ basis measurements, these unitary equivalence are given below:

$X=(H{)}^{†}ZH$$X = (H)^{\dagger} Z H$ $Y=(H{S}^{†}{)}^{†}Z(H{S}^{†})$$Y = (HS^{\dagger})^{\dagger} Z (HS^{\dagger})$

This feature can be exploited to create circuits for both $X$$X$ basis and $Y$$Y$ basis measurements from the standard $Z$$Z$ basis. It is done by applying the corresponding unitary operation to the prepared quantum state and then performing a measurement in $Z$$Z$ basis.

For example: In order to do a $Y$$Y$ basis measurement, first apply $H{S}^{†}$$HS^{\dagger}$ to the quantum state and then do a normal $Z$$Z$ basis measurement.

Similar to the one-qubit case, all multi-qubit Pauli-measurements can be written as a generalization of above.

For example: For two-qubit $X$$X$ basis measurement, $X\otimes X=(H\otimes H)[Z\otimes Z](H\otimes H)$$X\otimes X = (H\otimes H)[Z\otimes Z] (H\otimes H)$.

 

Constructing Measurement Circuits

xxxxxxxxxx
 
 
 
 
 
 
 
def measure_circuit(basis, num_qubits):
 
 
    """
 
    Generate measurement circuit according to the Pauli observable
 
    string provided. 
 

 
    Args:
 
    basis (str): String representation of tensor products of Pauli 
 
    observables
 
    num_qubits (int): Number of qubits in the circuit
 

 
    Returns:
 
    measure_qc (QuantumCircuit): Measurement Circuit for the 
 
    corresponding basis.
 
    
 
    """
 

 
    basis_qb = re.findall(r'[A-Za-z]|-?\d+\.\d+|\d+', basis)
 
    basis = basis_qb[0:][::2] 
 
    qubit = list(map(int, basis_qb[1:][::2]))
 

 
    measure_qc = QuantumCircuit(num_qubits, num_qubits)
 
    for base, qb in zip(basis, qubit):
 
      if base == 'I' or base == 'Z':
 
        pass
 
      elif base == 'X':
 
        measure_qc.h(qb)
 
      elif base == 'Y':
 
        measure_qc.sdg(qb)
 
        measure_qc.h(qb)
 
      else:
 
        raise ValueError("Wrong Basis provided")
 
    
 
    for qb in range(num_qubits):
 
        measure_qc.measure(qb, qb)
 

 
    return measure_qc
 

 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
def calculate_expecation_val(circuit, basis, shots=2048, backend='qasm_simulator'):
 
 
    """
 
    Calculate expectation value for the measurement of a circuit in a 
 
    given basis.
 

 
    Args:
 
    circuit (QuantumCircuit): Circuit using which expectation value 
 
    will be calculated for a given basis.
 
    basis (str): String representation of tensor products of Pauli 
 
    observables.
 
    shots (int): Number of times measurements needed to be done for calculating
 
    probability.
 
    backend (str): Backend for running the circuit.
 

 
    Returns:
 
    exp (float): Expectation value for the measurement of a circuit 
 
    in a given basis.
 

 
    """
 
    exp_circuit = circuit + measure_circuit(basis, circuit.num_qubits)
 

 
    result = execute(exp_circuit, backend=Aer.get_backend(backend),
 
                    shots=shots).result()
 
    
 
    exp = 0.0
 
    for key, counts in result.get_counts().items():
 
        exp += (-1)**(int(key[0])+int(key[1])) * counts
 

 
    return exp/shots
 

 
 
 
 
 
 

 

Examples

Measurements in XX, YY, ZZ basis

xxxxxxxxxx
 
 
 
 
 
 
measure_circuit('X0X1', 2).draw()
 
 
 
 
 
 

     ┌───┐┌─┐   
q_0: ┤ H ├┤M├───
     ├───┤└╥┘┌─┐
q_1: ┤ H ├─╫─┤M├
     └───┘ ║ └╥┘
c: 2/══════╩══╩═
           0  1 

 

 

xxxxxxxxxx
 
 
 
 
 
 
measure_circuit('Y0Y1', 2).draw()
 
 
 
 
 
 

     ┌─────┐┌───┐┌─┐   
q_0: ┤ SDG ├┤ H ├┤M├───
     ├─────┤├───┤└╥┘┌─┐
q_1: ┤ SDG ├┤ H ├─╫─┤M├
     └─────┘└───┘ ║ └╥┘
c: 2/═════════════╩══╩═
                  0  1 

 

 

xxxxxxxxxx
 
 
 
 
 
 
measure_circuit('Z0Z1', 2).draw()
 
 
 
 
 
 

     ┌─┐   
q_0: ┤M├───
     └╥┘┌─┐
q_1: ─╫─┤M├
      ║ └╥┘
c: 2/═╩══╩═
      0  1 

 

Measurements in XY, YZ, ZX basis (To show flexibility of the function)

xxxxxxxxxx
 
 
 
 
 
 
measure_circuit('X0Y1', 2).draw()
 
 
 
 
 
 

      ┌───┐      ┌─┐   
q_0: ─┤ H ├──────┤M├───
     ┌┴───┴┐┌───┐└╥┘┌─┐
q_1: ┤ SDG ├┤ H ├─╫─┤M├
     └─────┘└───┘ ║ └╥┘
c: 2/═════════════╩══╩═
                  0  1 

 

 

xxxxxxxxxx
 
 
 
 
 
 
measure_circuit('Y0Z1', 2).draw()
 
 
 
 
 
 

     ┌─────┐┌───┐┌─┐
q_0: ┤ SDG ├┤ H ├┤M├
     └─┬─┬─┘└───┘└╥┘
q_1: ──┤M├────────╫─
       └╥┘        ║ 
c: 2/═══╩═════════╩═
        1         0 

 

 

xxxxxxxxxx
 
 
 
 
 
 
measure_circuit('Z0X1', 2).draw()
 
 
 
 
 
 

          ┌─┐   
q_0: ─────┤M├───
     ┌───┐└╥┘┌─┐
q_1: ┤ H ├─╫─┤M├
     └───┘ ║ └╥┘
c: 2/══════╩══╩═
           0  1 

 

Determining the Ansatz

 

Designing Ansatze

Ansatze are simply a parameterized quantum circuits (PQC), which play an essential role in the performance of many variational hybrid quantum-classical (HQC) algorithms. Major challenge while designing an asatz is to choose an effective template circuit that well represents the solution space while maintaining a low circuit depth and number of parameters. Here, we make a choice of two ansatze, one randomly and another inspired from the given hint.

 

Ansatz 1 (Random Choice)

xxxxxxxxxx
 
 
 
 
 
 
 
def ansatz1(params, num_qubits): 
 
 
    """
 
    Generate an templated ansatz with given parameters
 

 
    Args:
 
    params (array[float]): Parameters to initialize the parameterized unitary.
 
    num_qubits (int): Number of qubits in the circuit.
 
    
 
    Returns:
 
    ansatz (QuantumCircuit): Generated ansatz circuit
 

 
    """
 

 
    ansatz = QuantumCircuit(num_qubits, num_qubits)
 
    params = params.reshape(2,2)
 
    
 
    for idx in range(num_qubits):
 
      ansatz.rx(params[0][idx], idx)
 
    
 
    ansatz.cnot(0, 1)
 

 
    for idx in range(num_qubits):
 
      ansatz.rz(params[1][idx], idx)
 

 
    return ansatz
 

 
ansatz1(np.random.uniform(-np.pi, np.pi, (2,2)), 2).draw()
 
 
 
 
 
 

     ┌─────────────┐     ┌──────────────┐
q_0: ┤ RX(0.74906) ├──■──┤ RZ(-0.10279) ├
     ├─────────────┤┌─┴─┐├─────────────┬┘
q_1: ┤ RX(0.40493) ├┤ X ├┤ RZ(-2.3803) ├─
     └─────────────┘└───┘└─────────────┘ 
c: 2/════════════════════════════════════
                                         

 

Ansatz 2 (From Hint)

xxxxxxxxxx
 
 
 
 
 
 
 
def ansatz2(params, num_qubits): 
 
 
    """
 
    Generate an templated ansatz with given parameters
 

 
    Args:
 
    params (array[float]): Parameters to initialize the parameterized unitary.
 
    num_qubits (int): Number of qubits in the circuit.
 
    
 
    Returns:
 
    ansatz (QuantumCircuit): Generated ansatz circuit
 

 
    """
 
    params = np.array(params).reshape(1)
 
    ansatz = QuantumCircuit(num_qubits, num_qubits)
 
    ansatz.h(0)
 
    ansatz.cx(0, 1)
 
    ansatz.rx(params[0], 0)
 

 
    return ansatz
 

 
ansatz2([np.pi], 2).draw()
 
 
 
 
 
 

     ┌───┐     ┌────────┐
q_0: ┤ H ├──■──┤ RX(pi) ├
     └───┘┌─┴─┐└────────┘
q_1: ─────┤ X ├──────────
          └───┘          
c: 2/════════════════════
                         

 

Running Simulations

Variational Quantum Eigensolver (VQE) is based on Rayleigh-Ritz variational principle. To perform VQE, the first step is to enocde the problem into a Hermitian matrix $M$$M$, whose expectation value w.r.t a trial wave function $|\psi (\theta )⟩$$|\psi(\theta)\rangle$ which is yielded by an ansatz $U(\theta )$$U(\theta)$, i.e., nothing but a unitary paramterized by $\{\overrightarrow{\theta }\}$$\{\vec{\theta}\}$.

Here, we optimize these parameters $\{\overrightarrow{\theta }\}$$\{\vec{\theta}\}$ using scipy.optimize library based on a cost function $C(\theta )$$C(\theta)$, which is nothing but the calculated expectation value of Hamiltonian, $H_1$$H_1$ for $|\psi (\theta )⟩$$|\psi(\theta)\rangle$: $C(\theta )=⟨\psi (\theta )|H_1|\psi (\theta )⟩$$C(\theta) = \langle \psi (\theta) | H_{1} | \psi (\theta) \rangle$

xxxxxxxxxx
 
 
 
 
 
 
 
def vqe(params, meas_basis, coeffs, circuit, num_qubits, shots=2048):
 
 
    """
 
        Return the calculated energy scalar for a given asnatz and 
 
        decomposed Hamiltoian. 
 

 
        Args:
 
        params (matrix(np.array)): Parameters for initializing the ansatz.
 
        meas_basis (list[str]): String representation of measurement basis, i.e.
 
        the decomposed pauli term with their corresponding qubits.
 
        coeffs (vector(np.array)): Coefficients for the decomposed Pauli Term
 
        circuit (QuantumCircuit): Template for Ansatz circuit 
 
        num_qubits (int): Number of qubits in the given asatze.
 
        shots (int): Number of shots to get the probability distribution.
 

 
        Return:
 
        energy (float): Expectation value of the Hamiltonian whose 
 
        decomposition was provided.
 
    """
 
    
 
    N = num_qubits
 
    circuit = circuit(params, num_qubits)
 
    energy = 0
 

 
    for basis, coeff in zip(meas_basis, coeffs):
 
      if basis.count('I') != N:
 
        energy += coeff*calculate_expecation_val(circuit, basis, shots)
 

 
    energy += 0.5
 

 
    return energy
 
 
 
 
 
 

 

Ansatz 1 (Random Choice)

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 

 
params = np.random.uniform(-np.pi, np.pi, (2,2))
 
func = ft.partial(vqe, meas_basis=a, coeffs=b, circuit=ansatz1, 
 
                  num_qubits=2, shots=2048)
 

 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
res = sp.optimize.minimize(func, params, method='Powell')
 
 
res
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
   direc: array([[1., 0., 0., 0.],
 
 
       [0., 1., 0., 0.],
 
       [0., 0., 1., 0.],
 
       [0., 0., 0., 1.]])
 
     fun: array(-1.)
 
 message: 'Optimization terminated successfully.'
 
    nfev: 190
 
     nit: 3
 
  status: 0
 
 success: True
 
       x: array([ 4.7130171 , -3.15044759, -1.59558296, -0.02529809])
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
assert(float(res.fun) == min(np.linalg.eig(H1)[0])) #Sanity Check
 
 
 
 
 
 

 

Visualizing the Result

After looking at multiple results' res.x, we realize that first and second parameters always have the values $\pm n\pi /2$$\pm n \pi/2$ and $\pm \pi$$\pm\pi$ respectively. So, while keeping them fixed, we do a parameters search-scan in the range ${\theta }_{3,4}\in [0,2\pi ]$$\theta_{3,4} \in [0,2\pi]$ and visualize how they affect our expectation value.

xxxxxxxxxx
 
 
 
 
 
 
 
def energy_expectation(x, y):
 
 
  """ Returns meshgrid values for plotting """
 
  energy = np.zeros(x.shape)
 
  for idx, thetas in enumerate(x):
 
    for ind, theta1 in enumerate(thetas):
 
      params = np.array([np.pi/2, np.pi, theta1, y[idx][ind]])
 
      energy[idx][ind] = func(params)
 
  return energy
 

 
theta1 = np.linspace(0.0, 2*np.pi, 200)
 
theta2 = np.linspace(0.0, 2*np.pi, 200)
 

 
X, Y = np.meshgrid(theta1, theta2)
 
Z = energy_expectation(X, Y)
 

 
fig = plt.figure()
 
ax = plt.axes(projection='3d')
 
ax.contour3D(X, Y, Z, 50, cmap='summer')
 
ax.set_xlabel('theta_1')
 
ax.set_ylabel('theta_2')
 
ax.set_zlabel('Expectation of H_1');
 
plt.show()
 

 
 
 
 
 
 

Ansatz 2 (From Hint)

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 
params = np.random.uniform(-np.pi, np.pi, 1)
 
func = ft.partial(vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, shots=2048)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
res = sp.optimize.minimize(func, params, method='Powell')
 
 
res
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
   direc: array([[1.]])
 
 
     fun: array(-1.)
 
 message: 'Optimization terminated successfully.'
 
    nfev: 26
 
     nit: 2
 
  status: 0
 
 success: True
 
       x: array([3.13590442])
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
assert(float(res.fun) == min(np.linalg.eig(H1)[0])) #Sanity Check
 
 
 
 
 
 

 

Visualizing the Result

We scan over the paramter $\theta \in [0,2\pi ]$$\theta \in [0,2\pi]$, to visualize the change in expectation value for the ansatz2, which is inspired from the hint given.

xxxxxxxxxx
 
 
 
 
 
 
 
thetas = np.linspace(0.0, 2*np.pi, 500)
 
 
energy = np.zeros(len(thetas))
 

 
for idx, theta in enumerate(thetas):
 
    energy[idx] = func([theta])    
 

 
fig = plt.figure()
 
plt.plot(thetas, energy)
 
plt.title('Expectation value w.r.t angle theta')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Theta (radians)')
 
indices = [idx for idx, x in enumerate(energy) if x <= -1.0]
 
xmin = thetas[indices[len(indices)//2]]
 
ymin = energy.min()
 
text= "Theta={:.3f}, \n Energy={:.3f}".format(xmin, ymin)
 
plt.plot(xmin,ymin,'x')
 
plt.annotate(text, xy=(xmin+0.75, ymin))
 
fig.show()
 

 
 
 
 
 
 

 

Number of Measurements (or Shots)

In order to obtain results from a quantum computer, one needs to perform sampling. This allows us to obtain the probability distribution, and hence know the most probable results from the measurement. Here, we see the change of expectation value of $H_1$$H_1$ with the number of shots, i.e., number of measurements for ansatz 2.

xxxxxxxxxx
 
 
 
 
 
 
 
thetas = np.linspace(0.0, 2*np.pi, 500)
 
 
shots = [128, 256, 512, 1024, 2048]
 
shot_energy = np.zeros((len(shots), len(thetas)))
 

 
for ind, shot in enumerate(shots):
 
    for idx, theta in enumerate(thetas):
 
        func = ft.partial(vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                      num_qubits=2, shots=shot)
 
        shot_energy[ind][idx] = func([theta])    
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
plt.figure()
 
 
for ind, shot in enumerate(shots):
 
    plt.plot(thetas, shot_energy[ind], label='Shots=' + str(shot))
 
    
 
plt.title('Expectation value w.r.t angle theta and shots')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Theta (radians)')
 
plt.legend()
 
plt.show()
 
 
 
 
 
 

 

Performance of Classical Optimizers

Successful application of hybrid-quantum classical algorithms, with the classical step involving an optimizer, on current hardware, requires the classical optimizer to be noise-aware. In order to study this, we study the change of expectation value of $H_1$$H_1$ with different classical optimizers for ansatz 2, first in absence of noise, then later in the presence of noise.

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 
params = np.random.uniform(-np.pi, np.pi, 1)
 
func = ft.partial(vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, shots=2048)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
optimizers = ['Nelder-Mead', 'L-BFGS-B', 'TNC', 'Powell', 'COBYLA', 'SLSQP']
 
 
opt_val = []
 
opt_param = []
 
opt_feval = []
 
opt_iter = []
 
init_params = params
 

 
for opt in optimizers:
 

 
    params = init_params
 
    res = sp.optimize.minimize(func, params, method=opt)
 

 
    opt_val.append(float(res.fun))
 
    opt_param.append(float(res.x))
 
    opt_feval.append(int(res.nfev))
 
    try:
 
        opt_iter.append(int(res.nit))
 
    except:
 
        opt_iter.append(1)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers, opt_val, color=['black', 'red', 'green', 'blue', 'yellow', 'orange'])
 
plt.title('Minimum expectation value w.r.t different optimizers')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers[:-1], opt_param[:-1], color=['black', 'red', 'green', 'blue', 'yellow'])
 
plt.title('Optimal parameter w.r.t different optimizers')
 
plt.ylabel('Value of theta')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers[:-1], opt_feval[:-1], color=['black', 'red', 'green', 'blue', 'yellow'])
 
plt.title('Function evaluations w.r.t different optimizers')
 
plt.ylabel('Number of function evaluations$')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers[:-1], opt_iter[:-1], color=['black', 'red', 'green', 'blue', 'yellow'])
 
plt.title('Number of iterations w.r.t different optimizers')
 
plt.ylabel('Number of iterations')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

Noisy Simulation

Computational power of NISQ processors suffers a lot from the limited capabilities of their physical qubits. This is essentially due to present of decoherence, limited connectivity, and absence of error-correction. Hence, it is essential to see how our ansatz would perform on a real device. Here, we replicate the noise model from a real device, and then test the perfrormance of our VQE module.

xxxxxxxxxx
 
 
 
 
 
 
 
# Noise Model Preparation
 
 
provider = IBMQ.load_account()
 
noisy_backend = provider.get_backend('ibmq_vigo')
 
coupling_map = noisy_backend.configuration().coupling_map
 
noise_model = providers.aer.noise.NoiseModel.from_backend(noisy_backend)
 
basis_gates = noise_model.basis_gates
 
 
 
 
 
 

Qubit coupling map of this real device is the following:

xxxxxxxxxx
 
 
 
 
 
 
qiskit.transpiler.CouplingMap(noisy_backend.configuration().coupling_map).draw()
 
 
 
 
 
 

 

 

xxxxxxxxxx
 
 
 
 
 
 
noise_model
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
NoiseModel:
 
 
  Basis gates: ['cx', 'id', 'u2', 'u3']
 
  Instructions with noise: ['id', 'measure', 'cx', 'u2', 'u3']
 
  Qubits with noise: [0, 1, 2, 3, 4]
 
  Specific qubit errors: [('id', [0]), ('id', [1]), ('id', [2]), ('id', [3]), ('id', [4]), ('u2', [0]), ('u2', [1]), ('u2', [2]), ('u2', [3]), ('u2', [4]), ('u3', [0]), ('u3', [1]), ('u3', [2]), ('u3', [3]), ('u3', [4]), ('cx', [0, 1]), ('cx', [1, 0]), ('cx', [1, 2]), ('cx', [1, 3]), ('cx', [2, 1]), ('cx', [3, 1]), ('cx', [3, 4]), ('cx', [4, 3]), ('measure', [0]), ('measure', [1]), ('measure', [2]), ('measure', [3]), ('measure', [4])]
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 

 
 
def calculate_noisy_expecation_val(circuit, basis, mitigated=False, shots=2048, backend='qasm_simulator'):
 
    """
 
    Calculate expectation value for the measurement of a circuit in a 
 
    given basis in presence of noise.
 

 
    Args:
 
    circuit (QuantumCircuit): Circuit using which expectation value 
 
    will be calculated for a given basis.
 
    basis (str): String representation of tensor products of Pauli 
 
    observables.
 
    mitigated (bool): Whether mitigation has to be performed or not.
 
    shots (int): Number of times measurements needed to be done for calculating
 
    probability.
 
    backend (str): Backend for running the circuit.
 

 
    Returns:
 
    exp (float): Expectation value for the measurement of a circuit 
 
    in a given basis.
 

 
    """
 
    exp_circuit = circuit + measure_circuit(basis, circuit.num_qubits)
 
    
 
    result = execute(exp_circuit, backend=Aer.get_backend(backend),
 
                     shots=shots, coupling_map=coupling_map,
 
                     basis_gates=basis_gates,
 
                     noise_model=noise_model).result()
 

 
    if mitigated:
 
        meas_calibs, state_labels = complete_meas_cal(qr=exp_circuit.qregs[0], cr=exp_circuit.cregs[0])
 
        cal_results = qiskit.execute(meas_calibs, backend=Aer.get_backend(backend),
 
                                     coupling_map=coupling_map, basis_gates=basis_gates,
 
                                     shots=shots, noise_model=noise_model).result()
 
        meas_fitter = CompleteMeasFitter(cal_results, state_labels)
 
        meas_filter = meas_fitter.filter
 
        result = meas_filter.apply(result)
 
        
 
    exp = 0.0
 
    for key, counts in result.get_counts().items():
 
        exp += (-1)**(int(key[0])+int(key[1])) * counts
 

 
    return exp/shots
 

 

 
def noisy_vqe(params, meas_basis, coeffs, circuit, num_qubits, mitigated=False, shots=2048):
 
    """
 
        Return the calculated energy scalar for a given asnatz and 
 
        decomposed Hamiltoian in presence of noise. 
 

 
        Args:
 
        params (matrix(np.array)): Parameters for initializing the ansatz.
 
        meas_basis (list[str]): String representation of measurement basis, i.e.
 
        the decomposed pauli term with their corresponding qubits.
 
        coeffs (vector(np.array)): Coefficients for the decomposed Pauli Term
 
        circuit (QuantumCircuit): Template for Ansatz circuit 
 
        num_qubits (int): Number of qubits in the given asatze.
 
        mitigated (bool): Whether mitigation has to be performed or not.
 
        shots (int): Number of shots to get the probability distribution.
 

 
        Return:
 
        energy (float): Expectation value of the Hamiltonian whose 
 
        decomposition was provided.
 
    """
 
    
 
    N = num_qubits
 
    circuit = circuit(params, num_qubits)
 
    energy = 0
 

 
    for basis, coeff in zip(meas_basis, coeffs):
 
      if basis.count('I') != N:
 
        energy += coeff*calculate_noisy_expecation_val(circuit, basis, mitigated, shots)
 

 
    energy += 0.5
 

 
    return energy
 
 
 
 
 
 

 

Ansatz 1 (Random Choice)

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 

 
params = np.random.uniform(-np.pi, np.pi, (2,2))
 
func = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz1, 
 
                  num_qubits=2, mitigated=False, shots=8192)
 

 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
res = sp.optimize.minimize(func, params, method='Powell')
 
 
res
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
   direc: array([[1., 0., 0., 0.],
 
 
       [0., 1., 0., 0.],
 
       [0., 0., 1., 0.],
 
       [0., 0., 0., 1.]])
 
     fun: array(-0.87451172)
 
 message: 'Optimization terminated successfully.'
 
    nfev: 266
 
     nit: 3
 
  status: 0
 
 success: True
 
       x: array([1.64212384, 3.14620703, 2.15383843, 0.53103981])
 
 
 
 
 
 

 

 

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 

 
params = np.random.uniform(-np.pi, np.pi, (2,2))
 
func = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz1, 
 
                  num_qubits=2, mitigated=True, shots=8192)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
res = sp.optimize.minimize(func, params, method='Powell')
 
 
res
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
   direc: array([[1., 0., 0., 0.],
 
 
       [0., 1., 0., 0.],
 
       [0., 0., 1., 0.],
 
       [0., 0., 0., 1.]])
 
     fun: -0.9619728377818637
 
 message: 'Optimization terminated successfully.'
 
    nfev: 257
 
     nit: 3
 
  status: 0
 
 success: True
 
       x: array([ 1.58409174, -3.11397559, -0.50929574, -1.96419365])
 
 
 
 
 
 

 

Ansatz 2 (From Hint)

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 
params = np.random.uniform(-np.pi, np.pi, 1)
 
func = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, mitigated=False, shots=8192)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
res = sp.optimize.minimize(func, params, method='Powell')
 
 
res
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
   direc: array([[-0.1127984]])
 
 
     fun: array(-0.88122559)
 
 message: 'Optimization terminated successfully.'
 
    nfev: 88
 
     nit: 5
 
  status: 0
 
 success: True
 
       x: array([3.17376448])
 
 
 
 
 
 

 

 

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 
params = np.random.uniform(-np.pi, np.pi, 1)
 
func = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, mitigated=True, shots=8192)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
res = sp.optimize.minimize(func, params, method='Powell')
 
 
res
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
   direc: array([[3.01441921e-09]])
 
 
     fun: -0.9702586386528804
 
 message: 'Optimization terminated successfully.'
 
    nfev: 120
 
     nit: 5
 
  status: 0
 
 success: True
 
       x: array([-3.12601157])
 
 
 
 
 
 

 

Visualizing the Result (Noisy v/s Mitigated v/s Ideal)

xxxxxxxxxx
 
 
 
 
 
 
 
thetas = np.linspace(0.0, 2*np.pi, 500)
 
 
energy1 = np.zeros(len(thetas))
 
energy2 = np.zeros(len(thetas))
 
func1 = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, mitigated=False, shots=8192)
 
func2 = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, mitigated=True, shots=8192)
 

 
for idx, theta in enumerate(thetas):
 
    energy1[idx] = func1([theta])    
 
    energy2[idx] = func2([theta])    
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
plt.plot(thetas, energy, label='simulation')
 
 
plt.plot(thetas, energy1, label='noisy')
 
plt.plot(thetas, energy2, label='mitigated')
 
plt.title('Expectation value w.r.t angle theta')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Theta (radians)')
 
plt.legend()
 
plt.show()
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
<IPython.core.display.Javascript object>
 
 
 
 
 
 

 

 

Number of Measurements (or Shots)

Let's see the change of expectation value of $H_1$$H_1$ with the number of shots for ansatz 2.

xxxxxxxxxx
 
 
 
 
 
 
 
thetas = np.linspace(0.0, 2*np.pi, 500)
 
 
shots = [128, 256, 512, 1024, 2048]
 
noisy_shot_energy = np.zeros((len(shots), len(thetas)))
 

 
for ind, shot in enumerate(shots):
 
    for idx, theta in enumerate(thetas):
 
        func = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                      num_qubits=2, shots=shot)
 
        noisy_shot_energy[ind][idx] = func([theta])    
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
plt.figure()
 
 
for ind, shot in enumerate(shots):
 
    plt.plot(thetas, noisy_shot_energy[ind], label='Shots=' + str(shot))
 
    
 
plt.title('Expectation value w.r.t angle theta and shots')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Theta (radians)')
 
plt.legend()
 
plt.show()
 
 
 
 
 
 

 

Peformance of Classical Optimizers

Let's see the change of expectation value of $H_1$$H_1$ with different classical optimizers for ansatz 2.

xxxxxxxxxx
 
 
 
 
 
 
 
a, b, c = decompose_ham_to_pauli(H1)
 
 
params = np.random.uniform(-np.pi, np.pi, 1)
 
func = ft.partial(noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, mitigated=True, shots=8192)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
optimizers = ['Nelder-Mead', 'L-BFGS-B', 'TNC', 'Powell', 'COBYLA', 'SLSQP']
 
 
opt_val = []
 
opt_param = []
 
opt_feval = []
 
opt_iter = []
 

 
init_params = params
 

 
for opt in optimizers:
 

 
    params = init_params
 
    res = sp.optimize.minimize(func, params, method=opt)
 

 
    opt_val.append(float(res.fun))
 
    opt_param.append(float(res.x))
 
    opt_feval.append(int(res.nfev))
 
    try:
 
        opt_iter.append(int(res.nit))
 
    except:
 
        opt_iter.append(1)
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers, opt_val, color=['black', 'red', 'green', 'blue', 'yellow', 'orange'])
 
plt.title('Minimum expectation value w.r.t different optimizers')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers[:-1], opt_param[:-1], color=['black', 'red', 'green', 'blue', 'yellow'])
 
plt.title('Optimal parameter w.r.t different optimizers')
 
plt.ylabel('Value of theta')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers[:-1], opt_feval[:-1], color=['black', 'red', 'green', 'blue', 'yellow'])
 
plt.title('Function evaluations w.r.t different optimizers')
 
plt.ylabel('Number of function evaluations')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
fig = plt.figure()
 
 
plt.bar(optimizers[:-1], opt_iter[:-1], color=['black', 'red', 'green', 'blue', 'yellow'])
 
plt.title('Number of iterations w.r.t different optimizers')
 
plt.ylabel('Number of iterations')
 
plt.xlabel('Optimizer')
 
plt.xticks(rotation=5)
 
fig.show()
 
 
 
 
 
 

 

Comparing the performance of VQE in presence of errors due to Readout, Thermal Relaxation, Amplitude-Phase damping and Depolarization

Readout

Describes classical readout errors

xxxxxxxxxx
 
 
 
 
 
 
 
# Measurement miss-assignement probabilities
 
 
p0given1 = 0.1
 
p1given0 = 0.05
 
noise_readout = ReadoutError([[1 - p1given0, p1given0], [p0given1, 1 - p0given1]])
 
noise_readout
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
ReadoutError([[0.95 0.05]
 
 
 [0.1  0.9 ]])
 
 
 
 
 
 

 

Thermal Relaxation

Single qubit thermal relaxation error is characterized by relaxation time constants $T_1$$T_1$, $T_2$$T_2$, and the gate time $t$$t$.

xxxxxxxxxx
 
 
 
 
 
 
 
num_qubits = 2
 
 

 
T1s = np.random.normal(50e3, 10e3, num_qubits) # Sampled from normal distribution
 
T2s = np.random.normal(80e3, 20e3, num_qubits)  # Sampled from normal distribution
 
T2s = np.array([min(T2s[j], 2 * T1s[j]) for j in range(num_qubits)]) # Ensuring T2s < 2*T1s
 

 
# Instruction times (in nanoseconds)
 
time_u1 = 0   # virtual gate
 
time_u2 = 50  # (single X90 pulse)
 
time_u3 = 100 # (two X90 pulses)
 
time_cx = 300
 
time_measure = 1000 # 1 microsecond
 

 
# Add depolarizing error to all single qubit u1, u2, u3, cx, measure gates
 
errors_thermal_u1  = [thermal_relaxation_error(t1, t2, time_u1)
 
              for t1, t2 in zip(T1s, T2s)]
 
errors_thermal_u2  = [thermal_relaxation_error(t1, t2, time_u2)
 
              for t1, t2 in zip(T1s, T2s)]
 
errors_thermal_u3  = [thermal_relaxation_error(t1, t2, time_u3)
 
              for t1, t2 in zip(T1s, T2s)]
 
errors_thermal_cx = [[thermal_relaxation_error(t1a, t2a, time_cx).expand(
 
             thermal_relaxation_error(t1b, t2b, time_cx))
 
              for t1a, t2a in zip(T1s, T2s)]
 
               for t1b, t2b in zip(T1s, T2s)]
 
errors_thermal_measure = [thermal_relaxation_error(t1, t2, time_measure)
 
                  for t1, t2 in zip(T1s, T2s)]
 

 
# Add errors to noise model
 
noise_thermal = NoiseModel()
 
for j in range(num_qubits):
 
    noise_thermal.add_quantum_error(errors_thermal_u1[j], "u1", [j])
 
    noise_thermal.add_quantum_error(errors_thermal_u2[j], "u2", [j])
 
    noise_thermal.add_quantum_error(errors_thermal_u3[j], "u3", [j])
 
    for k in range(num_qubits):
 
        noise_thermal.add_quantum_error(errors_thermal_cx[j][k], "cx", [j, k])
 
    noise_thermal.add_quantum_error(errors_thermal_measure[j], "measure", [j])
 

 
noise_thermal
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
NoiseModel:
 
 
  Basis gates: ['cx', 'id', 'u2', 'u3']
 
  Instructions with noise: ['u2', 'u3', 'cx', 'measure']
 
  Qubits with noise: [0, 1]
 
  Specific qubit errors: [('u2', [0]), ('u2', [1]), ('u3', [0]), ('u3', [1]), ('cx', [0, 0]), ('cx', [0, 1]), ('cx', [1, 0]), ('cx', [1, 1]), ('measure', [0]), ('measure', [1])]
 
 
 
 
 
 

 

Amplitude-Phase Damping

Single-qubit generalized combined phase and amplitude damping error is given by an amplitude damping parameter $\lambda$$\lambda$, and a phase damping parameter $\gamma$$\gamma$.

xxxxxxxxxx
 
 
 
 
 
 
 
lamb = 0.05
 
 
gamma = 0.1
 

 
errors_phase_amplitude_1 =  phase_amplitude_damping_error(lamb, gamma)
 
errors_phase_amplitude_2 =  phase_amplitude_damping_error(lamb, gamma).expand(phase_amplitude_damping_error(lamb, gamma))
 

 
# Add errors to noise model
 
noise_phase_amplitude = NoiseModel()
 
noise_phase_amplitude.add_all_qubit_quantum_error(errors_phase_amplitude_1, ['u1', 'u2', 'u3', 'measure'])
 
noise_phase_amplitude.add_all_qubit_quantum_error(errors_phase_amplitude_2, ['cx'])
 
noise_phase_amplitude
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
NoiseModel:
 
 
  Basis gates: ['cx', 'id', 'u1', 'u2', 'u3']
 
  Instructions with noise: ['u2', 'u3', 'u1', 'cx', 'measure']
 
  All-qubits errors: ['u1', 'u2', 'u3', 'measure', 'cx']
 
 
 
 
 
 

 

Depolarization

N-qubit depolarizing error is given by a depolarization probability $p$$p$.

xxxxxxxxxx
 
 
 
 
 
 
 
# Add depolarizing error to all single qubit u1, u2, u3, cx, measure gates
 
 
error_depol_1 = depolarizing_error(0.05, 1)
 
error_depol_2 = depolarizing_error(0.1, 2)
 
error_depol_3 = depolarizing_error(0.1, 1)
 

 
# Add errors to noise model
 
noise_depol = NoiseModel()
 
noise_depol.add_all_qubit_quantum_error(error_depol_1, ['u1', 'u2', 'u3'])
 
noise_depol.add_all_qubit_quantum_error(error_depol_2, ['cx'])
 
noise_depol.add_all_qubit_quantum_error(error_depol_3, ['measure'])
 
noise_depol
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
NoiseModel:
 
 
  Basis gates: ['cx', 'id', 'u1', 'u2', 'u3']
 
  Instructions with noise: ['u2', 'u3', 'u1', 'cx', 'measure']
 
  All-qubits errors: ['u1', 'u2', 'u3', 'cx', 'measure']
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
def calculate_restricted_noisy_expecation_val(circuit, basis, noise_model, shots=2048, backend='qasm_simulator'):
 
 
    """
 
    Calculate expectation value for the measurement of a circuit in a 
 
    given basis in presence of a particular kind of noise.
 

 
    Args:
 
    circuit (QuantumCircuit): Circuit using which expectation value 
 
    will be calculated for a given basis.
 
    basis (str): String representation of tensor products of Pauli 
 
    observables.
 
    noise_model (NoiseModel): Noise Model for execution
 
    shots (int): Number of times measurements needed to be done for calculating
 
    probability.
 
    backend (str): Backend for running the circuit.
 

 
    Returns:
 
    exp (float): Expectation value for the measurement of a circuit 
 
    in a given basis.
 

 
    """
 
    exp_circuit = circuit + measure_circuit(basis, circuit.num_qubits)
 
    
 
    result = execute(exp_circuit, backend=Aer.get_backend(backend),
 
                     shots=shots, noise_model=noise_model).result()
 
        
 
    exp = 0.0
 
    for key, counts in result.get_counts().items():
 
        exp += (-1)**(int(key[0])+int(key[1])) * counts
 

 
    return exp/shots
 

 

 
def restricted_noisy_vqe(params, meas_basis, coeffs, circuit, num_qubits, noise_model, shots=2048):
 
    """
 
        Return the calculated energy scalar for a given asnatz and 
 
        decomposed Hamiltoian in presence of a particular kind of noise.
 

 
        Args:
 
        params (matrix(np.array)): Parameters for initializing the ansatz.
 
        meas_basis (list[str]): String representation of measurement basis, i.e.
 
        the decomposed pauli term with their corresponding qubits.
 
        coeffs (vector(np.array)): Coefficients for the decomposed Pauli Term
 
        circuit (QuantumCircuit): Template for Ansatz circuit 
 
        num_qubits (int): Number of qubits in the given asatze.
 
        noise_model (NoiseModel): Noise Model for execution
 
        shots (int): Number of shots to get the probability distribution.
 

 
        Return:
 
        energy (float): Expectation value of the Hamiltonian whose 
 
        decomposition was provided.
 
    """
 
    
 
    N = num_qubits
 
    circuit = circuit(params, num_qubits)
 
    energy = 0
 

 
    for basis, coeff in zip(meas_basis, coeffs):
 
      if basis.count('I') != N:
 
        energy += coeff*calculate_restricted_noisy_expecation_val(circuit, basis, noise_model, shots)
 

 
    energy += 0.5
 

 
    return energy
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
noise_models = [noise_readout, noise_thermal, noise_phase_amplitude, noise_depol]
 
 

 
thetas = np.linspace(0.0, 2*np.pi, 500)
 
noisy_energy = np.zeros((len(noise_models), len(thetas)))
 

 
for ind, model in enumerate(noise_models):
 
    func3 = ft.partial(restricted_noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, noise_model=model, shots=4096)
 
    for idx, theta in enumerate(thetas):
 
        noisy_energy[ind][idx] = func3([theta])
 

 

 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
noise_models_name = ['Readout Error (0.1, 0.05)', 'Thermal Relaxation (5e-5, 8e-5)', 
 
 
                     'Amplitude-Phase Damping (0.05, 0.1)', 'Depolarization (0.05)']
 

 
plt.figure()
 
for ind, model in enumerate(noise_models):
 
    plt.plot(thetas, noisy_energy[ind], label=noise_models_name[ind])
 
    
 
plt.title('Expectation value w.r.t angle theta and shots')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Theta (radians)')
 
plt.legend()
 
plt.show()
 
 
 
 
 
 

We see from the above plot that the presence of noise due to depolarization has affected our VQE results the most. This can be explained via its definiton that it is the uniform contraction of the Bloch sphere, parameterized by $p$$p$. So, we now do a quick comparision between the effects of depolarization on single qubit gates v/s two-qubit gates v/s measure gates.

xxxxxxxxxx
 
 
 
 
 
 
 
# Add depolarizing error to all single qubit u1, u2, u3, cx, measure gates
 
 
error_depol_1 = depolarizing_error(0.05, 1)
 
error_depol_2 = depolarizing_error(0.1, 2)
 
error_depol_3 = depolarizing_error(0.1, 1)
 

 
# Add errors to noise model
 
noise_depol_single = NoiseModel()
 
noise_depol_two = NoiseModel()
 
noise_depol_measure = NoiseModel()
 
noise_depol_single.add_all_qubit_quantum_error(error_depol_1, ['u1', 'u2', 'u3'])
 
noise_depol_two.add_all_qubit_quantum_error(error_depol_2, ['cx'])
 
noise_depol_measure.add_all_qubit_quantum_error(error_depol_3, ['measure'])
 
noise_depol_single, noise_depol_two, noise_depol_measure
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
(NoiseModel:
 
 
   Basis gates: ['cx', 'id', 'u1', 'u2', 'u3']
 
   Instructions with noise: ['u2', 'u3', 'u1']
 
   All-qubits errors: ['u1', 'u2', 'u3'],
 
 NoiseModel:
 
   Basis gates: ['cx', 'id', 'u3']
 
   Instructions with noise: ['cx']
 
   All-qubits errors: ['cx'],
 
 NoiseModel:
 
   Basis gates: ['cx', 'id', 'u3']
 
   Instructions with noise: ['measure']
 
   All-qubits errors: ['measure'])
 
 
 
 
 
 

 

xxxxxxxxxx
 
 
 
 
 
 
 
depol_noise_models = [noise_depol_single, noise_depol_two, noise_depol_measue]
 
 

 
thetas = np.linspace(0.0, 2*np.pi, 500)
 
depol_noisy_energy = np.zeros((len(depol_noise_models), len(thetas)))
 

 
for ind, model in enumerate(depol_noise_models):
 
    func3 = ft.partial(restricted_noisy_vqe, meas_basis=a, coeffs=b, circuit=ansatz2, 
 
                  num_qubits=2, noise_model=model, shots=4096)
 
    for idx, theta in enumerate(thetas):
 
        depol_noisy_energy[ind][idx] = func3([theta])
 

 

 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
 
depol_noise_models_name = ['One-Qubit Gates Depolarization', ' Two-Qubit Gates Depolarization', 
 
 
                           'Measure Gate Depolarization']
 

 
plt.figure()
 
for ind, model in enumerate(depol_noise_models):
 
    plt.plot(thetas, depol_noisy_energy[ind], label=depol_noise_models_name[ind])
 
    
 
plt.title('Expectation value w.r.t angle theta and shots')
 
plt.ylabel('Expectation value of $H_1$')
 
plt.xlabel('Theta (radians)')
 
plt.legend()
 
plt.show()
 
 
 
 
 
 

 

Excited States of the Hamiltonian

 

One of the way to find the $k^{th}$$k^{th}$ exicted states of a Hamiltonian would be to update the Hamiltonain $H$$H$ to $H_k$$H_{k}$ as follows: $H_k=H+\sum _{i=0}^{k-1}{\beta }_i|i⟩⟨i|$$H_{k} = H + \sum_{i=0}^{k-1}\beta_{i}|i\rangle\langle i|$ For example: To find the first excited state, we find $H_1=H+{\beta }_0|g⟩⟨g|$$H_1 = H + \beta_0 |g\rangle\langle g|$, where $g$$g$ is the ground state of the Hamiltonian. This works because of spectral decomposition.

This can be used to update our cost function, $C_k$$C_{k}$, as follows: $C({\theta }_k)=⟨\psi ({\theta }_k)|H|\psi ({\theta }_k)⟩+\sum _{i=0}^{k-1}{\beta }_i|⟨\psi ({\theta }_k)|\psi ({\theta }_i)⟩{|}^2$$C(\theta_{k}) = \langle \psi(\theta_k)| H | \psi (\theta_k) \rangle + \sum_{i=0}^{k-1} \beta_{i} |\langle \psi(\theta_k) | \psi(\theta_{i})\rangle|^{2}$

We'd already calculated the first term in VQE, but to calculate the next summation term, we refer to [2]. It suggests us to rewrite the overlap term $|⟨\psi ({\theta }_k)|\psi ({\theta }_i)⟩{|}^2$$|\langle \psi(\theta_k)|\psi(\theta_{i})\rangle|^{2}$ as $|⟨00|U({\theta }_k{)}^{†}U({\theta }_i)|00⟩{|}^2$$|\langle00|U(\theta_{k})^{\dagger}U(\theta_{i})|00\rangle |^{2}$. Therefore, we can prepare the state $U({\theta }_k{)}^{†}U({\theta }_i)|00⟩$$U(\theta_{k})^{\dagger}U(\theta_{i})|00\rangle$ using the trial state preparation circuit for current state and $i^{th}$$i^{th}$ previously-computed state.

This technique is known as Variational Quantum Deflation and requires the same number of qubits as Variational Quantum Eigensolver and around twice the circuit depth.

xxxxxxxxxx
 
 
 
 
 
 
 
def vqd_cost(thetak, thetai, func, circuit, num_qubits, beta, shots=2048):
 
 
    """  Returns cost for VQD """
 
    
 
    energy = func(thetak)
 

 
    for idx, theta in enumerate(thetai):
 
        ansatz = circuit(theta, num_qubits) + circuit(thetak, num_qubits).inverse()
 
        ansatz.measure_all()
 
        backend = Aer.get_backend('qasm_simulator')
 
        result = execute(ansatz, backend, shots=shots).result()
 
        energy += beta[idx]*result.get_counts().get('00 00', 0)/shots
 

 
    return energy
 

 
def vqd(hamiltonian, circuit, num_qubits, params_shape, beta):
 
    """ 
 
        Runs Variational Quantum Deflation algorithm on a given Hamiltonian
 
        to give its excited eigenvalues.
 

 
        Args: 
 
        hamiltonian (matrix(ndarray)): Hamiltonian for 
 
        circuit (Quantum Circuit): ansatz template
 
        num_qubits (int): number of qubits to make an ansatz
 
        params_shape (tuple): shape tuple for parameters in ansatz
 
        beta (int): beta coefficent for VQD
 

 
        Returns:
 
        energy (vector(ndarray)): Calculated eigen energies for the Hamiltonian. 
 
    """
 

 
    a, b, c = decompose_ham_to_pauli(hamiltonian)
 
    initial_params  = np.random.uniform(-np.pi, np.pi, params_shape)
 
    func  = ft.partial(vqe, meas_basis=a, coeffs=b, circuit=circuit, 
 
                  num_qubits=num_qubits, shots=2048)
 
    res1 = sp.optimize.minimize(func, initial_params, method='Powell')
 
    energies = [float(res1.fun)]
 
    thetas = [np.reshape(res1.x, params_shape)]
 
    if len(beta) != np.shape(hamiltonian)[0]-1:
 
        raise ValueError('Length mismatch! Provide sufficient amount of beta_{i}')
 

 
    for eigst in range(np.shape(hamiltonian)[0]-1):
 
        cost = ft.partial(vqd_cost, thetai=thetas, func=func, circuit=circuit, 
 
                  num_qubits=num_qubits, beta=beta, shots=2048)
 
        res2 = sp.optimize.minimize(cost, initial_params, method='Powell')
 
        if res2.success:
 
            thetas.append(np.reshape(res2.x, params_shape))
 
            energies.append(float(res2.fun))
 
        else:
 
            raise ValueError('Minimization was not successful.')
 

 
    return energies
 

 
 
 
 
 
 

 

Ansatz 1

xxxxxxxxxx
 
 
 
 
 
 
 
excited_energies = vqd(H1, ansatz1, 2, (2,2), [3.5, 1e-1, 1e-3])
 
 
excited_energies
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
[-1.0, 0.9921875, 0.994189453125, 1.0003808593749999]
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
assert(np.allclose(np.asarray(excited_energies),  np.sort(np.linalg.eig(H1)[0]), atol=1e-2))
 
 
 
 
 
 

 

Ansatz 2

xxxxxxxxxx
 
 
 
 
 
 
 
excited_energies = vqd(H1, ansatz2, 2, 1, [2.7, 1e-2, 1e-3])
 
 
excited_energies
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
[-1.0, 0.99365234375, 1.0012158203125, 1.0074306640625]
 
 
 
 
 
 

xxxxxxxxxx
 
 
 
 
 
 
assert(np.allclose(np.asarray(excited_energies),  np.sort(np.linalg.eig(H1)[0]), atol=1e-2))
 
 
 
 
 
 

 

Interpreting the Results

1. Both of our ansatz produced accurate results despite of ansatz 2 being more expressible than ansatz 1. One reason for this can be that the the ansatz 2 expressibility was covered the interested solution subspace. Another reason could be due to its more entangling capability.
2. Convergence of VQE depends upon both the type of optimizer (gradient-based or point-search) and the choice of ansatz. In general more the number of parameters, and the depth of ansatz, more difficult will be the convergence.
3. Number of shots did matter for non-optimal values of theta in ansatz 2. However, at optimal values of theta, less shots also gave good results. For noisy simulations, higher number of shots were required to get agreeable results.
4. COBYLA and Powell performed far better than others in both noisy and ideal conditions. Next one in line was Nelder-Mead. Gradient-based optimizers gave much poor results. In future, it would nice to see results from specialized optimizers such as rotosolve and rotoselect, and also optimizers that calculate gradients via parameter-shift rule.
5. Ansatz 2 performed better than ansatz 1 under Noisy simulation. However, mitigation improved the result for both of them.
6. Out of all the all the noise models tested, we saw VQE still gave the (deviated) minima at the optimal value of theta, and performed worst for errors due to damping and depolarization.
7. Even with depolarization, presence of errors in measurement gates gave much poor result than two-qubit gates which in turn gave worse results than single-qubit gates.
8. VQD was able to perform better as we adjusted the values of ${\beta }_i$$\beta_{i}$ based on the difference between the eigen-energies. In general, values of these ${\beta }_i$$\beta_{i}$ should be larger than the largest difference between the consecutive eigen-energies.