Repository: quantumlib/OpenFermion
Pull request #506 (2019-02-19 17:13:43)

Wedge product code between two tensors
The two tensors are the moments for a particular sequence of
fermionic ladder operators.

Symmetric vs anitysmmetric could be an easy extension to this code.


Patch:
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index 07a5565..8f6454c 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -76,6 +76,10 @@ from ._unitary_cc import (uccsd_convert_amplitude_format,
                           uccsd_singlet_get_packed_amplitudes,
                           uccsd_singlet_paramsize)
 
+from ._wedge_product import (generate_parity_permutations,
+                             wedge)
+
+
 # Imports out of alphabetical order to avoid circular dependency.
 from ._jellium_hf_state import hartree_fock_state_jellium
 
diff --git a/src/openfermion/utils/_wedge_product.py b/src/openfermion/utils/_wedge_product.py
new file mode 100644
index 0000000..65a6fb5
--- /dev/null
+++ b/src/openfermion/utils/_wedge_product.py
@@ -0,0 +1,143 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+"""This module contains methods to lift tensors to higher spaces
+through the Grassmann wedge product."""
+from typing import Sequence, Union
+from itertools import product
+from math import factorial
+import copy
+import numpy
+
+# numpy einsum alphabet
+EINSUM_CHARS = 'abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ'
+
+
+def generate_parity_permutations(seq):
+    """
+    Generates the permutations and sign of a sequence by constructing a tree
+    where the nth level contains all n-permutations of m (n < m) objects.
+
+    At the last level where n == m all permutations are generated.  The sign
+    is kept updated by determining where the next number is inserted into the
+    current leaf's set of numbers.
+
+    Example:
+         Constructing the permutations of the sequence [A, B, C] constructs the
+         following tree:
+
+         [[(A, +1)], [(AB, +1), (BA, -1)], [(ABC, +1), (ACB, -1), (CAB, +1],
+         [(BAC, -1), (BCA, +1), (CBA, -1)]]
+
+    Args:
+        seq: a sequence of a string to provide permutations
+
+    returns:
+        a permutation list with the elements of the seq permuted and a sign
+        associated with the permutation.
+    """
+    if isinstance(seq, str):
+        seq = [x for x in seq]
+
+    indices = seq[1:]
+    permutations = [([seq[0]], 1)]
+    while indices:
+        index_to_inject = indices.pop(0)
+
+        new_permutations = []  # permutations in the tree
+        for perm in permutations:
+            # now loop over positions to insert
+            for put_index in range(len(perm[0]) + 1):
+                new_index_list = copy.deepcopy(perm[0])
+                # insert new object starting at end of the list
+                new_index_list.insert(len(perm[0]) - put_index, index_to_inject)
+
+                new_permutations.append(
+                    (new_index_list, perm[1] * (-1) ** (put_index)))
+
+        permutations = new_permutations
+
+    return permutations
+
+
+def wedge(left_tensor, right_tensor, left_index_ranks, right_index_ranks):
+    """
+    Implement the wedge product between left_tensor and right_tensor
+
+    The wedge product is defined as
+
+    .. math::
+        \\begin{align}
+        a_{j_{1}, j_{2}, ...,j_{p}}^{i_{1}, i_{2}, ..., i_{p}} \\wedge
+        b_{j_{p+1}, j_{p+2}, ..., j_{N}}^{i_{p+1}, i_{p + 2}, ..., i_{N}} =
+        \\left(\\frac{1}{N!}\\right)^{2} = \\sum_{\\pi, \\sigma}\\epsilon(\\pi)
+        \\epsilon(\\sigma)\\pi \\sigma
+        a_{j_{1}, j_{2}, ...,j_{p}}^{i_{1}, i_{2}, ..., i_{p}}
+        b_{j_{p+1}, j_{p+2}, ..., j_{N}}^{i_{p+1}, i_{p + 2}, ..., i_{N}}
+        \\end{align}
+
+    The top indices are those that transform contravariently.  The bottom
+    indices transform covariently.
+
+    The tensor storage convention for marginals follows the OpenFermion
+    convention. tpdm[i, j, k, l] = <i^ j^ k l>,
+    rtensor[u1, u2, u3, d1] = <u1^ u2^ u3^ d1>
+
+    Args:
+        left_tensor: left tensor to wedge product
+        right_tensor: right tensor to wedge product
+        left_index_ranks: tuple of number of indices that transform
+                          contravariently and covariently
+        right_index_ranks: tuple of number of indices that transform
+                           contravariently and covariently
+    Returns:
+        new tensor constructed as the wedge product of the left_tensor and
+        right_tensor
+    """
+    if left_tensor.ndim != sum(left_index_ranks):
+        raise IndexError(
+            "n_tensor shape is not consistent with the input n_index rank")
+    if right_tensor.ndim != sum(right_index_ranks):
+        raise IndexError(
+            "n_tensor shape is not consistent with the input n_index rank")
+    # assign upper and lower indices for n_tensor
+    total_upper = left_index_ranks[0] + right_index_ranks[0]
+    total_lower = left_index_ranks[1] + right_index_ranks[1]
+    upper_characters = EINSUM_CHARS[:total_upper]
+    lower_characters = EINSUM_CHARS[total_upper:total_upper + total_lower]
+    new_tensor = numpy.zeros(left_tensor.shape + right_tensor.shape,
+                             dtype=complex)
+    ordered_einsum_string = upper_characters + lower_characters
+
+    for upper_order_parities, lower_order_parities in product(
+            generate_parity_permutations(upper_characters),
+            generate_parity_permutations(lower_characters[::-1])):
+        # we reverse the order in the lower_chars so because
+        # <a^ b^ c d> = D_{dc}^{ab} in this code.
+        n_upper_einsum_chars = upper_order_parities[0][:left_index_ranks[0]]
+        m_upper_einsum_chars = upper_order_parities[0][left_index_ranks[0]:]
+        n_lower_einsum_chars = lower_order_parities[0][:left_index_ranks[1]][::-1]
+        m_lower_einsum_chars = lower_order_parities[0][left_index_ranks[1]:][::-1]
+
+        n_string = "".join(n_upper_einsum_chars + n_lower_einsum_chars)
+        m_string = "".join(m_upper_einsum_chars + m_lower_einsum_chars)
+
+        # we are doing lots of extra += operations but with the benefit of not
+        # having to write a python loop over the entire new_tensor object.
+        new_tensor += upper_order_parities[1] * lower_order_parities[1] * \
+                      numpy.einsum('{},{}->{}'.format(n_string, m_string,
+                                                      ordered_einsum_string),
+                                   left_tensor, right_tensor)
+
+    new_tensor /= factorial(total_upper) * factorial(total_lower)
+
+    return new_tensor


