Shortcuts



Note

Run interactively in jupyter notebook.

Warm-starting quantum optimization

Introduction

Optimization problems with integer variables or constraints are often hard to solve. For example, the Quadratic Unconstrained Binary Optimization (QUBO) problem, i.e.

\begin{align} \min_{x\in\{0,1\}^n}x^T\Sigma x + \mu^Tx, \end{align}

is NP-Hard. Here, \(\Sigma\) is an \(n\times n\) matrix and \(x\) is a vector of \(n\) binary variables. Note that we could have added the linear term \(\mu\) to the diagonal as \(x_i^2=x_i\) for \(x_i\in\{0, 1\}\). While QUBOs are hard to solve there exists many ways to relax them to problems that are easier to solve. For example, if \(\Sigma\) is semi-definite positive the QUBO can be relaxed and results in a convex Quadratic Program

\begin{align} \min_{x\in[0,1]^n}x^T\Sigma x, \end{align}

which becomes easy to solve as \(x\) now represents \(n\) continuous variables bound to the range \([0, 1]\). Such relaxations can be leveraged to warm-start quantum optimization algorithms as shown in [1].

References

[1] D. J. Egger, J Marecek, S. Woerner, Warm-starting quantum optimization, arXiv:2009.10095

[1]:
import numpy as np
import copy

# Problem modelling imports
from docplex.mp.model import Model

# Qiskit imports
from qiskit import BasicAer
from qiskit.utils import QuantumInstance
from qiskit.algorithms import QAOA, NumPyMinimumEigensolver
from qiskit.utils.algorithm_globals import algorithm_globals
from qiskit_optimization.algorithms import MinimumEigenOptimizer, CplexOptimizer
from qiskit_optimization import QuadraticProgram
from qiskit_optimization.problems.variable import VarType
from qiskit_optimization.converters.quadratic_program_to_qubo import QuadraticProgramToQubo

Preliminaries: relaxing QUBOs

First, we show how to relax a QUBO built with a semi-definite positive matrix to obtain an easy-to-solve QP.

[2]:
def create_problem(mu: np.array, sigma: np.array, total: int = 3) -> QuadraticProgram:
    """ Solve the quadratic program using docplex. """

    mdl = Model()
    x = [mdl.binary_var('x%s' % i) for i in range(len(sigma))]

    objective = mdl.sum([mu[i]*x[i] for i in range(len(mu))])
    objective -= 2*mdl.sum([sigma[i, j]*x[i]*x[j] for i in range(len(mu)) for j in range(len(mu))])
    mdl.maximize(objective)
    cost = mdl.sum(x)
    mdl.add_constraint(cost == total)

    qp = QuadraticProgram()
    qp.from_docplex(mdl)

    return qp

def relax_problem(problem) -> QuadraticProgram:
    """Change all variables to continuous."""
    relaxed_problem = copy.deepcopy(problem)
    for variable in relaxed_problem.variables:
        variable.vartype = VarType.CONTINUOUS

    return relaxed_problem

For this example, we use a positive semi-definite matrix \(\Sigma\) and a linear term \(\mu\) as defined below.

[3]:
mu = np.array([3.418,  2.0913,  6.2415,  4.4436, 10.892 ,  3.4051])
sigma = np.array([[ 1.07978412,  0.00768914,  0.11227606, -0.06842969, -0.01016793, -0.00839765],
                  [ 0.00768914,  0.10922887, -0.03043424, -0.0020045 ,  0.00670929, 0.0147937 ],
                  [ 0.11227606, -0.03043424,  0.985353  ,  0.02307313, -0.05249785, 0.00904119],
                  [-0.06842969, -0.0020045 ,  0.02307313,  0.6043817 ,  0.03740115, -0.00945322],
                  [-0.01016793,  0.00670929, -0.05249785,  0.03740115,  0.79839634, 0.07616951],
                  [-0.00839765,  0.0147937 ,  0.00904119, -0.00945322,  0.07616951, 1.08464544]])

Using DOCPLEX we build a model with binary variables.

[4]:
qubo = create_problem(mu, sigma)
qubo
[4]:
\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: docplex_model1

Maximize
 obj: 3.418000000000 x0 + 2.091300000000 x1 + 6.241500000000 x2
      + 4.443600000000 x3 + 10.892000000000 x4 + 3.405100000000 x5 + [
      - 4.319136480000 x0^2 - 0.061513120000 x0*x1 - 0.898208480000 x0*x2
      + 0.547437520000 x0*x3 + 0.081343440000 x0*x4 + 0.067181200000 x0*x5
      - 0.436915480000 x1^2 + 0.243473920000 x1*x2 + 0.016036000000 x1*x3
      - 0.053674320000 x1*x4 - 0.118349600000 x1*x5 - 3.941412000000 x2^2
      - 0.184585040000 x2*x3 + 0.419982800000 x2*x4 - 0.072329520000 x2*x5
      - 2.417526800000 x3^2 - 0.299209200000 x3*x4 + 0.075625760000 x3*x5
      - 3.193585360000 x4^2 - 0.609356080000 x4*x5 - 4.338581760000 x5^2 ]/2
