Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 20 commits into from
Aug 25, 2021

Conversation

Emieeel
Copy link
Contributor

@Emieeel Emieeel commented Apr 30, 2021

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.

@google-cla google-cla bot added the cla: yes The Google CLA has been signed label Apr 30, 2021
(1/4 * two_body_integrals[p,q,p,q]))

htildepq = np.zeros(one_body_integrals.shape)
for p in range(n_orb):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ncrubin will know this better than I, but for the purposes of speed could this possibly be improved by using np.einsum?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because of the restriction on the sums I didn't know how to employ einsum here. There could be a way to do that, it would certainly make it faster... In the calculations of the paper I actually used numba http://numba.pydata.org/ to speed it up, which helped a lot, but I didn't include it here

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why not antisymmeterize the integrals and then do the summation?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also it seems like there is no spin-function information.

from openfermion import get_one_norm, MolecularData, jordan_wigner
from openfermion.config import DATA_DIRECTORY

class test_get_one_norm(unittest.TestCase):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: I think we're moving slightly away from unittest and towards pytest? I have no personal preference here, just wanting to flag.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm fine with everything :) I just mimicked some other test scripts

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes please move over to pytest

@Emieeel Emieeel requested a review from obriente May 6, 2021 11:09
@obriente
Copy link
Collaborator

obriente commented May 6, 2021

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
Copy link
Collaborator

ncrubin commented May 6, 2021

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

(1/4 * two_body_integrals[p,q,p,q]))

htildepq = np.zeros(one_body_integrals.shape)
for p in range(n_orb):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also it seems like there is no spin-function information.

return one_norm


def _get_one_norm_woconst(one_body_integrals, two_body_integrals):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

still needs a doc string

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given my above comments, I think this function should be exposed to the user, in which case the docstring from 'get_one_norm_int' can be copied down with minor adjustments.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, in quantum computing 1-norms should almost always be reported without the constant term. I'd suggest calling this function 'get_one_norm', and making it clear in the docs that the constant term is not added. ( @ncrubin , let me know if you have worries about this.)

n_orb = one_body_integrals.shape[0]

htildepq = np.zeros(one_body_integrals.shape)
for p in range(n_orb):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

n_orbs or spin-orbitals?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is all in spatial orbital basis

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So can you add this as a doc string?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Last thing. See above. Just document input with "parameters" like the other functions since this is now a public function.


one_norm = np.sum(np.absolute(htildepq))

for p in range(n_orb):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

use itertools product or convert to contraction over antisymmeterized integrals.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Converted to contraction over antisymmetrized integrals. Thank you for the suggestion!

from openfermion import MolecularData


def get_one_norm(mol_or_int, return_constant=True):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Be explicit with inputs. Use default None and allow the user to specify. This partially documents the function.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It will also give you freedom to calculate norms for 1-body ops only. or 2-body only. Is there a difficulty extending this to arbitrary interactions?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm, this is only valid in the first place for two-body tensors which have real orbital symmetries. I'm not sure how to extend this, as an n-body interaction will not be a molecular system. There is probably a way to write an arbitrary n-body Fermionic Hamiltonian in Majorana representation in terms of the original coefficients though? (I just wouldn't know how)

@Emieeel Emieeel requested a review from ncrubin May 20, 2021 09:26
@ncrubin
Copy link
Collaborator

ncrubin commented Jun 21, 2021

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
Copy link
Contributor Author

Emieeel commented Jun 30, 2021

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
Copy link
Contributor Author

Emieeel commented Jun 30, 2021

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

if isinstance(mol_or_int, MolecularData):
return _get_one_norm_mol(mol_or_int, no_constant=no_constant)
else:
if no_constant is not None:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this tree - is the idea that no_constant should be either True, False, or None? Also, do we plan the default 'None' to return the constant, or not? I'd recommend setting the default to False or True, and then having a single if statement.

(Also, it seems one of your issues is that you're not covering all three branches of this in your test case).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm that is what it was before. I thought this is what Nick meant with "Use default None and allow the user to specify", but I'm probably wrong.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Heh, fair enough. @ncrubin , can you comment here? I think that a 'None' default is less explanatory as it doesn't tell the user whether the constant will be by default returned or not. Is there any issue with using True/False as a default?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @Emieeel , just chatted with @ncrubin about this, and I think that the original comment was misunderstood. The point is that instead of taking a tuple as input with all the different data, you should split it up into the different parts, i.e.

def get_one_norm_int(constant, one_body_integrals, two_body_integrals, no_constant=True),

and make a second function 'get_one_norm_mol' which takes a molecular Hamiltonian.

Sorry for the delay in sorting this out!

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could also do a 'get_one_norm' and 'get_one_norm_mol' as a convenience function.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm sorry it took so long, I forgot a bit about it; Is it fine now? I defined two functions: get_one_norm_int and get_one_norm_mol, with no_constant=True.


def _get_one_norm_mol(molecule, no_constant=None):
"""Compute one_norm for a MolecularData class"""
if no_constant is not None:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same issue here as above.

@obriente
Copy link
Collaborator

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!

@Emieeel Emieeel requested a review from obriente August 11, 2021 11:53
obriente
obriente previously approved these changes Aug 11, 2021
Copy link
Collaborator

@obriente obriente left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM. I notice that Nick has two outstanding comments regarding docstrings - as the functions are private and the public functions are clearly documented then I think this is ok. Would suggest merging as-is to prevent this dragging on any longer.

Copy link
Collaborator

@ncrubin ncrubin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The get_one_norm_mol needs to still be changed. Overloading the input like that should not be done. just make a separate function.

"""Compute 1-norm given molecular integrals"""
n_orb = one_body_integrals.shape[0]

htilde = constant
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be htilde = abs(constant), right? Otherwise we can wind up with negative 1 norms...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No it is actually correct like this. The reason is that you have in principle two (or three if you consider an active space) contributions to the constant term in majorana operators: The nuclear repulsion (always positive), the mean field energy of the frozen orbitals if you consider a CAS (always negative) and the ((1 / 2 * g_[ppqq]) - (1 / 4 *g_[pqpq])) term from rearranging the majorana operators. It is important to actually take the absolute value after counting them up together as they can cancel eachother. I take the absolute value of htilde on the line below, so the 1-norm can never be negative.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah ok, thanks for the clarification!

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no I think this goes back to the comment I made in the "issue" about this. It depends on how you are defining and accessing the LCU. There is a way to represent the LCU such that what @Emieeel is saying is true. And it is a clever trick that he's pointed out there. But under the simplest way to simulate the LCU, you don't get that subtraction from using the Majorana operators.

@obriente
Copy link
Collaborator

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

@obriente
Copy link
Collaborator

@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
Copy link
Collaborator

ncrubin commented Aug 24, 2021

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

Copy link
Collaborator

@ncrubin ncrubin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@Emieeel Thanks for the changes. There are only two remaining things that are small. Please document the functions that you have now made public. This can either reference another function or you can duplicate to docstring.

Lastly, please add your name and institutions to the contributor list for the package! This is a very useful addition to the code.

molecule.two_body_integrals)


def get_one_norm_mol_woconst(molecule):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay this is the very last thing and then this is done. It seems that now that we aren't dispatching these functions need to be properly documented. You You can just reference get_one_norm_int or specify that the parameters molecule: molecularData class.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think its fixed now. Thank you for being precise Nick, it was a learning experience for me :)

n_orb = one_body_integrals.shape[0]

htildepq = np.zeros(one_body_integrals.shape)
for p in range(n_orb):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Last thing. See above. Just document input with "parameters" like the other functions since this is now a public function.

@PabloAMC
Copy link

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 Emieeel requested a review from ncrubin August 25, 2021 11:18
@ncrubin ncrubin merged commit 8ef0c8b into quantumlib:master Aug 25, 2021
@Emieeel
Copy link
Contributor Author

Emieeel commented Sep 13, 2021

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
Copy link

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
Copy link
Contributor Author

Emieeel commented Sep 21, 2021

@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
Copy link

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
Copy link
Contributor Author

Emieeel commented Sep 21, 2021

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

@PabloAMC
Copy link

Ah, of course, you're totally right

@PabloAMC
Copy link

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
Copy link
Contributor Author

Emieeel commented Sep 21, 2021

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
Copy link

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
Copy link
Contributor Author

Emieeel commented Sep 21, 2021

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
Copy link

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 !

ncrubin added a commit to ncrubin/OpenFermion that referenced this pull request Jul 25, 2022
… from molecular integrals (quantumlib#725)

* Added get_one_norm with tests

* format incremental applied

* reference to the paper

* changed from unittest to pytest

* Now pytest works...

* My bad... fixed formatting now

* anti-symmetric integral implementation and better documentation

* Specified input must be RHF or ROHF Hamiltonian

* I think everything is good now?

* no overloading functions...

* type correction

* Documented the functions properly

* Added Type hints to the functions

Co-authored-by: Nicholas Rubin <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
cla: yes The Google CLA has been signed
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants