7. Introducción.#
En este apartado vamos a desarrollar las denominadas funciones de penalización que junto con el valor comunmente llamado \(\gamma\) o coeficiente de lagrange son las que nos van a permitir introducir en los programas tipo QUBO o Ising las restricciones del problema planteado.
De una forma muy resumida, se puede decir que una función de penalización es aquella función que devuelve sus valores mínimos cuando el valor de las variables cumplen con las condiciones de la restricción planteada. Por lo tanto el valor de esta función es la penalización que añadiremos cuando no se cumplen las restricciones planteadas en el problema, haciendo por lo tanto que la solución para estos últimos casos sea menos óptima.
La idea de las funciones de penalización es la siguiente: supongamos que queremos minimizar una función \(f\) sujeta a cierta restricciones, es decir se tiene el siguiente problema:
Entonces podemos intentar transformar \(f\) en otra función \(g\) de manera que los valores que no cumplen la restricción tomarán valores “grandes” en \(g\). Entonces de una forma mucho más perfilada, se puede definir la función de penalización de la siguiente manera:
Definición: Sean \(P\) (esta va a ser la función de penalización) y \(f\) funciones de \(\{0,1\}^N \) en \(\mathbb{R}\), r una restricción e \(y=argmin\{f(z):z\,cumple\,la\,restricción\,r \}\), es decir, y es el punto donde se alcanza el mínimo de \(f\) pero restringido a los valores que cumplen con la restricción r. Diremos que \(P\) es una función de penalización para la restricción r si para los valores de x que no satisfacen r, se cumple que \(f(x)+p(x)>f(y)\) y además \(p(z)=0\) para los valores z que cumplen con la restricción r.
De esta manera, añadiendo mediante suma y multiplicada por el coeficiente de lagrange dicha función de penalización a la función \(f\) de coste, se tendrá que las mejores soluciones serán aquellas que cumplen con la restricción.
7.1. Construcción de una función de penalización.#
En este apartado vamos a ver un ejemplo de cómo podemos construir una función de penalización, para ello supongamos que tenemos la siguiente restricción: \(x_1+x_2\le 1\).
Aunque la anterior restricción es una restricción de desigualdad y para ello existe un método general de construcción, como se verá en un apartado posterior titulado restricciones de desigualdad, aquí vamos a ver cómo construir esta función de penalización de una forma diferente y más eficiente, ya que de esta manera se evita el incremento del número de variables.
La idea que subyace en este tipo de construcción es otorgar el valor de cero si se cumple la restricción y 1 si no se cumple. Es decir mínima energía si cumplimos restricción y mayor energía si no se cumple. Entonces para ello podemos construir para este ejemplo la siguiente tabla de la verdad que contiene el cruce de cada uno de los valores de las variable \(x_1\) y \(x_2\), y en la tercer fila el valor 0 ó 1 dependiendo si se cumple o no la restricción.
x1 |
x2 |
x1+x2<=1? |
|---|---|---|
0 |
0 |
0 |
0 |
1 |
0 |
1 |
0 |
0 |
1 |
1 |
1 |
La función de penalización que vamos buscando tiene el siguiente formato:
Entonces vamos a calcular los valores de a,b,c y d de tal manera que al sustituir \(x_1\) y \(x_2\) por los correspondientes valores de las dos primeras columnas de cada fila de la tabla de la verdad anterior, se obtenga el valor de la tercera columna de dicha tabla de la verdad. Es decir de esta manera vamos a tener un sistema fácil de resolver de cuatro ecuaciones con cuatro incógnitas.
Si tomamos el par de valores de x1 y x2 (0,0) se tendrá que d=0
Si tomamos el par de valores de x1 y x2 (0,1) se tendrá que b=0
Si tomamos el par de valores de x1 y x2 (1,0) se tendrá que a=0
Si tomamos el par de valores de x1 y x2 (1,1) se tendrá que c=1
En consecuencia, en este caso la función de penalización será \(x_1\cdot x_2\)
En la siguiente tablas se muestran las funciones de penalización utilizadas con mayor frecuencia, y su deducción se puede obtener fácilmente siguiendo los pasos marcados en el párrafo anterior:
Restricción |
penalty function |
|---|---|
x+y<=1 |
x.y |
x+y>=1 |
1-x-y+x.y |
x+y=1 |
1-x-y+2xy |
x<=y |
x-xy |
x+y+z<=1 |
xy+xz+yz |
x=y |
x+y-2xy |
A continuación se muestra una función que permite calcular funciones de penalización de una forma rápida y cómoda
import itertools
import numpy as np
from math import factorial
from numpy.linalg import matrix_rank
def penalties(expresion, n = 2):
"""
Con esta función se van a obtener los coeficientes de una ecuación de penatys, utilizada en QUBO
Parámetros:
expresión: es la expresión que se debe avaluar para definir los penaltys
n: Puede ser 2 ´0 3. Indica el número de variables a tener en cuenta
"""
# Definimos esta función, para obtener el factorial de un número, es decir Cn,m
def comb(n,m):
resul=factorial(n)/(factorial(m)*factorial(n-m))
return resul
# definimos el número de combinaciones que se pueden dar
x1 = [0,1]
A = np.array([p for p in itertools.product(x1, repeat=n)])
# Calculamos los términos independientes
independiente= np.zeros(A.shape[0])
for h in range(len(independiente)):
if n == 2:
x = A[h,0]
y = A[h,1]
independiente[h]=1-eval(expresion)
elif n == 3:
x = A[h,0]
y = A[h,1]
z = A[h,2]
independiente[h]=1-eval(expresion)
# Calculamos los productos cruzados o cuadráticos,
#es decir xy, xz, yz
cruzados = np.zeros((A.shape[0],int(comb(A.shape[1],2))))
for i in range(A.shape[0]):
if n == 2:
combinaciones = list(itertools.combinations([0,1],2))
cruzados[i] = A[i,combinaciones[0][0]]*A[i,combinaciones[0][1]]
if n == 3:
combinaciones = list(itertools.combinations([0,1,2],2))
for j in range(len(combinaciones)):
cruzados[i,j] = A[i,combinaciones[j][0]]*A[i,combinaciones[j][1]]
# Ahora obtenemos la matriz B de la Expresión Bx=y
B = np.append(A, cruzados, axis=1)
B = np.append(B, np.ones((A.shape[0],1)),axis=1)
if n == 3:
# borramos la última fila para resolver el sistema, ya que es c.l. del resto
B = np.delete(B, B.shape[0]-1, 0)
independiente = independiente[:-1]
# obtenemos el resultado
resul = np.linalg.solve(B, independiente)
if n==2:
return "{}*x+({}*y)+({}*xy)+({})".format(resul[0],resul[1],resul[2],resul[3])
elif n==3:
return "{}*x+({}*y)+({}*z)+({}*xy)+({}*xz)+({}*yz)+({})".format(resul[0],resul[1],resul[2],resul[3],resul[4],resul[5],resul[6])
La utilización de la función anterior es muy fácil y a continuación se muestran algunos ejemplos, con la finalidad de que sirvan de ejemplo para entender el uso de la misma.
penalties("x + y == z", n = 3)
'1.0*x+(1.0*y)+(1.0*z)+(-1.0*xy)+(-2.0*xz)+(-2.0*yz)+(0.0)'
penalties("x+y==2",n=2)
'0.0*x+(0.0*y)+(-1.0*xy)+(1.0)'
# Esta es la puerta NOT
penalties("x!=y",n=2)
'-1.0*x+(-1.0*y)+(2.0*xy)+(1.0)'
Dependiendo del tipo de restricción existen también otros procedimientos de calcular estas funciones de penalización que se pasan a desarrollar en los siguientes apartados.
7.2. Restricciones de igualdad#
Supongamos que tenemos una restrcción de igualdad de la forma \(a=p(x_1,...,x_n)\), siendo p un polinomio de grado 1 en las variables \(x_1,...,x_n\). Se cumple entonces que la función \(f=:\left(a-p(x_{1},...,x-n)\right)^{2}\) toma su valor mínimo 0 en los puntos para los cuales p=a. De esta manera, encontrar \(x_1,...,x_n\) tales que \(p(x_1,...,x_n)=a\) coincide con encontrar los valores \(x_1,...,x_n\) donde \(f\) alcanza el mínimo. Por lo tanto la función \(f\) anterior se corresponde con una función de penalización para la restricción lineal \(a=p(x_1,...,x_n)\). Una vez que se consigue esto, sólo queda multiplicar la función de penalización \(f\) por un valor \(\gamma\) que llamaremos coeficiente de lagrange, de manera que se garantice que las mejores soluciones son aquellas que satisfacen la restricción lineal.
En consecuencia, supongamos que queremos minimizar la función \(g(x_1,...,x_n)\) sujeto a la condición lineal de que \(f(x_1,...x_n)=0\), siendo f un polinomio de grado 1. Entonces se tomaría como función de coste la siguiente función:
donde \(\gamma\) es el denominado coeficiente de Lagrange. Si \(\gamma\) tiene un valor suficientemente alto, tendremos entonces que cuando no se satisfaga la restricción de la función \(f\), la función de coste tomará valores más altos que en aquellos puntos donde sí se cumpla la restricción correspondiente a la función \(f\).
Se darán más detalles sobre un correcto ajuste de los coeficientes de Lagrange en la sección denominada coeficientes de Lagrange un poco más adelante.
Una vez que hemos explicado cómo construir funciones de penalización que garanticen que se cumplen restricciones lineales de igualdad, explicaremos como garantizar que se cumplen las restricciones de desigualdad.
7.3. Restricciones de desigualdad.#
Otras de las posibles restricciones que pueden aparecer al modelizar un problema en un formato de tipo QUBO son las restricciones de desigualdad. En estos casos se convierten en restricciones de igualdad introduciendo variables de holgura (o slack variables).
Como en las modelizaciones QUBO las variables deben ser binarias, estas variables deben introducirse expresándose de manera binaria. Veámoslo con el siguiente ejemplo.
Supongamos que tenemos la restricción \(f(x_1,...,x_n)\le b\) y sea m el mínimo de \(f\). Calculemos el número \(na:=(log_2(b-m))+1\). Podemos convertir entonces la restricción de desigualdad anterior en otra de igualdad de la siguiente manera: \(f(x_1,...,x_n)+\sum_{j=0}^{na}2^j\cdot a_j=b\).
Este método tiene el inconveniente de que aumenta el número de variables necesarias en el modelo por cada restricción de desigualdad, añadiendo gran cantidad de cúbits de cara a su programación.
La solución óptima consiste en encontrar una función de penalización \(P\) que tome el mismo valor cuando se cumpla la restricción, y valores grandes cuando no. De esta manera evitaríamos introducir variables auxiliares extra. Sin embargo, encontrar estas funciones de penalización suele ser en ocasiones un trabajo bastante complicado.
No obstante lo anterior, vamos a desarrollar con mayor detalle, como pasar de una restricción de desigualdad a otra de igualdad, en los casos en los que la restricción tiene el siguiente formato (que suele ser muy común en este tipo de problemas):
Para simplificar los siguientes cálculos, denotamos \(Y_i=m_ix_i\), quedando la expresión anterior de la siguiente manera:
Transformamos la inecuación anterior en una igualdad, de la siguiente manera:
Y ahora pasamos todo a un miembro y elevamos al cuadrado:
La variable s es una variable de holgura que nos permite pasar de una inecuación a una ecuación. Desarrollando la expresión anterior se tiene lo siguiente:
Finalmente y desarrollando el primer sumatorio se tiene la siguiente expresión:
Rehaciendo la expresión \(Y_i=m_ix_i\), introducimos de nuevo nuestras variables. Pero en vez de usar la notación \(x_i\) usaremos \(q_i\), denotando que nuestras variables ya son cubits. Por lo tanto quedaría lo siguiente.
Y ahora simplificando sumatorios, se obtiene lo siguiente:
Para terminar el proceso, debemos representar la variable de holgura \(s\) como una variable compuesta por cubits. Estas variables auxiliares tienen la propiedad de sólo poder obtener valores 0 ó 1 por lo que se hace la siguiente igualdad
Donde \(s_k\) es cada cubit que representa la variable auxiliar. De esta forma introducimos las variables en las ecuaciones y minimizamos el so de cubits a los estrictamente necesarios. Sustituyendo esto último en la ecuación anterior y simplificando los sumatorios resultantes queda la ecuación final que buscamos:
De la expresión anterior, modificando los valores de n, el vector de valores \(m\) y el valor \(d\), podemos obtener la función de penalización de muchas de las restricciones que implican desigualdad, ahorrando de esta forma tiempo en el cálculo de cada función y pudiendo implementar una función general que genere estas penalizaciones o restricciones de desigualdad.
7.4. Problema de las N-reinas.#
Introducimos en este apartado, el problema de las N-reinas pues a él nos referiremos en el apartado posterior. Este problema consiste en lo siguiente:
Dado un tablero de ajedrez de tamaño NxN y disponer N reinas de tal manera que las mismas no se puedan atacar, es decir dispuesta una reina en una casilla del tablero, no puede haber otra en una casilla de la horizontal, ni en la vertical ni en diagonal.
La explicación del problema la podemos ver en el siguiente vídeo: https://youtu.be/ZGKSoTUaphU
La explicación del código se puede ver en : https://youtu.be/3wVwpPzoyDs
Nuestra función a minimizar será:
\(\begin{equation} f(x) = -\sum_i^n \sum_j^n x_{i,j} + \lambda\sum_{i_1}^n \sum_{i_2}^n \sum_{i_3}^n \sum_{i_4}^n J_{i_1,i_2,i_3,i_4}x_{i_1,i_2}x_{i_3,i_4} \end{equation}\)
donde \(x_{i,j}\) valdrá 1 si en la casilla (i,j) hay una dama y 0 si no.
Daremos además a \(J_{i_1,i_2,i_3,i_4}\) el valor 1 si las casillas \((i_1,i_2)\) y \((i_3, i_4)\) están conectadas y 0 en otro caso.
qubovert lo encontramos en https://qubovert.readthedocs.io/en/latest/.
El código que se utiliza para resolver el problema lo ponemos a continuación.
import qubovert
#size = 4
size = 6
lagrange = size ** 2
# Creamos las variables de nuestro modelo
Q = qubovert.QUBO()
for i in range(size):
for j in range(size):
Q.create_var(f"x_{i}_{j}") # Creamos las variables x_i y x_j
# Añadimos el primer bloque de la función objetivo
for i in range(size):
for j in range(size):
Q[(f"x_{i}_{j}",)] = -1
# Incluimos las restricciones finales
for i1 in range(size):
for i2 in range(size):
for i3 in range(size):
for i4 in range(size):
if i1 == i3 or i2 == i4 or i1 - i3 == i2 - i4 or i1 - i3 == i4 - i2:
if not (i1 == i3 and i2 == i4):
Q[(f"x_{i1}_{i2}", f"x_{i3}_{i4}")] = lagrange
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[5], line 1
----> 1 import qubovert
3 #size = 4
4 size = 6
ModuleNotFoundError: No module named 'qubovert'
Ahora ya tenemos definido nuestro qubo. Vale -1 para los coeficientes individuales y gamma = 16 (\(=2^N=2^4=16\)) cuando los cuadrados están conectados
# Lo hacemos para que dwave lo admita
dwave_dic = {}
for i in Q:
if len(i) == 1:
dwave_dic[(i[0],i[0])] = Q[i]
else:
dwave_dic[i] = Q[i]
from neal import SimulatedAnnealingSampler
#from dwave.system import DWaveSampler, EmbeddingComposite
n_samples = 2000 # número de veces que ejecutamos el sistema
sampler = SimulatedAnnealingSampler()
#sampler = EmbeddingComposite(DWaveSampler())
sampleset = sampler.sample_qubo(dwave_dic, num_reads = n_samples)
solution = sampleset.first.sample
for i in range(size):
for j in range(size):
if solution[f"x_{i}_{j}"] == 0:
print("O", end = " ")
else:
print("X", end = " ")
print()
7.5. Coeficientes de Lagrange#
Una de las tareas más importantes para realizar una buena modelización QUBO es el ajuste de los coeficientes de Lagrange. En muchos artículos que abordan problemas mediante el algoritmo Quantum Annealing se proponen diversas modelizaciones sin hacer referencia sobre los valores que deben tener los coeficientes de Lagrange.
Una buena modelización con unos malos coeficientes de Lagrange provocará que el algoritmo de atemperamiento cuántico no arroje buenos resultados. En muchas modelizaciones el ajuste de estos coeficientes resulta ser una tarea complicada y delicada.
Si los coeficientes de Lagrange no son lo suficientemente grandes, aparecerán soluciones que no cumplen ciertas restricciones.
Si por el contrario toman valores demasiado grandes, el algoritmo se atascará en soluciones que garanticen que se cumplen ciertas restricciones y no conseguirá buenos valores de la función objetivo o en otras restricciones.
Podemos apreciar este hecho con el código correspondiente al problema de las N-reinas que encontramos en el anterior apartado. Si tomamos en el código lagrange = size ^2 (es decir, la variable de nombre lagrange se obtiene al elevar al cuadrado la variable de nombre size, que son las reinas a tener en cuenta), tendremos que el algoritmo deja de encontrar soluciones para el problema a partir de \(N \ge 11\).
Si por el contrario tomamos en el código lagrange = 2, el algoritmo proporciona soluciones para valores de N tan grandes como 100. El caso óptimo en el que el algoritmo de Quantum Annealing proporciona resultados realmente buenos se produce cuando el óptimo está suficientemente alejado del resto de valores y no existen puntos para los cuales la función toma valores grandes muy alejados del resto.
Vamos a tratar de mostrar algunos resultados que nos ayuden a estimar los coeficientes de Lagrange. El problema de estimar los coeficientes de Lagrange no es un problema novedoso del Quantum Annealing. En el artículo “Penalty and partitioning techniques to improve performance of QUBO solvers” , se realiza un estudio sobre cómo estimar estos coeficientes cuando se está realizando una modelización QUBO previamente a aplicar annealing clásico.
Dicho artículo propone unas condiciones que ha de cumplir el coeficiente de Lagrange cuando el algoritmo heurístico con el que se va a resolver el problema es de tipo 1-flip. Los algoritmos heurísticos de tipo 1-flip son aquellos en los que, para pasar de una posible solución a una nueva, simplemente se produce un cambio en una de las variables que forma el modelo. Esta condición conduce a un caso más sencillo para el cual el coeficiente de Lagrange debe cumplir condiciones menos restrictivas. Sin embargo, algunos de los resultados de los que se habla son similares a los que vamos a manejar ahora.
Tratemos de generalizar lo expuesto en el artículo citado a estimar el valor que debe tomar \(\lambda\) para cada penalización en el caso del Quantum Annealing.
Supongamos que tenemos una función de coste \(h(x_1,...,x_n)=g(x_1,...,x_n)+\lambda f(x_1,...,x_n)\), donde \(g\) es la función que queremos minimizar y \(f\) representa una restricción que queremos que se satisfaga.
Sea \(y\in \{0,1\}^n\) y \(x\in \{y:f(y)=0 \}\), es decir los puntos donde se hace cero la función que representa la restricción. Supongamos que \(x_0\) es el valor que minimiza \(g\) y además cumple que \(f(x_0)=0\), es decir seria la solución de nuestro problema. Entonces necesitaremos que se cumpla para todo y lo siguiente:
Si tomamos como \(p=min\{f(y):f(y)>0, y\in \{0,1\}^n\}\) tendremos que para todo y:
por lo que si garantizamos que \(min\, g(y)+\lambda p >g(x_0)\) tendremos que se cumplirá lo que buscamos. Como en muchos casos calcular \(g(x_0)\) será equivalente a resolver el problema, es suficiente con tomar
siendo x uno de los puntos donde se cumple que f(x)=0.
Todos los mínimos que aparecen en el párrafo anterior se corresponden con el mínimo de cierta función en un conjunto de valores finito, por lo que existen, son finitos y se cumple que \(p \ne 0\).
Si no fuese sencillo estimar un valor de g(x), podría sustituirse por \(\underset{y}{max}\,g(y)\).
Estimar de manera exacta el valor p se correspondería con resolver el denominado problema de suma 0, cuya complejidad computacional es \(2^n\) . Sin embargo, si la precisión de los coeficientes de f se corresponde con \(10^{-m}\), podría simplemente aproximarse p por el valor \(a10^{-m}\) con \(a\in(0,1)\).
A pesar de esto, como comentamos anteriormente, \(\lambda\) no debe ser excesivamente grande ya que ocasionará que la función a minimizar presente mínimos locales con valores mucho menores que los de su alrededor.
Por lo tanto, trataremos de ajustar \(\lambda\) de manera que sea lo menor posible tratando de satisfacer la desigualdad
Apliquemos lo que se acaba de exponer el problema de la N-Reinas visto en el anterior apartado. En este caso la función de coste es la siguiente: