<jupyter_start><jupyter_text>MARKDOWN EQUATIONS LACK FULL SUPPORT ON GITHUB'S PREVIEWER. PLEASE CHOOSE ONE OF THE FOLLOWING OPTIONS.1. Edit and run on [Google Colab](https://colab.research.google.com/github/alexander-soare/A-Framework-Agnostic-VQE-Tutorial/blob/master/02_Tutorial.ipynb)    [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/alexander-soare/A-Framework-Agnostic-VQE-Tutorial/blob/master/02_Tutorial.ipynb)2. [Preview on Jupyter nbviewer](https://nbviewer.jupyter.org/github/alexander-soare/A-Framework-Agnostic-VQE-Tutorial/blob/master/02_Tutorial.ipynb) Please check out [Section 1]('./01_Prereading.md') in case you haven't already.  2.1 Our toy HamiltonianOkay so we're ready to dive in with a toy VQE example. To keep things very simple we'll stick with a 1-qubit system for now. The Hamiltonian is a 2x2 matrix, and one of the other prerequisites for VQE which I haven't mentioned yet, is that it needs to be expressible as a sum of products of Pauli matrices (shown below for reference).$$I \equiv \sigma_0 = \biggl(\begin{array} \ 1 & 0 \\ 0 & 1 \end{array}\biggr),\quad \sigma_x = \biggl(\begin{array} \ 0 & 1 \\ 1 & 0 \end{array}\biggr),\quad \sigma_y = \biggl(\begin{array} \ 0 & -i \\ i & 0 \end{array}\biggr),\quad \sigma_z = \biggl(\begin{array} \ 1 & 0 \\ 0 & -1 \end{array}\biggr)\tag{1}$$In the multi-qubit case we need to express the Hamiltonian as a sum of tensor products of Pauli matrices. It turns out that any Hamiltonian can be written in such a format, and that furthermore, the number of terms in such an expansion can be polynomial in the size of the system for some practical problems of interest (from the [VQE paper](https://arxiv.org/pdf/1304.3061.pdf)). That's important, because the purpose of VQE would be defeated if the prerequisites were exponentially hard to achieve.So an example for a single qubit (where it only makes sense to go to first order) might be:$$H = 2\sigma_0 + \sigma_x + 3\sigma_z \tag{2}$$Notice that we can also include the identity Pauli operator $I \equiv \sigma_0$ in the expansion.Right so now we can prepare a toy Hamiltonian. Let's code..._**IN PRACTICE** we wouldn't be the ones defining the Hamiltonian. It's already implicit in some quantum system we want to learn about._<jupyter_code>import cmath  # for complex math operations
import numpy as np  # for matrix multiplication
import random

# some boiler plate work
# feel free to skip reading this class, or just skim it
class HamiltonianBuilder():
    """
    Helper class for generating Hamiltonians from sums of single-qubit
    pauli matrices.
    """
    def __init__(self, coefficients_={}):
        # define standard single qubit puali matrices
        identity = np.eye(2)
        sigma_x = np.array([[0, 1], [1, 0]])
        sigma_y = np.array([[0, -1j], [1j, 0]])
        sigma_z = np.array([[1, 0], [0, -1]])
        self.pauli_matrices = {
            'i': identity,
            'x': sigma_x,
            'y': sigma_y,
            'z': sigma_z
        }
        # initialise random coefficents
        self._zero_coefficients()
        self.reset_coefficients(coefficients_)
        # build the hamiltonian
        self._build()

    def reset_coefficients(self, coefficients_={}):
        """
        coefficients_ is a dictionary expected in the form
        {'i': float, 'x': float, 'y': float, 'z': float}
        not all keys need to be present
        if coefficients_ is not provided, a random set is initialised
        """
        self._zero_coefficients()
        if len(coefficients_) <= 0:
            self._init_random_coefficients()
        else:
            assert all(key in self.coefficients.keys() for key in coefficients_.keys()), \
                "Coefficients dict incorrectly formed"
            self.coefficients.update(coefficients_)
        self._build()

    def _zero_coefficients(self):
        self.coefficients = {'i': 0, 'x': 0, 'y': 0, 'z': 0}

    def _build(self):
        """
        build the hamiltonian
        """
        # initialise zeros matrix for hamiltonian
        self.hamiltonian = np.zeros(shape=(2, 2)).astype('complex128')
        # run through coefficients dict and build hamiltonian additively
        for key in self.coefficients:
            self.hamiltonian += self.coefficients[key] * \
                self.pauli_matrices[key]

    @property
    def matrix(self):
        return self.hamiltonian

    def _init_random_coefficients(self):
        for k, v in self.coefficients.items():
            self.coefficients[k] = random.randint(-3, 3)
