Note

This page was generated from docs/tutorials/10_warm_start_qaoa.ipynb.

# 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
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)
mdl.add_constraint(cost == total)
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
%qiskit_copyright
```

### Version Information

Qiskit Software | Version |
---|---|

`qiskit-terra` | 0.18.3 |

`qiskit-aer` | 0.9.0 |

`qiskit-optimization` | 0.2.3 |

System information | |

Python | 3.8.12 (default, Sep 13 2021, 08:28:12) [GCC 9.3.0] |

OS | Linux |

CPUs | 2 |

Memory (Gb) | 6.790924072265625 |

Wed Oct 06 16:50:44 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.

```
[ ]:
```

```
```