Huffman Ternary Tree#

This notebook shows how to generate a Huffman-code Ternary Tree as described by:

Li, Q. S., Liu, H. Y., Wang, Q., Wu, Y. C., & Guo, G. P. (2025). Huffman-Code-based Ternary Tree Transformation. Chinese Physics Letters.

import ferrmion as fr
from ferrmion.encode.ternary_tree_node import TTNode
import json
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt

folder = Path.cwd().joinpath(Path("../../../python/tests/"))
with open(folder.joinpath("./data/h2o_sto-3g.json"), "rb") as file:
    data = json.load(file)
ones = np.array(data["ones"])
twos = np.array(data["twos"])

Majorana Frequency#

First we need to find the weighting of majorana operators, “frequency” in the language of huffman codes.

For some hamiltonian \(H= \sum_l h_l\), such as the Molecular Structure Hamiltonian: $\(H = \sum_{i,j} h_{ij}a^{\dagger}_i a_j + \sum_{i,j,k,l} h_{ijkl}a^{\dagger}_i a^{\dagger}_j a_k a_l \)$

We define the frequency of a majorana operator as:

\[F_j = \sum_l I_l(j)|h_l|\]

\(I_l(j)=1\) iff \(m_j \in H_l\)

n_modes = ones.shape[0]
majorana_freq = np.zeros(n_modes)

for i in range(n_modes):
    for j in range(n_modes):
        val = np.abs(ones[i, j])
        positions = {i, j}
        for p in positions:
            majorana_freq[p] += val

for i in range(n_modes):
    for j in range(n_modes):
        for k in range(n_modes):
            for l in range(n_modes):
                val = np.abs(twos[i, j, k, l])
                positions = {i, j, k, l}
                for p in positions:
                    majorana_freq[p] += val
majorana_freq = majorana_freq.repeat(2)

Build the Tree#

Now we can build the tree, using the huffman compression algorithm.

