.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tn_simulation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_tn_simulation.py: Usage and application of tensor network simulations ===================================================== In this example, we will illustrate the usage and application of TN simulation of MBQC. Firstly, let's import the relevant modules: .. GENERATED FROM PYTHON SOURCE LINES 11-24 .. code-block:: default from functools import reduce import cotengra as ctg import matplotlib.pyplot as plt import networkx as nx import numpy as np import opt_einsum as oe import quimb.tensor as qtn from scipy.optimize import minimize from graphix import Circuit .. GENERATED FROM PYTHON SOURCE LINES 25-28 This application will be for the QAOA (Quantum Approximate Optimization Algorithm), which is an algorithm that can be used for example to solve combinatorical optimization problems. For this example a max cut problem will be solved on a star-like graph, so we can easier validate the results. .. GENERATED FROM PYTHON SOURCE LINES 30-31 Let's start with defining a helper function for buidling the circuit. .. GENERATED FROM PYTHON SOURCE LINES 31-43 .. code-block:: default def ansatz(circuit, n, gamma, beta, iterations): for j in range(0, iterations): for i in range(1, n): circuit.cnot(i, 0) circuit.rz(0, gamma[j]) circuit.cnot(i, 0) for i in range(0, n): circuit.rx(i, beta[j]) .. GENERATED FROM PYTHON SOURCE LINES 44-45 Let's look at how the quantum circuit is going to be built. .. GENERATED FROM PYTHON SOURCE LINES 45-60 .. code-block:: default n = 5 # This will result in a graph that would be too large for statevector simulation. iterations = 2 # Define the number of iterations in the quantum circuit. gamma = 2 * np.pi * np.random.rand(iterations) beta = 2 * np.pi * np.random.rand(iterations) # Define the circuit. circuit = Circuit(n) ansatz(circuit, n, gamma, beta, iterations) # Transpile Circuit into pattern as it is needed for creating the TN. pattern = circuit.transpile(opt=True).pattern # Optimizing according to standardization algorithm of graphix. pattern.standardize() pattern.shift_signals() .. rst-class:: sphx-glr-script-out .. code-block:: none {0: set(), 5: set(), 7: set(), 6: {0, 7}, 8: {5}, 9: {0, 6, 7}, 10: {8, 5}, 12: set(), 11: {0, 6, 7, 9, 12}, 13: {8, 10, 5}, 14: {0, 6, 7, 9, 11, 12}, 15: {8, 10, 5, 13}, 17: set(), 16: {0, 6, 7, 9, 11, 12, 14, 17}, 18: {5, 8, 10, 13, 15}, 19: {0, 6, 7, 9, 11, 12, 14, 16, 17}, 20: {5, 8, 10, 13, 15, 18}, 22: set(), 21: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 22}, 23: {5, 8, 10, 13, 15, 18, 20}, 24: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22}, 25: {5, 8, 10, 13, 15, 18, 20, 23}, 1: {6, 7}, 27: set(), 2: {11, 12}, 29: set(), 3: {16, 17}, 31: set(), 4: {21, 22}, 33: set(), 26: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24}, 35: {5, 8, 10, 13, 15, 18, 20, 23, 25, 27}, 37: set(), 36: {0, 37, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26}, 38: {35, 5, 8, 10, 13, 15, 18, 20, 23, 25}, 39: {0, 36, 37, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26}, 40: {35, 5, 38, 8, 10, 13, 15, 18, 20, 23, 25, 29}, 42: set(), 41: {0, 36, 37, 6, 39, 7, 9, 42, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26}, 43: {35, 5, 38, 40, 8, 10, 13, 15, 18, 20, 23, 25}, 44: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26, 36, 37, 39, 41, 42}, 45: {35, 5, 38, 40, 8, 10, 43, 13, 15, 18, 20, 23, 25, 31}, 47: set(), 46: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26, 36, 37, 39, 41, 42, 44, 47}, 48: {35, 5, 38, 40, 8, 10, 43, 45, 13, 15, 18, 20, 23, 25}, 49: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26, 36, 37, 39, 41, 42, 44, 46, 47}, 50: {33, 35, 5, 38, 40, 8, 10, 43, 45, 13, 15, 48, 18, 20, 23, 25}, 52: set(), 51: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26, 36, 37, 39, 41, 42, 44, 46, 47, 49, 52}, 53: {35, 5, 38, 40, 8, 10, 43, 45, 13, 15, 48, 50, 18, 20, 23, 25}, 54: {0, 6, 7, 9, 11, 12, 14, 16, 17, 19, 21, 22, 24, 26, 36, 37, 39, 41, 42, 44, 46, 47, 49, 51, 52}, 55: {35, 5, 38, 40, 8, 10, 43, 45, 13, 15, 48, 50, 18, 20, 53, 23, 25}, 28: {1, 36, 37, 6, 7}, 57: {27}, 30: {2, 11, 12, 41, 42}, 59: {29}, 32: {3, 16, 17, 46, 47}, 61: {31}, 34: {4, 21, 22, 51, 52}, 63: {33}} .. GENERATED FROM PYTHON SOURCE LINES 61-62 Print some properties of the graph. .. GENERATED FROM PYTHON SOURCE LINES 62-67 .. code-block:: default nodes, edges = pattern.get_graph() print(f"Number of nodes: {len(nodes)}") print(f"Number of edges: {len(edges)}") .. rst-class:: sphx-glr-script-out .. code-block:: none Number of nodes: 65 Number of edges: 76 .. GENERATED FROM PYTHON SOURCE LINES 68-69 Optimizing by performing Pauli measurements in the pattern using efficient stabilizer simulator. .. GENERATED FROM PYTHON SOURCE LINES 69-72 .. code-block:: default pattern.perform_pauli_measurements(use_rustworkx=True) .. GENERATED FROM PYTHON SOURCE LINES 73-76 Simulate using the TN backend of graphix, which will return an MBQCTensorNet object. The graph_prep argument is optional, but with 'parallel' the TensorNetworkBackend will prepeare the graph state faster. .. GENERATED FROM PYTHON SOURCE LINES 76-81 .. code-block:: default mbqc_tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="parallel") sv = mbqc_tn.to_statevector().flatten() print("Statevector after the simulation:", sv) .. rst-class:: sphx-glr-script-out .. code-block:: none Statevector after the simulation: [-0.25675858-0.09162405j 0.09480265-0.05174428j 0.09480265-0.05174428j -0.07759243+0.04417738j 0.09480265-0.05174428j -0.07759243+0.04417738j -0.07759243+0.04417738j -0.00359709+0.14241443j 0.09480265-0.05174428j -0.07759243+0.04417738j -0.07759243+0.04417738j -0.00359709+0.14241443j -0.07759243+0.04417738j -0.00359709+0.14241443j -0.00359709+0.14241443j 0.16638679-0.47151125j 0.16638679-0.47151125j -0.00359709+0.14241443j -0.00359709+0.14241443j -0.07759243+0.04417738j -0.00359709+0.14241443j -0.07759243+0.04417738j -0.07759243+0.04417738j 0.09480265-0.05174428j -0.00359709+0.14241443j -0.07759243+0.04417738j -0.07759243+0.04417738j 0.09480265-0.05174428j -0.07759243+0.04417738j 0.09480265-0.05174428j 0.09480265-0.05174428j -0.25675858-0.09162405j] .. GENERATED FROM PYTHON SOURCE LINES 82-85 Let's explore what is really happening, how the TN is being constructed. The tensor network is created using the graph structure, so from the list of nodes as well as the edges. The graph must be preprocessed before the constraction of the TN itself. .. GENERATED FROM PYTHON SOURCE LINES 87-92 The goal is to represent quantum states, so for every node a list is created, which stores the index labels for each dimension for that given node. Its length will be one larger than the number of edges of the node. This is due to the fact that the first entry of the list represents the "dangling" index, which corresponds to the physical index of the qubit (i.e., the index that represents the local Hilbert space of the qubit). The following entries in the list are then correspond to neighbouring tensors, and can be contracted with them. For additional details and visualization visit: https://quimb.readthedocs.io/en/latest/tensor-1d.html. .. GENERATED FROM PYTHON SOURCE LINES 94-96 Let's take a closer look at an MPS tensor (left plot) and an MPS tensor network that consists of two MPS tensors (right plot). By the network on the right the middle index is shared between the two tensors, essentially allowing for contraction between them by summing over it. .. GENERATED FROM PYTHON SOURCE LINES 96-112 .. code-block:: default t = qtn.rand_tensor([2], "a") fig, ax = plt.subplots(1, 2) t.draw( ax=ax[0], title="MPS tensor", legend=False, ) t1 = qtn.rand_tensor([2, 2], ["a", "b"]) t2 = qtn.rand_tensor([1, 2], ["b", "c"]) t1.add_tag("T1") t2.add_tag("T2") t = qtn.TensorNetwork([t1, t2]) t.draw(ax=ax[1], title="MPS", legend=False, color=["T1", "T2"]) plt.show() .. image-sg:: /gallery/images/sphx_glr_tn_simulation_001.png :alt: MPS tensor, MPS :srcset: /gallery/images/sphx_glr_tn_simulation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 113-120 Additionally, the type of edges are also stored, in a binary valued list for each node. These are used to construct the tensor itself. From each node in the graph a tensor is constructed, which has a dimension that is exactly one larger than its neighbour count. The tensor is described using two outer products, for which the list from above is used, that describes the edges for every node. For additional information on TN construction please refer to: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.76.052315 . Section III A provides further information on Matrix Porduct States and section III C gives an example using a 1-D cluster state. In section IV novel resource states are explored, where parts A, B can be used for getting a deeper understanding. .. GENERATED FROM PYTHON SOURCE LINES 122-124 Let's also plot the resulting tensor network (notice that there are five dangling edges, which is exactly the number of qubits that were defined in the quantum circuit). Here open means that the node has a dangling index, but the "dangling" edge itself is not drawn, rather the tensors are color coded. .. GENERATED FROM PYTHON SOURCE LINES 124-138 .. code-block:: default fig, ax = plt.subplots(figsize=(13, 10)) color = ["Open", "Close"] mbqc_tn.draw( ax=ax, color=color, show_tags=False, show_inds=False, iterations=100, k=0.01, ) plt.show() .. image-sg:: /gallery/images/sphx_glr_tn_simulation_002.png :alt: tn simulation :srcset: /gallery/images/sphx_glr_tn_simulation_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 139-140 Let's calculate the measuring probability corresponding to the first basis state. .. GENERATED FROM PYTHON SOURCE LINES 140-144 .. code-block:: default value = mbqc_tn.get_basis_amplitude(0) print(f"Probability for {0} is {value}") .. rst-class:: sphx-glr-script-out .. code-block:: none Probability for 0 is 0.07431993491418348 .. GENERATED FROM PYTHON SOURCE LINES 145-147 It is also possible to change the path contraction algorithm. Let's explore that too and define a custom optimizer for contraction, that we can use later. .. GENERATED FROM PYTHON SOURCE LINES 147-154 .. code-block:: default opt = ctg.HyperOptimizer( minimize="combo", reconf_opts={}, progbar=True, ) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/graphix/envs/v0.2.16/lib/python3.10/site-packages/cotengra/hyperoptimizers/hyper.py:54: UserWarning: Couldn't find `optuna`, `baytune (btb)`, `chocolate`, `nevergrad` or `skopt` so will use completely random sampling in place of hyper-optimization. warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 155-157 Let's also calculate the expectation value for the measurement in the computational basis. The expectation value can be optiained using a function of graphix. .. GENERATED FROM PYTHON SOURCE LINES 157-165 .. code-block:: default pauli_z = np.array([[1, 0], [0, -1]]) identity = np.array([[1, 0], [0, 1]]) operator = reduce(np.kron, [pauli_z] * n) # Use the defined optimizer by setting the 'optimize' parameter. exp_val = mbqc_tn.expectation_value(operator, range(n), optimize=opt) print("Expectation value for Z^n: ", exp_val) .. rst-class:: sphx-glr-script-out .. code-block:: none Expectation value for Z^n: (6.33410273312767e-17+2.3949689821869273e-17j) .. GENERATED FROM PYTHON SOURCE LINES 166-169 If we want to find the solution for our initial max-cut problem, then we must deploy a classical minimizer too for an apropriate cost function. Create a cost function using the elements of graphix, which were already discussed above. .. GENERATED FROM PYTHON SOURCE LINES 169-188 .. code-block:: default def cost(params, n, ham, quantum_iter, slice_index, opt=None): circuit = Circuit(n) gamma = params[:slice_index] beta = params[slice_index:] ansatz(circuit, n, gamma, beta, quantum_iter) pattern = circuit.transpile(opt=True).pattern pattern.standardize() pattern.shift_signals() pattern.perform_pauli_measurements(use_rustworkx=True) mbqc_tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="parallel") exp_val = 0 for op in ham: exp_val += np.real(mbqc_tn.expectation_value(op, range(n), optimize=opt)) return exp_val .. GENERATED FROM PYTHON SOURCE LINES 189-190 We want to find the ground state energy for the Hamiltonian :math:`\hat{H} = \sum \hat{Z}_k + \sum \hat{Z}_i \hat{Z}_j` with i,j running over the edges. .. GENERATED FROM PYTHON SOURCE LINES 190-213 .. code-block:: default ham = [reduce(np.kron, [pauli_z] * n)] for i in range(1, n): op = [identity] * n op[0] = pauli_z op[i] = pauli_z op = reduce(np.kron, op) ham.append(op) # Use yet again another optimizer for path contraction. class MyOptimizer(oe.paths.PathOptimizer): def __call__(self, inputs, output, size_dict, memory_limit=None): return [(0, 1)] * (len(inputs) - 1) opt = MyOptimizer() # Define initial parameters, which will be optimized through running the algorithm. params = 2 * np.pi * np.random.rand(len(gamma) + len(beta)) # Run the classical optimizer and simulate the quantum circuit with TN backend. res = minimize(cost, params, args=(n, ham, iterations, len(gamma), opt), method="COBYLA") print(res.message) .. rst-class:: sphx-glr-script-out .. code-block:: none Optimization terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 214-215 Finally, run the circuit once again with the optimized parameters. .. GENERATED FROM PYTHON SOURCE LINES 215-223 .. code-block:: default circuit = Circuit(n) ansatz(circuit, n, res.x[: len(gamma)], res.x[len(gamma) :], iterations) pattern = circuit.transpile(opt=True).pattern pattern.standardize() pattern.shift_signals() mbqc_tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="parallel") .. GENERATED FROM PYTHON SOURCE LINES 224-225 Let's use the defined optimizer and find the most probable basis states. .. GENERATED FROM PYTHON SOURCE LINES 225-233 .. code-block:: default max_prob = 0 most_prob_state = 0 bars = [] for i in range(0, 2**n): value = mbqc_tn.get_basis_amplitude(i) bars.append(value) .. GENERATED FROM PYTHON SOURCE LINES 234-235 Plot the output. .. GENERATED FROM PYTHON SOURCE LINES 235-244 .. code-block:: default fig, ax = plt.subplots(figsize=(10, 5)) ax.bar(range(0, 2**n), bars, color="maroon", width=0.2) ax.set_xticks(range(0, 2**n)) ax.set_xlabel("States") ax.set_ylabel("Probabilites") ax.set_title("Measurement probabilities using the optimized parameters") plt.show() .. image-sg:: /gallery/images/sphx_glr_tn_simulation_003.png :alt: Measurement probabilities using the optimized parameters :srcset: /gallery/images/sphx_glr_tn_simulation_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 245-249 As we can see the most probable are 15 and 16 ( ``|11110>`` and ``|00001>`` because of bit ordering), which mean that splitting the graph so that node number 0 is in one set, and all other nodes in the other solves the max cut problem. This result is what we would expect from this star-like graph. .. GENERATED FROM PYTHON SOURCE LINES 251-253 The following illustration shows the starting graph on the left, and the graph with the resulting sets found on the right, where the nodes with different colours belong to different groups. .. GENERATED FROM PYTHON SOURCE LINES 253-264 .. code-block:: default fig, ax = plt.subplots(ncols=2, figsize=(8, 6)) ax = ax.flatten() g = nx.Graph() for i in range(1, n): g.add_edge(0, i) color = ["blue"] * n color[0] = "red" nx.draw(g, ax=ax[0], with_labels=True) nx.draw(g, ax=ax[1], node_color=color, with_labels=True) plt.show() .. image-sg:: /gallery/images/sphx_glr_tn_simulation_004.png :alt: tn simulation :srcset: /gallery/images/sphx_glr_tn_simulation_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 3 minutes 43.192 seconds) .. _sphx_glr_download_gallery_tn_simulation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tn_simulation.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tn_simulation.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_