Added functionality to compute one-norm of qubit/majorana Hamiltonian from molecular integrals

[Emieeel]

Added a script in functionals to compute the 1-norm lambda = sum_i |h_i| of a qubit/majorana Hamiltonian H = sum_i h_i P_i directly from the molecular integrals (either input the integrals directly or input a MolecularData class). This is significantly faster than doing an explicit Jordan-wigner/bravyi kitaev transformation, which can become very costly for large systems. See https://arxiv.org/abs/2103.14753 for more information on its usefulness and a derivation for the 1-norm in terms of the integrals.

[obriente]

This looks good to me. @ncrubin , do you have any issues before I merge? I left a question in the code about whether this can be sped up using einsum, but this sounds like it might be a lot of work and Emiel has this working on a hundred spin orbitals so I think it's good enough for now?

[ncrubin]

Does antisymmeterization not work to remove the if statements for p>q etc?

[ncrubin]

I think once @obriente pulls back in @Emieeel to answer 1 or 2 more questions we are good to go. I would also like the documents to say that this is explicitly for Hamiltonians generated from RHF and ROHF. I don't think this L1 norm is applicable to Hamiltonians built from UHF orbitals

[Emieeel]

I think once @obriente pulls back in @Emieeel to answer 1 or 2 more questions we are good to go. I would also like the documents to say that this is explicitly for Hamiltonians generated from RHF and ROHF. I don't think this L1 norm is applicable to Hamiltonians built from UHF orbitals

Thats true, my understanding was that openfermion only handled restricted molecular Hamiltonians, but maybe I'm wrong?

[Emieeel]

I don't understand -- why does the coverage check fail?

[obriente]

Seems that there are two of the branches of your handling of no_constant that you don't enter. I've left comments at the relevant places - I think the structure of the code there could be simplified a bit as well.

Thanks for getting onto this!

[obriente]

