From 32d35db9db0bfe6ce357995f599200f9be941cc2 Mon Sep 17 00:00:00 2001 From: Mark Fingerhuth Date: Tue, 27 Mar 2018 18:13:33 -0400 Subject: [PATCH 1/2] added XY mixer example to QAOA docs --- docs/qaoa.rst | 147 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 147 insertions(+) diff --git a/docs/qaoa.rst b/docs/qaoa.rst index 9d8d25f..c284486 100644 --- a/docs/qaoa.rst +++ b/docs/qaoa.rst @@ -408,6 +408,153 @@ As expected the approximate evolution becomes more accurate as the number of steps (\\(\\alpha\\)) is increased. For this simple model \\(\\alpha = 2\\) is sufficient to find the two degnerate cuts of the four ring graph. +XY mixer Hamiltonian Example +~~~~~~~~~~~~~~~~~~~~~~~ + +In this example we will use the important optional argument ``driver_operators`` to solve a simple toy problem +using the Ising QAOA wrapper. This example is loosely based on the work 'QAOA with hard and soft constraints' by Hadfield et. al (2017) +[`4 `_]. Here, the core idea is to implement hard constraints into the +driver/mixer Hamiltonian which ensures that the search takes place only in the feasible subspace. Just keep reading and it will get +clear what *hard constraints* and *feasible subspace* means. +Suppose we encode the directions right, left and up into 3 bits each: + +.. image:: ising_qaoa/3_directions.png + :align: center + :scale: 75% + +All of these three bit strings have one thing in common: a Hamming weight of one. (basically the number of +1's in the bit string [`5 `_]). Let's say we want to use our +Ising QAOA to find self-avoiding walks (SAW) with three moves. Due to the way we assigned bit strings to +directions we know that every *feasible* solution for a walk with three moves must have a Hamming weight +of three! Hence, our *hard constraint* is a constant Hamming weight of one for each triplet of bits. +Consequently, the bit string \\( 110 \\, 111 \\, 011 \\) is not a feasible solution since its Hamming weight is +too large and the three triplets carry no meaning in the context of our problem! Thus, we can define the +*feasible subspace* as the set of all *feasible* bit strings. +For example, a valid SAW with three moves (a feasible solution) could look like this: + +.. image:: ising_qaoa/valid_saw.png + :align: center + :scale: 75% + +The goal in this example is to create a driver/mixer Hamiltonian which preserves the Hamming weight of each bit triplet. +The \\( SWAP_{i,j} \\) operation is a great building block since it simply swaps the values of two bits which +of course preserves Hamming weight. In our case, we use the following simplified definition of +the \\( SWAP_{i,j} \\) operation (see [`4 `_]): + +$$ SWAP_{i,j} = \\frac{1}{2} ( X_{i}X_{j} + Y_{i}Y_{j} ) $$ + +The mixer should mix amplitudes between feasible solutions which means it should swap around the bits within each triplet. +This can be achieved with the following driver/mixer Hamiltonian: + +$$ \\mathbf{H}_{M} = \\sum_{k \\in \\{ 0,3,6\\}} \\sum_{i=k}^{k+1} \\sum_{i+1}^{k+2} SWAP_{i,j} $$ + +Let's first do our imports: + +.. code-block:: python + + from grove.alpha.arbitrary_state.arbitrary_state import create_arbitrary_state + import numpy as np + import pyquil.api as api + import itertools + from functools import reduce + import operator + +In pyQuil code we can generate this driver/mixer Hamiltonian like this: + +.. code-block:: python + + mixer_operators = [] + for k in [0,3,6]: + l = list(range(k,k+3)) + ij_pairs = [ (a,b) for a in l for b in l[l.index(a):] if a!=b and b-a<=2] + for i, j in ij_pairs: + mixer_operators.append(PauliSum([PauliTerm("X", i, 0.5) * PauliTerm("X", j)])) + mixer_operators.append(PauliSum([PauliTerm("Y", i, 0.5) * PauliTerm("Y", j)])) + +So far so good but now let's look at this walk: + +.. image:: ising_qaoa/valid_sol_not_SAW.png + :align: center + :scale: 75% + +This walk is a feasible solution since it satisfies the hard constraint of having an overall Hamming weight of three. +Yet, it is not a valid SAW. In order to ensure that we are getting valid SAWs we will introduce *soft constraints* into +the cost Hamiltonian. The only soft constraint that we have to enforce is that there shouldn't be a move to the right +followed by a move to the left or vice versa. In order to evaluate if the \\( j \\)-th move went to the right or left, let's define the following the functions: + +$$ m^{j}_{x^{+}} = (1-Z_{3j})Z_{3j+1}$$ +$$ m^{j}_{x^{-}} = (1-Z_{3j+1})Z_{3j}$$ + +Let's define these functions in pyQuil: + +.. code-block:: python + + m_xplus = lambda j: (1-PauliTerm("Z", 3*j))*PauliTerm("Z", 3*j+1) + m_xminus = lambda j: (1-PauliTerm("Z", 3*j+1))*PauliTerm("Z", 3*j) + +From this we can construct our cost Hamiltonian: + +$$ \\mathbf{H}_{C} = \\lambda_{penalty} \\sum_{j=0}^{1} (m^{j}_{x^{-}} \\land m^{j+1}_{x^{+}}) + (m^{j}_{x^{+}} \\land m^{j+1}_{x^{-}}) $$ + +We can generate this with: + +.. code-block:: python + + penalty = 5 + cost_operators = [] + for j in range(2): + term = (m_xminus(j)*m_xplus(j+1)) + (m_xplus(j)*m_xminus(j+1)) + cost_operators.append(penalty*term) + +Initial state!!! + +.. code-block:: python + + from grove.alpha.arbitrary_state.arbitrary_state import create_arbitrary_state + import numpy as np + import pyquil.api as api + import itertools + from functools import reduce + import operator + + connection = api.QVMConnection() + + moves = 3 + zero = np.array([[1,0],]) + one = np.array([[0,1],]) + + tensor = lambda variables: reduce(np.kron,variables)[0] + + probabilities = np.zeros((2**(3*moves),1)).T[0] + + for elements in itertools.product([(one, zero, zero), (zero, one, zero), (zero, zero, one)], repeat=moves): + all_qubits = reduce(operator.add, elements) + probabilities += (1/np.sqrt(moves*3))*tensor(all_qubits) + + initial_state = create_arbitrary_state(probabilities) + +Now we're ready to construct the QAOA instance! + +.. code-block:: python + + n_nodes = 9 # 3 moves with 3 qubits each + num_steps = 2 + qaoa_inst = QAOA(connection, n_nodes, num_steps, cost_operators, driver_operators, initial_state) + +Let's evaluate the wavefunction: + +.. code-block:: python + + t = np.hstack((betas, gammas)) + param_prog = inst.get_parameterized_program() + prog = param_prog(t) + wf = qvm_connection.wavefunction(prog) + wf = wf.amplitudes + + for state_index in range(2**inst.n_qubits): + print(inst.states[state_index], np.conj(wf[state_index])*wf[state_index]) + + Source Code Docs ---------------- From 0efd89e802aff6f2c508b58f955eb85bd24b4481 Mon Sep 17 00:00:00 2001 From: Mark Fingerhuth Date: Tue, 27 Mar 2018 19:46:16 -0400 Subject: [PATCH 2/2] XY mixer doc addition to QAOA completed --- docs/qaoa.rst | 213 +++++++++++++++++++++---------------- docs/qaoa/3_directions.png | Bin 0 -> 4252 bytes docs/qaoa/invalid_saw.png | Bin 0 -> 2451 bytes docs/qaoa/valid_saw.png | Bin 0 -> 2319 bytes 4 files changed, 122 insertions(+), 91 deletions(-) create mode 100644 docs/qaoa/3_directions.png create mode 100644 docs/qaoa/invalid_saw.png create mode 100644 docs/qaoa/valid_saw.png diff --git a/docs/qaoa.rst b/docs/qaoa.rst index c284486..eb6fb10 100644 --- a/docs/qaoa.rst +++ b/docs/qaoa.rst @@ -105,7 +105,7 @@ algorithm for finding "a 'good' solution to an optimization problem" [`1 `_, `2 `_]. What's with the name? For a given NP-Hard problem an approximate algorithm is a -polynomial-time algorithm that solves every instance of the problem with some +polynomial-time algorithm that solves every instance of the problem with some guaranteed quality in expectation. The value of merit is the ratio between the quality of the polynomial time solution and the quality of the true solution. @@ -154,8 +154,8 @@ barbell graph there are two equal weight partitionings that correspond to a maximum cut (the right two partitonings)--i.e. cutting the barbell in half. One can denote which set \\( S \\) or \\( \\overline{S} -\\) a node is in with either a \\(0\\) or a \\(1\\), respectively, in a bit -string of length \\( N \\). The four partitionings of the barbell graph listed +\\) a node is in with either a \\(0\\) or a \\(1\\), respectively, in a bit +string of length \\( N \\). The four partitionings of the barbell graph listed above are, \\(\\{ 00, 11, 01, 10 \\} \\)---where the left most bit is node \\(A\\) and the right most bit is node \\(B\\). The bit string representation makes it easy to represent a particular partition of the graph. Each bit @@ -165,7 +165,7 @@ For any graph, the bit string representations of the node partitionings are alwa length \\(N\\). The total number of partitionings grows as \\(2^{N}\\). For example, a square ring graph -.. image:: qaoa/square.png +.. image:: qaoa/square.png :scale: 55% :align: center @@ -193,9 +193,9 @@ particular quality. For example, a famous polynomial time algorithm is the randomized partitioning approach. One simply iterates over the nodes of the graph and flips a coin. If the coin is heads the node is in \\( S \\), if tails the node is in \\( \\overline{S} \\). The quality of the random -assignment algorithm is at least 50 percent of the maximum cut. -For a coin-flip process the probability of an edge being in the cut is 50\%. -Therefore, the expectation value of a cut produced by random assignment can be +assignment algorithm is at least 50 percent of the maximum cut. +For a coin-flip process the probability of an edge being in the cut is 50\%. +Therefore, the expectation value of a cut produced by random assignment can be written as follows: $$\\sum_{e \\in E} w_{e} \\cdot \\mathrm{Pr}(e \\in \\mathrm{cut}) = \\frac{1}{2} \\sum_{e \\in E}w_{e}$$ @@ -210,7 +210,7 @@ give cuts of expected value at least 0.87856 times the maximum cut [`3 Quantum Approximate Optimization ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -One can think of the bit strings (or set of bit strings) that correspond to the +One can think of the bit strings (or set of bit strings) that correspond to the maximum cut on a graph as the ground state of a Hamiltonian encoding the cost function. The form of this Hamiltonian can be determined by constructing the classical function that returns a 1 (or the weight of the edge) if the edge spans two-nodes in different sets, or 0 if the nodes are in the same set. @@ -219,13 +219,13 @@ C_{ij} = \\frac{1}{2}(1 - z_{i}z_{j}) \\end{align} \\( z_{i}\\) or \\(z_{j}\\) is \\(+1\\) if node \\(i\\) or node \\(j\\) is in \\(S\\) or \\(-1\\) if node \\(i\\) or node \\(j\\) is in \\(\\overline{S}\\). The total cost is the -sum of all \\( (i ,j) \\) node pairs that form the edge set of the graph. -This suggests that for MAX-CUT the Hamiltonian that encodes the problem is +sum of all \\( (i ,j) \\) node pairs that form the edge set of the graph. +This suggests that for MAX-CUT the Hamiltonian that encodes the problem is $$\\sum_{ij}\\frac{1}{2}(\\mathbf{I} - \\sigma_{i}^{z}\\sigma_{j}^{z})$$ where the sum is over \\( (i,j) \\) node pairs that form the edges of the graph. -The quantum-approximate-optimization-algorithm relies on the fact that we can +The quantum-approximate-optimization-algorithm relies on the fact that we can prepare something approximating the ground state of this Hamiltonian and -perform a measurement on that state. Performing a measurement on the \\(N\\)-body +perform a measurement on that state. Performing a measurement on the \\(N\\)-body quantum state returns the bit string corresponding to the maximum cut with high probability. @@ -252,9 +252,9 @@ QAOA identifies the ground state of the MAXCUT Hamiltonian by evolving from a reference state. This reference state is the ground state of a Hamiltonian that couples all \\( 2^{N} \\) states that form the basis of the cost Hamiltonian---i.e. the diagonal basis for cost function. -For MAX-CUT this is the \\(Z\\) computational basis. +For MAX-CUT this is the \\(Z\\) computational basis. -The evolution between the ground state of the reference Hamiltonian +The evolution between the ground state of the reference Hamiltonian and the ground state of the MAXCUT Hamiltonian can be generated by an interpolation between the two operators \\begin{align} @@ -292,8 +292,8 @@ eigenvectors of the \\(\\sigma_{x}\\) operator (\\(\\mid + \\end{align} The reference state is easily generated by performing a Hadamard gate on each -qubit--assuming the initial state of the system is all zeros. The Quil code -generating this state is +qubit--assuming the initial state of the system is all zeros. The Quil code +generating this state is .. code-block:: c @@ -306,7 +306,7 @@ pyQAOA requires the user to input how many slices (approximate steps) for the evolution between the reference and MAXCUT Hamiltonian. The algorithm then variationally determines the parameters for the rotations (denoted \\(\\beta\\) and -\\(\\gamma\\)) +\\(\\gamma\\)) using the quantum-variational-eigensolver method [`4 `_][`5 `_] that maximizes the cost function. @@ -325,11 +325,11 @@ where \\begin{align} U(\\hat{H}_{\\mathrm{ref}}, \\beta_{i}) = e^{-i \\hat{H}_{\\mathrm{ref}} \\beta_{i}} \\end{align} -and +and \\begin{align} U(\\hat{H}_{\\mathrm{MAXCUT}}, \\gamma_{i}) = e^{-i \\hat{H}_{\\mathrm{MAXCUT}} \\gamma_{i}} \\end{align} -\\( U(\\hat{H}_{\\mathrm{ref}}, \\beta_{i}) \\) and \\( U(\\hat{H}_{\\mathrm{MAXCUT}}, \\gamma_{i})\\) can be expressed as a short quantum circuit. +\\( U(\\hat{H}_{\\mathrm{ref}}, \\beta_{i}) \\) and \\( U(\\hat{H}_{\\mathrm{MAXCUT}}, \\gamma_{i})\\) can be expressed as a short quantum circuit. For the \\(U(\\hat{H}_{\\mathrm{ref}}, \\beta_{i})\\) term (or mixing term) all operators in the sum commute and thus can be split into a product of @@ -347,7 +347,7 @@ e^{-i\\hat{H}_{\\mathrm{ref}} \\beta_{i}} = \\prod_{n = H 1 RZ(beta_i) 1 H 1 - + Of course, if RX is in the natural gate set for the quantum-processor this Quil is compiled into a set of RX rotations. The Quil code for the cost function @@ -359,21 +359,21 @@ looks like this: .. code-block:: c X 0 - PHASE(gamma{i}/2) 0 + PHASE(gamma{i}/2) 0 X 0 PHASE(gamma{i}/2) 0 CNOT 0 1 RZ(gamma{i}) 1 CNOT 0 1 -Executing the Quil code will generate the +Executing the Quil code will generate the \\( \\mid + \\rangle_{1}\\otimes\\mid + \\rangle_{0}\\) state and perform the evolution with selected \\(\\beta\\) and \\(\\gamma\\) angles. \\begin{align} \\mid \\beta, \\gamma \\rangle = e^{-i \\hat{H}_{\\mathrm{ref}} \\beta_{1}}e^{-i \\hat{H}_{\\mathrm{MAXCUT}} \\gamma_{1}}e^{-i \\hat{H}_{\\mathrm{ref}} \\beta_{0}}e^{-i \\hat{H}_{\\mathrm{MAXCUT}} \\gamma_{0}} \\mid + \\rangle_{N-1,...,0} \\end{align} In order to indentify the set of \\(\\beta\\) and \\(\\gamma\\) angles that -maximize the objective function +maximize the objective function \\begin{align} \\mathrm{Cost} = \\langle \\beta, \\gamma \\mid \\hat{H}_{\\mathrm{MAXCUT}} \\mid \\beta, \\gamma \\rangle @@ -383,10 +383,10 @@ pyQAOA leverages the classical-quantum hybrid approach known as the quantum-variational-eigensolver[`4 `_][`5 `_]. The quantum processor is used to prepare a state through a polynomial number of operations which is -then used to evaluate the cost. Evaluating the cost +then used to evaluate the cost. Evaluating the cost (\\( \\langle \\beta, \\gamma \\mid \\hat{H}_{\\mathrm{MAXCUT}} \\mid \\beta, \\gamma \\rangle\\)) requires -many preparations and measurements to generate enough samples to accurately -construct the distribution. The classical computer then generates a new set of +many preparations and measurements to generate enough samples to accurately +construct the distribution. The classical computer then generates a new set of parameters (\\( \\beta, \\gamma\\)) for maximizing the cost function. .. image:: qaoa/VQE.png @@ -404,35 +404,33 @@ state with \\(\\beta, \\gamma\\) and sampling. :scale: 55% The probability distributions above are for the four ring graph discussed earlier. -As expected the approximate evolution becomes more accurate as the number of +As expected the approximate evolution becomes more accurate as the number of steps (\\(\\alpha\\)) is increased. For this simple model \\(\\alpha = 2\\) is sufficient to find the two degnerate cuts of the four ring graph. -XY mixer Hamiltonian Example +XY Mixer Hamiltonian Example ~~~~~~~~~~~~~~~~~~~~~~~ -In this example we will use the important optional argument ``driver_operators`` to solve a simple toy problem -using the Ising QAOA wrapper. This example is loosely based on the work 'QAOA with hard and soft constraints' by Hadfield et. al (2017) -[`4 `_]. Here, the core idea is to implement hard constraints into the +In this example we will use the optional argument ``ref_hamiltonian`` to solve a simple toy problem +using the QAOA. This example is loosely based on the work 'QAOA with hard and soft constraints' by Hadfield et. al (2017) +[`4 `_]. The core idea here is to implement hard constraints into the driver/mixer Hamiltonian which ensures that the search takes place only in the feasible subspace. Just keep reading and it will get -clear what *hard constraints* and *feasible subspace* means. -Suppose we encode the directions right, left and up into 3 bits each: +clear what *hard constraints* and *feasible subspace* means. Suppose we encode the directions right, left and up into 3 bits each: -.. image:: ising_qaoa/3_directions.png +.. image:: qaoa/3_directions.png :align: center :scale: 75% -All of these three bit strings have one thing in common: a Hamming weight of one. (basically the number of -1's in the bit string [`5 `_]). Let's say we want to use our -Ising QAOA to find self-avoiding walks (SAW) with three moves. Due to the way we assigned bit strings to -directions we know that every *feasible* solution for a walk with three moves must have a Hamming weight -of three! Hence, our *hard constraint* is a constant Hamming weight of one for each triplet of bits. -Consequently, the bit string \\( 110 \\, 111 \\, 011 \\) is not a feasible solution since its Hamming weight is -too large and the three triplets carry no meaning in the context of our problem! Thus, we can define the -*feasible subspace* as the set of all *feasible* bit strings. -For example, a valid SAW with three moves (a feasible solution) could look like this: - -.. image:: ising_qaoa/valid_saw.png +All of these three bit strings have one thing in common: a Hamming weight of one (basically the number of +1's in the bit string [`5 `_]). Let's say we want to find +self-avoiding walks (SAW) with two moves. Due to the way we assigned bit strings to directions we know that +every *feasible* solution for a walk with three moves must have a Hamming weight of three! Hence, our *hard constraint* +is a constant Hamming weight of one for each triplet of bits. Consequently, the bit string \\( 110 \\, 111 \\, 011 \\) +is not a feasible solution since its Hamming weight is too large and the three triplets carry no meaning in the context +of our problem! Thus, we can define the *feasible subspace* as the set of all *feasible* bit strings. For example, a valid +SAW with two moves (a feasible solution) could look like this: + +.. image:: qaoa/valid_saw.png :align: center :scale: 75% @@ -446,51 +444,60 @@ $$ SWAP_{i,j} = \\frac{1}{2} ( X_{i}X_{j} + Y_{i}Y_{j} ) $$ The mixer should mix amplitudes between feasible solutions which means it should swap around the bits within each triplet. This can be achieved with the following driver/mixer Hamiltonian: -$$ \\mathbf{H}_{M} = \\sum_{k \\in \\{ 0,3,6\\}} \\sum_{i=k}^{k+1} \\sum_{i+1}^{k+2} SWAP_{i,j} $$ +$$ \\mathbf{H}_{M} = \\sum_{k \\in \\{ 0,3\\}} \\sum_{i=k}^{k+1} \\sum_{i+1}^{k+2} SWAP_{i,j} $$ -Let's first do our imports: +Let's first do our imports and open a QVMConnection: .. code-block:: python - from grove.alpha.arbitrary_state.arbitrary_state import create_arbitrary_state - import numpy as np - import pyquil.api as api import itertools - from functools import reduce import operator + from functools import reduce + + import numpy as np + import pyquil.api as api + from grove.pyqaoa.qaoa import QAOA + from pyquil.paulis import PauliSum, PauliTerm + from pyquil.gates import * + from pyquil.quil import Program + from grove.alpha.arbitrary_state.arbitrary_state import create_arbitrary_state + + connection = api.QVMConnection() In pyQuil code we can generate this driver/mixer Hamiltonian like this: .. code-block:: python + moves = 2 mixer_operators = [] - for k in [0,3,6]: - l = list(range(k,k+3)) - ij_pairs = [ (a,b) for a in l for b in l[l.index(a):] if a!=b and b-a<=2] + for k in range(0, 3*moves, 3): + l = list(range(k, k+3)) + ij_pairs = [ (a, b) for a in l for b in l[l.index(a):] if a != b and b - a <= 2] for i, j in ij_pairs: mixer_operators.append(PauliSum([PauliTerm("X", i, 0.5) * PauliTerm("X", j)])) mixer_operators.append(PauliSum([PauliTerm("Y", i, 0.5) * PauliTerm("Y", j)])) So far so good but now let's look at this walk: -.. image:: ising_qaoa/valid_sol_not_SAW.png +.. image:: qaoa/invalid_saw.png :align: center :scale: 75% This walk is a feasible solution since it satisfies the hard constraint of having an overall Hamming weight of three. -Yet, it is not a valid SAW. In order to ensure that we are getting valid SAWs we will introduce *soft constraints* into -the cost Hamiltonian. The only soft constraint that we have to enforce is that there shouldn't be a move to the right -followed by a move to the left or vice versa. In order to evaluate if the \\( j \\)-th move went to the right or left, let's define the following the functions: +Yet, it is not a valid SAW since it intersects with itself. In order to ensure that we are getting valid SAWs we will +introduce *soft constraints* into the cost Hamiltonian. The only soft constraint that we have to enforce is that there +shouldn't be a move to the right followed by a move to the left or vice versa. In order to evaluate if the \\(j\\)-th +move went to the right or left, let's define the following functions: -$$ m^{j}_{x^{+}} = (1-Z_{3j})Z_{3j+1}$$ -$$ m^{j}_{x^{-}} = (1-Z_{3j+1})Z_{3j}$$ +$$ m^{j}_{x^{+}} = (1 - Z_{3j})Z_{3j+1}$$ +$$ m^{j}_{x^{-}} = (1 - Z_{3j+1})Z_{3j}$$ -Let's define these functions in pyQuil: +We can define these functions in pyQuil: .. code-block:: python - m_xplus = lambda j: (1-PauliTerm("Z", 3*j))*PauliTerm("Z", 3*j+1) - m_xminus = lambda j: (1-PauliTerm("Z", 3*j+1))*PauliTerm("Z", 3*j) + m_xplus = lambda j: (1 - PauliTerm("Z", 3 * j))*PauliTerm("Z", 3 * j + 1) + m_xminus = lambda j: (1 - PauliTerm("Z", 3 * j + 1))*PauliTerm("Z", 3 * j) From this we can construct our cost Hamiltonian: @@ -500,60 +507,84 @@ We can generate this with: .. code-block:: python - penalty = 5 + penalty = 10 cost_operators = [] - for j in range(2): - term = (m_xminus(j)*m_xplus(j+1)) + (m_xplus(j)*m_xminus(j+1)) - cost_operators.append(penalty*term) + for j in range(moves - 1): + term = (m_xminus(j) * m_xplus(j + 1)) + (m_xplus(j) * m_xminus(j + 1)) + cost_operators.append(penalty * term) -Initial state!!! +In vanilla QAOA with the standard X mixer we start from a superposition over all states by applying a Hadamard gate to each qubit. But in +this case, we would like to start from the superposition of strings that constitute valid SAWs with two moves: -.. code-block:: python +$$ \\mid\\psi> = \\frac{1}{\\sqrt{9}} \\big( \\mid001001> + \\mid001010> + \\mid001100> + ... + \\mid100100> \\big) $$ - from grove.alpha.arbitrary_state.arbitrary_state import create_arbitrary_state - import numpy as np - import pyquil.api as api - import itertools - from functools import reduce - import operator +The following code makes use of grove's ``create_arbitrary_state`` function in order to generate a circuit that will generate +this initial superposition: - connection = api.QVMConnection() +.. code-block:: python - moves = 3 - zero = np.array([[1,0],]) - one = np.array([[0,1],]) + zero = np.array([[1, 0],]) + one = np.array([[0, 1],]) tensor = lambda variables: reduce(np.kron,variables)[0] - probabilities = np.zeros((2**(3*moves),1)).T[0] + probabilities = np.zeros((2 ** (3 * moves),1)).T[0] + quantum_states = itertools.product([(one, zero, zero), + (zero, one, zero), + (zero, zero, one)], repeat=moves) - for elements in itertools.product([(one, zero, zero), (zero, one, zero), (zero, zero, one)], repeat=moves): - all_qubits = reduce(operator.add, elements) - probabilities += (1/np.sqrt(moves*3))*tensor(all_qubits) + for qubits in quantum_states: + all_qubits = reduce(operator.add, qubits) + probabilities += (1 / np.sqrt(moves * 3))*tensor(all_qubits) initial_state = create_arbitrary_state(probabilities) -Now we're ready to construct the QAOA instance! +Now we're ready to construct the QAOA instance by providing it our custom SWAP mixer and the quantum circuit that generates the initial state: .. code-block:: python - n_nodes = 9 # 3 moves with 3 qubits each - num_steps = 2 - qaoa_inst = QAOA(connection, n_nodes, num_steps, cost_operators, driver_operators, initial_state) + n_nodes = 3 * moves # moves with 3 qubits each + num_steps = 6 + inst = QAOA(qvm=connection, n_qubits=n_nodes, steps=num_steps, cost_ham=cost_operators, + ref_hamiltonian=mixer_operators, driver_ref=initial_state, store_basis=True) Let's evaluate the wavefunction: .. code-block:: python + betas, gammas = inst.get_angles() t = np.hstack((betas, gammas)) param_prog = inst.get_parameterized_program() prog = param_prog(t) - wf = qvm_connection.wavefunction(prog) + wf = connection.wavefunction(prog) wf = wf.amplitudes - for state_index in range(2**inst.n_qubits): - print(inst.states[state_index], np.conj(wf[state_index])*wf[state_index]) - + for state_index in range(2 ** inst.n_qubits): + print(inst.states[state_index], np.conj(wf[state_index]) * wf[state_index]) + +This yields a long list of probabilities (most of which are \\(0\\)). Here, I am only printing +a small subset of them:: + + 000000 (4.814824860968226e-35+0j) + [...] + 001001 (0.25374558548209253+0j) + 001010 (4.4094400090988824e-05+0j) + 001011 (3.759565398379776e-32+0j) + 001100 (0.1075389948215445+0j) + [....] + 010010 (0.10332997853766406+0j) + 010011 (2.1304324480120207e-32+0j) + 010100 (0.038396483138820936+0j) + [...] + 100001 (0.10753899482154447+0j) + 100010 (0.038396483138820985+0j) + 100011 (6.39803705190637e-32+0j) + 100100 (0.35096529125930287+0j) + [...] + 111111 (4.333342374871216e-34+0j) + +As you can see, the quantum states, \\(010100\\) and \\(100010\\), which violate the soft constraints +end up with the lowest probabilities whilst the other solutions range between 10-35% probability. Source Code Docs ---------------- diff --git a/docs/qaoa/3_directions.png b/docs/qaoa/3_directions.png new file mode 100644 index 0000000000000000000000000000000000000000..3702de369c1413a264c3d7db1c5d12f128863427 GIT binary patch literal 4252 zcmcIoXH*l)woV9LN~nS~=^X(9iC_>2EkFRlBTYiU5RphPA}x_n0)o=TplB?h!XbbZ zsiBFWzy&Fxqewf1suV+e@!j|1z5D;oT2uC{y}tSO%&c$r8loe@ikJHgHvj z0RVs^%yB#?JG13}n>4`eSi{V0TsWB}mh(m`)8-1cb_)Xl_&WboV4gak1ane6+`={7 zImj>6zJk)HmJyLvqgv%c~b;}vznnEARlGnTyUfaQ#3GZYp( z#ltcKSZCMb?~NV@?aNV&C_rtQM{HBvo~JYhRxo0eOWFc7!JGWY@m!l%@OKE{<}apH zqsx8OngD2Fhl5104&`FkSoNcuuR@eIStz~pzO7Ryg<$enA*V)L;OGcdlNv(B;8ux& zbhPVd`HV1P;9(v9C;kQg9MHiq;v7K!%HdQxGFL&(>=$YRS3(RN6%xck?T5D?*$lH0 zs}(7_6cPw4{Vs?jFfubD9`2dY>8>}%f@F{EWJQqWU@{op!Hm0@v}yay=|O3dgcNDs z*dTx<3`FW`+pyt*?@#7gyAYOtPF-%1h*0icujATDth|Bc70PH+ZP5yim2y7&7qJ&` zcEw-rZP#(z0Gev98tLFeoC+@9|1irXCZ2o3U1q(dKdmbBP#v!0{jk)9RD<71+~yfV z{!9+KO%*JY=>DfJveyA4xhFmlZ)WYH*3VjUc$Y6`%r9LS={D~+0~Eb`m)An)|XU<0xs3!xi5*(SSAe`jY1g2b38Zp*pzkZgsirk-}! zX}LcuDAM!*vH&7{Bzi!(YB3hWZ?7~*Ul+?E9=4*^{KSU8gN7q$W*gbwMzC)WZ6g+3#mg0 zBSGHzYw5i zn5>?y+e;lW!>SR})LXJn0{(Qi*ta95oU#o50w1L0>yhmwFz+`_EO{1M&Ifbhh)5tj z;0LX;q2st3KV`ry_6<8#HwLRz2;ie+4x*@|J24JZ1vbm|5mr)I@5QFn=*c;W4<8BG zNT!w_FIEMqvj~{|;tza%aHi^9ot~U54DDYMzY`ZVS(X$+OnbKeAhlY8Qd@$U9g)|~ zDor)yC>geKA%K$#O&NPfokI0yoZrySRjt1Q#_}zc)$k1_?dVWG(h*5I2~3{VenQXh z?a)j|a5||bgFb%}&C8>H675n^!!?-5-AM4Z=4<4m8W%bC+jBvPUE8DVliAEw4I{^yueWpsz@j~}6+L}lj6hi1hrZdr!0FK#6K zbjH=T*(K%S5gZ$p!yR0)JYRYs&-eVG4jc*W=14$t#-2yD>f@h0KEqU->?zZ@Q#T;N ztOnAr*qcOYPt5m5d!w7B%E08H^Q^A0=4X=YOr!+J?re#di}YHJbrLvl?yMZ%m0j|m z;RyDGVv$LV+D>%AU?2Y+hT;1LUNyhOdmm(n2q_((J@0vPnmSaJ~gu#{6K z+QK6BqJmm`3d|zcj4T&x?E+I}npSsimj(I@lS7%}#$WUAx&LrU2~{Ht3yOR**A#K$ zCZ=sNFW{>`h=Tvugb$}(`F{qnQc>%Dlgd;^!h_8TM_UWwKY=ih6V9XrRrnv4kbl#9 zUX-Si#Ju{BdshED-v16}ZCuF*w{5U_c2j#3k?t#C#;N1zA8Cu`p+K~`U0KRBRxsMy zd^8P}bL_WxQyn-56%2rxfOvo%$1z)zAoey}94`w3==Cx1FxZ3rkZ<1+h+Jj42K?#U zGQ%3r%A+nuT1E;t8@%gu^k$)4qL`-pa*4AtA4(pYa12#hH-#Z#&>{i}F~^yp_dQr7 zA9suyJi7DpH-m!`kKcFiyE$iVr=nbEn19@^O3k1xU99S=)VtS*N22zKG3Bj1MeySS z^?CdY98(&>`hx4~0$&~HVyFo|I(-~}14q+0z^zUvDK%=nsl5Y5M_#Oa9WJ{aR2T27 zj&)+Ri(Lcm)PT$i4;=9kj*hQ+-QLiR&t%crRk}^(qmEPMX(zL6MYP^$4B6A(W9^}H z1Jt_8lq(~K09P0z@fVb?w5P*oTUK%@HJxIx8Ea`5MiCu3u;n3x?gr;V(yj;svBl=2P9F=N9f)ioZw2eW|95RiNy)?C%Z5 zW)l>r@?B-xQtnWF2)XlAwad08yp{>$4F}P8uHtozu)G*jfO3VVOm?Wi_Cvl z`u?2P!&q3LwViY$v{t6G&Kb50AWBVWZ(5b+>-R3!S_o1vp@m10%}?n~4k{|e8(sVC zx~H?^ptq>2u7RuC?z3rymxp^3fC)1{gkW1IwYZZi^7i}c;tC^3uFAPaLhe9SX({~& z?F8BpdCxZ{!L0e@O z)gq$22680XYaK$k_DeqNCoN*&B;&pWPg1Xg+J7t+#YYG9qlGZ~V(RxGij zHfsM`>l12Y;yxf>Y5y&XYf>&gO_zrzVaRF-rqvT3DBPl6!>gb9;gCE3Wo}9Q(&iwL zE{gU3@F=G8?}Kn}6n|`Dg#SnOg}P7D5ywF##w{t`2;cK1hs8XTcL&6^e=2r zy<9JPaka7`_-&lDNnpNl!1RL}To`V?GON7r2FK$u7p;zAFX)UFB^H)YA|+9wBd{uc z9pC)d_3T<_o~XAy+mgFz>6d5~p=;(A=t}4QoUBMI)|P^~04eOJle0r(KbL=XKJ1L0 zBhj2=FmtTlCK{iomq&Z&X>)aD`F{mvjl|y8Ywg}=APe@oSR)Fm#Tzc;VAm&ix7;8<+(~>j7ZY7RRLC(Uid}1~-0On*UF)sLjI!Gd zd&X_cgmxFmxjXmb)Ez$CUPOm9EW?W{NOiPX>|S4W5BW>td!DV7m@@(3h=Q5hxusu% zd)+<|)>18}`mIYJ9Ng9IDu^jlKw7M{TM9rk?5}v7wyi~coB?O3EjzRioiQ^VZ|rPA z1^kGgL(ntq&c8vK8`ic3&Y;2_xa#ZqhY}2Qw8r*aEAqa*1~>8Vk>_0|QY{%tZr(DO zi!C{pkJrBs&pr*AqR32IwEtRV|C&jpt~ral32eZXPgJA`?sV&iTlcr14vRw5@q6#BeuL!2(QM&ub8|Wt?O$W^@ilfUeELklXJo@ zmC$$<-ABFoepi&}l6BK=UhLHvil9CzZ+b9L$`l1_y~B9mbc#eBikYpo`hz+Uk#uSD z0e8?U>U@qlNw;sxatL4kQ&%M5(r9za(E|}cquOM#le+Fg83hF6jZ`0uQfY=Lg?+X zX(={eZ4}Ah;PUuA&%LjuuGZ+3FpxZz;qLKWzYA_Qa0HjzDJzhxfqJY%QsZ2MMsa z(+BT)kQbb0!JN)~D7fs5qv3_FvAq74ARqpn7ZkMdiFA4tr)$9t?#RIL;ub2*a*aBX#X^$HYV4^EvK}Iaiy&&bxB-VZ6qUq;zn641pViG0bpZ+fcP)$|> zAn+vcOu6C|gG+&(z$)PVz>G45X+w+y-UwUHe@~G#$^ClSc;@AfX@>@K+d%XxEJ^oFuM$4 zbHq6GTKP7x8NH^Hj>9hC#O90KT1cZf_Ieyx2`mF9HeKrUn1igz&B#(@>2swp2A?C_ zAUfAxWQ8sUUfE>X!!Qw8ifc_q(?3UHUf^d~jRy088_4dTwC&Y_OnIxI1FaGtQV{zBrB3)dlcvd-Nb@HqWPG5V#Zh8Amm#AJdU9^)ukv9Q*RsCduD6=mG8{pTXUNW-Lo7 ze?hPd`HU9>zXkR)y1uPwuK2Bs9$+2W2XZHxtz1OwZ$0!OE3^{WPNg`RkX2!{Hvt>|5^$1wgh^C^&}glr-AE%ccZ6uETRwf?|RhVN^yMXMQaaX12O_^ zM5cHFdMf{)&l>dN>m^w1HVAl z1ssE7i9kt+{*!^vWZ*5p(dZQCb`<^iXa-($fzyGRC}!s$z#4KoV|5O27mCNNp<}R{ z(Q(=|$4zFS%uOqR8_|yU&A_WN$aFODJ>Urx9x#G_h^v8DWZ;)KW~0M@SE3M*-FV|* zKAMu>IFzDgFTnv>ca!|Y#StjtASm&`5O6u82+Hq9KW<#S?wlAt-v#~^!DCTMxt7!T zBiL&%@V=Bh^TtY|$1RQ0bfB|OPTT#!Lr59zVZEn_SDfiUziG1@#W+2PH$reV86zgI zF--b_PXJ4R^MGHPx<`=N48#324w}{SG)l7w9kq1=xvoOi>rAAaL%1QMZ%)BCZ#YVG z038Kh1>8e6rTLg?_gB1@xfh)nb?P?~?RHoht^=l+y17Vsocb@OkP%S+6sPguffxDe zJ5l(dLvvXKuhn&7FGm4GbATU_UTM0Jf;vCOIRP|FaG#id$~*gXD)1&0B<0`>q*t0= zU`M^p9`src!WzT-?IAR$8JZ0QTauZ$hx~St@_hPQbA!`EXbwfrTn(LOOp~D@i_)BH z+T0L9ceSDWKnj_PMi!-cmud5Z2)g-(Zfyw7O~n7EAOqmgeI%irOhq%K8Zxa>=I#n} zgL9MKfNV-*I~|pv@Hj(vR|=VmMi!;9KIN$~w2#;6hG}#cx^exVokaU9cKjL?Fs$Js z+(T&2aoI53XUeBV(B%P#{~+){1l^yJNe8_m+Bs{elL}@8_)!OJo&Oc~6 z9fxhEe59{#rYR@kWoc~NP^GccuxDeK*Kk)5r$iBs5@}@e(-bLJk+2hsJqhS7&7=rk zaTd$ZtJhY9DUBr>N^k>5hzS%o#Z4>Ie?!`~6KvMfBng=8xb*nyCYW-PZJ=RO_V5Z* znkdIN2zGXrwrx@r#kZAV6+g+tXqC)wO{_vC;boP~W97n>W-vsMFgW{9!{n2aEC^X; z0xgxdBW$1jmLmxEl{Dz98w{aK%_~f422Adi+JoRV1Fq zoh8mnKyQ;eZS&Px$JXy8$MM>3E5ekWAh1(ggHY?jq$X68mtc zzO#t85HB<3zxnDOG3A%X@ae+cuGoPV0}mohY3vP7h`E7X{MiyiI~wf3{S&NXiOzu4 z0lcm_=S&!>IFC}^# z!L9<|6G4{;9R5AH0~CbgaWiJEqqT#x;&6v4FOQ%*3H_T3C>@*Mrx`(;nxG|!s0e`b zO!@W%G9JVoxcd)hh0sjF4a^!KnSIG)D9v$!iPRtBVGmbM@HEkowJCx}lyHw`(-PMdvbFK;W}z_?Xp z%qP@^H+nDV+C+T*-t%L4_^;t$HV3hc=xpH6qbIja~ED~KrhrkA5fas6~(DKG&v>@fgX z3lxB#fG2^;K&Nt&s&I?|YLRPJ13U_x4RkCgsVt#ImAGQ%PqwxdV6(nRI*v?1ALe*_z6)z}3hD(y7y@1)k<;Ze&j%W-*u2~*PmLo@w96562 z$dMyQjvP61VUx#^BT)3DrCU=!@)#e-FG2d=sNR06HLH z2*1!Y9^)Tiks=+1ZomZKR^TP%Eys1mG9?DDF~D2Ezs-2924+w>Arsw@{)R;qPKs#b zP$a~12@;dryp0VpcwL0VM4Na$S3Iv1x&mK>G+IE?Xq^?&&V1xOAX6KcJMx-AOQAF` zi-Bs#c9TX2H-N2nlSmF`Ic66emM;zKs19vMf zAJH-n=t}&))M9?3J;5{ppCHMWrm16H1k)4AF<3xeeZ8dF_lE(WDc-Ye8=k5T`jX%P zp4*hoqL8*`1J59yVHD6FxtXqlX@~r_ddFQ5VABUbxcpzr7oFA}|7qQU`8AkF@%wX^ z!KM?EwygQn-iW(k9#ZGK-(WKU8K-(ewh96LP&Is-!KNA+Boz=BH_(za#1qU!bR^y^*+En@-%W8lD+R z*+(3~%*Cz~o$!OEH&UJfupa16h_GF$Y==3V!&NY&)Op^Gu$`i8=X)QTUIZJ9@BTiB zXjp$Abzw@y_~PIF$e6%{UC5XVWNdS6ikwUsRs35(CJtfvJm9!Apxa($v=V4eqrHkr zFvqL<>WGF{BNqkat$?)90U7*(JW<;t>^1MKQ*&vqwa8#K4C_@rO|H7_^+qefc4X`Y z^1xMrSujVb`fm{pdo!Nq$0dz{$c+H_-NIh;9%{5o2i@_xZn3c6tm=J|rxV!rW}{V0 zo6f~z6pVgJ+a_JPr^Jq=mnTi9pV(ov@}YTW{nELf9@}KP(l*YEWe(_UIHw8PCR8w* zw9WFrZaA%*qdU9D!afU5Cg{|rdoApBf6`1$=*p{6)ei=J#Rop?TTl>(#s(z$_R|+UAYBEIM3?oJ`Q?v-ye99xzR= zuDrBOJ(vZf#|XVF8un}&XzIq2M4hT1W?`=vrRt1U>7Wk(=JdC)AE4?s1Nyko>orll zs(>Sl^{o_#CMLoBsOp0v8t#GY$U?e_^ETqU)WH_^n)i=Jt8`F6R=`1KOnoElb&8Jf zjaH@OJ5`Tbn;MAkrT;(9B$ygikIGu4zjs~QcENsp;=PI9hS0ofj8^GEFN}?iusx&r z*UQtF7FOY}nb~skjswhsS&9ENpN(u5nlMq---#emDs=njhS>GNU!ANWUc5-DSgPt* zhS-h7U+`RszowN2HXtcyLDH;*Ewo;x>aW=}zL@B}o-X+R9ZNmper*wFFKM(o<}1vC zzem~K>%b{Vv`2oT3&<4q)Bu|*z;t2+t#@5^y@ck*s>T>}Ze3mn#)y4Q+|3EhMyo=ulSfIZDg774d zQ;=+-x_BGJbphX*_VwE1Z|4+%XYnTyRlxbcpT+O?6Yj}dK;P7Dzyrura1%`<@oET9 z1bilbIHp&;Af5uUeBS(fQ_`WKGz~3W0bGZ;&?k_eQm?bux^y-r%%*AZ*j6OSt5JRS zJJGWi?TKv4f{?a`A_Deo_Y^tlptb^Vhf7C66hDXX@XD`=o=HuL=v+_eEX{e@GHem4jws9ikfa{SQv4ibCC%(sr`OB>aQryYhYh;2d!xtlDHKNNpM(*})1 p?~;S80eSazC$a?bF^2z(;Xl(NNEXA}^56gf002ovPDHLkV1n*pTQdLv literal 0 HcmV?d00001