Refine HF solver to compute correct gradients

[soranjh]

Context:
The HF solver returns nan values when gradients are computed for any molecule larger than H2. The issue is related to the autograd failure in computing gradients of functions that perform matrix diagonalization for matrices with degenerate eigenvalues. See for example:

import autograd
import autograd.numpy as anp

s = anp.array([[ 1.00000000, 0.38675779, 0.58564814, 0.00000000, 0.00000000, -0.47793911],
               [ 0.38675779, 1.00000000, 0.24113657, 0.00000000, 0.00000000,  0.00000000],
               [ 0.58564814, 0.24113657, 1.00000000, 0.00000000, 0.00000000,  0.00000000],
               [ 0.00000000, 0.00000000, 0.00000000, 1.00000000, 0.00000000,  0.00000000],
               [ 0.00000000, 0.00000000, 0.00000000, 0.00000000, 1.00000000,  0.00000000],
               [-0.47793911, 0.00000000, 0.00000000, 0.00000000, 0.00000000,  1.00000000]])

def grad_fn():
    return lambda s: anp.linalg.eigh(s)[0][0]

print(autograd.grad(grad_fn())(s))

which returns

[[nan nan nan nan nan nan]
 [nan nan nan nan nan nan]
 [nan nan nan nan nan nan]
 [nan nan nan nan nan nan]
 [nan nan nan nan nan nan]
 [nan nan nan nan nan nan]]

The issue can be fixed by lifting the degeneracy through adding a small random value (in the order of 1.0e-12) to each of the diagonal elements of the matrix given above with

s = s + anp.diag(anp.random.rand(len(s)) * 1.0e-12)

This returns the gradients as

[[ 0.53101783 -0.15440008 -0.3480074   0.          0.          0.32262883]
 [-0.15440008  0.04489376  0.10118751  0.          0.         -0.09380837]
 [-0.3480074   0.10118751  0.22806983  0.          0.         -0.21143775]
 [ 0.          0.          0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.          0.          0.        ]
 [ 0.32262883 -0.09380837 -0.21143775  0.          0.          0.19601858]]

Description of the Change:
In the generate_scf, the overlap matrix is modified such that a small random perturbation in the order of 1.0e-12, is added to its diagonal elements.

Benefits:
This change allows computing gradient values for any given molecule.

Possible Drawbacks:
The diagonal elements of the overlap matrix are all 1 since the atomic orbitals are normalised. We have assumed that adding a small random value to these elements does not affect the accuracy of the results returned by the HF solver. However, this is not a fundamental fix to the problem of autograd failing to compute the gradients when eigh is applied to a matrix with degenerate eigenvalues.

Related GitHub Issues:
No related issues.

[github-actions]

Hello. You may have forgotten to update the changelog!
Please edit doc/releases/changelog-dev.md with:

• A one-to-two sentence description of the change. You may include a small working example for new features.
• A link back to this PR.
• Your name (or GitHub username) in the contributors section.

[codecov]

Codecov Report

Merging #1780 (afb5fbf) into qchem_hf (55a8d41) will increase coverage by 0.00%.
The diff coverage is 100.00%.

@@            Coverage Diff            @@
##           qchem_hf    #1780   +/-   ##
=========================================
  Coverage     98.95%   98.95%           
=========================================
  Files           215      215           
  Lines         16175    16181    +6     
=========================================
+ Hits          16006    16012    +6     
  Misses          169      169           

Impacted Files	Coverage Δ
pennylane/hf/hamiltonian.py	100.00% <100.00%> (ø)
pennylane/hf/hartree_fock.py	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 55a8d41...afb5fbf. Read the comment docs.

[josh146]

@soranjh is there a deeper technical issue with computing gradients with degenerate eigenvalues? Or is it simply a matter of Autograd not supporting it (but other frameworks do?)

[soranjh]

@soranjh is there a deeper technical issue with computing gradients with degenerate eigenvalues? Or is it simply a matter of Autograd not supporting it (but other frameworks do?)

@josh146, I tried jax for the same example given above and it worked. I believe it is just autograd that has this issue.

[josh146]

@josh146, I tried jax for the same example given above and it worked. I believe it is just autograd that has this issue.

Oh, that's good to know! This means there is the option for us to simply patch/fix the autograd.numpy.linalg.eigh function to get the correct behaviour

[rmoyard]

It looks like a good temporary fix for making it into the release, but we should as @josh146 mentioned probably rewrite the autograd function in qml.math afterwards :)

[josh146]

Thanks @soranjh, looks good from my end!

[josh146]

minor question: does adding in a constant * I term cause any redundant quantum evaluations 🤔 E.g., is PL evaluating <psi|I|psi> even though it doesn't need to?

Not an issue for this PR, but if this is the case, this is a bug we should fix in a separate PR!

[soranjh]

I have modified computation of the Hamiltonian such that all coefficients and ops are computed first and the PL Hamiltonian is computed once at the end using h = qml.Hamiltonian(coeffs, ops, simplify=True). This avoids multiple PL Hamiltonian additions and makes the process much faster and also fixes the issue of having redundant terms.

[ixfoduap]

Looks great!
Indeed, let's revisit this in the future for a cleaner fix

[agran2018]

Looks good.