Ctrl+K

ADAPT con el Modelo de Lipkin

May 06, 2024 | 3 min read

ADAPT con el Modelo de Lipkin#

Hide code cell source 
import scipy.optimize
import scipy.sparse
import numpy as np
import random
import qibo
from qibo import hamiltonians as ham, models as m, gates as g
from qibo.symbols import X, Y, Z, I

El Modelo de Lipkin#

Un problema abierto en física nuclear es la determinación de la energía del estado fundamental para núcleos atómicos de toda la carta de Segrè. El motivo de que esté abierto es, principalmente, que los modelos que se utilizan habitualmente no son resolubles analíticamente y muy complicados de resolver numéricamente [4], [15], debido a la escala exponencial de los grados de libertad con el tamaño del núcleo. Es por ello que algoritmos como los VQE pueden jugar un papel importante, ya que pueden obtener dicha energía del estado fundamental y ser implementados en un ordenador cuántico, donde el número de estados computacionales que pueden describir los grados de libertad del núcleo también escala exponencialmente [20], [19], [6].

No obstante, el modelo que vamos a resolver aquí, llamado el modelo de Lipkin-Meshkov-Glick [12] o simplemente de Lipkin, no es más que un preámbulo, pues este sí es resoluble analíticamente (como el caso del Modelo de Ising 1D y 2D).

En este modelo consideramos un sistema fermiónico de n partículas que pueden ocupar dos niveles cuya separación en energía es un gap de magnitud \(\varepsilon\) y ambos tienen degeneración \(\Omega = n\). Este modelo presenta una interacción a dos cuerpos constante de magnitud \(V\). El Hamiltoniano del modelo de Lipkin viene dado por:

\[ \mathcal{H} = \frac{1}{2} \varepsilon \sum_{i=1}^n Z_i - \frac{1}{2} V \sum_{i < j} (X_i X_j - Y_i Y_j) \]

Un ejemplo gráfico de este modelo puede observarse en la siguiente imagen:

docs/pictures/Lipkin-model.png

Fig. 2 Representación esquemática del modelo de Lipkin.#

# Lipkin hamiltonian
# eps es el gap, V es la interacción a dos cuerpos
def lipkin_ham_qibo(n, eps, V):

    sym_Z = sum(Z(i) for i in range(n))
    sym_XY = sum(X(i) * X(j) - Y(i) * Y(j) for i in range(n) for j in range(i+1, n))
    sym_H = .5*eps*sym_Z -.5*V *sym_XY

    H = ham.SymbolicHamiltonian(sym_H)

    return np.matrix(H.matrix)

ADAPT-VQE#

Para este problema vamos a explorar un nuevo tipo de ansatz utilizable en problemas VQE: ADAPT [9]. Se trata de un sistema muy similar a los ansätze tipo UCC: Tenemos un conjunto de operadores \(\{\mathcal{A}_i (\theta_i)\}\) con el que podemos codificar nuestros fenómenos físicos, como excitaciones o desexcitaciones dobles. En un problema tipo UCC, obtendríamos el operador unitario \(U = \exp\left[\sum_i \mathcal{A}_i\right]\) y haríamos una descomposición de Trotter-Suzuki.

Sin embargo, en este caso lo que hacemos es definir nuestro conjunto de operadores posibles de forma que \([\mathcal{A}_i, \mathcal{A}_j] = 0\). Con esto, conseguimos que

\[ U = \exp\left[\sum_i \mathcal{A}_i\right] = \prod_i e^{\mathcal{A}_i} \]

Lo que hacemos entonces es aplicar los operadores de forma sucesiva. Para ello, comenzamos calculando el gradiente de todos los posibles operadores de nuestro conjunto imponendo \(\theta_i = 0\). Esto se puede hacer de forma analítica, imponiendo \(| \psi^{(0)} \rangle\) como nuestro estado de referencia y calculando

\[ \langle \psi^{(0)} | [\mathcal{H}, \mathcal{A}_i (0)] | \psi^{(0)} \rangle \]

tomando \(\mathcal{A^{(1)}}\) como el operador que maximiza dicha cantidad.

Entonces, calculamos nuestro nuevo estado como \(|\psi^{(1)} (\theta) \rangle = e^{\mathcal{A}^{(1)} (\theta)} |\psi^{(0)}\rangle\). Con este estado podemos calcular la energía de la forma usual, \(E^{(1)} = \langle \psi^{(1)} (\theta) | \mathcal{H} | \psi^{(1)} (\theta) \rangle\), y aquí ya usamos un proceso de minimización clásico (BFGS, SPSA, Powell…) para obtener el valor óptimo de \(\theta\).

Repetimos este procedimiento de forma iterativa, seleccionando un nuevo operador \(\mathcal{A}^{(2)}\) maximizando esta vez \(\langle \psi^{(1)} (\theta_{\text{opt}})| [\mathcal{H}, \mathcal{A}_i (0)] | \psi^{(1)} (\theta_{\text{opt}}) \rangle \), calculando \(|\psi^{(2)} (\theta) \rangle = e^{\mathcal{A}^{(2)} (\theta)} |\psi^{(1)}\rangle\) y minimizando \(E^{(2)} = \langle \psi^{(2)} (\theta) | \mathcal{H} | \psi^{(2)} (\theta) \rangle\); y así sucesivamente hasta converger al estado fundamental (o a un mínimo local). Cabe destacar que en esta última nueva minimización obtenemos un conjunto de parámetros óptimos que vuelve a ajustar los parámetros de iteraciones anteriores.

Al final de este proceso, el circuito que habríamos ejecutado, visto desde la perspectiva del sistema global, sería

\[ |g\rangle \approx e^{\mathcal{A}^{(n)}(\theta_n)} e^{\mathcal{A}^{(n-1)}(\theta_{n-1})}\dots e^{\mathcal{A}^{(2)}(\theta_2)}e^{\mathcal{A}^{(1)}(\theta_1)} |\psi^{(0)} \rangle \]

El punto clave de este formalismo está en el conjunto \(\{\mathcal{A}_i\}\), ya que debe conformar un Conjunto Completo de Operadores que Conmutan, es decir, debemos considerar operadores que nos permitan obtener cualquier estado posible de nuestro modelo, barriendo por completo el espacio de Hilbert.

# Construimos todos los operadores necesarios para el modelo de Lipkin
# Esto devuelve todos los operadores a dos cuerpos para un sistema de n niveles
# Esto tiene que ser un CONJUNTO COMPLETO (en el sentido algebraico) para cualquier configuración del modelo de Lipkin
def oper_pool(n):
    import numpy as np
    from qibo import matrices as m

    oper_mats=[]
    oper_labs=[]

    for p in range(0,n):
        oper_mats.append(X(p).full_matrix(n))
        oper_labs.append('X'+str(p))
        oper_mats.append(Y(p).full_matrix(n))
        oper_labs.append('Y'+str(p))
        oper_mats.append(Z(p).full_matrix(n))
        oper_labs.append('Z'+str(p))


        for q in range(p+1,n):

            oper_mats.append(0.5*(np.dot(X(p).full_matrix(n),X(q).full_matrix(n)) + np.dot(Y(p).full_matrix(n),Y(q).full_matrix(n))))
            oper_labs.append('0.5 ('+'X'+str(p)+' '+'X'+str(q)+' '+'+'+' '+'Y'+str(p)+' '+'Y'+str(q)+')')

            oper_mats.append(0.5*(np.dot(X(p).full_matrix(n),X(q).full_matrix(n)) - np.dot(Y(p).full_matrix(n),Y(q).full_matrix(n))))
            oper_labs.append('0.5 ('+'X'+str(p)+' '+'X'+str(q)+' '+'-'+' '+'Y'+str(p)+' '+'Y'+str(q)+')')

            oper_mats.append(0.5*(np.dot(X(p).full_matrix(n),Y(q).full_matrix(n)) + np.dot(Y(p).full_matrix(n),X(q).full_matrix(n))))
            oper_labs.append('0.5 ('+'X'+str(p)+' '+'Y'+str(q)+' '+'+'+' '+'Y'+str(p)+' '+'X'+str(q)+')')

            oper_mats.append(0.5*(np.dot(Y(p).full_matrix(n),X(q).full_matrix(n)) - np.dot(X(p).full_matrix(n),Y(q).full_matrix(n))))
            oper_labs.append('0.5 ('+'Y'+str(p)+' '+'X'+str(q)+' '+'-'+' '+'X'+str(p)+' '+'Y'+str(q)+')')

            oper_mats.append(np.dot(X(p).full_matrix(n),Z(q).full_matrix(n)))
            oper_labs.append('X'+str(p)+' '+'Z'+str(q) )

            oper_mats.append(np.dot(Y(p).full_matrix(n),Z(q).full_matrix(n)))
            oper_labs.append('Y'+str(p)+' '+'Z'+str(q) )

            oper_mats.append(np.dot(Z(p).full_matrix(n),Z(q).full_matrix(n)))
            oper_labs.append('Z'+str(p)+' '+'Z'+str(q) )


    return [oper_labs, oper_mats]

# Obtenemos una lista de Pauli-strings y sus representaciones matriciales






# En cada iteración, obtenemos el operador que nos da un mayor gradiente
# Evaluamos conmutadores tipo <state|[H, A(0)]|state>
# Devuelve el operador óptimo (con máximo gradiente) y el valor de dicho gradiente
# Con verbose=True nos da el valor de todos los posibles gradientes
def optimal_op(state, ham, ops, mats, verbose=False):

    #Evaluate commutators
    opt_op = []
    com_pre = 0.
    for m in range(len(ops)):
        Om, Am = ops[m], mats[m]
        com_mat = ham.dot(Am)-Am.dot(ham)
        com_val = -1j*state.T*com_mat*state
        if verbose:
         if (abs(com_val) > 1e-6): print(m, Om, com_val)
        if abs(com_val) > (abs(com_pre)+1e-6):
            opt_op.append(Om)
            com_pre = com_val

    return [opt_op, abs(com_pre)]





# Aquí actualizamos el estado, haciendo |psi(i+1)> = e^{A(i+1)(theta)} |psi(i)>
def update_state(state,mats,pars):

    

    assert len(mats) == len(pars)
    update_state = state[:]
    update_state = np.matrix(update_state).T

    for k in reversed(range(0, len(mats))):
        update_state = scipy.sparse.linalg.expm_multiply(1j*pars[k]*mats[k], update_state)

    #print("in update state: ", pars)

    return update_state


# Computamos el valor esperado <psi|H|psi> de forma que podamos minimizar después en función del valor del parámetro
def ener(pars, state, ham, mats): 

    state_it = update_state(state,mats,pars)

    ener = state_it.conj().T * ham * state_it

    assert abs(ener.imag) < 1e-6

    return ener.real


# Utilizando otro Hamiltoniano en lugar de Lipkin, seguimos teniendo una implementación apropiada de ADAPT

Vamos a presentar un caso práctico utilizando un sistema de 4 partículas con unos valores concretos para \(\varepsilon\) y \(V\).

####

# Parámetros del Hamiltoniano de Lipkin
n = 4
eps, V = 0.42, 0.482

# Parámetros del proceso de optimización
max_iter = 300
options = {'disp': False, 'gtol': 1e-10}


# Construimos el hamiltoniano
hamiltonian = lipkin_ham_qibo(n, eps, V)
eigvals, eigvecs = scipy.linalg.eigh(hamiltonian)
print("Min eigenvalue from diagonalization: ", min(eigvals))
print(" ")

# Construimos el conjunto de operadores
[ops_pool, mats_pool] = oper_pool(n)

# Estado de referencia
initial_state = np.zeros(2**n,dtype='complex128')
mbeners = hamiltonian.diagonal()
index_lowest_energy = np.argmin(mbeners)
state = initial_state[:]
state[index_lowest_energy]=1.
ener_new = np.min(mbeners)

# Inicializamos los arrays de soluciones
pars = []
sol_mat = []
sol_ops = []
opt_mat = []
opt_ops = []



print("Start of iterative procedure")
print("============================")
print(" ")
print(" ")
print("it  ener")
print("===========")
state = np.matrix(state).T
for it in range(max_iter):

    

    [op_it_array, toler]     = optimal_op(np.matrix(state), hamiltonian, ops_pool, mats_pool, verbose=0)
    if (len(op_it_array) > 0 ): op_it = op_it_array[-1]

    ener_prev = ener_new

    if (abs(toler) < 1e-6): # Convergencia en norma (todos los gradientes son 0)
        print(" ")
        print("============================")
        print("Calculations have converged in norm!")
        print("============================")
        break

    op_it_mat          = mats_pool[ops_pool.index(op_it)]
    opt_ops.append(op_it)
    opt_mat.append(op_it_mat)
    opt_pars = np.append(pars, 0.) #np.insert(pars,0,0.)


    opt_sol = scipy.optimize.minimize(ener, opt_pars, args=(initial_state, hamiltonian, opt_mat), options=options, method='BFGS')

   # opt_sol = scipy.optimize.minimize(ener, opt_pars, args=(initial_state, hamiltonian, opt_mat), method='Nelder-Mead')


    ener_new = opt_sol.fun

    if (abs(ener_new - ener_prev) < 1e-6): # Convergencia en energía (alcanzado un mínimo, aunque podría ser local)
        print(" ")
        print("============================")
        print("Calculations have converged in energy!")
        print("============================")
        break



    pars = opt_sol.x
    sol_ops.append(op_it)
    sol_mat.append(op_it_mat)

    state = update_state(initial_state,sol_mat,pars)

    print(it, opt_sol.fun, toler, pars)


print(" ")
print(" ")
print("Parameters and operators")
print("------------------------")
for k in range(len(sol_ops)):
 print(k, pars[k], sol_ops[k])


# Estado fundamental obtenido (notación sparse)
gs = scipy.sparse.csc_matrix(update_state(initial_state,sol_mat,pars))

print(" ")

print("Ground state: ")
print(gs)

[Qibo 0.1.12.dev0|INFO|2024-06-13 07:33:06]: Using tensorflow backend on /device:CPU:0

[Qibo 0.1.12.dev0|WARNING|2024-06-13 07:33:06]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.

Min eigenvalue from diagonalization:  -1.8690874778886049
 
Start of iterative procedure
============================
 
 
it  ener
===========
0 -1.0593152586948005 [[0.964]] [0.42701317]
1 -1.278630517389601 [[0.964]] [0.42701317 0.42701317]
2 -1.4399699907734906 [[0.63330259]] [0.42701317 0.42701316 0.32300631]
3 -1.6013094641573797 [[0.79865131]] [0.42701316 0.42701316 0.32300631 0.32300632]

4 -1.7351984710229948 [[0.50568615]] [0.42701317 0.42701317 0.3230063  0.32300629 0.27094707]
5 -1.869087477888605 [[0.73484309]] [0.42701317 0.4270132  0.3230063  0.3230063  0.27094706 0.27094704]
 
============================
Calculations have converged in norm!
============================
 
 
Parameters and operators
------------------------
0 0.42701316857130084 0.5 (X0 Y1 + Y0 X1)
1 0.4270132009161558 0.5 (X2 Y3 + Y2 X3)
2 0.3230063008194932 0.5 (X0 Y2 + Y0 X2)
3 0.32300629767953715 0.5 (X1 Y3 + Y1 X3)
4 0.27094705584500045 0.5 (X0 Y3 + Y0 X3)
5 0.27094704405041853 0.5 (X1 Y2 + Y1 X2)
 
Ground state: 
  (0, 0)	(0.27529140908499256+0j)
  (3, 0)	(0.2578798480867337+0j)
  (5, 0)	(0.25787985637875865+0j)
  (6, 0)	(0.25787983306969847+0j)
  (9, 0)	(0.2578798212751166+0j)
  (10, 0)	(0.2578798532388026+0j)
  (12, 0)	(0.2578798804315887+0j)
  (15, 0)	(0.7247085909150072+0j)
Licencia Creative Commons

License: Licencia Creative Commons Atribución-CompartirIgual 4.0 Internacional.

This work has been financially supported by the Ministry for Digital Transformation and of Civil Service of the Spanish Government through the QUANTUM ENIA project call - Quantum Spain project, and by the European Union through the Recovery, Transformation and Resilience Plan - NextGenerationEU within the framework of the Digital Spain 2026 Agenda.

https://quantumspain-project.es/wp-content/uploads/2024/02/Banner-QS_GOB_v2.png
Contents