# now with the boiler plate work out of the way let's build a Hamiltonian
h_builder = HamiltonianBuilder()

# we can initialise a random hamiltonian (already done from object initialisation)
print(h_builder.coefficients) # print the coefficients of the pauli terms
print('\n')
print(h_builder.matrix) # and the resulting Hamiltonian

# or we can specify coefficients
# lets go with the ones in equation (2) for now
coefficients = {'i': 2, 'x': 1, 'z': 3}
h_builder.reset_coefficients(coefficients)
print(h_builder.coefficients)
print('\n')
print(h_builder.matrix)<jupyter_output>{'i': 2, 'x': 1, 'y': 0, 'z': 3}


[[ 5.+0.j  1.+0.j]
 [ 1.+0.j -1.+0.j]]<jupyter_text>2.2 Energy and energy eigenvaluesNow recall from [Section 1](01_Prereading.md) that we're measuring the energy of a quantum state with respect to our Hamiltonian. To jog your memory on that, just consider some eigenstate $|E⟩$, then by definition:$$H|E⟩ = E|E⟩ \tag{3}$$where $E$ is the associated eigenvalue (or equivalently "energy"). So by doing an inner product between the result and the original state (this is called getting the expectation value of observable $H$) we get:$$⟨H⟩_E \equiv ⟨E|H|E⟩ = E⟨E|E⟩ = E \tag{4}$$This extends to arbitrary states (not just eigenstates) because they can be expresses as weighted sums of the eigenstates:$$|\Psi⟩ = \sum_{i} c_i|E_i⟩ \tag{5}$$where $i$ enumerates eignenstates and $c_i$ are complex coefficients. It takes a couple of lines of applying linearity and the principle that distinct eignenstates are orthogonal (meaning their inner product is zero) to get:$$⟨H⟩_{\Psi} = \sum_{i} |c_i|^2E_i \tag{6}$$_Pause... While we're getting all carried away with doing the math, let's just remind ourselves of the physical significance of (6). Remember that the coefficients in a wavefunction expansion really represent the probability that when you measure a qubit, you get the corresponding eigenstate as an outcome. This equation is calculating a weighted sum or "expectation value" of the energy distribution, which is the same as the expectation value in probability theory. Unpause..._Anyway, going back to equation (2) and by the linearity of (4) we can express the expectation value of $H$ in terms of weighted expectation values of its expansion. Here I'll just show that for our toy example:$$⟨H⟩_{\Psi} = 2⟨\sigma_0⟩_{\Psi} + ⟨\sigma_x⟩_{\Psi} + 3⟨\sigma_z⟩_{\Psi} \tag{7}$$where the subscript $|\Psi⟩$ means this equation applies to arbitrary states._Pause... Why am I leading you this way? And why have I been making such a big deal out of Pauli matrices? ... Because the Pauli matrices represent operations we know how to apply with real hardware in real quantum systems. You'll soon see where it comes into play. Unpause..._Just because we're going to the trouble of doing this whole tutorial without a framework, let's verify that equation (7) is true for our single-qubit toy system. To do that we'll first need to create a Python class to represent a qubit. We can write an arbitrary state as:$$|\Psi⟩ = c_1|0⟩ + c_2|1⟩ \tag{8}$$where $|0⟩$ and $|1⟩$ are our computational basis states and the $c$ are complex coeffiecients. Let's code..._**IN PRACTICE** we can't just write down an arbitrary state and do stuff with it. We have to prepare it somehow. More on that soon._<jupyter_code># lets create our Qubit class here
import numpy as np
import random
import cmath
import math

class Qubit():
    def __init__(self, coefficents=(0, 0)):
        self.comp_0 = np.array([1, 0]) # |0⟩
        self.comp_1 = np.array([0, 1]) # |0⟩
        self.reinitialise_state(coefficents)

    # set the coefficients
    def reinitialise_state(self, coefficents=(0, 0)):
        c0, c1 = coefficents
        if c0 == 0 and c1 == 0:
            c0 = complex(random.random(), random.random())
            c1 = complex(random.random(), random.random())
        self.c0, self.c1 = self._normalise_coefficients(c0, c1)

    # since |c_0|^2 and |c_1|^2 are probabilities, they need to sum to 1
    def _normalise_coefficients(self, c0, c1):
        denom = math.sqrt(c0*np.conj(c0) + c1*np.conj(c1))
        return c0/denom, c1/denom

    @property
    def state_vector(self):
        return np.array([self.c0, self.c1])

    def pprint_state(self):
        """
        pretty prints the state
        also prints it in polar form (https://en.wikipedia.org/wiki/Bloch_sphere)
        """
        r_0, theta_0 = cmath.polar(self.c0)
        r_1, theta_1 = cmath.polar(self.c1)
        dphase = theta_1 - theta_0
        print(f'({self.c0:.2f})|0⟩ + ({self.c1:.2f})|1⟩ ≡ {r_0:.2f}|0⟩ + {r_1:.2f}⋅e^({dphase/math.pi:.2f}iπ)⋅|1⟩')

# we can initialise a random state (already done from object initialisation)
qubit = Qubit()
qubit.pprint_state()
print('\n')

# or we can initialise a state of our choosing
qubit.reinitialise_state((0.5, 0.5))
qubit.pprint_state()<jupyter_output>(0.40+0.47j)|0⟩ + (0.57+0.55j)|1⟩ ≡ 0.61|0⟩ + 0.79⋅e^(-0.03iπ)⋅|1⟩


(0.71)|0⟩ + (0.71)|1⟩ ≡ 0.71|0⟩ + 0.71⋅e^(0.00iπ)⋅|1⟩<jupyter_text>Now that we have a qubit we can just do matrix multiplication with our hamiltonian from before to compute the expectation value.<jupyter_code># If we we get the energy of our state w.r.t. our Hamiltonian it's
energy_hamiltonian = np.matmul(
    np.conj(qubit.state_vector), # complex conjugate
    np.matmul(h_builder.matrix, qubit.state_vector)
)

# notice that the value is purely real as we expect for a Hamiltonian
#  (remember: Hamiltonians are hermitian, and so are Pauli matrices,
#   eigenvalues of hermitian matrices are purely real)
print('Energy w.r.t. toy hamiltonian:', f'{energy_hamiltonian:.2f}')<jupyter_output>Energy w.r.t. toy hamiltonian: 3.00+0.00j<jupyter_text>We can also use (7) to do the same as above, but for each term in the pauli matrix expansion, then sum the result.<jupyter_code># first define our pauli terms
identity = np.eye(2)
sigma_x = np.array([[0, 1], [1, 0]])
sigma_y = np.array([[0, -1j], [1j, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
terms = [2 * identity, sigma_x, 3 * sigma_z]

# now we might as well write a function to get the energy
def analytically_measure_observable(observable, state_vector):
    energy = np.matmul(
        np.conj(state_vector),
        np.matmul(observable, state_vector)
    )
    return energy

# get the energies
energies = [analytically_measure_observable(term, qubit.state_vector) for term in terms]

print('Energy w.r.t. pauli matrix expansion:',
    f'{np.sum(energies):.2f}')<jupyter_output>Energy w.r.t. pauli matrix expansion: 3.00<jupyter_text>Let's also find the smallest eigenvalue of the Hamiltonian by solving the characteristic equation:$$H - EI = 0 \tag{9}$$<jupyter_code># we'll let numpy do the work for us
eigenvalues = np.linalg.eigvals(h_builder.matrix)
min_eigenvalue = np.min(eigenvalues)
print(f'Energy of ground state for single-qubit example is {np.real(min_eigenvalue):.3f}')
print('\n')

# oh and while we're at it, lets prove a point
# lets see how long that would take to run for a 10 qubit system
# play with higher numbers at your own risk!
import time
dim = 2**10 # 10 qubits
# not a hermitian matrix, but good enough for our demo
big_h = np.array([complex(*doublet) for doublet in np.random.random((dim**2, 2))]).reshape(dim, dim)
start = time.time()
eigenvalues = np.linalg.eigvals(big_h)
min_eigenvalue = np.min(eigenvalues)
print(f'Energy of ground sate for multi-qubit example is {np.real(min_eigenvalue):.3f} ' + \
        f'and it took {time.time() - start:.3f} s to solve')<jupyter_output>Energy of ground state for single-qubit example is -1.162


Energy of ground sate for multi-qubit example is -13.082 and it took 6.652 s to solve<jupyter_text>2.3 State preparation and ansatzesNow, with real quantum systems in the lab, usually we can only reliably initialise one of few states, and then we need to apply quantum operations to reach some other desired state. In our toy model let's say we can only reliably initialise to $|0⟩$. Then we need some mechanism to get us to some arbitrary state $|\Psi⟩$. Now I'll introduce the term "ansatz", which is just a fancy term for the control knobs we can use to reach $|\Psi⟩$.For example, consider the Bloch sphere with the qubit initialised to $|0⟩$.If your ansatz is $Ry$ (all possible rotations around the $y$-axis) then we can only reach anywhere on the longitude crossing the $x$ and $z$ axes. But if you combine this with $Rx$ now we can reach all possible states._**IN PRACTICE** the number of parameters needed to fully describe a state scales exponentially as the number of qubits. On the other hand, we do not want to have to apply an exponentially increasing number of gates to prepare our states, as that would defeat the purpose of VQE. Finding a good ansatz that's likely to cover the ground state we're looking for is a non-trivial task. And now you know why the concept of an ansatz deserves its own name!_We need to modify our Qubit class such that we're only able to initialise to $|0⟩$. We'll also add a method for applying operator to the qubit, as we'll need to do so to prepare a state.<jupyter_code># subclassing our Qubit class as a shortcut to avoid rewriting code
class Qubit2(Qubit):
    def __init__(self):
        self.comp_0 = np.array([1, 0])
        self.comp_1 = np.array([0, 1])
        self.reinitialise_state()

    def reinitialise_state(self):
        # now we can only initialise to |0⟩
        self.c0, self.c1 = self._normalise_coefficients(1, 0)

    def _apply_operation(self, operator):
        """
        operator is the matrix representation of the operator to be applied
        """
        self.c0, self.c1 = np.matmul(operator, self.state_vector)

    def apply_operations(self, operators):
        """
        operators is a list of matrix representations of operators to be applied
        they will be applied from left to right (opposite of mathematical notation)
        """
        for operator in operators:
            self._apply_operation(operator)

qubit = Qubit2()
qubit.pprint_state()<jupyter_output>(1.00)|0⟩ + (0.00)|1⟩ ≡ 1.00|0⟩ + 0.00⋅e^(0.00iπ)⋅|1⟩<jupyter_text>Now we need to define parametrised $R_x$ and $R_y$ rotations. In matrix form they look like this:$$R_x(\theta) = \biggl(\begin{array} \ \cos{(\theta_x/2)} & -i\sin{(\theta_x/2)} \\ -i\sin{(\theta_x/2)} & \cos{(\theta_x/2)} \end{array}\biggr),\quad R_y(\theta) = \biggl(\begin{array} \ \cos{(\theta_y/2)} & -\sin{(\theta_y/2)} \\ \sin{(\theta_y/2)} & \cos{(\theta_y/2)} \end{array}\biggr)\tag{10}$$where $\theta_x$ and $\theta_y$ are the angles of rotation about the $x$- and $y$- axes respectively. The following [figure](https://www.researchgate.net/figure/Example-universal-set-of-quantum-gates-consisting-of-three-single-qubit-rotation-gates_fig3_327671865) will help you get your bearings (we don't need the z-rotation for our toy example):![](./figures/rotations.png)Let's code...<jupyter_code>def get_rx(theta):
    rx = np.array([
        [math.cos(theta/2), complex(0, -math.sin(theta/2))],
        [complex(0, -math.sin(theta/2)), math.cos(theta/2)]
    ])
    return rx

def get_ry(theta):
    rx = np.array([
        [math.cos(theta/2), -math.sin(theta/2)],
        [math.sin(theta/2), math.cos(theta/2)]
    ])
    return rx

# let's test it out
# a rotation of pi/4 by y followed by a rotation of -pi/2
#  by x should get us on the positive diagonal of the x-y plane
#  (use the bloch sphere to trace that out)
# so that would be an equal superposition of |0⟩ and |1⟩
#  with a relative phase of pi/4
ry = get_ry(math.pi/4)
rx = get_rx(-math.pi/2)
qubit.reinitialise_state() # just in case
qubit.apply_operations([ry, rx])
qubit.pprint_state()<jupyter_output>(0.65+0.27j)|0⟩ + (0.27+0.65j)|1⟩ ≡ 0.71|0⟩ + 0.71⋅e^(0.25iπ)⋅|1⟩<jupyter_text>Check out the polar representation in the output of the cell above. As we expected, it has equal amplitudes in $|0⟩$ and $|1⟩$, and the relative phase is $\pi/4$.  2.4 Energy and energy eigenvalues (again)_Pause... "But we already did this in section 2.2!" Yes but that's how we would solve for energy eigenvalues analytically. But we need to think about how we would actually do the measurements of the different Pauli observables in a real quantum device. Unpause..._Let's bring up our previous equation for getting the energy associated with an arbitrary state ${\Psi}$:$$⟨H⟩_{\Psi} = 2⟨\sigma_0⟩_{\Psi} + ⟨\sigma_x⟩_{\Psi} + 3⟨\sigma_z⟩_{\Psi} \tag{7}$$And now let's translate that into lab talk. So we've prepared a state ${|\Psi⟩}$ and we want to measure the expectation values on the RHS of the equation in order to effectively simulate measuring the expectation value of $H$. Starting with $2⟨\sigma_0⟩_{\Psi}$This one's easy. The math already tells us that $⟨\sigma_0⟩_{\Psi} = 1$ irrespective of $|\Psi⟩$, so the whole term will just add a constant $2$ to our result. Moving onto $3⟨\sigma_z⟩_{\Psi}$$⟨\sigma_z⟩_{\Psi}$ just represents measuring $|\Psi⟩$ in our computational basis (remember that $|0⟩$ and $|1⟩$ lie on the $z$-axis of the Bloch sphere).We also know that we can't measure twice, because upon one measurement, we collapse the wavefunction to one of ${|0⟩}$ or ${|1⟩}$. In practice, we can only measure one of the two, and we'd have to prepare the same state and do the measurement many times in order get an estimate for the true energy eigenstate (just like flipping a loaded coin many times to determine an estimate for it's true probability of giving heads)._Pause... Maybe you want to know what a real experiment looks like. When I used to work in a lab, we would measure the state of trapped ion qubits by pointing a laser with a precise wavelength at them. Ions in the ${|1⟩}$ state would scatter the photons and we would capture that with a detector, whereas ions in the ${|0⟩}$ state would not interact with the photons. Unpause..._So say we do 1000 experiments, and we measure ${|0⟩}$ 600 times, and ${|1⟩}$ 400 times. Then the resulting estimate for $⟨\sigma_z⟩_{\Psi}$ is: $600/1000 \times (+1) + 400/1000 \times (-1) = 0.2$. More generally, we can write that if we take $N$ measurements, we get:$$⟨H⟩_\Psi = \frac{n_0}{N}\times(+1) + \frac{n_1}{N}\times(-1) = \frac{n_0 - n_1}{N} \tag{11}$$where $n_1$ and $n_2$ are the number of times we measure ${|0⟩}$ and ${|1⟩}$ respectively. Also, from statistics, our precision will vary with $\sqrt{N}$ - in other words, we'll need to set up our experiment to do polynomially more measurements if we need more precision. And finally $⟨\sigma_x⟩_{\Psi}$In many quantum devices, we can only measure directly in some axes (usually the z-axis, or "computational basis"). So given that we can't measure $⟨\sigma_x⟩_{\Psi}$ directly, we first rotate by $-\pi/2$ around the $y$-axis then do a normal measurement in $z$ as before. What we actually did was the equivalent of rotating our co-ordinate system about a fixed state-vector. As an answer in the comments section of this [article](https://www.mustythoughts.com/variational-quantum-eigensolver-explained) puts it:> The rationale is to change the basis of the measurement. If you want to measure the length of a table, but you have tape which allows only for measuring the height, you can't do that – unless you rotate the whole table._Pause... It's not part of our toy example, but if we needed to do the same for $⟨\sigma_y⟩$ we'd have to rotate by $+\pi/2$ around the $x$-axis. I invite you to pull up a picture of a Bloch sphere and think about it. Unpause..._So going back to the last bit of code we wrote, where we prepared a state:$$(0.65+0.27i)|0⟩ + (0.27+0.65i)|1⟩ ≡ 0.71|0⟩ + 0.71\cdot e^{0.25i\pi}|1⟩ \tag{12}$$we can think about what our terms would give.- $2⟨\sigma_0⟩_{\Psi}$ will give 2- $3⟨\sigma_z⟩_{\Psi}$ will give something near 0 because we have an equal superposition so we'd measure $+1$ about as often as $-1$ (and so they would cancel each other out in (11)).- For $⟨\sigma_x⟩_{\Psi}$ we first need to mentally rotate by $-\pi/2$ around the $y$-axis, which will put us on the positive diagonal of the $y-z$ plane. So we should get something between 0 and 1 as we are in the northern Hemisphere and more than half of our measurements will collapse to ${|0⟩}$.So summing that all up should get something in between 2 and 3. Now let's code...<jupyter_code># let's modify our Qubit class so that we an do measurements in
# the computational basis
class Qubit3(Qubit2):
    def __init__(self):
        super(Qubit3, self).__init__()

    def measure_state(self):
        """
        simulate a random measurement outcome
        returns +1 for |0⟩ or -1 for |1⟩
        """
        # probability of |0⟩
        prob_0 = self.c0 * np.conj(self.c0)
        if random.random() < prob_0:
            # collapse to |0⟩ and return +1
            self.c0 = complex(1,0)
            self.c1 = complex(0,0)
            return 1
        else:
            # collapse to |1⟩ and return -1
            self.c0 = complex(0,0)
            self.c1 = complex(1,0)
            return -1

# let's define a function to simulate multiple rounds of state prep
#  and measurement
def estimate_energy(qubit, state_prep_operators, basis='z', n_exp=100):
    """
    qubit is a Qubit object
    state_prep_operators is a list of operators to be applied to get from |0⟩
     to |Ψ⟩
    basis is one of 'x', 'y', or 'z'
    n_exp is number of experiments
    """
    outcomes = []
    # prepare change of basis operator (in order to measure in x or y)
    if basis == 'x':
        change_basis_operator = get_ry(-math.pi/2)
    elif basis == 'y':
        change_basis_operator = get_rx(math.pi/2)
    elif basis != 'z':
        raise UserWarning("Please provide one of 'x', 'y', or 'z' for basis")
    for i in range(n_exp):
        qubit.reinitialise_state() # initialize to |0⟩
        qubit.apply_operations(state_prep_operators) # prepare state
        if basis == 'x' or basis == 'y':
            # change basis before measurement
            qubit.apply_operations([change_basis_operator]) 
        outcomes.append(qubit.measure_state())
    return np.mean(outcomes)

# now first let's just prepare our state once and get the analytical result
#  with the function we defined in section 2.2
# I'll repeat some local variable initialisation here for convenience

# set up our toy hamiltonian
h_builder = HamiltonianBuilder()
coefficients = {'i': 2, 'x': 1, 'z': 3}
h_builder.reset_coefficients(coefficients)

# prepare our example qubit state
qubit = Qubit3()
ansatz_gates = [get_ry(math.pi/4), get_rx(-math.pi/2)]
qubit.apply_operations(ansatz_gates)
analytical_result = analytically_measure_observable(h_builder.matrix, qubit.state_vector)
print(f"Analytical result is {np.real(analytical_result):.2f}")
print('\n')

# Okay, and now we can try it by simulating a real(ish) quantum experiment. How exciting!
print("Running experiment...")
# 2 from our first term with sigma_0
exp_result = 2
# then the 1*sigma_x term
result = estimate_energy(qubit, ansatz_gates, basis='x', n_exp=10000)
print(f"  sigma_x contribution is {result:.2f}")
exp_result += result
# then the 3*sigma_z term
result = estimate_energy(qubit, ansatz_gates, basis='z', n_exp=10000)
print(f"  sigma_z contribution is {3*result:.2f}")
exp_result += 3*result

print(f"Experimental result is {exp_result:.2f}")

# notice that if you run this multiple times you will get different experimental results each time<jupyter_output>Analytical result is 2.71


Running experiment...
  sigma_x contribution is 0.70
  sigma_z contribution is -0.00
Experimental result is 2.70<jupyter_text>2.5 Classical optimisationLet's just take stock of what we've got so far.- We can initialise a qubit.- We can apply a parametrised ansatz to get any arbitrary state.- We can measure the energy of this state with respect to a Hamiltonian by breaking it down into pauli observables and running measurement experiments on each one.Now we need a classical optimisation loop which will iteratively suggest qubit-states to trial. It's not my intention to go into the topic of optimisation here. Nevertheless, it may be worth your time to look into it as the choice of algorithm is intricately tied into our understanding of the dynamics of NISQ devices (one good spot would be to read this [recent paper](https://arxiv.org/pdf/2004.03004.pdf)).So this tutorial is not _entirely_ framework agnostic in that I'll be using the the optimisation functionality of the SciPy library. And we'll be specifically using Powell's method for optimisation. See here for [SciPy documentation](https://docs.scipy.org/doc/scipy/reference/optimize.minimize-powell.htmloptimize-minimize-powell), and here for a [neat explanation I found on YouTube](https://youtu.be/ep4C_R7dUMs?t=1168)._**IN PRACTICE**: [Limited-memory BFGS](https://en.wikipedia.org/wiki/Limited-memory_BFGS) (L-BFGS) is the optimisation algorithm typically used in real VQE. Powell's method isn't really used, but I've chosen it because it works well for our toy example (and to be honest, I couldn't get L-BFGS to work here. It's on my backlog, but a PR would be welcome if you can figure it out!)_We'll now take a moment to get the hang of using the SciPy optimizer. Let's code..._Pause... by now you may have noticed I'm totally inconsistent with using the "ise" vs "ize" suffix. That's what happens when you come from any english-speaking country other than the US and you read a lot of papers. I don't think I'll ever settle on one or the other. Unpause..._<jupyter_code>import scipy.optimize as optimize

# let's test our optimisation method on a 2d paraboloid
#  centered at (0.5, 2)
def paraboloid(x):
    """
    v - 1d, length 2 array representing a point in the domain of the paraboloid
    """
    return (x[0]-0.5)**2 + (x[1]-2.)**2

result = optimize.minimize(paraboloid,
                        np.random.rand(2) * 10, # randomly select initial point 
                        method='Powell',
                        tol=1e-3) # desired error tolerance tells our search when to stop

print(f"Result {result.x} found in {result.nit} steps")<jupyter_output>Result [0.5 2. ] found in 3 steps<jupyter_text>2.6 Let's run VQE!We're approaching the moment of truth! Here we'll take the code from section 2.4 one step further. We'll write a class which takes Hamiltonian expansion coefficients and parameters for our ansatz, then runs the necessary experiments to return an expectation value for that Hamiltonian. Let's code...<jupyter_code># in writing this class we'll copy over a bunch of function definitions from earlier
#  the `evaluate` function is the only new code of interest here
class QuantumExperiment():
    def __init__(self, h_builder, qubit, default_n_exp=1000):
        """
        h builder is an instance of the HamiltonianBuilder class
        qubit is an instance of the Qubit3 class
        default_n_exp is the default number of experiments to run 
         per evaluation of each pauli observable
        """
        self.h_builder = h_builder
        self.qubit = qubit
        self.default_n_exp = default_n_exp
        # we'll also keep a record of results
        self.results = []
    
    def clear_results():
        self.results = []

    def evaluate(self, ansatz_params):
        """
        ansatz_params is a 1d, length 2 list of the ansatz parameters
         it should be specified as [theta_y, theta_x]
        """
        # get a list of gates that we need to perform state preparation
        theta_y, theta_x = ansatz_params
        ansatz_gates = [self._get_ry(theta_y), self._get_rx(theta_x)]
        
        # outer experiment loop - looping through observable evaluations
        exp_result = 0
        # add in the constant energy associated with sigma_0 if it's there
        exp_result += self.h_builder.coefficients.get('i', 0)
        # now go through all the non-trivial pauli observables
        for key in '<KEY>':
            if key in self.h_builder.coefficients and self.h_builder.coefficients[key] > 0:
                weight = self.h_builder.coefficients[key]
                exp_result += weight * self.estimate_energy(
                    self.qubit, ansatz_gates, basis=key, n_exp=self.default_n_exp)
        self.results.append(exp_result)
        return exp_result

    # copying this over from earlier
    def estimate_energy(self, qubit, state_prep_operators, basis='z', n_exp=100):
        """
        qubit is a Qubit object
        state_prep_operators is a list of operators to be applied to get from |0⟩
            to |Ψ⟩
        basis is one of 'x', 'y', or 'z'
        n_exp is number of experiments
        """
        outcomes = []
        # prepare change of basis operator
        if basis == 'x':
            change_basis_operator = get_ry(-math.pi/2)
        elif basis == 'y':
            change_basis_operator = get_rx(math.pi/2)
        elif basis != 'z':
            raise UserWarning("Please provide one of 'x', 'y', or 'z' for basis")
        for i in range(n_exp):
            qubit.reinitialise_state()
            qubit.apply_operations(state_prep_operators) # prepare state
            if basis == 'x' or basis == 'y':
                # change basis before measurement
                qubit.apply_operations([change_basis_operator]) 
            outcomes.append(qubit.measure_state())
        return np.mean(outcomes)

    # copying this over from earlier
    def _get_rx(self, theta):
        """
        define r_x matrix representation
        """
        rx = np.array([
            [math.cos(theta/2), complex(0, -math.sin(theta/2))],
            [complex(0, -math.sin(theta/2)), math.cos(theta/2)]
        ])
        return rx

    # copying this over from earlier
    def _get_ry(self, theta):
        """
        define r_y matrix representation
        """
        rx = np.array([
            [math.cos(theta/2), -math.sin(theta/2)],
            [math.sin(theta/2), math.cos(theta/2)]
        ])
        return rx


# set up our toy hamiltonian
h_builder = HamiltonianBuilder()
coefficients = {'i': 2, 'x': 1, 'z': 3}
h_builder.reset_coefficients(coefficients)

# prepare a qubit
qubit = Qubit3()

# get an instance of our QuantumExperiment class
q_exp = QuantumExperiment(h_builder, qubit, 5000)

# use a code snippet from earlier to get the analytical result
#  for the minimum eigenvalue
eigenvalues = np.linalg.eigvals(h_builder.matrix)
min_eigenvalue = np.min(eigenvalues)
print(f'Analytical result is {np.real(min_eigenvalue):.3f}')
print('\n')

# now the moment of truth... let's run VQE!
# # this could take 30~120 s with n_exp = 5000
result = optimize.minimize(q_exp.evaluate,
                        (np.random.rand(2) - 0.5) * 4*math.pi, # randomly select initial point 
                        method='Powell',
                        tol=1e-3) # desired error tolerance tells our search when to stop

print(f"Experimental result is {result.fun:.3f} found in {result.nit} steps")<jupyter_output>Analytical result is -1.162


Experimental result is -1.172 found in 3 steps
{'fraction_non_alphanumeric': 0.10717339667458432, 'fraction_numerical': 0.021060965954077592, 'mean_word_length': 3.427990464170523, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 12, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '455518', 'license': 'Apache-2.0', 'n_tokens_mistral': 10244, 'n_tokens_neox': 9770, 'n_words': 3774, 'path': '02_Tutorial.ipynb', 'repo_name': 'alexander-soare/A-Framework-Agnostic-VQE-Tutorial', 'size': '31581'}