Note

# 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]:

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.translators import from_docplex_mp


## 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)

qp = from_docplex_mp(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

Minimize
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.012055025682685
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.17524995761801482,
1.480388816398847e-07,
0.9709053264087595,
0.7384168677494172,
0.9999999916475085,
0.14438904470168706]


## 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]:


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]:


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.5%',
'011010: value: 15.744, probability: 2.2%',
'001011: value: 14.671, probability: 1.4%',
'101010: value: 14.626, probability: 1.2%',
'010110: value: 14.234, probability: 0.0%',
'100110: value: 13.953, probability: 1.8%',
'000111: value: 13.349, probability: 2.6%',
'110010: value: 12.410, probability: 2.1%',
'010011: value: 12.013, probability: 0.4%',
'100011: value: 11.559, probability: 1.7%']

[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 behavior 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: 48.4%',
'001011: value: 14.671, probability: 3.4%',
'101010: value: 14.626, probability: 4.5%',
'100110: value: 13.953, probability: 0.6%',
'000111: value: 13.349, probability: 0.4%',
'100011: value: 11.559, probability: 0.1%']

[23]:

import qiskit.tools.jupyter
%qiskit_version_table


### Version Information

Qiskit SoftwareVersion
qiskit-terra0.18.3
qiskit-aer0.9.0
qiskit-optimization0.2.3
System information
Python3.8.12 (default, Sep 13 2021, 08:28:12) [GCC 9.3.0]
OSLinux
CPUs2
Memory (Gb)6.790924072265625
Wed Oct 06 16:50:44 2021 UTC

### This code is a part of Qiskit

[ ]: