Skip to content
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

Add a modified simplify function for building Hamiltonians #2103

Merged
merged 13 commits into from
Jan 28, 2022
8 changes: 8 additions & 0 deletions doc/releases/changelog-dev.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,14 @@

<h3>New features since last release</h3>

* Added a modified version of the `simplify` function to the `hf` module.
[(#2103)](https://github.com/PennyLaneAI/pennylane/pull/2103)

This function combines redundant terms in a Hamiltonian and eliminates terms with a coefficient
smaller than a cutoff value. The new function makes construction of molecular Hamiltonians more
efficient. For LiH, as an example, the time to construct the Hamiltonian is reduced roughly by a
factor of 20.

* The JAX interface now supports evaluating vector-valued QNodes. Vector-valued
QNodes include those with:
* `qml.probs`;
Expand Down
53 changes: 51 additions & 2 deletions pennylane/hf/hamiltonian.py
Original file line number Diff line number Diff line change
Expand Up @@ -247,7 +247,8 @@ def hamiltonian(*args):
for n, t in enumerate(h_ferm[1]):

if len(t) == 0:
coeffs = np.array([h_ferm[0][n]])
coeffs = anp.array([0.0])
coeffs = coeffs + np.array([h_ferm[0][n]])
soranjh marked this conversation as resolved.
Show resolved Hide resolved
ops = ops + [qml.Identity(0)]

elif len(t) == 2:
soranjh marked this conversation as resolved.
Show resolved Hide resolved
Expand Down Expand Up @@ -282,7 +283,7 @@ def hamiltonian(*args):
coeffs = np.concatenate([coeffs, np.array(op[0]) * h_ferm[0][n]])
ops = ops + op[1]

h = qml.Hamiltonian(coeffs, ops, simplify=True)
h = simplify(qml.Hamiltonian(coeffs, ops), cutoff=cutoff)

return h

Expand Down Expand Up @@ -355,6 +356,54 @@ def _generate_qubit_operator(op):
return c, o


def simplify(h, cutoff=1.0e-12):
r"""Add together identical terms in the Hamiltonian.

The Hamiltonian terms with identical Pauli words are added together and eliminated if the
overall coefficient is smaller than a cutoff value.

Args:
h (Hamiltonian): PennyLane Hamiltonian
cutoff (float): cutoff value for discarding the negligible terms

Returns:
Hamiltonian: Simplified PennyLane Hamiltonian

**Example**

>>> c = np.array([0.5, 0.5])
soranjh marked this conversation as resolved.
Show resolved Hide resolved
>>> h = qml.Hamiltonian(c, [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)])
>>> print(simplify(h))
(1.0) [X0 Y1]
"""
wiremap = dict(zip(h.wires, range(len(h.wires) + 1)))

c = []
o = []
for i, op in enumerate(h.terms[1]):
soranjh marked this conversation as resolved.
Show resolved Hide resolved
op = qml.operation.Tensor(op).prune()
op = qml.grouping.pauli_word_to_string(op, wire_map=wiremap)
if op not in o:
c.append(h.terms[0][i])
o.append(op)
else:
c[o.index(op)] += h.terms[0][i]

coeffs = []
ops = []
nonzero_ind = np.argwhere(abs(np.array(c)) > cutoff).flatten()
soranjh marked this conversation as resolved.
Show resolved Hide resolved
for i in nonzero_ind:
coeffs.append(c[i])
ops.append(qml.grouping.string_to_pauli_word(o[i], wire_map=wiremap))

try:
coeffs = qml.math.stack(coeffs)
except ValueError:
soranjh marked this conversation as resolved.
Show resolved Hide resolved
pass

return qml.Hamiltonian(coeffs, ops)


Copy link
Contributor

Choose a reason for hiding this comment

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

It would be great to reproduce this pattern in the Hamiltonian class. At the moment, H.simplify() is a method that manipulates the internal data of H, which is dangerous. We'd like it to return a new operator. Another option is qml.simplify(H), which could be cool because other ops may be simplified in some way too.

I don't think this relates to this PR though, which is just about building Hamiltonians in a special way, right?

Copy link
Contributor

Choose a reason for hiding this comment

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

This PR is literally adding a function that simplifies any input Hamiltonian. As far as I can tell, it should be possible to reproduce in the Hamiltonian class, or as you suggest, having qml.simplify(H). I really like this latter suggestion actually

Copy link
Contributor Author

@soranjh soranjh Jan 25, 2022

Choose a reason for hiding this comment

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

Hi @mariaschuld, @ixfoduap and @josh146. Here are my thoughts on the simplify function.

This is a general functions that does not necessary belong to qml.hf. However, the way it is written currently is based on one important point: convert Pauli words to strings and manipulate the strings with python. For example, instead of looping over the words and using compare, one can simply use
if o in ['XYZ', ...]
However, I am not sure we want to manipulate Hamiltonians/operators by turning them to strings. Another philosophy could be generating the binary representation of the words and Hamiltonians and manipulate those objects numerically to, for example, compare or multiply them. I think we need to decide first what approach we want to follow for operator manipulation/arithmetic in general. This is the only reason that I hided _simplify inside hf.hamiltonian for now.

So, I want your thoughts on these options:

  1. Convert hf.hamiltonian._simplify to qml.simplify(H).
  2. Wait until we have an ADR on how we really want to efficiently work with operators/Hamiltonian in PL and then create qml.simplify(H) accordingly.
  3. Change hf.hamiltonian._simplify to hf.hamiltonian.simplify.

Copy link
Member

Choose a reason for hiding this comment

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

@soranjh I think these are very valid points, still a lot needs to be figured out in terms of how to best implement comparison/simplification long-term.

Having said that, it's not quite clear when this research will be done/implemented. Since your version already is an improvement over what is already in PL, I am very much in favour of making it more general now, and continuing to improve it iteratively down the line.

In other words: since this is already a net benefit over what PL users have access to, I don't think there is much disadvantage to making this generally available in the short term?

Copy link
Contributor

@mariaschuld mariaschuld Jan 26, 2022

Choose a reason for hiding this comment

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

@soranjh @josh146 I would suggest to go for 2 in this PR - this is the least work for you now, and it incurs only 100% non-user-facing changes in the future!

Then, echoing Josh, I would suggest to soon make a PR to call the private function (hacky for now) in Hamiltonian.simplify if possible, or to otherwise improve that method with the insights gained here. We could combine it with making this Hamiltonian method return a new operator, but this breaks a few user-facing things so it is not trivial.

Lastly, we win more time for the ADR of how to do things properly.

Copy link
Member

Choose a reason for hiding this comment

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

@mariaschuld makes the following proposal, which I think makes a lot of sense:

  1. Continuing working on this PR

  2. Once this PR is merged, improve the existing Hamiltonian.simplify() method with this code

  3. Finally, we can use what we learnt to (a) inform the ADR, and (b) help us decide any final UI/design decisions.

def _pauli_mult(p1, p2):
r"""Return the result of multiplication between two tensor products of Pauli operators.

Expand Down
58 changes: 6 additions & 52 deletions pennylane/hf/tapering.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,11 +17,12 @@
# pylint: disable=unnecessary-lambda

import functools
import numpy

import autograd.numpy as anp
import numpy
import pennylane as qml
from pennylane import numpy as np
from pennylane.hf.hamiltonian import _generate_qubit_operator
from pennylane.hf.hamiltonian import _generate_qubit_operator, simplify


def _binary_matrix(terms, num_qubits, wire_map=None):
Expand Down Expand Up @@ -193,7 +194,7 @@ def get_generators(nullspace, num_qubits):
x, z = op
tau @= pauli_map[f"{x}{z}"](idx)

ham = qml.Hamiltonian([1.0], [tau], simplify=True)
ham = simplify(qml.Hamiltonian([1.0], [tau]))
generators.append(ham)

return generators
Expand Down Expand Up @@ -328,54 +329,7 @@ def _observable_mult(obs_a, obs_b):
o.append(op)
c.append(phase * obs_a.terms[0][i] * obs_b.terms[0][j])

return _simplify(qml.Hamiltonian(qml.math.stack(c), o))


def _simplify(h, cutoff=1.0e-12):
r"""Add together identical terms in the Hamiltonian.

The Hamiltonian terms with identical Pauli words are added together and eliminated if the
overall coefficient is zero.

Args:
h (Hamiltonian): PennyLane Hamiltonian
cutoff (float): cutoff value for discarding the negligible terms

Returns:
Hamiltonian: Simplified PennyLane Hamiltonian

**Example**

>>> c = np.array([0.5, 0.5])
>>> h = qml.Hamiltonian(c, [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)])
>>> print(_simplify(h))
(1.0) [X0 Y1]
"""
s = []
wiremap = dict(zip(h.wires, range(len(h.wires) + 1)))

for term in h.terms[1]:
term = qml.operation.Tensor(term).prune()
s.append(qml.grouping.pauli_word_to_string(term, wire_map=wiremap))

o = list(set(s))
c = [0.0] * len(o)
for i, item in enumerate(s):
c[o.index(item)] += h.terms[0][i]
c = qml.math.stack(c)

coeffs = []
ops = []
nonzero_ind = np.argwhere(abs(c) > cutoff).flatten()
for i in nonzero_ind:
coeffs.append(c[i])
ops.append(qml.grouping.string_to_pauli_word(o[i], wire_map=wiremap))
try:
coeffs = qml.math.stack(coeffs)
except ValueError:
pass

return qml.Hamiltonian(coeffs, ops)
return simplify(qml.Hamiltonian(qml.math.stack(c), o))


def clifford(generators, paulix_ops):
Expand Down Expand Up @@ -478,7 +432,7 @@ def transform_hamiltonian(h, generators, paulix_ops, paulix_sector):
c = anp.multiply(val, h.terms[0])
c = qml.math.stack(c)

return _simplify(qml.Hamiltonian(c, o))
return simplify(qml.Hamiltonian(c, o))


def optimal_sector(qubit_op, generators, active_electrons):
Expand Down
47 changes: 47 additions & 0 deletions tests/hf/test_hamiltonians.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
_generate_qubit_operator,
_pauli_mult,
_return_pauli,
simplify,
generate_electron_integrals,
generate_fermionic_hamiltonian,
generate_hamiltonian,
Expand Down Expand Up @@ -379,3 +380,49 @@ def circuit(*args):
grad_finitediff = (e_2 - e_1) / 0.0002

assert np.allclose(grad_autograd[0][0], grad_finitediff)


@pytest.mark.parametrize(
("hamiltonian", "result"),
[
(
qml.Hamiltonian(
np.array([0.5, 0.5]), [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)]
),
qml.Hamiltonian(np.array([1.0]), [qml.PauliX(0) @ qml.PauliY(1)]),
),
(
qml.Hamiltonian(
np.array([0.5, -0.5]),
[qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)],
),
qml.Hamiltonian([], []),
),
(
qml.Hamiltonian(
np.array([0.0, -0.5]),
[qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)],
),
qml.Hamiltonian(np.array([-0.5]), [qml.PauliX(0) @ qml.PauliZ(1)]),
),
(
qml.Hamiltonian(
np.array([0.25, 0.25, 0.25, -0.25]),
[
qml.PauliX(0) @ qml.PauliY(1),
qml.PauliX(0) @ qml.PauliZ(1),
qml.PauliX(0) @ qml.PauliY(1),
qml.PauliX(0) @ qml.PauliY(1),
],
),
qml.Hamiltonian(
np.array([0.25, 0.25]),
[qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)],
),
),
],
)
def test_simplify(hamiltonian, result):
r"""Test that simplify returns the correct hamiltonian."""
h = simplify(hamiltonian)
assert h.compare(result)
52 changes: 3 additions & 49 deletions tests/hf/test_tapering.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,16 +15,16 @@
Unit tests for functions needed for qubit tapering.
"""
import functools
import scipy
import pytest

import pennylane as qml
import pytest
import scipy
from pennylane import numpy as np
from pennylane.hf.tapering import (
_binary_matrix,
_kernel,
_observable_mult,
_reduced_row_echelon,
_simplify,
clifford,
generate_paulis,
generate_symmetries,
Expand Down Expand Up @@ -358,52 +358,6 @@ def test_observable_mult(obs_a, obs_b, result):
assert o.compare(result)


@pytest.mark.parametrize(
("hamiltonian", "result"),
[
(
qml.Hamiltonian(
np.array([0.5, 0.5]), [qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)]
),
qml.Hamiltonian(np.array([1.0]), [qml.PauliX(0) @ qml.PauliY(1)]),
),
(
qml.Hamiltonian(
np.array([0.5, -0.5]),
[qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliY(1)],
),
qml.Hamiltonian([], []),
),
(
qml.Hamiltonian(
np.array([0.0, -0.5]),
[qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)],
),
qml.Hamiltonian(np.array([-0.5]), [qml.PauliX(0) @ qml.PauliZ(1)]),
),
(
qml.Hamiltonian(
np.array([0.25, 0.25, 0.25, -0.25]),
[
qml.PauliX(0) @ qml.PauliY(1),
qml.PauliX(0) @ qml.PauliZ(1),
qml.PauliX(0) @ qml.PauliY(1),
qml.PauliX(0) @ qml.PauliY(1),
],
),
qml.Hamiltonian(
np.array([0.25, 0.25]),
[qml.PauliX(0) @ qml.PauliY(1), qml.PauliX(0) @ qml.PauliZ(1)],
),
),
],
)
def test_simplify(hamiltonian, result):
r"""Test that simplify returns the correct hamiltonian."""
h = _simplify(hamiltonian)
assert h.compare(result)


@pytest.mark.parametrize(
("generator", "paulix_ops", "result"),
[
Expand Down