Nota

Esta página fue generada a partir de docs/tutorials/03_minimum_eigen_optimizer.ipynb.

Optimizador Propio Mínimo

Introducción

Una clase interesante de problemas de optimización que debe abordar la computación cuántica son los problemas de Optimización Binaria Cuadrática sin Restricciones (Quadratic Unconstrained Binary Optimization, QUBO). Encontrar la solución a un QUBO es equivalente a encontrar el estado fundamental de un Hamiltoniano de Ising correspondiente, que es un problema importante no solo en la optimización, sino también en química cuántica y física. Para esta traducción, las variables binarias que toman valores en \(\{0, 1\}\) se reemplazan por variables de espín que toman valores en \(\{-1, +1\}\), lo que permite reemplazar las variables de espín resultantes por matrices Pauli Z y, por lo tanto, un Hamiltoniano de Ising. Para más detalles sobre este mapeo nos referimos a [1].

Qiskit proporciona conversión automática de un QuadraticProgram adecuado a un Hamiltoniano de Ising, que luego permite aprovechar todas las implementaciones de SamplingMinimumEigensolver, como

• SamplingVQE,

• QAOA, o

• NumpyMinimumEigensolver (método clásico exacto).

Nota 1: MinimumEigenOptimizer no es compatible con qiskit.algorithms.minimum_eigensolver.VQE. Pero en su lugar se puede usar qiskit.algorithms.minimum_eigensolver.SamplingVQE.

Nota 2: MinimumEigenOptimizer puede usar NumpyMinimumEigensolver como un caso de excepción, aunque hereda de MinimumEigensolver (no de SamplingMinimumEigensolver).

Qiskit Optimization proporciona una clase MinimumEigenOptimizer, que envuelve la traducción a un Hamiltoniano de Ising (en Qiskit Terra también llamado Operator), la llamada a un MinimumEigensolver y la traducción de los resultados de vuelta a un OptimizationResult.

A continuación, primero ilustramos la conversión de un QuadraticProgram a un Operator y luego mostramos cómo usar el MinimumEigenOptimizer con diferentes MinimumEigensolvers para resolver un QuadraticProgram dado. Los algoritmos en Qiskit intentan convertir automáticamente un problema determinado a la clase de problema soportada si es posible, por ejemplo, el MinimumEigenOptimizer traducirá automáticamente las variables enteras a variables binarias o agregará restricciones de igualdad lineal como un término de penalización cuadrática al objetivo. Se debe mencionar que se generará un QiskitOptimizationError si se intenta la conversión de un programa cuadrático con variables enteras.

La profundidad del circuito de QAOA podría ser incrementada con el tamaño del problema, lo que podría ser prohibitivo para los dispositivos cuánticos a corto plazo. Un posible método alternativo es el QAOA recursivo, como se introduce en la referencia [2]. Qiskit generaliza este concepto a RecursiveMinimumEigenOptimizer, el cual se introduce al final de este tutorial.

Referencias

[1] A. Lucas, Ising formulations of many NP problems, Front. Phys., 12 (2014).

[2] S. Bravyi, A. Kliesch, R. Koenig, E. Tang, Obstacles to State Preparation and Variational Optimization from Symmetry Protection, arXiv preprint arXiv:1910.08980 (2019).

Convertir un QUBO a un Operator

from qiskit.utils import algorithm_globals
from qiskit.algorithms.minimum_eigensolvers import QAOA, NumPyMinimumEigensolver
from qiskit.algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler
from qiskit_optimization.algorithms import (
    MinimumEigenOptimizer,
    RecursiveMinimumEigenOptimizer,
    SolutionSample,
    OptimizationResultStatus,
)
from qiskit_optimization import QuadraticProgram
from qiskit.visualization import plot_histogram
from typing import List, Tuple
import numpy as np

# create a QUBO
qubo = QuadraticProgram()
qubo.binary_var("x")
qubo.binary_var("y")
qubo.binary_var("z")
qubo.minimize(linear=[1, -2, 3], quadratic={("x", "y"): 1, ("x", "z"): -1, ("y", "z"): 2})
print(qubo.prettyprint())

Problem name:

Minimize
  x*y - x*z + 2*y*z + x - 2*y + 3*z

Subject to
  No constraints

  Binary variables (3)
    x y z

A continuación traducimos este QUBO en un operador de Ising. Esto no solo resulta en un Operator, sino también en un desplazamiento constante que se debe tener en cuenta para cambiar el valor resultante.

op, offset = qubo.to_ising()
print("offset: {}".format(offset))
print("operator:")
print(op)

offset: 1.5
operator:
-0.5 * IIZ
+ 0.25 * IZI
- 1.75 * ZII
+ 0.25 * IZZ
- 0.25 * ZIZ
+ 0.5 * ZZI

A veces, un QuadraticProgram también podría ser dado directamente en forma de Operator. Para tales casos, Qiskit también proporciona un traductor de un Operator de vuelta a un QuadraticProgram, que ilustramos a continuación.

qp = QuadraticProgram()
qp.from_ising(op, offset, linear=True)
print(qp.prettyprint())

Problem name:

Minimize
  x0*x1 - x0*x2 + 2*x1*x2 + x0 - 2*x1 + 3*x2

Subject to
  No constraints

  Binary variables (3)
    x0 x1 x2

Este traductor permite, por ejemplo, traducir un Operator a un QuadraticProgram y luego resolver el problema con otros algoritmos que no se basan en la representación del Hamiltoniano de Ising, como el GroverOptimizer.

Resolver un QUBO con el MinimumEigenOptimizer

Empezamos inicializando el MinimumEigensolver que queremos usar.

algorithm_globals.random_seed = 10598
qaoa_mes = QAOA(sampler=Sampler(), optimizer=COBYLA(), initial_point=[0.0, 0.0])
exact_mes = NumPyMinimumEigensolver()

Luego, usamos el MinimumEigensolver para crear un MinimumEigenOptimizer.

qaoa = MinimumEigenOptimizer(qaoa_mes)  # using QAOA
exact = MinimumEigenOptimizer(exact_mes)  # using the exact classical numpy minimum eigen solver

Primero usamos el MinimumEigenOptimizer basado en el NumPyMinimumEigensolver exacto clásico para obtener la solución de referencia óptima para este pequeño ejemplo.

exact_result = exact.solve(qubo)
print(exact_result.prettyprint())

objective function value: -2.0
variable values: x=0.0, y=1.0, z=0.0
status: SUCCESS

A continuación aplicamos el MinimumEigenOptimizer basado en QAOA al mismo problema.

qaoa_result = qaoa.solve(qubo)
print(qaoa_result.prettyprint())

objective function value: -2.0
variable values: x=0.0, y=1.0, z=0.0
status: SUCCESS

Análisis de Muestras

OptimizationResult proporciona información útil en forma de SolutionSamples (aquí indicado como muestras). Cada SolutionSample contiene información sobre los valores de entrada (x), el valor de la función objetivo correspondiente (fval), la fracción de muestras correspondiente a esa entrada (probability) y la solución status (SUCCESS, FAILURE, INFEASIBLE). Múltiples muestras correspondientes a la misma entrada se consolidan en una sola SolutionSample (siendo su atributo probability la fracción agregada de muestras representada por esa SolutionSample).

print("variable order:", [var.name for var in qaoa_result.variables])
for s in qaoa_result.samples:
    print(s)

variable order: ['x', 'y', 'z']
SolutionSample(x=array([0., 1., 0.]), fval=-2.0, probability=0.4410306097905226, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([0., 0., 0.]), fval=0.0, probability=0.22763693649873265, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([1., 1., 0.]), fval=0.0, probability=0.14136368254300133, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([1., 0., 0.]), fval=1.0, probability=0.12574358779906872, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([0., 0., 1.]), fval=3.0, probability=0.020510231887331747, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([1., 0., 1.]), fval=3.0, probability=0.030444770051099967, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([0., 1., 1.]), fval=3.0, probability=0.012349858838771063, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([1., 1., 1.]), fval=4.0, probability=0.0009203225914718031, status=<OptimizationResultStatus.SUCCESS: 0>)

También es posible que deseemos filtrar muestras de acuerdo con su estado o probabilidades.

def get_filtered_samples(
    samples: List[SolutionSample],
    threshold: float = 0,
    allowed_status: Tuple[OptimizationResultStatus] = (OptimizationResultStatus.SUCCESS,),
):
    res = []
    for s in samples:
        if s.status in allowed_status and s.probability > threshold:
            res.append(s)

    return res

filtered_samples = get_filtered_samples(
    qaoa_result.samples, threshold=0.005, allowed_status=(OptimizationResultStatus.SUCCESS,)
)
for s in filtered_samples:
    print(s)

SolutionSample(x=array([0., 1., 0.]), fval=-2.0, probability=0.4410306097905226, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([0., 0., 0.]), fval=0.0, probability=0.22763693649873265, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([1., 1., 0.]), fval=0.0, probability=0.14136368254300133, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([1., 0., 0.]), fval=1.0, probability=0.12574358779906872, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([0., 0., 1.]), fval=3.0, probability=0.020510231887331747, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([1., 0., 1.]), fval=3.0, probability=0.030444770051099967, status=<OptimizationResultStatus.SUCCESS: 0>)
SolutionSample(x=array([0., 1., 1.]), fval=3.0, probability=0.012349858838771063, status=<OptimizationResultStatus.SUCCESS: 0>)

Si queremos obtener una mejor perspectiva de los resultados, la estadística es de gran ayuda, tanto con respecto a los valores de la función objetivo como a sus respectivas probabilidades. Por lo tanto, la media y la desviación estándar son los elementos básicos para comprender los resultados.

fvals = [s.fval for s in qaoa_result.samples]
probabilities = [s.probability for s in qaoa_result.samples]

np.mean(fvals)

1.5

np.std(fvals)

1.9364916731037085

Finalmente, a pesar de todo el procesamiento numérico, la visualización suele ser el mejor enfoque de análisis temprano.

samples_for_plot = {
    " ".join(f"{qaoa_result.variables[i].name}={int(v)}" for i, v in enumerate(s.x)): s.probability
    for s in filtered_samples
}
samples_for_plot

{'x=0 y=1 z=0': 0.4410306097905226,
 'x=0 y=0 z=0': 0.22763693649873265,
 'x=1 y=1 z=0': 0.14136368254300133,
 'x=1 y=0 z=0': 0.12574358779906872,
 'x=0 y=0 z=1': 0.020510231887331747,
 'x=1 y=0 z=1': 0.030444770051099967,
 'x=0 y=1 z=1': 0.012349858838771063}

plot_histogram(samples_for_plot)

RecursiveMinimumEigenOptimizer

El RecursiveMinimumEigenOptimizer toma un MinimumEigenOptimizer como entrada y aplica el esquema de optimización recursiva para reducir el tamaño del problema una variable a la vez. Una vez que el tamaño del problema intermedio generado está por debajo de un umbral dado (min_num_vars), el RecursiveMinimumEigenOptimizer usa otro solucionador (min_num_vars_optimizer), por ejemplo, un solucionador clásico exacto como CPLEX o el MinimumEigenOptimizer basado en el NumPyMinimumEigensolver.

A continuación, mostramos cómo usar el RecursiveMinimumEigenOptimizer usando los dos MinimumEigenOptimizers presentados anteriormente.

En primer lugar, construimos el RecursiveMinimumEigenOptimizer tal que reduce el tamaño del problema de 3 variables a 1 variable y luego utiliza el solucionador exacto para la última variable. Entonces llamamos a solve para optimizar el problema considerado.

rqaoa = RecursiveMinimumEigenOptimizer(qaoa, min_num_vars=1, min_num_vars_optimizer=exact)

rqaoa_result = rqaoa.solve(qubo)
print(rqaoa_result.prettyprint())

objective function value: -2.0
variable values: x=0.0, y=1.0, z=0.0
status: SUCCESS

filtered_samples = get_filtered_samples(
    rqaoa_result.samples, threshold=0.005, allowed_status=(OptimizationResultStatus.SUCCESS,)
)

samples_for_plot = {
    " ".join(f"{rqaoa_result.variables[i].name}={int(v)}" for i, v in enumerate(s.x)): s.probability
    for s in filtered_samples
}
samples_for_plot

{'x=0 y=1 z=0': 1.0}

plot_histogram(samples_for_plot)

import qiskit.tools.jupyter

%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
qiskit-terra	0.23.0
qiskit-aer	0.11.1
qiskit-optimization	0.5.0
qiskit-machine-learning	0.6.0
System information
Python version	3.9.15
Python compiler	Clang 14.0.0 (clang-1400.0.29.102)
Python build	main, Oct 11 2022 22:27:25
OS	Darwin
CPUs	4
Memory (Gb)	16.0
Mon Dec 05 22:42:36 2022 JST

This code is a part of Qiskit

© Copyright IBM 2017, 2022.

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.