(For the record, if you click on the details of the coverage check, and scroll to the bottom it will tell you the lines in question. There's a bunch of other text, but all the important stuff is at the bottom.)

[obriente]

Damn, just realised I was commenting on the wrong branch. Let me change to the new one.

[obriente]

@Emieeel , thanks for the patience. I left one comment which I think is a genuine bug, please ignore the others.

@ncrubin , please confirm that this now fits your expectations.

[ncrubin]

This will not be merged unless the commit authors actually document the functions. Spatial integrals is still not referenced.

[PabloAMC]

Hi @Emieeel, awesome and super useful job! By the way, would there be a possibility to make this also return capital Lambda (the maximum value of the amplitudes in the Jordan-Wigner operator) and Gamma (the number of non-zero terms)? I think it is pretty straightforward and it is used in some old phase estimation methods. Doing it myself would be easy but a bit dirty if I intend to publish the code. Thanks again!

[Emieeel]

Hi @PabloAMC, I didn't see you you're comment up to now unfortunately. I think it is definitely possible to return both Lambda and Gamma from spatial integrals with a similar code. The easiest way would just be during the summation to check every term whether its bigger than the current maximum, and keep count how many terms there are (if they are non-zero).

[PabloAMC]

Hi @Emieeel, I have made a small Jupiter notebook collab to try to compute lambda using the function you made and also from the corresponding Jordan-Wigner operator, but I think I don't get the same result (roughly double in the former vs the latter). I've also computed Gamma (the number of non-zero terms) and I think I got it wrong too. Capital Lambda (not shown) is right. The link is this one. Do you know what I may not be getting right for each of these two parameters? Thanks!

[Emieeel]

@PabloAMC I had a look at your notebook, and I have a few corrections. First of all -- and I guess this is why it is important that the function was properly documented -- you're working with spin-orbital integrals, while the input of the 1-norm code is spatial integrals. Then you'll see you get the right 1-norm. Then, computing Lambda is quite easy, as its just the maximum of everything we sum up to the one_norm.

Computing Gamma is a bit more tedious. If you look at https://doi.org/10.1103/PhysRevResearch.3.033127 Eq. (C13), the Hamiltonian is expressed in terms of Majorana operators. Taking into account that we sum over spin, we have 2# of 1bdy terms. For the 2bdy terms, the antisymmetric integrals are summed with p>r, s>q and a sum over spin, so we have 1/4 2 # of antisym terms. In trying to compute Gamma here, we find that there is another symmetry for the very last term that doesn't impact the 1-norm at all, but is important for Gamma. That is: switching all together pq \sigma <-> rs \tau (It doesn't matter for the 1-norm as g{pqrs} = g{rspq}). This results in a factor of 1/2 for the Gamma, so we have another 1/2 2 * # 2bdy terms. Also, don't forget to put a tolerance to set numerical errors to zero (OpenFermion has a standard tolerance of 1e-8 I believe).

I put my corrected code in this notebook, where I compute the 1-norm, Lambda and Gamma for O2.

[PabloAMC]

Quick question on the above:

that is: switching all together pq \sigma <-> rs \tau (It doesn't matter for the 1-norm as g{pqrs} = g{rspq}). This results in a factor of 1/2 for the Gamma, so we have another 1/2 2 # 2bdy terms.

Why are we multiplying by 2? I thought we had already taken into account the sum over spin

[Emieeel]

We still have the sum over sigma != tau in the last term

[PabloAMC]

Ah, of course, you're totally right

[PabloAMC]

Final question, if I may: lets say I have defined a molecule_geometry and a grid and run

spinless = False

V_dual = dual_basis_potential(grid = grid, spinless = spinless)
U_dual = dual_basis_external_potential(grid = grid, geometry = molecule_geometry, spinless = spinless)
T_primal = plane_wave_kinetic(grid, spinless = spinless)

lambda_value_V = V_dual.induced_norm() - V_dual.constant
lambda_value_U = U_dual.induced_norm() - U_dual.constant
lambda_value_T = T_primal.induced_norm() - T_primal.constant

However, if I run

JW_op = jordan_wigner(V_dual)
l = abs(np.array(list(JW_op.terms.values())))
lambda_V = np.sum(np.abs(l[1:]))

JW_op = jordan_wigner(U_dual)
l = abs(np.array(list(JW_op.terms.values())))
lambda_U = 2*np.sum(np.abs(l[:]))

JW_op = jordan_wigner(T_primal)
l = abs(np.array(list(JW_op.terms.values())))
lambda_T = 2*np.sum(np.abs(l[1:]))

I get different values. In the case of the T term and U terms I get them right, but in V I don't no matter how I do it. If you want to take a look I've done an example here. Thanks!

[Emieeel]

i'm not sure what you're calculating with the grid here, but it seems that the V_dual, U_dual and T_primal have a different induced norm before and after a jordan-wigner transformation. This is not weird as indeed some coefficients can cancel each other. If you calculate JW_op.induced_norm() you get the 1-norm of the QubitOperator. Also, the constant terms are zero here so watch out with l[1:] -- convert it to l[:] if the constant is zero.

[PabloAMC]

You're right, we should carry out some calculations similar to the ones you did in the article to match the Jordan-Wigner operator one-norm. However, what is clear is that

lambda_value_V = (V_dual.induced_norm() - V_dual.constant)/2
lambda_value_U = (U_dual.induced_norm() - U_dual.constant)/2
lambda_value_T = (T_primal.induced_norm() - T_primal.constant)/2

should be an upper bound to the actual JW one-norm, isn't it? I could use

JW_op.induced_norm()

but the problem is that the bottleneck is in computing JW_op

[Emieeel]

It is an upper bound, yes. If you express analytically the form of the V_dual, U_dual and T_primal in terms of (unique) Majorana operators, you can find out what the 1-norm is without actually doing the JW transform.

[PabloAMC]

I think I know the reason why both U_dual and T_primal give correct results, because they are one_body. Anyway, I think this is what I was looking for?

maj = openfermion.transforms.get_majorana_operator(V_dual)
val = maj.terms.values()
l = np.abs(np.array(list(val)))
one_norm = sum(l[1:])

and I get the same as if I do

JW_op = jordan_wigner(V_dual)
l_V = abs(np.array(list(JW_op.terms.values())))
lambda_V = np.sum(np.abs(l_V[1:]))

but much quicker. So thanks a lot @Emieeel !