From f26798363eb67b38b584c65e0a0aa3fbf92c39c3 Mon Sep 17 00:00:00 2001 From: mat70593 <matic.petric@fokus.fraunhofer.de> Date: Tue, 6 Feb 2024 11:30:21 +0100 Subject: [PATCH 1/3] Fixed the solve_QUBO function to return the optimal_solution --- .../MaxCut.rst | 339 ------------------ .../QUBO.rst | 170 --------- .../source/reference/Examples/Shor.rst | 110 ------ src/qrisp/qaoa/problems/QUBO.py | 21 +- 4 files changed, 9 insertions(+), 631 deletions(-) delete mode 100644 documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst delete mode 100644 documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst delete mode 100644 documentation/source/reference/Examples/Shor.rst diff --git a/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst b/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst deleted file mode 100644 index cce390f4..00000000 --- a/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst +++ /dev/null @@ -1,339 +0,0 @@ -.. _MaxCutQAOA: - -MaxCut QAOA Implementation -========================== - -Welcome to our Quantum Approximate Optimization Algorithm (QAOA) implementation tutorial for the MaxCut problem! This tutorial is designed to provide you with a hands-on understanding of how to apply QAOA using Qrisp to solve the MaxCut optimization problem. - -To refresh your memory, `Maximum Cut, or MaxCut for short <https://en.wikipedia.org/wiki/Maximum_cut>`_ involves dividing the nodes of a graph into two groups such that the number of edges between the groups is maximized. It has applications in various fields, including computer science, operations research, and statistical physics. - -In this interactive tutorial, we’ll first tackle the MaxCut problem to build a solid QAOA foundation. You'll be guided to learn how to formulate the problem for QAOA, implement the actual algorithm, and interpret the results. - -Whether you’re a quantum computing novice or an experienced algorithm developer checking out this remarkable new framework, you may have observed that virtually every other quantum programming framework and their quantum grandmother has implemented QAOA for MaxCut in a similar manner (authors used irony, it was super effective). We'll go one step further and showcase on why using the Qrisp QAOA module and its :ref:`QAOAProblem class <QAOAProblem>` **is** a no-brainer when it comes to simplicity of implementation as well as the speed of the optimization itself. - -Your very first QAOA layer --------------------------- - -Let's start slowly with deepening the understanding of QAOA at its core - one layer of the algorithm and see how it looks like. We start with importing necessary quantum gates, as well as `networkx`, which is a package for creating and manipulating graphs. In this block of code you also define the angle parameters $\gamma$ and $\beta$: -:: - from qrisp import QuantumVariable, h, barrier, rz, rx, cx - import networkx as nx - from sympy import Symbol - - G = nx.Graph() - G.add_edges_from([[0,3],[0,4],[1,3],[1,4],[2,3],[2,4]]) - #nx.draw(G, with_labels=True) - N = G.number_of_nodes() - -Feel free to play around with this block of code changing the graph through the number of nodes, the edges connecting them, and print the result by uncommenting ``nx.draw(G, with_labels=True)``. - -Next thing on our list is defining the MaxCut cost operator and the mixer. Taking a peek in the appendices section of Hadfield's paper, we see, the MaxCut phase separator we have to implement is defined as -$$U_P=e^{-i\\gamma_pH_P}=\e^{-i\\gamma_p\\sum_{ij\\in E}(-Z_iZ_j)}.$$ Instead of using the `ZZ` gate, we will perform the same operation using only elementary gates. -:: - gamma = Symbol('γ') - - def apply_cost_operator(qv, gamma): - for pair in list(G.edges()): - cx(qv[pair[0]], qv[pair[1]]) - rz(2 * gamma, qv[pair[1]]) - cx(qv[pair[0]], qv[pair[1]]) - barrier(qv) - return qv - -The first line starts a loop over all edges in the graph ``G``. Each edge is a pair of nodes, which we refer to as ``pair[0]`` and ``pair[1]``. We then apply a controlled-X (CNOT) gate to the pair of qubits before applying a rotation around the Z-axis by an angle of $2\gamma$ to the second qubit in the pair. At the end we apply another CNOT gate to the pair of qubits. - -As suggested in the table of our QAOA theoretical overview, we are using the ``rx`` gate as the mixer. -:: - beta = Symbol('β') - - def apply_mixer(qv, beta): - for i in range(0, N): - rx(2 * beta, qv[i]) - barrier(qv) - return qv - -With these two functions at our disposal we can now implement our first QAOA layer and `print(qv.qs)` it. -:: - qv_1 = QuantumVariable(N) - h(qv_1) - apply_cost_operator(qv_1, gamma) - apply_mixer(qv_1, beta) - print(qv_1.qs) - -Applying hadamard gates on ``qv`` prepares the initial state for our system - the superposition state $|s\rangle$ . - -If you have executed all the preceding cells correctly, you should have been able to print and examine the circuit of the initial QAOA layer. You may notice a pattern in the phase separator section of the circuit, which links two qubits corresponding to a specific pair of nodes in `G` that are connected by an edge (commencing flashback: `G.add_edges_from([[0,3],[0,4],[1,3],[1,4],[2,3],[2,4]])`). The use of barriers between such pairs facilitates a more intuitive understanding of the phase separator unitary's inner workings. - -We’ve used the parameters $\gamma$ and $\beta$ symbolically, without assigning specific numerical values to them. In the next section we will build on top of the first layer, assign values to the angle parameters and iteratively optimize their value with each additional layer (SPOILER ALERT). - -As we add more layers to our QAOA circuit, the values of $\gamma$ and $\beta$ are refined, leading us closer to the optimal solution. This iterative improvement is the core of what makes QAOA such a powerful tool. - -MaxCut QAOA implementation --------------------------- - -Putting focuss on the problem at hand (MaxCut if you're just scrolling through and not paying that much attention) it is reasonable to ask ourselves what will be the objective that we want to optimize. - -As the name suggests, we are looking for a cut going through as much edges as possible. It's therefore crucial to count the amount of cuts we cut through (I'll cut it out now) by checking if the nodes ``i`` and ``j`` belong to different groups. If they are, the edge $(i,j)$ is cut and added to the total. -:: - def maxcut_obj(x): - cut = 0 - for i, j in G.edges(): - if x[i] != x[j]: - cut -= 1 - return cut - -This optimization objective is important for the last building blocks described in the theoretical overview we haven't mentioned yet - the cost function. The cost function is important for keeping track of the energy of the system: -:: - def maxcut_cost_funct(meas_res): - energy = 0 - for meas, p in meas_res.items(): - obj_for_meas = maxcut_obj(meas) - energy += obj_for_meas * p - return energy - -We loop over the measurement in ``meas_res``. Each solution is represented by a string of 1 and 0 in ``meas`` with probability ``p`` keeping score in how likely a particular solution appears. - -For each such solution the cost is calculated using the ``maxcut_obj`` returning the number of cut edges. At the end, ``maxcut_cost_funct`` calculates the average “quality” of a set of solutions to the MaxCut problem. A lower average energy means that, on average, the solutions are better - they cut more edges. - -This is now nearly all the building blocks we need in order to run QAOA and start optimizing. Well, after we add more layers to our algorithm, that is. -:: - p = 3 - - def apply_p_layers(qv, beta, gamma): - assert(len(beta) == len(gamma)) - p = len(beta) - h(qv) - for i in range(p): - apply_cost_operator(qv, gamma[i]) - apply_mixer(qv, beta[i]) - barrier(qv) - return qv - -With our $p$-layered algorithm in hand there is one last piece of the puzzle: the objective function, which we need to generate in order to calculate the average energy of the resulting solutions. Oh, and it features the angle parameters $\gamma$ and $\beta$, which we combine in one list ``theta``. The solutions we get depend on the values of the angle parameters. -:: - - def quantum_objective(theta): - qv_p = QuantumVariable(N) - beta = theta[:p] - gamma = theta[p:] - qv = apply_p_layers(qv_p,beta, gamma) - results = qv.get_measurement() - return maxcut_cost_funct(results) - - -We can finally finish combine all the pieces of the mosaic and have it appraised by a professional - a classical optimizer. Running it online might take a second so we would like to encourage you to download Qrisp and try out for yourself! - -We first reset the initial state and randomize it in a list where the first $p$ values correspond to $\beta$ and the second one to $\gamma$. It has finally come the time to inclute the COBYLA optimiser, which will return the optimal parameters that minimze our objective function. This is done using the ``minimize`` method, which adjusts the parameters iteratively until it finds the best ones. - -Then, once we have our optimal parameters, we apply QAOA one more time to get our final set of solutions (``counts``), before finding the best solution by looking for the one with the lowest energy. Finally, we visualize our solution by coloring the nodes of our graph according to which group they belong to in the best solution. -:: - import numpy as np - from scipy.optimize import minimize - from operator import itemgetter - - init_point = np.pi * np.random.rand(2 * p) - - res_sample = minimize(quantum_objective, init_point, method='COBYLA', options={'maxiter':50}) - - optimal_theta = res_sample['x'] - qv_p = QuantumVariable(N) - qv = apply_p_layers(qv_p, optimal_theta[:p], optimal_theta[p:]) - counts = qv_p.get_measurement() - - best_cut, best_solution = min([(maxcut_obj(x),x) for x in counts.keys()], key=itemgetter(0)) - print(f"Best string: {best_solution} with cut: {-best_cut}") - - colors = ['r' if best_solution[node] == '0' else 'b' for node in G] - nx.draw(G,node_color = colors, pos=nx.bipartite_layout(G, [0,1,2])) - -And voila! We just solved the MaxCut problem using a relatively straightforward brute force approach. While this method can be effective, it’s not the most efficient or elegant way to tackle this problem, and other problem instances with more complex objectives, mixers, or both. - -That’s where our QAOAProblem class comes in. This class simplifies the process of running QAOA by introducing modularity, without sacrificing the generality of whichever problem instance we're dealing with. As we move forward in this tutorial, you’ll have the opportunity to try it out for yourself. - -QAOAProblem. Enough said. -------------------------- - -🎶 Enough chit chat, you ain't got all day - let's get to it: QAOAProblem, lead the way! 🎶 - -Taking the essential building blocks from `QAOAnsatz <https://arxiv.org/abs/1709.03489>`_ into account, we built the :ref:`QAOAProblem class <QAOA>` with modularity in mind. We gathered and improved the functions shown in the example above and created a powerfull architecture with which it's easy to implement various problem instances of varying complexity.In this tutorial we focus on using QAOAProblem, with extensive documentation being available :ref:`here <QAOA>`. - -We start with renaming our quantum argument `qv` to a more general `qarg` because more often than not we'll combine QuantumVariables into a QuantumArray to make implementations of other problem instances more efficient. -:: - qarg = QuantumVariable(len(G)) - - depth = 3 - -BEHOLD, THE POWER OF QRISP! -:: - from qrisp.qaoa import QAOAProblem, RX_mixer - - maxcut_instance = QAOAProblem(apply_cost_operator, RX_mixer, maxcut_cost_funct) - - res = maxcut_instance.run(qarg, depth, max_iter = 50) # runs the simulation - -And that's pretty much it, really. Apart from visualizing the results again. -:: - best_cut, best_solution = min([(maxcut_obj(x),x) for x in res.keys()], key=itemgetter(0)) - print(f"Best string: {best_solution} with cut: {-best_cut}") - - res_str = list(res.keys())[0] - print("QAOA solution: ", res_str) - best_cut, best_solution = (maxcut_obj(res_str),res_str) - - colors = ['r' if best_solution[node] == '0' else 'b' for node in G] - nx.draw(G,node_color = colors) - -After thorough comparison certainly made after running both approaches for different graph topologies you are in position to be the judge regarding which approach is qrispier. - -.. _benchmark_maxcut: -Benchmarking the performance ----------------------------- -Lucky for you, the tutorial is not over since there is one more important functionality we would like to show you - the :meth:`benchmark <qrisp.qaoa.QAOAProblem.benchmark>` method of our QAOA module. This method allows you to observe the performance of your simulation: - -- get insights about the approximation ratio for each run together with the average approximation ratio with its variance, -- obtain total circuit depth, -- obtain the overall quantum volume, -- and rank the results for different depths, amount of shots, as well as iterations. - -Let's show how easy benchmarking QAOA is using Qrisp: -:: - print('RND') - - benchmark_data = maxcut_instance.benchmark(qarg = QuantumVariable(len(G)), - depth_range = [1,2,3,4,5], - shot_range = [1000, 10000], - iter_range = [100, 200], - optimal_solution = '00011', - repetitions = 1, - init_type='random' - ) - - temp = benchmark_data.rank(print_res = True) - - _,rndFO=benchmark_data.evaluate() - - print('Approximation ratio: ',sum(rndFO)/len(rndFO)) - print("Variance: ",np.var(rndFO)) - -In the above example we obtain the list of results for the approximation ratio for various depths of QAOA ranging from 1 to 5. We run the benchmark for both 1000, and 10000 shots (the latter number is to avoid sampling of the solution); and for 100 and 200 iterations of the optimizer. Since QAOA is probabilistic every run is unique so it's important to run more repetitions for clearer insights. - -It is important to note that in order to do the benchmark, one has to already know the optimal solution in order to calculate the optimal energy, and through that, the approximation ratio. This is an example of the output for the benchmark above: -:: - RND - Rank approx_ratio Overall QV p QC depth QB count Shots Iterations - ============================================================================ - 0 0.990e+0 4.4000e+8 5 44 5 10000 200 - 1 0.977e+0 4.4000e+7 5 44 5 1000 200 - 2 0.974e+0 3.6000e+8 4 36 5 10000 200 - 3 0.956e+0 1.8000e+8 4 36 5 10000 100 - 4 0.945e+0 2.2000e+7 5 44 5 1000 100 - 5 0.940e+0 2.2000e+8 5 44 5 10000 100 - 6 0.940e+0 3.6000e+7 4 36 5 1000 200 - 7 0.871e+0 2.0000e+8 2 20 5 10000 200 - 8 0.830e+0 2.8000e+8 3 28 5 10000 200 - 9 0.817e+0 1.4000e+7 3 28 5 1000 100 - 10 0.815e+0 1.4000e+8 3 28 5 10000 100 - 11 0.775e+0 2.8000e+7 3 28 5 1000 200 - 12 0.774e+0 1.8000e+7 4 36 5 1000 100 - 13 0.765e+0 2.0000e+7 2 20 5 1000 200 - 14 0.764e+0 1.0000e+7 2 20 5 1000 100 - 15 0.736e+0 1.0000e+8 2 20 5 10000 100 - 16 0.718e+0 1.2000e+8 1 12 5 10000 200 - 17 0.705e+0 1.1000e+7 1 11 5 1000 200 - 18 0.556e+0 6.0000e+6 1 12 5 1000 100 - 19 0.555e+0 5.5000e+7 1 11 5 10000 100 - Approximation ratio: 0.8201699999999998 - Variance: 0.016184990766666633 - -As expected, on average the runs with higher depths yield better approximation ratio (with some outliers, of course). - -The ``print('RND')`` was used because since the 0.4 update we have also included the `TQA warm-start initialization <https://arxiv.org/abs/2101.05742>`_, which can be used within the :meth:`benchmark <qrisp.qaoa.QAOAProblem.benchmark>` method by setting ``init_type='tqa'``. Let's try benchmarking this approach as well: -:: - print('TQA') - maxcut_instance = maxcut_problem(G) - - benchmark_data = maxcut_instance.benchmark(qarg = QuantumVariable(len(G)), - depth_range = [1,2,3,4,5], - shot_range = [1000, 10000], - iter_range = [100, 200], - optimal_solution = '00011', - repetitions = 1, - init_type='tqa' - ) - - temp = benchmark_data.rank(print_res = True) - - _,tqaFO=benchmark_data.evaluate() - - print('Approximation ratio: ',sum(tqaFO)/len(tqaFO)) - print("Variance: ",np.var(tqaFO)) - -Again, since QAOA is probabilistic, every run returns different results. You can find our try below: -:: - TQA - Rank approx_ratio Overall QV p QC depth QB count Shots Iterations - ============================================================================ - 0 0.989e+0 4.4000e+8 5 44 5 10000 200 - 1 0.983e+0 2.2000e+8 5 44 5 10000 100 - 2 0.963e+0 1.8000e+8 4 36 5 10000 100 - 3 0.963e+0 3.6000e+8 4 36 5 10000 200 - 4 0.957e+0 4.4000e+7 5 44 5 1000 200 - 5 0.897e+0 1.4000e+8 3 28 5 10000 100 - 6 0.897e+0 2.8000e+8 3 28 5 10000 200 - 7 0.892e+0 2.8000e+7 3 28 5 1000 200 - 8 0.880e+0 1.0000e+8 2 20 5 10000 100 - 9 0.880e+0 2.0000e+8 2 20 5 10000 200 - 10 0.875e+0 1.4000e+7 3 28 5 1000 100 - 11 0.867e+0 1.8000e+7 4 36 5 1000 100 - 12 0.866e+0 1.0000e+7 2 20 5 1000 100 - 13 0.860e+0 2.0000e+7 2 20 5 1000 200 - 14 0.847e+0 2.2000e+7 5 44 5 1000 100 - 15 0.822e+0 3.6000e+7 4 36 5 1000 200 - 16 0.717e+0 6.0000e+7 1 12 5 10000 100 - 17 0.717e+0 1.2000e+8 1 12 5 10000 200 - 18 0.715e+0 5.5000e+6 1 11 5 1000 100 - 19 0.691e+0 1.1000e+7 1 11 5 1000 200 - Approximation ratio: 0.8639683333333332 - Variance: 0.007950404719444411 - -As we can see, the TQA initialization tends to return better approximation ratios more consistently. Feel free to play around with the :meth:`benchmark <qrisp.qaoa.QAOAProblem.benchmark>` method (or leave it running over night) to try and compare the two approaches further. - -Summary and motivation ----------------------- - -To shortly summarize, in order to implement QAOA with QAOAProblem one needs to specify the problem using the following recipe - -I. define **CLASSICAL COST FUNCTION** of the problem you want to implement: ``maxcut_cost_funct(counts)``, -II. define the **INITIAL STATE** if it is not the superposition, which is set by default, -III. define **COST OPERATOR aka PHASE SEPARATOR** (or use the ones specified in `From QAOA to QAOA <https://arxiv.org/abs/1709.03489>`_) like ``apply_cost_operator`` above, and -IV. select **MIXER** from the :ref:`assortment we provide and list here <MIXers>`. - -Let's condense all of the above, and implement QAOA for MaxCut one last time in one block of code -:: - from qrisp.qaoa import QuantumArray, QuantumVariable, QAOAProblem, maxcut_obj,create_maxcut_cl_cost_function,create_maxcut_cost_operator, RX_mixer - import networkx as nx - from operator import itemgetter - - G = nx.Graph() - G.add_edges_from([[0,3],[0,4],[1,3],[1,4],[2,3],[2,4]]) - - qarg = QuantumArray(qtype = QuantumVariable(1), shape = len(G)) - - depth = 5 - - maxcut_instance = QAOAProblem(create_maxcut_cost_operator(G), RX_mixer, create_maxcut_cl_cost_function(G)) - - res = maxcut_instance.run(qarg, depth, max_iter = 50) - - best_cut, best_solution = min([(maxcut_obj(x,G),x) for x in res.keys()], key=itemgetter(0)) - - res_str = list(res.keys())[0] - print("QAOA solution: ", res_str) - best_cut, best_solution = (maxcut_obj(res_str,G),res_str) - - colors = ['r' if best_solution[node] == '0' else 'b' for node in G] - nx.draw(G,node_color = colors, pos=nx.bipartite_layout(G, [0,1,2])) - -You've got to admit that this is pretty cool, clean and qrispy, right?! - -If you are still not convinced, we provide a more complex problem instance in the next tutorial while also showcasing and putting some unique functionalities of Qrisp to the test. Let's make this transition a little more dramatic by saying that it's time to put our money where out mouths are (this is funny, because Qrisp is open source) and put this recipe to the test tackling the problem class which generalizes MaxCut: :ref:`Max-$\\kappa$-Colorable Subgraph <MkCSQAOA>`. \ No newline at end of file diff --git a/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst b/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst deleted file mode 100644 index 4ff60392..00000000 --- a/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst +++ /dev/null @@ -1,170 +0,0 @@ -.. _QUBOQAOA: - -QUBO as a QAOAProblem Problem Instance -====================================== - -In this tutorial you’ll be guided through the process of defining a new phase separator to be used within the scope of the :ref:`Alternating Operator Ansatz <AOA>` focussed on solving various QUBO problems with only needing the QUBO matrix $Q$ as an input. -QUBO, or `Quadratic Unconstrained Binary Optimization <https://en.wikipedia.org/wiki/Quadratic_unconstrained_binary_optimization>`_, is a type of problem that involves optimizing a quadratic function of binary variables. - -After first translating the QUBO Hamiltonian $H_C$ from the binary basis $x_i$ to the state basis with the translation $x_i\rightarrow \frac{\mathbb{I}-Z_i}{2}$, we'll construct the phase separator unitary and run QAOA to solve the set partitioning problem. - -Not to get too ahead of ourselves with how to do it, let's first show how to solve a QUBO problem using Qrisp. To obtain the optimal solutions to a problem, we only need to know the QUBO matrix $Q$. In order to get familiar with that, we propose having a look at the `Tutorial on Formulating and Using QUBO Models <https://arxiv.org/abs/1811.11538>`_, which explains how to derive such a matrix $Q$. - -One can then simply minimize the cost function $\min C(x)=\min x^TQx$ for binary variables $x_i\in\{0,1\}$. This is done either by sending the QUBO matrix $Q$ to the annealer/quantum computer, which calculates and returns the bitstrings with corresponding values of $C(x)$. The lower the $C(x)$, the better the solution. In other words, we want to find the Hamiltonian $H_C|x\rangle=C(x)|x\rangle$. - -Let's borrow the QUBO matrix (explained in depth in the above mentioned `tutorial <https://arxiv.org/abs/1811.11538>`_) for the set partitioning problem. The QUBO matrix in that case is - -$$Q = \\begin{pmatrix}-17&10&10&10&0&20\\\\10&-18&10&10&10&20\\\\10&10&-29&10&20&20\\\\10&10&10&-19&10&10\\\\0&10&20&10&-17&10\\\\20&20&20&10&10&-28\\end{pmatrix}$$ - -Usually, QUBO matrices are upper triangular (by convention) - this means that only elements above the diagonal (with the diagonal included) are not equal to zero. This is because the variables in QUBO problems are binary, which results in the identity $q_iq_j=q_jq_i$ for binary variables. Because of this, the elements above the diagonal are doubled compared to the symmetrized version of $Q$ we wrote down above. - -Of course it's easy to go from the conventional QUBO $Q_{\text{up}\Delta}$ formulation to the symmetrized $Q_\text{sym}$ QUBO: $Q_\text{sym}=\frac{1}{2}\big(Q_{\text{up}\Delta}+Q_{\text{up}\Delta}^T\big)$, and back: -:: - import numpy as np - Q_up_triang = 2 * np.triu(Q_sym) - np.diag(np.diag(Q_sym)) - - -Our implementation of solving QUBO problems using QAOA works for both the upper triangular, as well as symmetrized matrix conventions. As promised, we can now immediately solve the QUBO by calling the ``solve_QUBO`` function: -:: - from qrisp.qaoa import * - import numpy as np - - Q = np.array( - [ - [-17, 10, 10, 10, 0, 20], - [ 10, -18, 10, 10, 10, 20], - [ 10, 10, -29, 10, 20, 20], - [ 10, 10, 10, -19, 10, 10], - [ 0, 10, 20, 10, -17, 10], - [ 20, 20, 20, 10, 10, -28], - ] - ) - - solve_QUBO(Q, depth = 1, backend = def_backend, n_solutions=10) - -That's it! You can try running this block of code on the website using the Thebe interactivity integration, or run it on your own Qrispy environment. -We see that we obtain the 10 best solutions with corresponding bitsring of the binary variables. -You can test the statements above, and convert the symmetrized matrix into the upper triangular one, and run ``solve_QUBO`` again. - -We see, that we only needed to provide the matrix $Q$, specify the depth of the QAOA circuit that is running in the background, and the backend you want to run the algorithm on. Clean and qrispy! - -From QUBO to QAOA ------------------ - -Of course it's beneficial to not only know the "how", but also understand the "why". So let's dig in! - -To construct such a circuit with quantum gates and run it on a quantum computer, one has to translate between the basis of $x_i$ to Pauli gate operators. - -$$Z_i|x\\rangle=(-1)^{x_i}\\rightarrow\\frac{\\mathbb{1}-Z_i}{2}|x\\rangle=x_i$$ - -To find $H_C$ one can calculate - -$$\\begin{align}\\begin{split}H_C\&=\\sum_{i,j=1}^nQ_{ij}\\frac{\\mathbb{1}-Z_i}{2}\\frac{\\mathbb{1}-Z_j}{2}\\\\&=\\sum_{i,j=1}^n\\frac{Q_{ij}}{4}+\\sum_{i,j}^n\\frac{Q_{ij}}{4}Z_iZ_j\\\\&\\space\\space\\space\\space-\\sum_{i=1}^n\\bigg(\\sum_{j=1}^n\\frac{Q_{ij}}{4}\\bigg)Z_i-\\sum_{j=1}^n\\bigg(\\sum_{i=1}^n\\frac{Q_{ij}}{4}\\bigg)Z_j\\\\&\\space\\space\\space\\space\\end{split}\\end{align}$$ - -Swapping the indices $i$ and $j$ in the last sum, and using the identity $Z_iZ_i=\mathbb{1}$, we get - -$$\\begin{align}\\begin{split}H_C&=\\frac{1}{4}\\sum_{i\\neq j}Q_{ij}Z_iZ_j-\\frac{1}{4}\\sum_{i=1}^n\\sum_{j=1}^n(Q_{ij}+Q_{ji})Z_i&+\\frac{1}{4}\\sum_{i,j=1}^nQ_{ij}+\\frac{1}{4}\\sum_{i=1}^nQ_{ii}\\end{split}\\end{align}$$ - -Note that for each single $Z_i$ we sum the $i$-th row and the $i$-th column of the matrix $Q$. - - -For the cost operator $U_C$, which we feed into ``QAOAProblem``, we get - -$$\\begin{align}\\begin{split}U_C=e^{-i\\gamma H_C}=\&\\prod_{i,j=1}^nR_{Z_iZ_j}\\Bigg(\\frac{\\gamma}{2}Q_{ij}\\Bigg)\\times\\prod_{i=1}^nR_{Z_i}\\Bigg(-\\frac{\\gamma}{2}\\bigg(\\sum_{j=1}^n(Q_{ij}+Q_{ji})\\bigg)\\Bigg)\\\\&\\times\\exp\\Bigg(-\\frac{i\\gamma}{4}\\sum_{i,j=1}^nQ_{ij}-\\frac{i\\gamma}{4}\\sum_{i=1}^nQ_{ii}\\Bigg)\\end{split}\\end{align}$$ - -Here, the last factor correspods to a global phase. - -Let's translate this into a function: -:: - def create_QUBO_cost_operator(Q): - - def QUBO_cost_operator(qv, gamma): - - gphase(-gamma/4*(np.sum(Q)+np.trace(Q)),qv[0]) - for i in range(len(Q)): - rz(-gamma/2*(sum(Q[i])+sum(Q[:,i])), qv[i]) - for j in range(len(Q)): - if i != j and Q[i][j] != 0: - rzz(gamma/2*Q[i][j], qv[i], qv[j]) - return QUBO_cost_operator - -Like we did for :ref:`MaxCut <MaxCutQAOA>` and :ref:`M$\\kappa$CS <MkCSQAOA>` we also define the general QUBO objective function, the classical cost function, as well as construct the ``QUBOProblem`` blueprint bringing everything together. -:: - from qrisp import rzz, rz, gphase - import numpy as np - - def QUBO_obj(bitstring, Q): - x = np.array(list(bitstring), dtype=int) - cost = x.T @ Q @ x - return cost - - def create_QUBO_cl_cost_function(Q): - - def cl_cost_function(counts): - - def QUBO_obj(bitstring, Q): - x = np.array(list(bitstring), dtype=int) - cost = x.T @ Q @ x - return cost - - energy = 0 - for meas, meas_count in counts.items(): - obj_for_meas = QUBO_obj(meas,Q) - energy += obj_for_meas * meas_count - return energy - - return cl_cost_function - - def QUBO_problem(Q,init_type='random'): - - from qrisp.qaoa import QAOAProblem, RX_mixer - - return QAOAProblem(create_QUBO_cost_operator(Q), RX_mixer, create_QUBO_cl_cost_function(Q),init_type=init_type) - -That's it for the necessary ingredients you learned about in the :ref:`QAOA theory 101 section <QAOA101>`! Let's solve the set partitioning problem from above using this newly acquired information, and combine with how we already ran the QAOA algorithm using the :meth:`run <qrisp.qaoa.QAOAProblem.run>` method: - -- define the QUBO matrix $Q$, -- define the quantum argument ``qarg`` as a :ref:`QuantumArray <QuantumArray>` of :ref:`QuantumVariables <QuantumVariable>`, -- create the QUBO instance using ``QUBO_problem`` we defined above, -- run the algorithm using the :meth:`run <qrisp.qaoa.QAOAProblem.run>` method, and last but not least, -- examine the QAOA solutions with the highest probabilities for classical post processing: compute the cost functions, sort the solutions by their cost in ascending order, and print the solutions with their costs. - -These are exactly the pieces in the mosaic of code that ``solve_QUBO`` consists of and performs: -:: - from qrisp.default_backend import def_backend - from qrisp import QuantumVariable, QuantumArray - from operator import itemgetter - - Q = np.array( - [ - [-17, 20, 20, 20, 0, 40], - [ 0, -18, 20, 20, 20, 40], - [ 0, 0, -29, 20, 40, 40], - [ 0, 0, 0, -19, 20, 20], - [ 0, 0, 0, 0, -17, 20], - [ 0, 0, 0, 0, 0, -28], - ] - ) - - qarg = QuantumArray(qtype = QuantumVariable(1), shape = len(Q)) - - QUBO_instance = QUBO_problem(Q) - - depth = 1 - res = QUBO_instance.run(qarg, depth, mes_kwargs={"backend" : def_backend}, max_iter = 50) - - n_solutions = 10 - res = dict(list(res.items())[:n_solutions]) - - costs_and_solutions = [(QUBO_obj(bitstring, Q), bitstring) for bitstring in res.keys()] - - sorted_costs_and_solutions = sorted(costs_and_solutions, key=itemgetter(0))# - - for i in range(n_solutions): - print(f"Solution {i+1}: {sorted_costs_and_solutions[i][1]} with cost: {sorted_costs_and_solutions[i][0]}") - - -Now you are prepared to solve all QUBOs you derive and want to solve. On the other hand, if you would just like to play around instead, try out some QUBOs from this `list of QUBO formulations <https://blog.xa0.de/post/List-of-QUBO-formulations>`_. - - - diff --git a/documentation/source/reference/Examples/Shor.rst b/documentation/source/reference/Examples/Shor.rst deleted file mode 100644 index 061f9ac2..00000000 --- a/documentation/source/reference/Examples/Shor.rst +++ /dev/null @@ -1,110 +0,0 @@ -.. _ShorExample: - -Exploring Shor's algorithm -========================== - -Shor’s algorithm, named after mathematician Peter Shor, is a quantum algorithm for integer factorization. It’s a way to find the prime factors of a composite number. This algorithm is particularly famous because it is expected solve this problem exponentially faster than the most efficient known classical factoring algorithm, providing a good enough quantum computer. - -This, as you shall see in the example below, has significant implications for cryptography, as many encryption systems rely on the difficulty of factoring really large numbers. So far only numbers up to 15 were crackable utiliziing an implementation of Shor's algorithm - this is where Qrisp comes in. We were able to utilize it's high level architecture and take advantage of it's features to significantly increase this. If provided with enough quality qubits this is the implementation to use to break encryption. 🔥🏦🔥🚒 - -Factorizing numbers -------------------- - -Let's start with a simple example of using Shor's algorithm to factor a number. As stated above, up to our implementation, the highest factorized number was 15. So let's factor number 65! -:: - from qrisp.shor import shors_alg - shors_alg(65) - - -Try running the code on the website yourself and feel free to try the algorithm out factorizing different numbers! The result we obtain is 5, which *checks notes* is indeed one of the factors! - -As we will see in the next example, number 65 is easy to crack in terms of the private and public key pairings, which is used for encryption. However, the bacis principles of encryption remain the same even with using much greater numbers. - -A tale of encryption and decryption ------------------------------------ - -Imagine a scenario where two characters, Alice and Bob, are trying to exchange a secure message. They decide to use RSA encryption, a popular method that uses the product of two prime numbers as a key. In this case, they choose 5 and 13 as their private keys, and 7 as one of the public keys. -:: - from qrisp.shor import rsa_encrypt_string - rsa_encrypt_string(p = 5, q = 13, e = 7, message = "Qrisp is awesome!") - -Enter our detective, let's call him Gadget, who manages to intercept the encrypted message using his highly advanced encrypted-message-interceptor tool. He knows that Alice and Bob have used RSA encryption, but he doesn’t know the private keys they used. "Aaaargh, so close!", he thought. - -Luckily for the future of his career as a detective, he remembered that he has recently stumbled upon the website of Eclipse Qrisp where he read the :ref:`enlightening tutorial about Shor's algorithm <shor_tutorial>`. Albeit thinking the text in the tutorial is bordering science fiction, he still decided to give the implementation a go. - -His console read: -:: - intercepted_message = '01010000000101001010001100100110010010000101000010001101000010100011010101110011101000100100011100000100000100110111101000011000111110111111' - - from qrisp.shor import rsa_decrypt_string - rsa_decrypt_string(e = 7, N = 65, ciphertext = intercepted_message) - -He ran the command and simply smirked at the result and said "You've got that right, Alice and Bob... Well played!". - -*fin* - -New adder, no problem ---------------------- - -Stories like the one above are fun and exciting way to showcase the elegant approach of utilizing Eclipse Qrisp's high level structure. Learning from existing frameworks, however, it is also of utmost importance to ask ourselves the serious, hard hitting question of how to futureproof such an implementation. You've asked the question, we've got the answer - let's look under the hood and delve into the nitty-gritty! - -As elaborated on in the :ref:`Fault-Tolerant compilation tutorial <ft_compilation_shor>`, the Qrisp implementation of Shor's algorithm allows you to provide an arbitrary adder for the execution of the required arithmetic. With our Qrispy structure one can write ones own adder, or implement a shiny new one future research publications might bring, and test its performance claims. - -As of right now, the following list of adders have been pre-implemented: - -* The :meth:`fourier_adder <qrisp.fourier_adder>` (`paper <https://arxiv.org/abs/quant-ph/0008033>`_) requires minimal qubit overhead and has a very efficient :meth:`custom_control <qrisp.custom_control>` but uses a lot of parametized phase gates, which increases the T-depth. The low qubit count makes it suitable for simulation, which is why it is the default adder. - -* The :meth:`cucarro_adder <qrisp.cuccaro_adder>` (`paper <https://arxiv.org/abs/quant-ph/0410184>`_) also requires minimal qubits but no parametrized phase gates. It doesn't have a custom controlled version. - -* The :meth:`gidney_adder <qrisp.gidney_adder>` (`paper <https://arxiv.org/abs/1709.06648>`_) requires $n$ ancillae but uses the ``gidney`` Toffoli method described above, making it very fast in terms of T-depth but also economical in terms of T-count. - -* The :meth:`qcla <qrisp.qcla>` (`paper <https://arxiv.org/abs/2304.02921>`_) requires quite a lot of ancillae but has only logarithmic scaling when it comes to T-depth. It is faster than the Gidney adder for any input size larger than 7. - -Using a diffent adder is as easy as adding an ``inpl_adder`` keyword to the :ref:`QuantumModulus <QuantumModulus>` variable. Literally! - -Let's provide an example of benchmarking the :meth:`gidney_adder <qrisp.gidney_adder>` and compare it to the :meth:`qcla <qrisp.qcla>` on the operation most relevant for Shor's algorithm: Controlled modular in-place multiplication. - -:: - - from qrisp import * - N = 3295 - qg = QuantumModulus(N, inpl_adder = gidney_adder) - - ctrl_qbl = QuantumBool() - - with control(ctrl_qbl): - qg *= 953 - - gate_speed = lambda op : t_depth_indicator(op, epsilon = 2**-10) - - qc = qg.qs.compile(gate_speed = gate_speed, compile_mcm = True) - print(qc.t_depth()) - # Yields 956 - print(qc.num_qubits()) - # Yields 79 - - -Now the :meth:`qcla <qrisp.qcla>`: - -:: - - qg = QuantumModulus(N, inpl_adder = qcla) - - ctrl_qbl = QuantumBool() - - with control(ctrl_qbl): - qg *= 10 - - qc = qg.qs.compile(workspace = 10, gate_speed = gate_speed, compile_mcm = True) - - print(qc.t_depth())s - # Yields 784 - print(qc.num_qubits()) - # Yields 88 - -We see that the T-depth is reduced by $\approx 20 \%$. Due to the logarithmic scaling of the adder, larger scales will profit even more! Note that we granted the compiler 10 qubits of :ref:`workspace <workspace>`, as this adder can profit a lot from this resource. - -The comparison analysis is intriguing on its own, but here we wanted to emphasize the simplicity of improving the performance of Shor's algorithm by the means of implementing possible new shiny adders with the least amount of headaches. Future 👏🏻 proven 👏🏻 - - - diff --git a/src/qrisp/qaoa/problems/QUBO.py b/src/qrisp/qaoa/problems/QUBO.py index 5a1d5651..08b89e1e 100644 --- a/src/qrisp/qaoa/problems/QUBO.py +++ b/src/qrisp/qaoa/problems/QUBO.py @@ -77,15 +77,7 @@ def create_QUBO_cost_operator(Q): cost_operator function. """ - #def QUBO_cost_operator(qv, gamma): - - # for i in range(len(Q)): - # rz(-0.5*2*gamma*(0.5*Q[i][i]+0.5*sum(Q[i])), qv[i]) - # for j in range(i+1, len(qv)): - # if Q[i][j] !=0: - # rzz(0.25*2*gamma*Q[i][j], qv[i], qv[j]) - #return QUBO_cost_operator - #new try + def QUBO_cost_operator(qv, gamma): gphase(-gamma/4*(np.sum(Q)+np.trace(Q)),qv[0]) @@ -141,8 +133,8 @@ def solve_QUBO(Q, depth, backend = None, n_solutions = 1, print_res = True): Returns ------- - None - The function prints the runtime of the QAOA algorithm and the ``n_solutions`` best solutions with their respective costs. + optimal_solution: tuple + The function returns the optimal solution as a tuple where the first element is the cost and the second element is the optimal bitstring. If print_res is set to True, the function prints the runtime of the QAOA algorithm and the ``n_solutions`` best solutions with their respective costs. """ @@ -168,7 +160,12 @@ def solve_QUBO(Q, depth, backend = None, n_solutions = 1, print_res = True): # Sort the solutions by their cost in ascending order sorted_costs_and_solutions = sorted(costs_and_solutions, key=itemgetter(0)) + optimal_solution = sorted_costs_and_solutions[0] + if print_res is True: # Get the top solutions and print them for i in range(n_solutions): - print(f"Solution {i+1}: {sorted_costs_and_solutions[i][1]} with cost: {sorted_costs_and_solutions[i][0]}") \ No newline at end of file + print(f"Solution {i+1}: {sorted_costs_and_solutions[i][1]} with cost: {sorted_costs_and_solutions[i][0]}") + + return optimal_solution + From f289829e2b98923496e2106ca9c4e507d65f91a4 Mon Sep 17 00:00:00 2001 From: mat70593 <matic.petric@fokus.fraunhofer.de> Date: Tue, 6 Feb 2024 11:36:16 +0100 Subject: [PATCH 2/3] Re-added the MarCut, QUBO, and Shor example .rst files --- .../MaxCut.rst | 339 ++++++++++++++++++ .../QUBO.rst | 170 +++++++++ .../source/reference/Examples/Shor.rst | 110 ++++++ 3 files changed, 619 insertions(+) create mode 100644 documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst create mode 100644 documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst create mode 100644 documentation/source/reference/Examples/Shor.rst diff --git a/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst b/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst new file mode 100644 index 00000000..cce390f4 --- /dev/null +++ b/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/MaxCut.rst @@ -0,0 +1,339 @@ +.. _MaxCutQAOA: + +MaxCut QAOA Implementation +========================== + +Welcome to our Quantum Approximate Optimization Algorithm (QAOA) implementation tutorial for the MaxCut problem! This tutorial is designed to provide you with a hands-on understanding of how to apply QAOA using Qrisp to solve the MaxCut optimization problem. + +To refresh your memory, `Maximum Cut, or MaxCut for short <https://en.wikipedia.org/wiki/Maximum_cut>`_ involves dividing the nodes of a graph into two groups such that the number of edges between the groups is maximized. It has applications in various fields, including computer science, operations research, and statistical physics. + +In this interactive tutorial, we’ll first tackle the MaxCut problem to build a solid QAOA foundation. You'll be guided to learn how to formulate the problem for QAOA, implement the actual algorithm, and interpret the results. + +Whether you’re a quantum computing novice or an experienced algorithm developer checking out this remarkable new framework, you may have observed that virtually every other quantum programming framework and their quantum grandmother has implemented QAOA for MaxCut in a similar manner (authors used irony, it was super effective). We'll go one step further and showcase on why using the Qrisp QAOA module and its :ref:`QAOAProblem class <QAOAProblem>` **is** a no-brainer when it comes to simplicity of implementation as well as the speed of the optimization itself. + +Your very first QAOA layer +-------------------------- + +Let's start slowly with deepening the understanding of QAOA at its core - one layer of the algorithm and see how it looks like. We start with importing necessary quantum gates, as well as `networkx`, which is a package for creating and manipulating graphs. In this block of code you also define the angle parameters $\gamma$ and $\beta$: +:: + from qrisp import QuantumVariable, h, barrier, rz, rx, cx + import networkx as nx + from sympy import Symbol + + G = nx.Graph() + G.add_edges_from([[0,3],[0,4],[1,3],[1,4],[2,3],[2,4]]) + #nx.draw(G, with_labels=True) + N = G.number_of_nodes() + +Feel free to play around with this block of code changing the graph through the number of nodes, the edges connecting them, and print the result by uncommenting ``nx.draw(G, with_labels=True)``. + +Next thing on our list is defining the MaxCut cost operator and the mixer. Taking a peek in the appendices section of Hadfield's paper, we see, the MaxCut phase separator we have to implement is defined as +$$U_P=e^{-i\\gamma_pH_P}=\e^{-i\\gamma_p\\sum_{ij\\in E}(-Z_iZ_j)}.$$ Instead of using the `ZZ` gate, we will perform the same operation using only elementary gates. +:: + gamma = Symbol('γ') + + def apply_cost_operator(qv, gamma): + for pair in list(G.edges()): + cx(qv[pair[0]], qv[pair[1]]) + rz(2 * gamma, qv[pair[1]]) + cx(qv[pair[0]], qv[pair[1]]) + barrier(qv) + return qv + +The first line starts a loop over all edges in the graph ``G``. Each edge is a pair of nodes, which we refer to as ``pair[0]`` and ``pair[1]``. We then apply a controlled-X (CNOT) gate to the pair of qubits before applying a rotation around the Z-axis by an angle of $2\gamma$ to the second qubit in the pair. At the end we apply another CNOT gate to the pair of qubits. + +As suggested in the table of our QAOA theoretical overview, we are using the ``rx`` gate as the mixer. +:: + beta = Symbol('β') + + def apply_mixer(qv, beta): + for i in range(0, N): + rx(2 * beta, qv[i]) + barrier(qv) + return qv + +With these two functions at our disposal we can now implement our first QAOA layer and `print(qv.qs)` it. +:: + qv_1 = QuantumVariable(N) + h(qv_1) + apply_cost_operator(qv_1, gamma) + apply_mixer(qv_1, beta) + print(qv_1.qs) + +Applying hadamard gates on ``qv`` prepares the initial state for our system - the superposition state $|s\rangle$ . + +If you have executed all the preceding cells correctly, you should have been able to print and examine the circuit of the initial QAOA layer. You may notice a pattern in the phase separator section of the circuit, which links two qubits corresponding to a specific pair of nodes in `G` that are connected by an edge (commencing flashback: `G.add_edges_from([[0,3],[0,4],[1,3],[1,4],[2,3],[2,4]])`). The use of barriers between such pairs facilitates a more intuitive understanding of the phase separator unitary's inner workings. + +We’ve used the parameters $\gamma$ and $\beta$ symbolically, without assigning specific numerical values to them. In the next section we will build on top of the first layer, assign values to the angle parameters and iteratively optimize their value with each additional layer (SPOILER ALERT). + +As we add more layers to our QAOA circuit, the values of $\gamma$ and $\beta$ are refined, leading us closer to the optimal solution. This iterative improvement is the core of what makes QAOA such a powerful tool. + +MaxCut QAOA implementation +-------------------------- + +Putting focuss on the problem at hand (MaxCut if you're just scrolling through and not paying that much attention) it is reasonable to ask ourselves what will be the objective that we want to optimize. + +As the name suggests, we are looking for a cut going through as much edges as possible. It's therefore crucial to count the amount of cuts we cut through (I'll cut it out now) by checking if the nodes ``i`` and ``j`` belong to different groups. If they are, the edge $(i,j)$ is cut and added to the total. +:: + def maxcut_obj(x): + cut = 0 + for i, j in G.edges(): + if x[i] != x[j]: + cut -= 1 + return cut + +This optimization objective is important for the last building blocks described in the theoretical overview we haven't mentioned yet - the cost function. The cost function is important for keeping track of the energy of the system: +:: + def maxcut_cost_funct(meas_res): + energy = 0 + for meas, p in meas_res.items(): + obj_for_meas = maxcut_obj(meas) + energy += obj_for_meas * p + return energy + +We loop over the measurement in ``meas_res``. Each solution is represented by a string of 1 and 0 in ``meas`` with probability ``p`` keeping score in how likely a particular solution appears. + +For each such solution the cost is calculated using the ``maxcut_obj`` returning the number of cut edges. At the end, ``maxcut_cost_funct`` calculates the average “quality” of a set of solutions to the MaxCut problem. A lower average energy means that, on average, the solutions are better - they cut more edges. + +This is now nearly all the building blocks we need in order to run QAOA and start optimizing. Well, after we add more layers to our algorithm, that is. +:: + p = 3 + + def apply_p_layers(qv, beta, gamma): + assert(len(beta) == len(gamma)) + p = len(beta) + h(qv) + for i in range(p): + apply_cost_operator(qv, gamma[i]) + apply_mixer(qv, beta[i]) + barrier(qv) + return qv + +With our $p$-layered algorithm in hand there is one last piece of the puzzle: the objective function, which we need to generate in order to calculate the average energy of the resulting solutions. Oh, and it features the angle parameters $\gamma$ and $\beta$, which we combine in one list ``theta``. The solutions we get depend on the values of the angle parameters. +:: + + def quantum_objective(theta): + qv_p = QuantumVariable(N) + beta = theta[:p] + gamma = theta[p:] + qv = apply_p_layers(qv_p,beta, gamma) + results = qv.get_measurement() + return maxcut_cost_funct(results) + + +We can finally finish combine all the pieces of the mosaic and have it appraised by a professional - a classical optimizer. Running it online might take a second so we would like to encourage you to download Qrisp and try out for yourself! + +We first reset the initial state and randomize it in a list where the first $p$ values correspond to $\beta$ and the second one to $\gamma$. It has finally come the time to inclute the COBYLA optimiser, which will return the optimal parameters that minimze our objective function. This is done using the ``minimize`` method, which adjusts the parameters iteratively until it finds the best ones. + +Then, once we have our optimal parameters, we apply QAOA one more time to get our final set of solutions (``counts``), before finding the best solution by looking for the one with the lowest energy. Finally, we visualize our solution by coloring the nodes of our graph according to which group they belong to in the best solution. +:: + import numpy as np + from scipy.optimize import minimize + from operator import itemgetter + + init_point = np.pi * np.random.rand(2 * p) + + res_sample = minimize(quantum_objective, init_point, method='COBYLA', options={'maxiter':50}) + + optimal_theta = res_sample['x'] + qv_p = QuantumVariable(N) + qv = apply_p_layers(qv_p, optimal_theta[:p], optimal_theta[p:]) + counts = qv_p.get_measurement() + + best_cut, best_solution = min([(maxcut_obj(x),x) for x in counts.keys()], key=itemgetter(0)) + print(f"Best string: {best_solution} with cut: {-best_cut}") + + colors = ['r' if best_solution[node] == '0' else 'b' for node in G] + nx.draw(G,node_color = colors, pos=nx.bipartite_layout(G, [0,1,2])) + +And voila! We just solved the MaxCut problem using a relatively straightforward brute force approach. While this method can be effective, it’s not the most efficient or elegant way to tackle this problem, and other problem instances with more complex objectives, mixers, or both. + +That’s where our QAOAProblem class comes in. This class simplifies the process of running QAOA by introducing modularity, without sacrificing the generality of whichever problem instance we're dealing with. As we move forward in this tutorial, you’ll have the opportunity to try it out for yourself. + +QAOAProblem. Enough said. +------------------------- + +🎶 Enough chit chat, you ain't got all day - let's get to it: QAOAProblem, lead the way! 🎶 + +Taking the essential building blocks from `QAOAnsatz <https://arxiv.org/abs/1709.03489>`_ into account, we built the :ref:`QAOAProblem class <QAOA>` with modularity in mind. We gathered and improved the functions shown in the example above and created a powerfull architecture with which it's easy to implement various problem instances of varying complexity.In this tutorial we focus on using QAOAProblem, with extensive documentation being available :ref:`here <QAOA>`. + +We start with renaming our quantum argument `qv` to a more general `qarg` because more often than not we'll combine QuantumVariables into a QuantumArray to make implementations of other problem instances more efficient. +:: + qarg = QuantumVariable(len(G)) + + depth = 3 + +BEHOLD, THE POWER OF QRISP! +:: + from qrisp.qaoa import QAOAProblem, RX_mixer + + maxcut_instance = QAOAProblem(apply_cost_operator, RX_mixer, maxcut_cost_funct) + + res = maxcut_instance.run(qarg, depth, max_iter = 50) # runs the simulation + +And that's pretty much it, really. Apart from visualizing the results again. +:: + best_cut, best_solution = min([(maxcut_obj(x),x) for x in res.keys()], key=itemgetter(0)) + print(f"Best string: {best_solution} with cut: {-best_cut}") + + res_str = list(res.keys())[0] + print("QAOA solution: ", res_str) + best_cut, best_solution = (maxcut_obj(res_str),res_str) + + colors = ['r' if best_solution[node] == '0' else 'b' for node in G] + nx.draw(G,node_color = colors) + +After thorough comparison certainly made after running both approaches for different graph topologies you are in position to be the judge regarding which approach is qrispier. + +.. _benchmark_maxcut: +Benchmarking the performance +---------------------------- +Lucky for you, the tutorial is not over since there is one more important functionality we would like to show you - the :meth:`benchmark <qrisp.qaoa.QAOAProblem.benchmark>` method of our QAOA module. This method allows you to observe the performance of your simulation: + +- get insights about the approximation ratio for each run together with the average approximation ratio with its variance, +- obtain total circuit depth, +- obtain the overall quantum volume, +- and rank the results for different depths, amount of shots, as well as iterations. + +Let's show how easy benchmarking QAOA is using Qrisp: +:: + print('RND') + + benchmark_data = maxcut_instance.benchmark(qarg = QuantumVariable(len(G)), + depth_range = [1,2,3,4,5], + shot_range = [1000, 10000], + iter_range = [100, 200], + optimal_solution = '00011', + repetitions = 1, + init_type='random' + ) + + temp = benchmark_data.rank(print_res = True) + + _,rndFO=benchmark_data.evaluate() + + print('Approximation ratio: ',sum(rndFO)/len(rndFO)) + print("Variance: ",np.var(rndFO)) + +In the above example we obtain the list of results for the approximation ratio for various depths of QAOA ranging from 1 to 5. We run the benchmark for both 1000, and 10000 shots (the latter number is to avoid sampling of the solution); and for 100 and 200 iterations of the optimizer. Since QAOA is probabilistic every run is unique so it's important to run more repetitions for clearer insights. + +It is important to note that in order to do the benchmark, one has to already know the optimal solution in order to calculate the optimal energy, and through that, the approximation ratio. This is an example of the output for the benchmark above: +:: + RND + Rank approx_ratio Overall QV p QC depth QB count Shots Iterations + ============================================================================ + 0 0.990e+0 4.4000e+8 5 44 5 10000 200 + 1 0.977e+0 4.4000e+7 5 44 5 1000 200 + 2 0.974e+0 3.6000e+8 4 36 5 10000 200 + 3 0.956e+0 1.8000e+8 4 36 5 10000 100 + 4 0.945e+0 2.2000e+7 5 44 5 1000 100 + 5 0.940e+0 2.2000e+8 5 44 5 10000 100 + 6 0.940e+0 3.6000e+7 4 36 5 1000 200 + 7 0.871e+0 2.0000e+8 2 20 5 10000 200 + 8 0.830e+0 2.8000e+8 3 28 5 10000 200 + 9 0.817e+0 1.4000e+7 3 28 5 1000 100 + 10 0.815e+0 1.4000e+8 3 28 5 10000 100 + 11 0.775e+0 2.8000e+7 3 28 5 1000 200 + 12 0.774e+0 1.8000e+7 4 36 5 1000 100 + 13 0.765e+0 2.0000e+7 2 20 5 1000 200 + 14 0.764e+0 1.0000e+7 2 20 5 1000 100 + 15 0.736e+0 1.0000e+8 2 20 5 10000 100 + 16 0.718e+0 1.2000e+8 1 12 5 10000 200 + 17 0.705e+0 1.1000e+7 1 11 5 1000 200 + 18 0.556e+0 6.0000e+6 1 12 5 1000 100 + 19 0.555e+0 5.5000e+7 1 11 5 10000 100 + Approximation ratio: 0.8201699999999998 + Variance: 0.016184990766666633 + +As expected, on average the runs with higher depths yield better approximation ratio (with some outliers, of course). + +The ``print('RND')`` was used because since the 0.4 update we have also included the `TQA warm-start initialization <https://arxiv.org/abs/2101.05742>`_, which can be used within the :meth:`benchmark <qrisp.qaoa.QAOAProblem.benchmark>` method by setting ``init_type='tqa'``. Let's try benchmarking this approach as well: +:: + print('TQA') + maxcut_instance = maxcut_problem(G) + + benchmark_data = maxcut_instance.benchmark(qarg = QuantumVariable(len(G)), + depth_range = [1,2,3,4,5], + shot_range = [1000, 10000], + iter_range = [100, 200], + optimal_solution = '00011', + repetitions = 1, + init_type='tqa' + ) + + temp = benchmark_data.rank(print_res = True) + + _,tqaFO=benchmark_data.evaluate() + + print('Approximation ratio: ',sum(tqaFO)/len(tqaFO)) + print("Variance: ",np.var(tqaFO)) + +Again, since QAOA is probabilistic, every run returns different results. You can find our try below: +:: + TQA + Rank approx_ratio Overall QV p QC depth QB count Shots Iterations + ============================================================================ + 0 0.989e+0 4.4000e+8 5 44 5 10000 200 + 1 0.983e+0 2.2000e+8 5 44 5 10000 100 + 2 0.963e+0 1.8000e+8 4 36 5 10000 100 + 3 0.963e+0 3.6000e+8 4 36 5 10000 200 + 4 0.957e+0 4.4000e+7 5 44 5 1000 200 + 5 0.897e+0 1.4000e+8 3 28 5 10000 100 + 6 0.897e+0 2.8000e+8 3 28 5 10000 200 + 7 0.892e+0 2.8000e+7 3 28 5 1000 200 + 8 0.880e+0 1.0000e+8 2 20 5 10000 100 + 9 0.880e+0 2.0000e+8 2 20 5 10000 200 + 10 0.875e+0 1.4000e+7 3 28 5 1000 100 + 11 0.867e+0 1.8000e+7 4 36 5 1000 100 + 12 0.866e+0 1.0000e+7 2 20 5 1000 100 + 13 0.860e+0 2.0000e+7 2 20 5 1000 200 + 14 0.847e+0 2.2000e+7 5 44 5 1000 100 + 15 0.822e+0 3.6000e+7 4 36 5 1000 200 + 16 0.717e+0 6.0000e+7 1 12 5 10000 100 + 17 0.717e+0 1.2000e+8 1 12 5 10000 200 + 18 0.715e+0 5.5000e+6 1 11 5 1000 100 + 19 0.691e+0 1.1000e+7 1 11 5 1000 200 + Approximation ratio: 0.8639683333333332 + Variance: 0.007950404719444411 + +As we can see, the TQA initialization tends to return better approximation ratios more consistently. Feel free to play around with the :meth:`benchmark <qrisp.qaoa.QAOAProblem.benchmark>` method (or leave it running over night) to try and compare the two approaches further. + +Summary and motivation +---------------------- + +To shortly summarize, in order to implement QAOA with QAOAProblem one needs to specify the problem using the following recipe + +I. define **CLASSICAL COST FUNCTION** of the problem you want to implement: ``maxcut_cost_funct(counts)``, +II. define the **INITIAL STATE** if it is not the superposition, which is set by default, +III. define **COST OPERATOR aka PHASE SEPARATOR** (or use the ones specified in `From QAOA to QAOA <https://arxiv.org/abs/1709.03489>`_) like ``apply_cost_operator`` above, and +IV. select **MIXER** from the :ref:`assortment we provide and list here <MIXers>`. + +Let's condense all of the above, and implement QAOA for MaxCut one last time in one block of code +:: + from qrisp.qaoa import QuantumArray, QuantumVariable, QAOAProblem, maxcut_obj,create_maxcut_cl_cost_function,create_maxcut_cost_operator, RX_mixer + import networkx as nx + from operator import itemgetter + + G = nx.Graph() + G.add_edges_from([[0,3],[0,4],[1,3],[1,4],[2,3],[2,4]]) + + qarg = QuantumArray(qtype = QuantumVariable(1), shape = len(G)) + + depth = 5 + + maxcut_instance = QAOAProblem(create_maxcut_cost_operator(G), RX_mixer, create_maxcut_cl_cost_function(G)) + + res = maxcut_instance.run(qarg, depth, max_iter = 50) + + best_cut, best_solution = min([(maxcut_obj(x,G),x) for x in res.keys()], key=itemgetter(0)) + + res_str = list(res.keys())[0] + print("QAOA solution: ", res_str) + best_cut, best_solution = (maxcut_obj(res_str,G),res_str) + + colors = ['r' if best_solution[node] == '0' else 'b' for node in G] + nx.draw(G,node_color = colors, pos=nx.bipartite_layout(G, [0,1,2])) + +You've got to admit that this is pretty cool, clean and qrispy, right?! + +If you are still not convinced, we provide a more complex problem instance in the next tutorial while also showcasing and putting some unique functionalities of Qrisp to the test. Let's make this transition a little more dramatic by saying that it's time to put our money where out mouths are (this is funny, because Qrisp is open source) and put this recipe to the test tackling the problem class which generalizes MaxCut: :ref:`Max-$\\kappa$-Colorable Subgraph <MkCSQAOA>`. \ No newline at end of file diff --git a/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst b/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst new file mode 100644 index 00000000..4ff60392 --- /dev/null +++ b/documentation/source/general/tutorial/Quantum Alternating Operator Ansatz/QUBO.rst @@ -0,0 +1,170 @@ +.. _QUBOQAOA: + +QUBO as a QAOAProblem Problem Instance +====================================== + +In this tutorial you’ll be guided through the process of defining a new phase separator to be used within the scope of the :ref:`Alternating Operator Ansatz <AOA>` focussed on solving various QUBO problems with only needing the QUBO matrix $Q$ as an input. +QUBO, or `Quadratic Unconstrained Binary Optimization <https://en.wikipedia.org/wiki/Quadratic_unconstrained_binary_optimization>`_, is a type of problem that involves optimizing a quadratic function of binary variables. + +After first translating the QUBO Hamiltonian $H_C$ from the binary basis $x_i$ to the state basis with the translation $x_i\rightarrow \frac{\mathbb{I}-Z_i}{2}$, we'll construct the phase separator unitary and run QAOA to solve the set partitioning problem. + +Not to get too ahead of ourselves with how to do it, let's first show how to solve a QUBO problem using Qrisp. To obtain the optimal solutions to a problem, we only need to know the QUBO matrix $Q$. In order to get familiar with that, we propose having a look at the `Tutorial on Formulating and Using QUBO Models <https://arxiv.org/abs/1811.11538>`_, which explains how to derive such a matrix $Q$. + +One can then simply minimize the cost function $\min C(x)=\min x^TQx$ for binary variables $x_i\in\{0,1\}$. This is done either by sending the QUBO matrix $Q$ to the annealer/quantum computer, which calculates and returns the bitstrings with corresponding values of $C(x)$. The lower the $C(x)$, the better the solution. In other words, we want to find the Hamiltonian $H_C|x\rangle=C(x)|x\rangle$. + +Let's borrow the QUBO matrix (explained in depth in the above mentioned `tutorial <https://arxiv.org/abs/1811.11538>`_) for the set partitioning problem. The QUBO matrix in that case is + +$$Q = \\begin{pmatrix}-17&10&10&10&0&20\\\\10&-18&10&10&10&20\\\\10&10&-29&10&20&20\\\\10&10&10&-19&10&10\\\\0&10&20&10&-17&10\\\\20&20&20&10&10&-28\\end{pmatrix}$$ + +Usually, QUBO matrices are upper triangular (by convention) - this means that only elements above the diagonal (with the diagonal included) are not equal to zero. This is because the variables in QUBO problems are binary, which results in the identity $q_iq_j=q_jq_i$ for binary variables. Because of this, the elements above the diagonal are doubled compared to the symmetrized version of $Q$ we wrote down above. + +Of course it's easy to go from the conventional QUBO $Q_{\text{up}\Delta}$ formulation to the symmetrized $Q_\text{sym}$ QUBO: $Q_\text{sym}=\frac{1}{2}\big(Q_{\text{up}\Delta}+Q_{\text{up}\Delta}^T\big)$, and back: +:: + import numpy as np + Q_up_triang = 2 * np.triu(Q_sym) - np.diag(np.diag(Q_sym)) + + +Our implementation of solving QUBO problems using QAOA works for both the upper triangular, as well as symmetrized matrix conventions. As promised, we can now immediately solve the QUBO by calling the ``solve_QUBO`` function: +:: + from qrisp.qaoa import * + import numpy as np + + Q = np.array( + [ + [-17, 10, 10, 10, 0, 20], + [ 10, -18, 10, 10, 10, 20], + [ 10, 10, -29, 10, 20, 20], + [ 10, 10, 10, -19, 10, 10], + [ 0, 10, 20, 10, -17, 10], + [ 20, 20, 20, 10, 10, -28], + ] + ) + + solve_QUBO(Q, depth = 1, backend = def_backend, n_solutions=10) + +That's it! You can try running this block of code on the website using the Thebe interactivity integration, or run it on your own Qrispy environment. +We see that we obtain the 10 best solutions with corresponding bitsring of the binary variables. +You can test the statements above, and convert the symmetrized matrix into the upper triangular one, and run ``solve_QUBO`` again. + +We see, that we only needed to provide the matrix $Q$, specify the depth of the QAOA circuit that is running in the background, and the backend you want to run the algorithm on. Clean and qrispy! + +From QUBO to QAOA +----------------- + +Of course it's beneficial to not only know the "how", but also understand the "why". So let's dig in! + +To construct such a circuit with quantum gates and run it on a quantum computer, one has to translate between the basis of $x_i$ to Pauli gate operators. + +$$Z_i|x\\rangle=(-1)^{x_i}\\rightarrow\\frac{\\mathbb{1}-Z_i}{2}|x\\rangle=x_i$$ + +To find $H_C$ one can calculate + +$$\\begin{align}\\begin{split}H_C\&=\\sum_{i,j=1}^nQ_{ij}\\frac{\\mathbb{1}-Z_i}{2}\\frac{\\mathbb{1}-Z_j}{2}\\\\&=\\sum_{i,j=1}^n\\frac{Q_{ij}}{4}+\\sum_{i,j}^n\\frac{Q_{ij}}{4}Z_iZ_j\\\\&\\space\\space\\space\\space-\\sum_{i=1}^n\\bigg(\\sum_{j=1}^n\\frac{Q_{ij}}{4}\\bigg)Z_i-\\sum_{j=1}^n\\bigg(\\sum_{i=1}^n\\frac{Q_{ij}}{4}\\bigg)Z_j\\\\&\\space\\space\\space\\space\\end{split}\\end{align}$$ + +Swapping the indices $i$ and $j$ in the last sum, and using the identity $Z_iZ_i=\mathbb{1}$, we get + +$$\\begin{align}\\begin{split}H_C&=\\frac{1}{4}\\sum_{i\\neq j}Q_{ij}Z_iZ_j-\\frac{1}{4}\\sum_{i=1}^n\\sum_{j=1}^n(Q_{ij}+Q_{ji})Z_i&+\\frac{1}{4}\\sum_{i,j=1}^nQ_{ij}+\\frac{1}{4}\\sum_{i=1}^nQ_{ii}\\end{split}\\end{align}$$ + +Note that for each single $Z_i$ we sum the $i$-th row and the $i$-th column of the matrix $Q$. + + +For the cost operator $U_C$, which we feed into ``QAOAProblem``, we get + +$$\\begin{align}\\begin{split}U_C=e^{-i\\gamma H_C}=\&\\prod_{i,j=1}^nR_{Z_iZ_j}\\Bigg(\\frac{\\gamma}{2}Q_{ij}\\Bigg)\\times\\prod_{i=1}^nR_{Z_i}\\Bigg(-\\frac{\\gamma}{2}\\bigg(\\sum_{j=1}^n(Q_{ij}+Q_{ji})\\bigg)\\Bigg)\\\\&\\times\\exp\\Bigg(-\\frac{i\\gamma}{4}\\sum_{i,j=1}^nQ_{ij}-\\frac{i\\gamma}{4}\\sum_{i=1}^nQ_{ii}\\Bigg)\\end{split}\\end{align}$$ + +Here, the last factor correspods to a global phase. + +Let's translate this into a function: +:: + def create_QUBO_cost_operator(Q): + + def QUBO_cost_operator(qv, gamma): + + gphase(-gamma/4*(np.sum(Q)+np.trace(Q)),qv[0]) + for i in range(len(Q)): + rz(-gamma/2*(sum(Q[i])+sum(Q[:,i])), qv[i]) + for j in range(len(Q)): + if i != j and Q[i][j] != 0: + rzz(gamma/2*Q[i][j], qv[i], qv[j]) + return QUBO_cost_operator + +Like we did for :ref:`MaxCut <MaxCutQAOA>` and :ref:`M$\\kappa$CS <MkCSQAOA>` we also define the general QUBO objective function, the classical cost function, as well as construct the ``QUBOProblem`` blueprint bringing everything together. +:: + from qrisp import rzz, rz, gphase + import numpy as np + + def QUBO_obj(bitstring, Q): + x = np.array(list(bitstring), dtype=int) + cost = x.T @ Q @ x + return cost + + def create_QUBO_cl_cost_function(Q): + + def cl_cost_function(counts): + + def QUBO_obj(bitstring, Q): + x = np.array(list(bitstring), dtype=int) + cost = x.T @ Q @ x + return cost + + energy = 0 + for meas, meas_count in counts.items(): + obj_for_meas = QUBO_obj(meas,Q) + energy += obj_for_meas * meas_count + return energy + + return cl_cost_function + + def QUBO_problem(Q,init_type='random'): + + from qrisp.qaoa import QAOAProblem, RX_mixer + + return QAOAProblem(create_QUBO_cost_operator(Q), RX_mixer, create_QUBO_cl_cost_function(Q),init_type=init_type) + +That's it for the necessary ingredients you learned about in the :ref:`QAOA theory 101 section <QAOA101>`! Let's solve the set partitioning problem from above using this newly acquired information, and combine with how we already ran the QAOA algorithm using the :meth:`run <qrisp.qaoa.QAOAProblem.run>` method: + +- define the QUBO matrix $Q$, +- define the quantum argument ``qarg`` as a :ref:`QuantumArray <QuantumArray>` of :ref:`QuantumVariables <QuantumVariable>`, +- create the QUBO instance using ``QUBO_problem`` we defined above, +- run the algorithm using the :meth:`run <qrisp.qaoa.QAOAProblem.run>` method, and last but not least, +- examine the QAOA solutions with the highest probabilities for classical post processing: compute the cost functions, sort the solutions by their cost in ascending order, and print the solutions with their costs. + +These are exactly the pieces in the mosaic of code that ``solve_QUBO`` consists of and performs: +:: + from qrisp.default_backend import def_backend + from qrisp import QuantumVariable, QuantumArray + from operator import itemgetter + + Q = np.array( + [ + [-17, 20, 20, 20, 0, 40], + [ 0, -18, 20, 20, 20, 40], + [ 0, 0, -29, 20, 40, 40], + [ 0, 0, 0, -19, 20, 20], + [ 0, 0, 0, 0, -17, 20], + [ 0, 0, 0, 0, 0, -28], + ] + ) + + qarg = QuantumArray(qtype = QuantumVariable(1), shape = len(Q)) + + QUBO_instance = QUBO_problem(Q) + + depth = 1 + res = QUBO_instance.run(qarg, depth, mes_kwargs={"backend" : def_backend}, max_iter = 50) + + n_solutions = 10 + res = dict(list(res.items())[:n_solutions]) + + costs_and_solutions = [(QUBO_obj(bitstring, Q), bitstring) for bitstring in res.keys()] + + sorted_costs_and_solutions = sorted(costs_and_solutions, key=itemgetter(0))# + + for i in range(n_solutions): + print(f"Solution {i+1}: {sorted_costs_and_solutions[i][1]} with cost: {sorted_costs_and_solutions[i][0]}") + + +Now you are prepared to solve all QUBOs you derive and want to solve. On the other hand, if you would just like to play around instead, try out some QUBOs from this `list of QUBO formulations <https://blog.xa0.de/post/List-of-QUBO-formulations>`_. + + + diff --git a/documentation/source/reference/Examples/Shor.rst b/documentation/source/reference/Examples/Shor.rst new file mode 100644 index 00000000..061f9ac2 --- /dev/null +++ b/documentation/source/reference/Examples/Shor.rst @@ -0,0 +1,110 @@ +.. _ShorExample: + +Exploring Shor's algorithm +========================== + +Shor’s algorithm, named after mathematician Peter Shor, is a quantum algorithm for integer factorization. It’s a way to find the prime factors of a composite number. This algorithm is particularly famous because it is expected solve this problem exponentially faster than the most efficient known classical factoring algorithm, providing a good enough quantum computer. + +This, as you shall see in the example below, has significant implications for cryptography, as many encryption systems rely on the difficulty of factoring really large numbers. So far only numbers up to 15 were crackable utiliziing an implementation of Shor's algorithm - this is where Qrisp comes in. We were able to utilize it's high level architecture and take advantage of it's features to significantly increase this. If provided with enough quality qubits this is the implementation to use to break encryption. 🔥🏦🔥🚒 + +Factorizing numbers +------------------- + +Let's start with a simple example of using Shor's algorithm to factor a number. As stated above, up to our implementation, the highest factorized number was 15. So let's factor number 65! +:: + from qrisp.shor import shors_alg + shors_alg(65) + + +Try running the code on the website yourself and feel free to try the algorithm out factorizing different numbers! The result we obtain is 5, which *checks notes* is indeed one of the factors! + +As we will see in the next example, number 65 is easy to crack in terms of the private and public key pairings, which is used for encryption. However, the bacis principles of encryption remain the same even with using much greater numbers. + +A tale of encryption and decryption +----------------------------------- + +Imagine a scenario where two characters, Alice and Bob, are trying to exchange a secure message. They decide to use RSA encryption, a popular method that uses the product of two prime numbers as a key. In this case, they choose 5 and 13 as their private keys, and 7 as one of the public keys. +:: + from qrisp.shor import rsa_encrypt_string + rsa_encrypt_string(p = 5, q = 13, e = 7, message = "Qrisp is awesome!") + +Enter our detective, let's call him Gadget, who manages to intercept the encrypted message using his highly advanced encrypted-message-interceptor tool. He knows that Alice and Bob have used RSA encryption, but he doesn’t know the private keys they used. "Aaaargh, so close!", he thought. + +Luckily for the future of his career as a detective, he remembered that he has recently stumbled upon the website of Eclipse Qrisp where he read the :ref:`enlightening tutorial about Shor's algorithm <shor_tutorial>`. Albeit thinking the text in the tutorial is bordering science fiction, he still decided to give the implementation a go. + +His console read: +:: + intercepted_message = '01010000000101001010001100100110010010000101000010001101000010100011010101110011101000100100011100000100000100110111101000011000111110111111' + + from qrisp.shor import rsa_decrypt_string + rsa_decrypt_string(e = 7, N = 65, ciphertext = intercepted_message) + +He ran the command and simply smirked at the result and said "You've got that right, Alice and Bob... Well played!". + +*fin* + +New adder, no problem +--------------------- + +Stories like the one above are fun and exciting way to showcase the elegant approach of utilizing Eclipse Qrisp's high level structure. Learning from existing frameworks, however, it is also of utmost importance to ask ourselves the serious, hard hitting question of how to futureproof such an implementation. You've asked the question, we've got the answer - let's look under the hood and delve into the nitty-gritty! + +As elaborated on in the :ref:`Fault-Tolerant compilation tutorial <ft_compilation_shor>`, the Qrisp implementation of Shor's algorithm allows you to provide an arbitrary adder for the execution of the required arithmetic. With our Qrispy structure one can write ones own adder, or implement a shiny new one future research publications might bring, and test its performance claims. + +As of right now, the following list of adders have been pre-implemented: + +* The :meth:`fourier_adder <qrisp.fourier_adder>` (`paper <https://arxiv.org/abs/quant-ph/0008033>`_) requires minimal qubit overhead and has a very efficient :meth:`custom_control <qrisp.custom_control>` but uses a lot of parametized phase gates, which increases the T-depth. The low qubit count makes it suitable for simulation, which is why it is the default adder. + +* The :meth:`cucarro_adder <qrisp.cuccaro_adder>` (`paper <https://arxiv.org/abs/quant-ph/0410184>`_) also requires minimal qubits but no parametrized phase gates. It doesn't have a custom controlled version. + +* The :meth:`gidney_adder <qrisp.gidney_adder>` (`paper <https://arxiv.org/abs/1709.06648>`_) requires $n$ ancillae but uses the ``gidney`` Toffoli method described above, making it very fast in terms of T-depth but also economical in terms of T-count. + +* The :meth:`qcla <qrisp.qcla>` (`paper <https://arxiv.org/abs/2304.02921>`_) requires quite a lot of ancillae but has only logarithmic scaling when it comes to T-depth. It is faster than the Gidney adder for any input size larger than 7. + +Using a diffent adder is as easy as adding an ``inpl_adder`` keyword to the :ref:`QuantumModulus <QuantumModulus>` variable. Literally! + +Let's provide an example of benchmarking the :meth:`gidney_adder <qrisp.gidney_adder>` and compare it to the :meth:`qcla <qrisp.qcla>` on the operation most relevant for Shor's algorithm: Controlled modular in-place multiplication. + +:: + + from qrisp import * + N = 3295 + qg = QuantumModulus(N, inpl_adder = gidney_adder) + + ctrl_qbl = QuantumBool() + + with control(ctrl_qbl): + qg *= 953 + + gate_speed = lambda op : t_depth_indicator(op, epsilon = 2**-10) + + qc = qg.qs.compile(gate_speed = gate_speed, compile_mcm = True) + print(qc.t_depth()) + # Yields 956 + print(qc.num_qubits()) + # Yields 79 + + +Now the :meth:`qcla <qrisp.qcla>`: + +:: + + qg = QuantumModulus(N, inpl_adder = qcla) + + ctrl_qbl = QuantumBool() + + with control(ctrl_qbl): + qg *= 10 + + qc = qg.qs.compile(workspace = 10, gate_speed = gate_speed, compile_mcm = True) + + print(qc.t_depth())s + # Yields 784 + print(qc.num_qubits()) + # Yields 88 + +We see that the T-depth is reduced by $\approx 20 \%$. Due to the logarithmic scaling of the adder, larger scales will profit even more! Note that we granted the compiler 10 qubits of :ref:`workspace <workspace>`, as this adder can profit a lot from this resource. + +The comparison analysis is intriguing on its own, but here we wanted to emphasize the simplicity of improving the performance of Shor's algorithm by the means of implementing possible new shiny adders with the least amount of headaches. Future 👏🏻 proven 👏🏻 + + + From 2cc12756c28abfac1f88badddd637e31a6f3b0aa Mon Sep 17 00:00:00 2001 From: mat70593 <matic.petric@fokus.fraunhofer.de> Date: Tue, 6 Feb 2024 11:49:17 +0100 Subject: [PATCH 3/3] Added the max_iter=50 as a keyword in solve_QUBO for more customizability, also set print_res=False as default --- src/qrisp/qaoa/problems/QUBO.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/qrisp/qaoa/problems/QUBO.py b/src/qrisp/qaoa/problems/QUBO.py index 08b89e1e..27ac1584 100644 --- a/src/qrisp/qaoa/problems/QUBO.py +++ b/src/qrisp/qaoa/problems/QUBO.py @@ -111,7 +111,7 @@ def QUBO_problem(Q): return QAOAProblem(create_QUBO_cost_operator(Q), RX_mixer, create_QUBO_cl_cost_function(Q)) -def solve_QUBO(Q, depth, backend = None, n_solutions = 1, print_res = True): +def solve_QUBO(Q, depth, max_iter = 50, backend = None, n_solutions = 1, print_res = False): """ Solves a Quadratic Unconstrained Binary Optimization (QUBO) problem using the Quantum Approximate Optimization Algorithm (QAOA). The function imports the default backend from the 'qrisp.default_backend' module. @@ -126,6 +126,8 @@ def solve_QUBO(Q, depth, backend = None, n_solutions = 1, print_res = True): QUBO matrix to solve. depth : int The depth (amount of layers) of the QAOA circuit. + max_iter : int + The maximal amount of iterations of the COBYLA optimizer in the QAOA algorithm. backend : str The backend to be used for the quantum/annealing simulation. n_solutions : int @@ -151,7 +153,7 @@ def solve_QUBO(Q, depth, backend = None, n_solutions = 1, print_res = True): backend = backend # Run QAOA with given quantum arguments, depth, measurement keyword arguments and maximum iterations for optimization - res = QUBO_instance.run(qarg, depth, mes_kwargs={"backend" : backend}, max_iter = 50) # runs the simulation + res = QUBO_instance.run(qarg, depth, mes_kwargs={"backend" : backend}, max_iter = max_iter) # runs the simulation res = dict(list(res.items())[:n_solutions]) # Calculate the cost for each solution