nodes = {i: None for i in range(len(majorana_freq))}
weights = {i: j for i, j in enumerate(majorana_freq)}
n_ops = len(majorana_freq)
for i in range(n_ops // 2):
    parent_index = 2 * n_ops - 1 - i
    mins = sorted(weights.items(), key=lambda kv: (kv[1], kv[0]))[:3]

    parent = nodes.get(parent_index, TTNode(parent=None, qubit_label=i))

    for min, child_string in zip(mins, ["x", "y", "z"][: len(mins)]):
        possible_child = nodes[min[0]]
        if isinstance(possible_child, TTNode):
            parent.add_child(which_child=child_string, child_node=possible_child)

    new_weight = 0
    for index, weight in mins:
        new_weight += weight
        weights.pop(index)
        nodes.pop(index)

    nodes[parent_index] = parent
    weights[parent_index] = new_weight

assert len(nodes) == 1
root_node = [*nodes.values()][0]
huffman_tree = fr.TernaryTree(n_modes=n_modes, root_node=root_node)
huffman_tree.enumeration_scheme = huffman_tree.default_enumeration_scheme()
huffman_tree.string_pairs
{'': ('xzz', 'yzzz'),
 'x': ('xxz', 'xyz'),
 'y': ('yyzz', 'yxz'),
 'xx': ('xxx', 'xxy'),
 'xy': ('xyy', 'xyx'),
 'xz': ('xzx', 'xzy'),
 'yx': ('yxy', 'yxx'),
 'yy': ('yyx', 'yyyz'),
 'yz': ('yzyz', 'yzxz'),
 'yyy': ('yyyy', 'yyyx'),
 'yyz': ('yyzx', 'yyzy'),
 'yzx': ('yzxy', 'yzxx'),
 'yzy': ('yzyx', 'yzyy'),
 'yzz': ('yzzy', 'yzzx')}
huffman_tree.enumeration_scheme
{'': (0, 13),
 'x': (1, 11),
 'y': (2, 12),
 'xx': (3, 5),
 'xy': (4, 6),
 'xz': (5, 7),
 'yx': (6, 8),
 'yy': (7, 9),
 'yz': (8, 10),
 'yyy': (9, 0),
 'yyz': (10, 1),
 'yzx': (11, 2),
 'yzy': (12, 3),
 'yzz': (13, 4)}
from ferrmion.visualise import draw_tt

draw_tt(root_node)
../_images/afa180f2d8a9d40219f7d1bebf5fef96ccf87b8bdd2d8117aba17794a0277541.png

Redefining 2e frequency#

Now that we have a tree structure, we need to assign fermionic modes to pairs of Pauli-strings derived from the structure of the tree.

We should account for how the operators interact, since we’re aiming to find the minimal coefficient Pauli-weight.

To do this, we redefine the two-electron term frequency.

With a two electron term in the form \(h_l = m_i m_j m_k m_l\)

\[F_{ij}=\sum_l I_l(ij)|h_l|\]

\(I_l(ij) = 1\) if \(m_i\) and \(m_j\) are in positons \((0,1)\) or \((2,3)\).

To simplify, we take the terms with \(i<j\).

n_modes = ones.shape[0]
two_e_freq = np.zeros(ones.shape)
print(two_e_freq.shape)
for j in range(n_modes):
    for i in range(n_modes):
        for l in range(n_modes):
            for k in range(n_modes):
                val = np.abs(twos[i, j, k, l])
                two_e_freq[i, j] += val
                two_e_freq[k, l] += val
two_e_freq = np.kron(two_e_freq, np.array([[1, 1], [1, 1]]))
two_e_freq = np.triu(two_e_freq, k=1)
(14, 14)
plt.matshow(two_e_freq)
<matplotlib.image.AxesImage at 0x110ac97f0>
../_images/368936120a2215afd0b2a842770cf0dcc977523b20a8786acf9f456ceae86af0.png

Pairing Majorana operators to Pauli-strings#

If we want to enforce vaccum state preservation (we do!), we require that we map majorana operators to paired strings.

This means we only need to use the freqencies of \(i\in 0,2,4,6...N-1\) and let \(j=i+1\)

Sorting Modes#

vaccum_frequencies = np.diag(two_e_freq, k=1)
plt.matshow([vaccum_frequencies])
sorted_pairs = [(i, i + 1) for i in np.argsort(vaccum_frequencies)[::-1] if i % 2 == 0]
sorted_modes = [i[0] // 2 for i in sorted_pairs]
../_images/a50a4bd682f51b060845f9e99085fc80cad2526c9ec5f5dd5a276a488e4dc308.png

Sorting Operators#

from ferrmion.utils import pauli_to_symplectic, symplectic_product, find_pauli_weight
weights = {}
for index, pair in enumerate(huffman_tree.string_pairs.values()):
    left, right = pair
    left = huffman_tree.branch_pauli_map[left]
    right = huffman_tree.branch_pauli_map[right]

    weights[index] = {}
    left, _ = pauli_to_symplectic(left, 0)
    right, _ = pauli_to_symplectic(right, 0)
    pair_weight = find_pauli_weight(np.array([left])) + find_pauli_weight(
        np.array([right])
    )
    _, product = symplectic_product(left, right)
    product_weight = find_pauli_weight(np.array([product]))

    weights[index]["pair_weight"] = pair_weight
    weights[index]["prod_weight"] = product_weight

    operator_order = sorted(
        weights.items(), key=lambda kv: (kv[1]["prod_weight"], kv[1]["pair_weight"])
    )

    operator_order = [index for index, _ in operator_order]
huffman_tree.string_pairs
{'': ('xzz', 'yzzz'),
 'x': ('xxz', 'xyz'),
 'y': ('yyzz', 'yxz'),
 'xx': ('xxx', 'xxy'),
 'xy': ('xyy', 'xyx'),
 'xz': ('xzx', 'xzy'),
 'yx': ('yxy', 'yxx'),
 'yy': ('yyx', 'yyyz'),
 'yz': ('yzyz', 'yzxz'),
 'yyy': ('yyyy', 'yyyx'),
 'yyz': ('yyzx', 'yyzy'),
 'yzx': ('yzxy', 'yzxx'),
 'yzy': ('yzyx', 'yzyy'),
 'yzz': ('yzzy', 'yzzx')}
huffman_tree.enumeration_scheme
{'': (0, 13),
 'x': (1, 11),
 'y': (2, 12),
 'xx': (3, 5),
 'xy': (4, 6),
 'xz': (5, 7),
 'yx': (6, 8),
 'yy': (7, 9),
 'yz': (8, 10),
 'yyy': (9, 0),
 'yyz': (10, 1),
 'yzx': (11, 2),
 'yzy': (12, 3),
 'yzz': (13, 4)}
operator_order = sorted(
    weights.items(), key=lambda kv: (kv[1]["prod_weight"], kv[1]["pair_weight"])
)
operator_order = [index for index, _ in operator_order]

Assigning a mode-operator map#

mode_op_map = [0] * len(sorted_modes)
for operator_index, mode_index in enumerate(sorted_modes):
    mode_op_map[mode_index] = operator_order[operator_index]
huffman_tree.default_mode_op_map = mode_op_map
print(huffman_tree.enumeration_scheme)
huffman_tree.enumeration_scheme = huffman_tree.default_enumeration_scheme()
print(huffman_tree.enumeration_scheme)
{'': (0, 13), 'x': (1, 11), 'y': (2, 12), 'xx': (3, 5), 'xy': (4, 6), 'xz': (5, 7), 'yx': (6, 8), 'yy': (7, 9), 'yz': (8, 10), 'yyy': (9, 0), 'yyz': (10, 1), 'yzx': (11, 2), 'yzy': (12, 3), 'yzz': (13, 4)}
{'': (0, 13), 'x': (1, 11), 'y': (2, 12), 'xx': (3, 5), 'xy': (4, 6), 'xz': (5, 7), 'yx': (6, 8), 'yy': (7, 9), 'yz': (8, 10), 'yyy': (9, 0), 'yyz': (10, 1), 'yzx': (11, 2), 'yzy': (12, 3), 'yzz': (13, 4)}
huffman_tree.default_enumeration_scheme()
{'': (0, 13),
 'x': (1, 11),
 'y': (2, 12),
 'xx': (3, 5),
 'xy': (4, 6),
 'xz': (5, 7),
 'yx': (6, 8),
 'yy': (7, 9),
 'yz': (8, 10),
 'yyy': (9, 0),
 'yyz': (10, 1),
 'yzx': (11, 2),
 'yzy': (12, 3),
 'yzz': (13, 4)}

Coefficient Pauli Weight#

The coefficient Puali-weight of an operator is the sum of each terms Pauli-weight times the norm of that terms coefficient:

\[H=\sum_i c_i h_i\]
\[W_{CP} = \sum_i |c_i| W_{P}(h_i)\]

This weight can be important for quantum algorithms which rely on samplying hamiltonian terms or measurement results.

Now we can find out how well the Huffman-tree reduces the coefficient Pauli-weight.

from ferrmion.optimize import coefficient_pauli_weight

fham = fr.molecular_hamiltonian(ones, twos, constant_energy=0, physicist_notation=True)
qham = huffman_tree.encode(fham)
coefficient_pauli_weight(qham)
[np.float64(301.7202300406311)]
encoding_funcs = [fr.JordanWigner, fr.ParityEncoding, fr.BravyiKitaev, fr.JKMN]
for encoding in encoding_funcs:
    fham = fr.molecular_hamiltonian(ones, twos, constant_energy=0, physicist_notation=True)
    qham = encoding(fham.n_modes).encode(fham)
    pw = coefficient_pauli_weight(qham)
    print(encoding.__name__, pw)
JordanWigner [np.float64(191.34214293107775)]
ParityEncoding [np.float64(256.8019697995479)]
BravyiKitaev [np.float64(292.03224577227707)]
JKMN [np.float64(277.8229667626987)]

Inbuilt Function#

ferrmion provides a method to generate vaccum-preserving Huffman-code trees by providing one and two electron coefficents.

from ferrmion.optimize import huffman_ternary_tree

inbuilt_huffman = huffman_ternary_tree(ones, twos)
draw_tt(
    inbuilt_huffman,
    enumeration_scheme=inbuilt_huffman.enumeration_scheme,
    type="spaced",
)

inbuilt_ham = fr.molecular_hamiltonian(
    ones, twos, constant_energy=0, physicist_notation=True
)
coefficient_pauli_weight(inbuilt_huffman.encode(inbuilt_ham))
[np.float64(205.52406742880444)]
../_images/4df75bd0376b3d221ebe51014677d269e0b5d599db686f95166d8f45a1dd0054.png
inbuilt_huffman.default_enumeration_scheme()
{'': (0, 13),
 'x': (1, 11),
 'y': (2, 12),
 'xx': (3, 5),
 'xy': (4, 6),
 'xz': (5, 7),
 'yx': (6, 8),
 'yy': (7, 9),
 'yz': (8, 10),
 'yyy': (9, 0),
 'yyz': (10, 1),
 'yzx': (11, 2),
 'yzy': (12, 3),
 'yzz': (13, 4)}
inbuilt_huffman.as_dict() == huffman_tree.as_dict()
True