Test patch:
diff --git a/src/openfermion/utils/_wedge_product_test.py b/src/openfermion/utils/_wedge_product_test.py
new file mode 100644
index 0000000..a233b33
--- /dev/null
+++ b/src/openfermion/utils/_wedge_product_test.py
@@ -0,0 +1,145 @@
+import pytest
+import numpy
+from itertools import product
+from openfermion.utils._wedge_product import (wedge,
+                                              generate_parity_permutations)
+
+
+def test_upper_lower_index_size_check():
+    n_tensor = numpy.zeros((2, 2, 2))
+    m_tensor = numpy.zeros((2, 2))
+    with pytest.raises(IndexError):
+        wedge(n_tensor, m_tensor, (3, 1), (1, 1))
+
+    with pytest.raises(IndexError):
+        wedge(n_tensor, m_tensor, (4, 0), (1, 2))
+
+    with pytest.raises(IndexError):
+        wedge(n_tensor, m_tensor, (2, 1), (1, 2))
+
+
+def test_eye_wedge():
+    """
+    Test wedge product code through identity wedge
+
+    (p, s)(q, r) - (q, s)(p, r) - (p, r)(q, s) + (q, r)(p, s)
+    """
+    n_tensor = numpy.eye(2)
+    m_tensor = numpy.eye(2)
+    true_eye2 = numpy.zeros((2, 2, 2, 2))
+    for p, q, r, s in product(range(2), repeat=4):
+        true_eye2[p, q, r, s] += (n_tensor[p, s] * m_tensor[q, r] -
+                                  n_tensor[q, s] * m_tensor[p, r] -
+                                  n_tensor[p, r] * m_tensor[q, s] +
+                                  n_tensor[q, r] * m_tensor[p, s])
+    true_eye2 /= 2 ** 2
+    test_eye2 = wedge(n_tensor, m_tensor, (1, 1), (1, 1))
+    assert numpy.allclose(test_eye2, true_eye2)
+
+
+def test_random_wedge():
+    """
+    Test wedge product code through identity wedge
+
+    (p, s)(q, r) - (q, s)(p, r) - (p, r)(q, s) + (q, r)(p, s)
+    """
+    dim = 4
+    n_tensor = numpy.random.random((dim, dim))
+    m_tensor = numpy.random.random((dim, dim))
+    true_eye2 = numpy.zeros(tuple([dim] * 4))
+    for p, q, r, s in product(range(dim), repeat=4):
+        true_eye2[p, q, r, s] += (n_tensor[p, s] * m_tensor[q, r] -
+                                  n_tensor[q, s] * m_tensor[p, r] -
+                                  n_tensor[p, r] * m_tensor[q, s] +
+                                  n_tensor[q, r] * m_tensor[p, s])
+    true_eye2 /= 2 ** 2
+    test_eye2 = wedge(n_tensor, m_tensor, (1, 1), (1, 1))
+    assert numpy.allclose(test_eye2, true_eye2)
+
+
+def test_random_two_wedge():
+    """
+    Test wedge product code with random 4-tensor and a random 2-tensor
+
+    (a, b, c, d) wedge (e, f) -> (a, b, e, f, c, d)
+
+    """
+    dim = 3
+    n_tensor = numpy.random.random((dim, dim, dim, dim))
+    m_tensor = numpy.random.random((dim, dim))
+    true_tensor = numpy.zeros(tuple([dim] * 6))
+    for a, b, c, d, e, f in product(range(dim), repeat=6):
+        for u_perm, u_phase in generate_parity_permutations([a, b, c]):
+            for l_perm, l_phase in generate_parity_permutations([f, e, d]):
+                true_tensor[a, b, c, d, e, f] += u_phase * l_phase * n_tensor[
+                    u_perm[0], u_perm[1], l_perm[1], l_perm[0]] * m_tensor[
+                                                     u_perm[2], l_perm[2]]
+
+    true_tensor /= 6 ** 2
+    test_tensor = wedge(n_tensor, m_tensor, (2, 2), (1, 1))
+    assert numpy.allclose(test_tensor, true_tensor)
+
+
+def test_random_two_wedge_two():
+    """
+    Test wedge product code with a random 4-tensor and a random 4-tensor
+
+    (a, b, h, g) wedge (c, d, f, e) -> (a, b, c, d, e, f, g, h)
+
+    """
+    dim = 3
+    n_tensor = numpy.random.random((dim, dim, dim, dim))
+    m_tensor = numpy.random.random((dim, dim, dim, dim))
+    true_tensor = numpy.zeros(tuple([dim] * 8))
+    for a, b, c, d, e, f, g, h in product(range(dim), repeat=8):
+        for u_perm, u_phase in generate_parity_permutations([a, b, c, d]):
+            for l_perm, l_phase in generate_parity_permutations([h, g, f, e]):
+                # D_{hgfe}^{abcd} = C [D_{hg}^{ab}D_{fe}^{cd} +
+                # -1 * D_{hg}^{ab}D_{ef}^{cd} + ...]
+                true_tensor[a, b, c, d, e, f, g, h] += u_phase * l_phase * \
+                                                       n_tensor[
+                                                           u_perm[0], u_perm[1],
+                                                           l_perm[1], l_perm[
+                                                               0]] * m_tensor[
+                                                           u_perm[2], u_perm[3],
+                                                           l_perm[3], l_perm[2]]
+
+    true_tensor /= (4 * 3 * 2) ** 2
+    test_tensor = wedge(n_tensor, m_tensor, (2, 2), (2, 2))
+    assert numpy.allclose(test_tensor, true_tensor)
+
+
+def test_random_1_2_wedge_1_1():
+    dim = 3
+    n_tensor = numpy.random.random((dim, dim, dim))
+    m_tensor = numpy.random.random((dim, dim))
+    true_tensor = numpy.zeros(tuple([dim] * 5))
+    for a, b, c, d, e in product(range(dim), repeat=5):
+        for u_perm, u_phase in generate_parity_permutations([a, b]):
+            for l_perm, l_phase in generate_parity_permutations([e, d, c]):
+                # D_{edc}^{ab} = C[D_{ed}^{a}D_{c}^{b} - D_{ed}^{b}D_{c}^{a} -
+                # D_{ec}^{a}D_{d}^{b} +...]
+                true_tensor[a, b, c, d, e] += u_phase * l_phase * n_tensor[
+                    u_perm[0], l_perm[2], l_perm[1]] * m_tensor[
+                                                  u_perm[1], l_perm[0]]
+
+    true_tensor /= (2 * 3 * 2)
+    test_tensor = wedge(n_tensor, m_tensor, (1, 2), (1, 1))
+    assert numpy.allclose(test_tensor, true_tensor)
+
+
+def test_random_2_1_wedge_1_1():
+    dim = 3
+    n_tensor = numpy.random.random((dim, dim, dim))
+    m_tensor = numpy.random.random((dim, dim))
+    true_tensor = numpy.zeros(tuple([dim] * 5))
+    for a, b, c, d, e in product(range(dim), repeat=5):
+        for u_perm, u_phase in generate_parity_permutations([a, b, c]):
+            for l_perm, l_phase in generate_parity_permutations([e, d]):
+                true_tensor[a, b, c, d, e] += u_phase * l_phase * n_tensor[
+                    u_perm[0], u_perm[1], l_perm[0]] * m_tensor[
+                                                  u_perm[2], l_perm[1]]
+
+    true_tensor /= (2 * 3 * 2)
+    test_tensor = wedge(n_tensor, m_tensor, (2, 1), (1, 1))
+    assert numpy.allclose(test_tensor, true_tensor)