<jupyter_start><jupyter_text>State preparation with circuit optimizationWe want to create a circuit that produces the Bell state $\vert\Phi^+\rangle = \dfrac{\vert00\rangle + \vert11\rangle}{\sqrt 2}$. We already know that this state can be produced by a circuit containing a Hadamard gate on the first qubit followed by a CNOT gate [[1](https://en.wikipedia.org/wiki/Bell_stateCreating_Bell_states)], but we would like to replace the Hadamard gate with rotation gates.In the following we will use [PennyLane](https://pennylane.readthedocs.io/), a cross-platform Python library for quantum machine learning, because it abstracts away several implementation details (especially when it comes to using different quantum simulators and frameworks such as pyQuil/Forest). PennyLane ships with its own version of NumPy, enriched to include [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation) - which lets the gradients on the circuits be computed in a more efficient way using clever transformations. In the following, I will also highlight my contributions to the library: the `SquaredErrorLoss` class and the ongoing work on a Hamiltonian decomposition method.Please check the README file for instructions on how to prepare the environment and execute this notebook.  Importing the library<jupyter_code>import pennylane as qml

import requests


from pennylane import numpy as np
from pennylane.qnn.cost import SquaredErrorLoss<jupyter_output><empty_output><jupyter_text>Circuit definitionFirst of all we create the "ansatz" circuit, a base structure consisting of a list of gates applied to specific wires. By doing this we can reuse the same configuration in different circuits, for instance if all we want to change is the measurement.Since we want to tackle the general case, we put both kinds of the allowed rotation gates on both qubits; then, we also add a CNOT gate. We will therefore have 4 different parameters, each corresponding to a different angle as a parameter for each rotation gate.<jupyter_code>requests.post("http://localhost:5000/qvm")
def circuit(angles, **kwargs):
    qml.RX(angles[0], wires=0)
    qml.RY(angles[1], wires=0)

    qml.RX(angles[2], wires=1)
    qml.RY(angles[3], wires=1)

    qml.CNOT(wires=[0, 1])<jupyter_output><empty_output><jupyter_text>In order to create a circuit, we need to include a _device_ and a _measurement_. For starting out we can use PennyLane's "default qubit" device, which can simulate a real quantum device with any number of qubits; we also set `analytic=False` and `shots=1000` to simulate a more realistic device._Other devices can be used with the `pennylane-forest` plugin, which can be installed as explained in the README file. Once the plugin is installed and the QVM is running, any `forest.qvm` device can be uncommented and used in place of `default.qubit`. The `forest.qpu` device instead requires access to a real QPU._<jupyter_code># dev = qml.device('default.qubit', wires=2, analytic=False, shots=1000)

# For a more realistic QVM simulation (needs `pennylane-forest`, `qvm` and `quilc`, see README)
#dev = qml.device('forest.qvm', device='2q-qvm', shots=10, noisy=True)

# should be able to run with as is- just resource issues
# For Aspen-8 simulation (needs `pennylane-forest`, `qvm` and `quilc`, see README)
# dev = qml.device('forest.qvm', device='Aspen-8-qvm', shots=1, noisy=False)

# For a real QPU - start with 1 shot up by factor of 10 each time - goal to 1000
dev = qml.device('forest.qpu', device='Aspen-8', shots=1000)

# For a list of all the devices supported by the Forest SDK:
#
# from pyquil import list_quantum_computers
# sorted(list_quantum_computers())

# To check the capabilities of a device:
#
# dev.capabilities()<jupyter_output><empty_output><jupyter_text>We then create a list with one `PauliZ` _observable_ that we will use for the measurement.<jupyter_code>observables = [qml.PauliZ(0)]<jupyter_output><empty_output><jupyter_text>Now we can initialize some parameters (the gates' angles) and "execute" the circuit, using the expectation value of the observables as a measurement. The `measure` argument default for `qml.map` is already `'expval'`, so here it is added just for clarity.<jupyter_code>params = [0.2, 0.8, 0.4, 0.1]
qnode = qml.map(circuit, observables, dev, measure='expval')
print(qnode(params))<jupyter_output>[0.582]<jupyter_text>The problem of measuring only one qubit is that, in our case where we have an entangled state, such measurement would affect the other qubit as well by "assigning" it the same value. For this reason, we need to measure the 2 qubits together instead.We can use different combinations of observables:- combined observables: `[qml.PauliZ(0), qml.PauliZ(1)]` (returns two values);- tensor product observables: `[qml.PauliZ(0) @ qml.PauliZ(1)]` (returns one value);- generic Hermitian observables: `[qml.Hermitian(observable, wires=[0, 1])]` (returns one value).The latter is the most flexible, and since it is supported not only by [PennyLane](https://pennylane.readthedocs.io/en/stable/code/api/pennylane.Hermitian.html) but also by other frameworks (e.g. [Forest](http://docs.rigetti.com/en/stable/apidocs/pauli.htmlpyquil.paulis.PauliSum) supports sums of Pauli operators, and [Qiskit](https://qiskit.org/textbook/ch-gates/fun-matrices.htmlPauli-decomposition) explains how to do it), we can make use of it.We start with defining our target state $\vert\psi\rangle = \dfrac{1}{\sqrt 2} \left[\begin{matrix}1\\0\\0\\1\end{matrix}\right]$:<jupyter_code>psi = 1. / np.sqrt(2) * np.array([[1, 0, 0, 1]]).T<jupyter_output><empty_output><jupyter_text>Now we construct the observable as the outer product $O = \vert\psi\rangle \langle\psi\vert$ of the state with itself:<jupyter_code>obs = psi @ psi.conj().T
print(obs)<jupyter_output>[[0.5 0.  0.  0.5]
 [0.  0.  0.  0. ]
 [0.  0.  0.  0. ]
 [0.5 0.  0.  0.5]]<jupyter_text>Why do we "like" this matrix? Because, when used for measurements, it will "boost" the states that have close values for their first and last element. This matrix is actually the projector onto the \$\psi\$ state with eigenvalues 0 and 1.<jupyter_code>eigenvalues, eigenvectors = np.linalg.eig(obs)
print(eigenvalues)<jupyter_output>[1. 0. 0. 0.]<jupyter_text>The eigenstate corresponding to the eigenvalue 1 is our desired state:<jupyter_code>eig_1 = np.array([eigenvectors[:, 0]]).T
print(eig_1)
print(np.allclose(eig_1, psi))<jupyter_output>[[0.70710678]
 [0.        ]
 [0.        ]
 [0.70710678]]
True<jupyter_text>Therefore, we need to optimize its expectation value to be 1. In fact, we have that obviously $\langle\psi\vert O \vert\psi\rangle = 1$:<jupyter_code>print(psi.T @ obs @ psi)<jupyter_output>[[1.]]<jupyter_text>and a vector $\tilde\psi$ close to $\psi$ will have a measurement close to 1:<jupyter_code>psi_tilde = psi - 0.001

# For normalization
psi_tilde[1:3] = np.sqrt(0.5 - psi_tilde[0] ** 2)
assert np.allclose(np.linalg.norm(psi_tilde), 1.0)

print(psi_tilde.T @ obs @ psi_tilde)<jupyter_output>[[0.99717357]]<jupyter_text>Let's also verify whether the observable matrix is Hermitian by checking that it is equal to its own adjunct, i.e. $M = M^\dagger$:<jupyter_code>print(np.allclose(obs, obs.conj().T))<jupyter_output>True<jupyter_text>We can also check whether the matrix can be decomposed as a linear combination of Pauli operators tensor products:$O = \dfrac{1}{4}((I \otimes I) + (X \otimes X) - (Y \otimes Y) + (Z \otimes Z))$This step has to be performed manually, but my contribution for a decomposition method is [in progress](https://github.com/XanaduAI/pennylane/pull/671).<jupyter_code># Manually-defined Pauli operators. They can also be derived from PennyLane observables,
# e.g. PauliX = qml.PauliX(0).matrix

identity = np.eye(2)
pauliX = np.array([[0, 1], [1, 0]])
pauliY = np.array([[0, -1j], [1j, 0]])
pauliZ = np.array([[1, 0], [0, -1]])

decomp = 0.25 * np.sum([
    np.kron(identity, identity),
    np.kron(pauliX, pauliX),
    -np.kron(pauliY, pauliY),
    np.kron(pauliZ, pauliZ)
], axis=0)

# Is the decomposed matrix the same as the observable matrix we created before?
print(np.allclose(obs, decomp))<jupyter_output>True<jupyter_text>Now we are ready to run the circuit again using the new observable:<jupyter_code>observables = [qml.Hermitian(obs, wires=[0, 1])]

qnode = qml.map(circuit, observables, dev, measure='expval')
print(qnode(params))<jupyter_output>[0.579]<jupyter_text>OptimizationSince the circuit optimization has to be performed via gradient descent, we need first of all a good choice for the initial parameters. We can try a few different ones, including:- all parameters equal to 0;- all parameters randomly chosen;- parameters initialized to some chosen "sensible" defaults.<jupyter_code>params_init_method = 'chosen'

if params_init_method == 'zero':
    params = np.array([0.] * 4)
elif params_init_method == 'random':
    import random
    
    random.seed(0)
    params = np.random.normal(0., np.pi, 4)
elif params_init_method == 'chosen':
    params = np.array([np.pi / 4] * 4)
else:
    raise ValueError('{} initialization method does not exist'.format(weight_init_method))<jupyter_output><empty_output><jupyter_text>Then, we need to decide how many optimization steps we will run:<jupyter_code>steps = 1000<jupyter_output><empty_output><jupyter_text>We also create an additional array that will collect the learned parameters per round:<jupyter_code>params_history = np.zeros((4, steps))<jupyter_output><empty_output><jupyter_text>We define the cost function as the square of the difference between the value of the observable (which can be between 0 and 1) and 1; this means that the closer the observable gets to one, the smaller the cost becomes.My contribution is the `SquaredErrorLoss` class, which gives an easy way to calculate the loss given a target.<jupyter_code>loss = SquaredErrorLoss(circuit, observables, dev)

def cost(params):
    return loss(params, target=[1])<jupyter_output><empty_output><jupyter_text>Finally, we define the optimizer that we will use in the optimization process. We can start with the simplest `GradientDescentOptimizer`, and optionally try a more advanced optimizer such as the `AdamOptimizer`:<jupyter_code># opt = qml.GradientDescentOptimizer(stepsize=0.1)
opt = qml.AdamOptimizer(stepsize=0.1)<jupyter_output><empty_output><jupyter_text>We are ready to kick off the experiments! It might take a while with the default parameters :)<jupyter_code>%%time

for i in range(steps):
    params = opt.step(cost, params)

    if i == 0:
        print(f'\tCost after step {i:4d}: {cost(params)[0]: .7f} ({params})')
    elif (i + 1) % 50 == 0:
        print(f'\tCost after step {i+1:4d}: {cost(params)[0]: .7f} ({params})')

    params_history[:, i] = params

result = qml.map(circuit, observables, dev)(params)

print('Optimized parameters: {}'.format(params))
print('Result: {}'.format(result))

# The state cannot be seen when using a real device
try:
    print('Output state: {}'.format(np.round(dev.state, decimals=3)))
except NotImplementedError:
    print('Cannot see state when using device ' + dev.name)<jupyter_output>Cost after step    0:  0.1814760 ([0.68539817 0.88539815 0.68539818 0.6853982 ])
	Cost after step   50:  0.1482250 ([ 0.22836194  1.6914895  -0.01949193  0.05238423])
	Cost after step  100:  0.1592010 ([ 0.10961997  1.68426687 -0.03884338  0.06502656])
	Cost after step  150:  0.1095610 ([ 0.12908711  1.70756227 -0.03640683  0.08619756])
	Cost after step  200:  0.1128960 ([ 0.26982714  1.66307079 -0.07512116 -0.01178008])
	Cost after step  250:  0.1513210 ([0.21250233 1.67784563 0.01258101 0.19851932])<jupyter_text>EvaluationLet's check the results of the last step of the last optimization:<jupyter_code>print(params_history[:, -1])<jupyter_output>[0. 0. 0. 0.]<jupyter_text>We can see that each row shows values around $\left[\begin{matrix}0&\dfrac{\pi}{2}&0&0\end{matrix}\right]$, corresponding to the following rotations:- $R_x(\phi) = \left[\begin{matrix}cos(\phi/2)&-i sin(\phi/2)\\-i sin(\phi/2)&cos(\phi/2)\end{matrix}\right] \implies R_x(0) = \left[\begin{matrix}1&0\\0&1\end{matrix}\right]$- $R_y(\phi) = \left[\begin{matrix}cos(\phi/2)&-sin(\phi/2)\\sin(\phi/2)&cos(\phi/2)\end{matrix}\right] \implies R_y(0) = \left[\begin{matrix}1&0\\0&1\end{matrix}\right], R_y(\pi/2) = \dfrac{1}{\sqrt{2}}\left[\begin{matrix}1&-1\\1&1\end{matrix}\right]$ This means that only a $\dfrac{\pi}{2}$ Y-rotation should be applied on the first qubit, while the second qubit should be left unchanged.  VisualizationIn this section we will see the calculated values across runs for every parameter. We will plot the points using `matplotlib`.<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline
# TODO: can be simplified
fig, axs = plt.subplots(2, 2, figsize=(15,15))

axs[0, 0].plot(params_history[0, :])
axs[0, 0].set_title('RX(0)')

axs[0, 1].plot(params_history[1, :])
axs[0, 1].set_title('RY(0)')

axs[1, 0].plot(params_history[2, :])
axs[1, 0].set_title('RX(1)')    

axs[1, 1].plot(params_history[3, :])
axs[1, 1].set_title('RY(1)')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12241786721548899, 'fraction_numerical': 0.04276240567116396, 'mean_word_length': 5.102325581395349, 'pattern_counts': {'":': 0, '<': 78, '<?xml version=': 0, '>': 78, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, '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': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '22287', 'license': 'MIT', 'n_tokens_mistral': 4758, 'n_tokens_neox': 4325, 'n_words': 1422, 'path': 'qpu_results/state_preparation_circuit_optimization-1000_shots.ipynb', 'repo_name': 'nvitucci/notebook_qc', 'size': '13119'}