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.

To help to visualize the evolution we just performed, we can re-run our simulation, this time recording the full state at each timestep and then using the Schedule.draw_eigenvalues() method to plot how our state evolved versus the eigenvalues of the Hamiltonian over time:

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
)
../../_images/annealing_example_schedule_eigenvalues.svg ../../_images/annealing_example_schedule_eigenvalues_dark.svg

As we can see from this plot, the state started in the ground state of the mixing Hamiltonian (the bottom line, at the very left) and then evolved to the ground state of the problem Hamiltonian (at the very right). By moving slow enough we stayed in the ground state, with only a small percentage of our state leaking out to some of the excited states. A larger total time (and thus a slower evolution) would result in less state leaking out, at the cost of taking longer to run.

Lectura adicional