Optimización con recocido cuántico

En este tutorial, exploraremos cómo usar el Recocido Cuántico para resolver un problema de optimización simple usando QiliSDK.

Nota

Si aún no lo has hecho, puede ser útil revisar primero estos tutoriales: Quantum Basics, Quantum Circuits y Quantum Annealing.

El problema

Supongamos que tenemos cuatro personas y necesitamos formar dos equipos de dos personas cada uno. Tenemos información sobre lo bien que trabajan juntos cada par de personas, y queremos formar los equipos de modo que la «compatibilidad» total de los equipos sea máxima. La información sobre quién se lleva bien con quién es la siguiente:

	Alice	Bob	Carol	Dave
Alice	N/A	1	3	4
Bob	1	N/A	5	2
Carol	3	5	N/A	6
Dave	4	2	6	N/A

Esto significa que Alice preferiría trabajar con Dave, mientras que Bob preferiría trabajar con Carol, y así sucesivamente.

Matemáticamente, podemos representar nuestros dos equipos como variables binarias, donde \(x_i\) es 0 si la persona \(i\) está en el equipo 0, y 1 si está en el equipo 1. La compatibilidad del equipo 0 puede expresarse entonces como:

\[C_0 = \sum_{i=0,j=0}^3 p_{ij} (1 - x_i) (1 - x_j)\]

Donde \(p_{ij}\) es la puntuación de compatibilidad entre la persona \(i\) y la persona \(j\). Nótese que aquí cada término solo es distinto de cero si tanto \(x_i\) como \(x_j\) son 0, lo que significa que ambas personas están en el equipo 0. Mientras tanto, la compatibilidad del equipo 1 puede expresarse como:

\[C_1 = \sum_{i=0,j=0}^3 p_{ij} x_i x_j\]

Por lo tanto, queremos maximizar la compatibilidad total \(C = C_0 + C_1\), sujeto a la restricción de que el equipo 1 (y por tanto el equipo 0) tenga exactamente dos personas:

\[\max_{x_i \in \{0, 1\}} C\]

\[\text{sujeto a} \quad x_{Alice} + x_{Bob} + x_{Carol} + x_{Dave} = 2\]

Un problema con esta formulación es que tiene una restricción, lo que la hace más difícil de resolver usando algoritmos de optimización cuántica. Para solucionar esto, podemos usar algunas técnicas para convertirla en un problema QUBO (Optimización Binaria Cuadrática Sin Restricciones). En este caso, dado que solo tenemos una restricción de igualdad, podemos añadir una versión al cuadrado de esta a la función objetivo como término de penalización, de modo que el término de penalización sea cero cuando se satisface la restricción y positivo en caso contrario. Para más detalles sobre este proceso, consulta la referencia QUBO. Ahora podemos asumir que el problema se ha reformulado en el siguiente problema de optimización sin restricciones:

\[\min_{x_i \in \{0, 1\}} \sum_{i=0,j=0}^4 c_{ij} x_i x_j\]

Donde \(c_{ij}\) son coeficientes que codifican la función objetivo así como la restricción.

La solución

Ahora que tenemos nuestro problema formulado como un QUBO, podemos usar el recocido cuántico para encontrar una solución aproximada. El recocido cuántico es un algoritmo de optimización metaheurístico inspirado en el proceso de recocido en metalurgia, donde un material se calienta y luego se enfría lentamente para eliminar defectos y encontrar un estado de baja energía. En el recocido cuántico, comenzamos con un Hamiltoniano simple cuyo estado fundamental es fácil de preparar y luego lo evolucionamos lentamente hacia un Hamiltoniano más complejo que codifica el problema de optimización que queremos resolver. La esperanza es que si hacemos esto con suficiente lentitud, el sistema permanecerá en su estado fundamental durante toda la evolución y así terminará en el estado fundamental del Hamiltoniano final, que corresponde a la solución óptima de nuestro problema.

Asumiendo que hemos reformulado nuestro problema en un QUBO, podemos construir el Hamiltoniano del problema para el recocido cuántico, que viene dado por:

\[H_{prob} = \sum_{i,j} c_{ij} Z(i) Z(j)\]

Donde \(Z(i)\) es el operador de Pauli-Z que actúa sobre el cúbit \(i\), y \(c_{ij}\) son los coeficientes de nuestra formulación QUBO. Tenga en cuenta que estas variables \(Z(i)\) están relacionadas con las variables binarias de nuestro problema original mediante la transformación \(x_i = (1 - Z(i))/2\).

Nuestro Hamiltoniano de mezcla se elige típicamente como el Hamiltoniano de campo transverso, que viene dado por:

\[H_{mix} = - \sum_i X(i)\]

El Hamiltoniano dependiente del tiempo que evolucionamos está dado entonces por:

\[H(t) = A(t) H_{mix} + B(t) H_{prob}\]

Donde \(A(t)\) y \(B(t)\) son funciones que determinan cómo interpolamos entre el Hamiltoniano de mezcla y el Hamiltoniano del problema a lo largo del tiempo. Por simplicidad, aquí asumiremos que realizamos una interpolación lineal, de modo que \(A(t) = 1 - t/T\) y \(B(t) = t/T\), donde \(T\) es el tiempo total de recocido.

La implementación

Para simular la evolución de este Hamiltoniano dependiente del tiempo, primero necesitamos formular nuestro problema:

from qilisdk.core.model import QUBO, ObjectiveSense
from qilisdk.core.comparison import EQ
from qilisdk.core.variables import BinaryVariable

num_people = 4
vars = [BinaryVariable(f"x{i}") for i in range(num_people)]
preferences = [[0, 1, 3, 4],
              [1, 0, 5, 2],
              [3, 5, 0, 6],
              [4, 2, 6, 0]]

model = QUBO("team_formation_example")

team_1 = sum(
    preferences[i][j] * vars[i] * vars[j]
    for i in range(num_people)
    for j in range(i+1, num_people)
)

team_0 = sum(
    preferences[i][j] * (1 - vars[i]) * (1 - vars[j])
    for i in range(num_people)
    for j in range(i+1, num_people)
)

model.set_objective(
    team_0 + team_1,
    label="obj",
    sense=ObjectiveSense.MAXIMIZE
)

model.add_constraint(
    "team_size_constraint",
    EQ(sum(vars[i] for i in range(num_people)), 2),
    lagrange_multiplier=10
)

Luego usamos el modelo para formar nuestro Hamiltoniano del problema, y también podemos definir nuestro Hamiltoniano de mezcla y el esquema de nuestra evolución:

from qilisdk.analog import X
from qilisdk.analog import Schedule
from qilisdk.core import QTensor

problem_hamiltonian = model.to_hamiltonian()
mixer_hamiltonian = sum(-1.0 * X(i) for i in range(num_people))

También necesitamos definir el esquema: cómo cambian los coeficientes de los Hamiltonianos del problema y de mezcla a lo largo del tiempo. En este caso, simplemente usaremos una mezcla lineal sencilla:

schedule = Schedule.linear(mixer_hamiltonian, problem_hamiltonian, total_time=100.0, dt=0.1)

Comenzaremos en el estado fundamental de nuestro Hamiltoniano de mezcla, que es el estado de superposición uniforme sobre todas las posibles asignaciones de equipos:

initial_state = QTensor.uniform(num_people)

Finalmente, inicializamos nuestro simulador cuántico, ejecutamos la evolución y leemos los resultados:

from qilisdk.functionals import AnalogEvolution
from qilisdk.readout import Readout
from qilisdk.backends import QiliSim

backend = QiliSim()
evolution = AnalogEvolution(schedule=schedule, initial_state=initial_state)
readout = Readout().with_sampling(1000).with_expectation([problem_hamiltonian])
results = backend.execute(evolution, readout)
print("Results:", results)

Esto imprimirá algo parecido a lo siguiente:

Functional Results: [

    Sampling Results: (
            nshots=1000,
            samples={'0110': 513, '1001': 487}
    )

    Expectation Value Results: (
            expectation_values=[-8.999696720325046],
    )

]

Como puede verse en estos resultados, las muestras satisfacen la restricción (es decir, tienen exactamente dos 1’s). Más concretamente, las dos muestras «0110» y «1001» (correspondientes a Alice/Dave en un equipo y Bob/Carol en el otro) son las soluciones óptimas a nuestro problema, siendo 9 la puntuación de compatibilidad máxima obtenible.

Para ayudar a visualizar la evolución que acabamos de realizar, podemos volver a ejecutar nuestra simulación, registrando esta vez el estado completo en cada paso temporal y utilizando después el método Schedule.draw_eigenvalues() para representar cómo evolucionó nuestro estado frente a los valores propios del Hamiltoniano a lo largo del tiempo:

evolution = AnalogEvolution(
    schedule=schedule,
    initial_state=initial_state,
    store_intermediate_results=True
)
readout = Readout().with_state_tomography()
results = backend.execute(evolution, readout)

schedule.draw_eigenvalues(
    intermediate_states=results.get_intermediate_states(),
    show_overlaps=True
)

Como podemos ver en esta gráfica, el estado comenzó en el estado fundamental del Hamiltoniano de mezcla (la línea inferior, en el extremo izquierdo) y después evolucionó hasta el estado fundamental del Hamiltoniano del problema (en el extremo derecho). Al avanzar con suficiente lentitud nos mantuvimos en el estado fundamental, con solo un pequeño porcentaje de nuestro estado escapándose hacia algunos de los estados excitados. Un tiempo total mayor (y por tanto una evolución más lenta) haría que se escapara menos estado, a costa de tardar más en ejecutarse.

Lectura adicional

• QUBO

• Quantum Annealing