Using quimb to find the lowest energy spin configuration of Ising model with external fields

[saraphys]

I would like to use the following code mentioned in the "Bayesian Optimizing the Energy" example of quimb documents:

def energy(x): p = len(x) // 2 gammas = x[:p] betas = x[p:] circ = qtn.circ_qaoa(terms, p, gammas, betas)

ZZ = qu.pauli('Z') & qu.pauli('Z')
ens = [
    circ.local_expectation(weight * ZZ, edge, optimize=opt)
    for edge, weight in terms.items()
]

However, my system is an Ising chain with external field. So, I have ZZ terms as well as Z terms. If anyone can help me if it is possible to find the lowest energy configuration of such a system with quimb or not? Thanks

[jcmgray]

Hi @saraphys, there are lots of options here, maybe you could provide some extra information:

1. what geometry do you want to address? arbitrary?
2. does it have to be a quantum circuit ansatz or could it be a regular TN?
3. is it ZZ and Z terms? (usually more studied is the 'tranverse field' case of ZZ and X terms.

[saraphys]

Thank you for the reply.

1) It is a 1D Ising chain with ZZ interactions and external field in z direction with the following Hamiltonian: H=\sum{} J{ij} Z_i Z_j +\sum_i h_i Z_i. 2) I want to use a quantum circuit ansatz and then use tensor contraction to make the problem easier because I have a long Ising chain. 3) Yes, I have the Hamiltonian mentioned in the above. It includes ZZ and Z term. It is not transverse Ising chain. Thanks

[jcmgray]

So for 1D something like DMRG should work very well. But if it needs to be 'quantum':

The absolute simplest would just be to add the terms separately:

terms = {
    (i, j): J[i, j] * ZZ 
    for i, j in edges,
} | {
    i: h[i] * Z
    for i in sites
}

then the energy of each term is still exactly:

[
    circ.local_expectation(G, where, optimize=opt)
    for where, G in terms.items()
]

however because the single site terms are entirely contained within the lightcones of the two site terms you may want for efficiencies sake to combine then, LocalHamGen can actually do this for you:

H2 = {(i, j): J[i, j] * ZZ for i, j in edges}
H1 = {i: h[i] * Z for i in sites}
ham = qtn.LocalHamGen(H2=H2, H1=H1)

now ham.terms contains just the minimal set of two-body terms with the single body terms merged in.

[saraphys]

Thanks for your clear explanations. since the command circ_ex = qtn.circ_qaoa(terms, p, gammas, betas) does not accept diagonal elements, I had to define two terms:

terms1 = {
    (i, j): J[i][j]
    for i, j in G1.edges
}

terms2 = {
    (i, j): J[i][j]*ZZ 
    for i, j in G1.edges
}| {
   i: h[i]* Z
    for i in G2.nodes
}

to be able to run the following command:

def energy(x):
    p = len(x) // 2
    gammas = x[:p]
    betas = x[p:]
    circ = qtn.circ_qaoa(terms1, p, gammas, betas)
    ens = [
        circ.local_expectation(G, where, optimize=opt)
        for where, G in terms2.items()
    ]
    return sum(ens).real

Would you please let me know if it is right way to do that. Thanks for your great help

[jcmgray]

since the command "circ_ex = qtn.circ_qaoa(terms, p, gammas, betas)" does not accept diagonal elements, I had to define two terms:

By 'diagonal elements' do you mean single site terms? If so yes this constructor just depends on the 'edges' of the graph, i.e. the two site terms since these define the geometry, but once the ansatz circuit is constructed you can compute what you like with it -> So, as long as the QAOA ansatz circuit is what you want, your snippet looks good.

---

(by the way on github if you put code snippets in three backticks like this:

```python
<CODE>

then it renders nicely and more readably - have just edited your comment for this).

[saraphys]

Thank you very much for you reply.

Thanks to your help, I could find the best "gammas" and "betas" that minimize the energy of my considering system. Now, I need to find the final state obtained by substituting the "gammas" and "betas" in QAOA Hamiltonian and then find it in the computational basis.

In order to find the final state, I could not find any circuit action that leads to the final state of the circuit:

circ = qtn.circ_qaoa(terms1, p, gammas, betas)

In order to find the final state in the computational basis, based on the quimb documentation, I guess the best way is to use "circ.compute_marginal" function which its input is a chosen set of qubits where the other qubits are fixed. But I could not understand how to choose the sets.

Thanks a lot

[jcmgray]

Hi @saraphys, sorry to be slow getting to this over the holidays!

So in QAOA you don't need to find the full state as such (which much be very large and complicated) instead you want to sample from it, the bit-strings you get will be the candidate solutions.

The method is to call e.g.:

for b in circ.sample(10):
    print(b)

see the docs here - https://quimb.readthedocs.io/en/latest/tensor-circuit.html#Generate-unbiased-samples. This calls compute_marginal for you. There are approximate ways to do this as well (i.e. sample substrings locally only) but that is a more advanced topic.

[saraphys]

Thank you very much for the reply. I really appreciate your time and great help.

In QAOA, we used the following code:

terms1 = {
    (i, j): J[i][j]
    for i, j in G1.edges
}

terms2 = {
    (i, j): J[i][j]*ZZ 
    for i, j in G1.edges
}| {
   i: h[i]* Z
    for i in G2.nodes
}
def energy(x):
    p = len(x) // 2
    gammas = x[:p]
    betas = x[p:]
    circ = qtn.circ_qaoa(terms1, p, gammas, betas)
    ens = [
        circ.local_expectation(G, where, optimize=opt)
        for where, G in terms2.items()
    ]
    return sum(ens).real

The output of the code will be x (the optimal values of gammas and betas) as well as energy (x). In your mentioned command:

for b in circ.sample(10):
    print(b)

I do not know how to define the circ. The following circuit can not be used because the ZZ operator can not be included in terms1:

circ = qtn.circ_qaoa(terms1, p, gammas, betas)

Moreover, the following command that computes expectation value is also not applicable:

circ.local_expectation(G, where, optimize=opt)
        for where, G in terms2.items()

Would you please let me know how I can define circuit to be able to get sample from it.

Thanks a lot

[jcmgray]

I'm a bit confused, in an earlier comment it seemed like you had the energy optimization working? In which case you should just be able to sample from...

circ = qtn.circ_qaoa(terms1, p, gammas, betas)

Maybe is mixing G1 and G2 a typo? If I just use the same graph everything looks like it works.

The only potential unexpected behaviour is if you are trying to optimize a graph of terms which has completely disconnected nodes?

Note the code for circ_qaoa is very simple - https://quimb.readthedocs.io/en/latest/_modules/quimb/tensor/circuit_gen.html#circ_qaoa - you could just modify it however you want.

If you still can't get things working, maybe post a minimal working example with G1, G2 etc. defined and with the exact definition of the circuit you want and the objective function you want to minimize.