Subject To
 c0: x0 + x1 + x2 + x3 + x4 + x5 = 3

Bounds
 0 <= x0 <= 1
 0 <= x1 <= 1
 0 <= x2 <= 1
 0 <= x3 <= 1
 0 <= x4 <= 1
 0 <= x5 <= 1

Binaries
 x0 x1 x2 x3 x4 x5
End

Such binary problems are hard to deal with but can be solved if the problem instance is small enough. Our example above has as solution

[5]:
CplexOptimizer().solve(qubo)
[5]:
optimal function value: 16.7689322
optimal value: [0. 0. 1. 1. 1. 0.]
status: SUCCESS

We can create a relaxation of this problem in which the variables are no longer binary. Note that we use the QuadraticProgramToQubo converter to convert the constraint into a quadratic penalty term. We do this to remain consistent with the steps that the Qiskit optimization module applies internally.

[6]:
qp = relax_problem(QuadraticProgramToQubo().convert(qubo))
qp
[6]:
\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: docplex_model1

Maximize
 obj: 259.553391640000 x0 + 258.226691640000 x1 + 262.376891640000 x2
      + 260.578991640000 x3 + 267.027391640000 x4 + 259.540491640000 x5 + [
      - 89.697600360000 x0^2 - 170.818440880000 x0*x1 - 171.655136240000 x0*x2
      - 170.209490240000 x0*x3 - 170.675584320000 x0*x4 - 170.689746560000 x0*x5
      - 85.815379360000 x1^2 - 170.513453840000 x1*x2 - 170.740891760000 x1*x3
      - 170.810602080000 x1*x4 - 170.875277360000 x1*x5 - 89.319875880000 x2^2
      - 170.941512800000 x2*x3 - 170.336944960000 x2*x4 - 170.829257280000 x2*x5
      - 87.795990680000 x3^2 - 171.056136960000 x3*x4 - 170.681302000000 x3*x5
      - 88.572049240000 x4^2 - 171.366283840000 x4*x5 - 89.717045640000 x5^2 ]/2
      -384.203087460000
Subject To

Bounds
       x0 <= 1
       x1 <= 1
       x2 <= 1
       x3 <= 1
       x4 <= 1
       x5 <= 1
End

The solution of this continuous relaxation is different from the solution to the binary problem but can be used to warm-start a solver when dealing with the binary problem.

[7]:
sol = CplexOptimizer().solve(qp)
sol
[7]:
optimal function value: 17.012055025682628
optimal value: [1.75249958e-01 1.48038882e-07 9.70905326e-01 7.38416868e-01
 9.99999992e-01 1.44389045e-01]
status: SUCCESS
[8]:
c_stars = sol.samples[0].x
c_stars
[8]:
[0.17524995761801193,
 1.4803888163984616e-07,
 0.9709053264087679,
 0.7384168677494152,
 0.9999999916475085,
 0.14438904470168312]

QAOA

Here, we illustrate how to warm-start the quantum approximate optimization algorithm (QAOA) by leveraging the relaxed problem shown above.

Standard QAOA

First, we use standard QAOA to solve the QUBO. To do this, we convert the QUBO to Qiskit’s QuadraticProgram class (note that the resulting problem is still a binary problem).

[9]:
algorithm_globals.random_seed = 12345
quantum_instance = QuantumInstance(BasicAer.get_backend('statevector_simulator'),
                                   seed_simulator=algorithm_globals.random_seed,
                                   seed_transpiler=algorithm_globals.random_seed)
qaoa_mes = QAOA(quantum_instance=quantum_instance, initial_point=[0., 1.])
exact_mes = NumPyMinimumEigensolver()
[10]:
qaoa = MinimumEigenOptimizer(qaoa_mes)
[11]:
qaoa_result = qaoa.solve(qubo)
print(qaoa_result)
optimal function value: 16.768932200000002
optimal value: [0. 0. 1. 1. 1. 0.]
status: SUCCESS

Warm-start QAOA

Next, we compare this result to a warm-start QAOA in which we use the solution to the continuous relaxation of the problem. First, we create the initial state

\begin{align} |\phi^*\rangle=\bigotimes_{i=0}^{n-1}R_y(\theta_i)|0\rangle_n . \end{align}

which is given by applying \(R_y\) rotations with an angle \(\theta=2\arcsin(\sqrt{c^*_i})\) that depends on the solution to the relaxed problem. Here, \(c^*_i\) the value of variable \(i\) of the relaxed problem.

[12]:
from qiskit import QuantumCircuit

thetas = [2*np.arcsin(np.sqrt(c_star)) for c_star in c_stars]

init_qc = QuantumCircuit(len(sigma))
for idx, theta in enumerate(thetas):
    init_qc.ry(theta, idx)

init_qc.draw(output='mpl')
[12]:
../_images/tutorials_10_warm_start_qaoa_20_0.png

Next, we create the mixer operator for QAOA. When warm-starting QAOA we must ensure that the mixer operator has the initial state as ground state. We therefore chose the Hamiltonian

\begin{align} H_{M,i}^{(ws)}= \begin{pmatrix} 2c_i^*-1 & -2\sqrt{c_i^*(1-c_i^*)} \\ -2\sqrt{c_i^*(1-c_i^*)} & 1-2c_i^* \end{pmatrix} \end{align}

as mixer operator for qubit \(i\). Once multiplied by \(-i\beta\) and exponentiated this mixer produces the following mixer circuit.

[13]:
from qiskit.circuit import Parameter

beta = Parameter('β')

ws_mixer = QuantumCircuit(len(sigma))
for idx, theta in enumerate(thetas):
    ws_mixer.ry(-theta, idx)
    ws_mixer.rz(-2*beta, idx)
    ws_mixer.ry(theta, idx)

ws_mixer.draw(output='mpl')
[13]:
../_images/tutorials_10_warm_start_qaoa_22_0.png

The initial state and mixer operator can then be passed to QAOA.

[14]:
ws_qaoa_mes = QAOA(quantum_instance=quantum_instance,
                   initial_state=init_qc,
                   mixer=ws_mixer,
                   initial_point=[0., 1.])
[15]:
ws_qaoa = MinimumEigenOptimizer(ws_qaoa_mes)
[16]:
ws_qaoa_result = ws_qaoa.solve(qubo)
print(ws_qaoa_result)
optimal function value: 16.768932200000002
optimal value: [0. 0. 1. 1. 1. 0.]
status: SUCCESS

Analysis

Both results appear to give the same result. However, when we look at the underlying probability distribution we observe that the warm-start QAOA has a much higher probability of sampling the optimal solution.

[17]:
def format_qaoa_samples(samples, max_len: int = 10):
    qaoa_res = []
    for s in samples:
        if sum(s.x) == 3:
            qaoa_res.append((''.join([str(int(_)) for _ in s.x]), s.fval, s.probability))

    res = sorted(qaoa_res, key=lambda x: -x[1])[0:max_len]

    return [(_[0] + f': value: {_[1]:.3f}, probability: {1e2*_[2]:.1f}%') for _ in res]

format_qaoa_samples(qaoa_result.samples)
[17]:
['001110: value: 16.769, probability: 2.6%',
 '011010: value: 15.744, probability: 1.1%',
 '001011: value: 14.671, probability: 1.7%',
 '101010: value: 14.626, probability: 1.7%',
 '010110: value: 14.234, probability: 1.7%',
 '100110: value: 13.953, probability: 4.7%',
 '000111: value: 13.349, probability: 0.0%',
 '110010: value: 12.410, probability: 0.4%',
 '010011: value: 12.013, probability: 4.0%',
 '100011: value: 11.559, probability: 1.5%']
[18]:
format_qaoa_samples(ws_qaoa_result.samples)
[18]:
['001110: value: 16.769, probability: 86.3%',
 '001011: value: 14.671, probability: 0.3%',
 '101010: value: 14.626, probability: 0.8%',
 '100110: value: 13.953, probability: 1.7%',
 '000111: value: 13.349, probability: 0.1%',
 '100011: value: 11.559, probability: 0.0%']

Warm-start QAOA

The warm-start features above are available in the Qiskit optimization module as a single optimizer named WarmStartQAOAOptimizer which is illustrated below. This solver will solve a QUBO with a warm-start QAOA. It computes \(c^*\) by relaxing the problem. This behaviour is controlled by setting relax_for_pre_solver to True.

[19]:
from qiskit_optimization.algorithms import WarmStartQAOAOptimizer
[20]:
qaoa_mes = QAOA(quantum_instance=quantum_instance, initial_point=[0., 1.])
ws_qaoa = WarmStartQAOAOptimizer(pre_solver=CplexOptimizer(), relax_for_pre_solver=True,
                                 qaoa=qaoa_mes, epsilon=0.0)
[21]:
ws_result = ws_qaoa.solve(qubo)
ws_result
[21]:
optimal function value: 16.768932200000002
optimal value: [0. 0. 1. 1. 1. 0.]
status: SUCCESS
[22]:
format_qaoa_samples(ws_result.samples)
[22]:
['001110: value: 16.769, probability: 86.3%',
 '001011: value: 14.671, probability: 0.3%',
 '101010: value: 14.626, probability: 0.8%',
 '100110: value: 13.953, probability: 1.7%',
 '000111: value: 13.349, probability: 0.1%',
 '100011: value: 11.559, probability: 0.0%']
[23]:
import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit SoftwareVersion
QiskitNone
Terra0.17.0
Aer0.8.0
IgnisNone
AquaNone
IBM Q ProviderNone
System information
Python3.8.8 (default, Feb 19 2021, 19:42:00) [GCC 9.3.0]
OSLinux
CPUs2
Memory (Gb)6.791343688964844
Fri Apr 02 20:13:08 2021 UTC

This code is a part of Qiskit

© Copyright IBM 2017, 2021.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

[ ]: