In this tequila tutorial you can learn how to apply the Krylov method to approximate the ground state of a given Hamiltonian.

Theoretical introduction

Here we briefly introduce the main idea of the MultiReference Selected Quantum Krylov (MRSQK) method motivating why it is useful to have it implemented in Tequila.

MRSQK is a low-cost alternative to the quantum phase estimation algorithm that allows to generate the ground state of an hamiltonian as a linear combination of non-orthogonal Krylov basis states \(\mathcal{K}_s = \{\psi_\alpha, \ \alpha=1, 2, ... , N\}\). This basis is usually obtained via real time evolution from a reference set of states.

So, given an hamiltonian \(H\) and a Krylov basis \(\mathcal{K}_s\), the ground state of \(H\) can be written as:

\[|\Psi\rangle = \sum_\alpha c_\alpha|\psi_\alpha\rangle \quad .\]

The coefficients \(c_\alpha\) and the ground energy value \(E\) can be obtained by solving the following generalized eigenvalue problem: \[\mathbf{Hc} = \mathbf{Sc}E\] where the elements of the overlap matrix (\(\mathbf{S}\)) and Hamiltonian (\(\mathbf{H}\)) are \[S_{\alpha\beta} = \langle\psi_\alpha|\psi_\beta\rangle \quad ,\] \[H_{\alpha\beta} = \langle\psi_\alpha|H|\psi_\beta\rangle \quad .\]

Luckily, with tequila we can easily compute terms like \(S_{\alpha\beta}\) and \(H_{\alpha\beta}\) thanks to the tq.braket function. The krylov_method function uses it and allows to easily solve the generalized eigenvalue problem giving as output the energy \(E\) and the coefficients \(c_\alpha\).

Simple example

import tequila as tq
from tequila.apps.krylov import krylov_method
from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian
from tequila.tools.random_generators import make_random_circuit
import itertools as it
import numpy as np

Here we present a simple/trivial example in which we apply MRSQK. In order to do this, we create two quantum circuits randomly, \(|\psi\rangle\) and \(|\phi\rangle\), and we use these as Krylov basis.

np.random.seed(111)
n_krylov_states = 2

krylov_circs = [make_random_circuit(2, enable_controls=True) 
                for i in range(n_krylov_states)]

krylov_states = [tq.simulate(circ) for circ in krylov_circs]

1: Set the random seed for reproducibility
2: Number of Krylov states to generate
3: Create random quantum circuits, in this way it is very unlikely they will be orthogonal
4: Create the wavefunctions from the circuits

Then we build an Hamiltonian from these as follows: \[H = -|\psi\rangle\langle\psi|-|\phi\rangle\langle\psi|-|\psi\rangle\langle\phi|-|\phi\rangle\langle\phi|\] In this way we have an hermitian operator and we are sure that the Krylov space contains the ground state. This is a toy hamiltonian that has nothing to do with the Krylov method itself, it’s only needed to check the obtained states are the correct ones.

krylov_states_couples = list(it.product(krylov_states, repeat=2))


H = QubitHamiltonian()
for i, j in krylov_states_couples:
    H -= tq.paulis.KetBra(ket = i, bra = j)

1: Generate a list of all posible couples of Krylov states
2: Create a Hamiltonian from the obtained wavefunctions
3: Initialize an empty QubitHamiltonian object
4: For each couple of Krylov states, compute the braket and substract the term from the Hamiltonian

At this point we just need to call the krylov_method function, providing the Krylov circuits and the hamiltonian. It will build the matrices \(\mathbf{S_{\alpha\beta}}\) and \(\mathbf{H_{\alpha\beta}}\) and then return the ground energy \(E\) and the coefficients \(c_\alpha\):

kry_energies, kry_coefficients_matrix = krylov_method(krylov_circs, H)

kry_ground_energy = kry_energies[0]
kry_ground_coefficients = kry_coefficients_matrix[:,0]

1: Applying the Krylov method
2: Extract the ground state energy
3: Extract the coefficients

In order to check if the method gives meaningful solutions we can directly diagonalize the hamiltonian \(H\):

eigenvalues, eigenvectors = np.linalg.eigh(H.to_matrix())

1: Perform exact diagonalization of the Hamiltonian

As you can see below the ground energy and the ground states do correspond:

print('Ground State Energy Krylov: {:.4f}'.format(kry_ground_energy))      
print('Ground State Energy: {:.4f}'.format( eigenvalues[0]))

Ground State Energy Krylov: -1.6530
Ground State Energy: -1.6530

ground_state = tq.QubitWaveFunction()

for i in range(n_krylov_states):
    ground_state += kry_ground_coefficients[i]*krylov_states[i]

print(ground_state)

1: Initialize the ground state
2: Construct the ground state wavefunction by adding the scaled Krylov states
3: Print the ground state

+0.7338e^(-0.4020πi)|00> +0.5664i|10> -0.2799|11> +0.2497e^(+0.8983πi)|01>

wfn = tq.QubitWaveFunction.from_array(eigenvectors[:,0])
#print(eigenvectors[0])
print(wfn)

1: Create a QubitWaveFunction object from the array representing the first eigenvector

+0.7338|00> +0.2497e^(-0.6997πi)|01> +0.5664e^(+0.9020πi)|10> +0.2799e^(-0.5980πi)|11>

The obtained ground state seems different at first glance, but the states are actually identical due to equivalence up to a global phase. This depends on the employed simulator that decomposes gates up to a global phase. Using different ones may lead to different global phases, still having the same state.

We can easily check that these two are the same state by computing the fidelity between the two:

fidelity = abs(wfn.inner(ground_state.normalize()))**2
print(fidelity)

1: Compute the fidelity between the two states

1.0000000000000004