Fix multiset ordering bug in perturbation module

qiskit_dynamics/perturbation/dyson_magnus.py

@@ -804,13 +804,13 @@ def _get_q_term_list(complete_index_multisets: List[Multiset]) -> List:
 def _get_dyson_lmult_rule(
     complete_index_multisets: List[Multiset], perturbation_labels: Optional[List[Multiset]] = None
 ) -> List:
-    """Given a complete list of index multisets, return
-    the lmult rule in the format required for ``CustomProduct``.
-    Note, the generator :math:`G(t)` is encoded as index ``-1``, as
-    it will be prepended to the list of A matrices.
+    """Given a complete list of index multisets, return the lmult rule in the format required for
+    ``CustomProduct``. Note, the generator :math:`G(t)` is encoded as index ``-1``, as it will be
+    prepended to the list of A matrices. Similarly the index of the solution to the LMDE with
+    generator :math:`G(t)` is encoded as ``-1``.
 
-    While not required within the logic of this function, the input
-    should be canonically ordered according to ``_get_all_submultisets``.
+    While not required within the logic of this function, the input should be canonically ordered
+    according to ``_get_all_submultisets``.
 
     Args:
         complete_index_multisets: List of complete multisets.
@@ -834,8 +834,15 @@ def _get_dyson_lmult_rule(
     lmult_rule = [(np.array([1.0]), np.array([[-1, -1]]))]
 
     for term_idx, term in enumerate(complete_index_multisets):
-        if len(term) == 1:
-            lmult_rule.append((np.array([1.0, 1.0]), np.array([[-1, term_idx], [term_idx, -1]])))
+        if len(term) == 1 and term in perturbation_labels:
+            # if term not in perturbation_labels for len(term) == 1 the corresponding Dyson integral
+            # will always be 0
+            lmult_rule.append(
+                (
+                    np.array([1.0, 1.0]),
+                    np.array([[-1, term_idx], [perturbation_labels.index(term), -1]]),
+                )
+            )
         else:
             # self multiplied by base generator
             lmult_indices = [[-1, term_idx]]

qiskit_dynamics/perturbation/multiset_utils.py

@@ -15,7 +15,7 @@
 Utility functions for working with multisets.
 """
 
-from typing import List, Optional, Tuple
+from typing import List, Optional, Tuple, Iterable
 import itertools
 
 from multiset import Multiset
@@ -46,15 +46,51 @@ def _multiset_to_sorted_list(multiset: Multiset) -> List:
     return sorted_list
 
 
-def _sorted_multisets(multisets: List[Multiset]) -> List[Multiset]:
-    """Sort in non-decreasing order according to:
-
-    ms1 <= ms2 if len(ms1) < len(ms2), or if
-    len(ms1) == len(ms2) and if
-    str(_multiset_to_sorted_list(ms1)) <= str(_multiset_to_sorted_list(ms2)).
+class _MultisetSortKey:
+    """Dummy class for usage as a key when sorting Multiset instances. This assumes the elements
+    of the multisets can themselves be sorted.
     """
 
-    return sorted(multisets, key=lambda x: str(len(x)) + ", " + str(_multiset_to_sorted_list(x)))
+    __slots__ = ("multiset",)
+
+    def __init__(self, multiset: Multiset):
+        self.multiset = multiset
+
+    def __lt__(self, other: Multiset) -> bool:
+        """Implements an ordering on multisets.
+
+        This orders first according to length (the number of elements in each multiset). If ``self``
+        and ``other`` are the same length, ``self < other`` if, when written as fully expanded and
+        sorted lists, ``self < other`` in lexicographic ordering. E.g. it holds that ``Multiset({0:
+        2, 1: 1}) < Multiset({0: 1, 1: 2})``, as the list versions are ``x = [0, 0, 1]``, and ``y =
+        [0, 1, 1]``. Here ``x[0] == y[0]``, but ``x[1] < y[1]``, and hence ``x < y`` in this
+        ordering.
+        """
+        if len(self.multiset) < len(other.multiset):
+            return True
+
+        if len(other.multiset) < len(self.multiset):
+            return False
+
+        unique_entries = set(self.multiset.distinct_elements())
+        unique_entries.update(other.multiset.distinct_elements())
+        unique_entries = sorted(unique_entries)
+
+        for element in unique_entries:
+            self_count = self.multiset[element]
+            other_count = other.multiset[element]
+
+            if self_count != other_count:
+                return self_count > other_count
+
+        return False
+
+
+def _sorted_multisets(multisets: Iterable[Multiset]) -> List[Multiset]:
+    """Sort in non-decreasing order according to the ordering described in the dummy class
+    _MultisetSort.
+    """
+    return sorted(multisets, key=_MultisetSortKey)
 
 
 def _clean_multisets(multisets: List[Multiset]) -> List[Multiset]:

releasenotes/notes/multiset-order-bug-fix-1f1603ee1e230cba.yaml

@@ -0,0 +1,6 @@
+---
+fixes:
+  - |
+    Fixes a bug in the perturbation module with internal sorting of ``Multiset`` instances, which
+    caused incorrect computation of perturbation theory terms when ``>10`` perturbations are
+    present.

test/dynamics/perturbation/test_dyson_magnus.py

@@ -555,6 +555,31 @@ def test__get_dyson_lmult_rule_power_series_case1(self):
             perturbation_labels=perturbation_labels,
         )
 
+    def test__get_dyson_lmult_rule_power_series_case1_missing(self):
+        """Test _get_dyson_lmult_rule for higher order terms in the generator
+        decomposition case 1 where there is no corresponding perturbation_label term for
+        a desired expansion term.
+        """
+
+        expansion_labels = [[0], [1], [0, 1], [1, 1], [0, 1, 1]]
+        expansion_labels = [Multiset(label) for label in expansion_labels]
+        perturbation_labels = [[0], [0, 1], [1, 1]]
+        perturbation_labels = [Multiset(label) for label in perturbation_labels]
+        expected_lmult_rule = [
+            (np.ones(1, dtype=float), np.array([[-1, -1]])),
+            (np.ones(2, dtype=float), np.array([[-1, 0], [0, -1]])),
+            (np.ones(1, dtype=float), np.array([[-1, 1]])),
+            (np.ones(3, dtype=float), np.array([[-1, 2], [0, 1], [1, -1]])),
+            (np.ones(2, dtype=float), np.array([[-1, 3], [2, -1]])),
+            (np.ones(4, dtype=float), np.array([[-1, 4], [0, 3], [1, 1], [2, 0]])),
+        ]
+
+        self._test__get_dyson_lmult_rule(
+            expansion_labels,
+            expected_lmult_rule,
+            perturbation_labels=perturbation_labels,
+        )
+
     def test__get_dyson_lmult_rule_power_series_case2(self):
         """Test _get_dyson_lmult_rule for higher order terms in the generator
         decomposition case 2.
@@ -579,6 +604,36 @@ def test__get_dyson_lmult_rule_power_series_case2(self):
             perturbation_labels=perturbation_labels,
         )
 
+    def test__get_dyson_lmult_rule_power_series_case2_unordered(self):
+        """Test _get_dyson_lmult_rule for higher order terms in the generator
+        decomposition case 2, with the perturbation_labels being non-canonically ordered.
+
+        Note that the conversion of case2 to case2_unordered requires relabelling 1 <-> 3
+        in expected_lmult_rule, but also changing the ordering of the pairs in the matrix-product
+        description of expected_lmult_rule, as the rule is constructed by iterating through the
+        entries of perturbation_labels. The matrix-product description must be ordered by the first
+        index.
+        """
+
+        expansion_labels = [[0], [1], [0, 1], [1, 1], [0, 1, 1]]
+        expansion_labels = [Multiset(label) for label in expansion_labels]
+        perturbation_labels = [[0], [0, 1], [2], [1]]
+        perturbation_labels = [Multiset(label) for label in perturbation_labels]
+        expected_lmult_rule = [
+            (np.ones(1, dtype=float), np.array([[-1, -1]])),
+            (np.ones(2, dtype=float), np.array([[-1, 0], [0, -1]])),
+            (np.ones(2, dtype=float), np.array([[-1, 1], [3, -1]])),
+            (np.ones(4, dtype=float), np.array([[-1, 2], [0, 1], [1, -1], [3, 0]])),
+            (np.ones(2, dtype=float), np.array([[-1, 3], [3, 1]])),
+            (np.ones(4, dtype=float), np.array([[-1, 4], [0, 3], [1, 1], [3, 2]])),
+        ]
+
+        self._test__get_dyson_lmult_rule(
+            expansion_labels,
+            expected_lmult_rule,
+            perturbation_labels=perturbation_labels,
+        )
+
     def test__get_dyson_lmult_rule_case1(self):
         """Test __get_dyson_lmult_rule case 1."""
         expansion_labels = [[0], [1], [0, 1], [1, 1], [0, 1, 1]]
@@ -652,11 +707,10 @@ def test__get_dyson_lmult_rule_case2(self):
         self._test__get_dyson_lmult_rule(expansion_labels, expected_lmult_rule)
 
     def _test__get_dyson_lmult_rule(
-        self, complete_symmetric_dyson_term_list, expected, perturbation_labels=None
+        self, complete_index_multisets, expected, perturbation_labels=None
     ):
-        """Run a test case for _get_symmetric_dyson_mult_rules."""
-        lmult_rule = _get_dyson_lmult_rule(complete_symmetric_dyson_term_list, perturbation_labels)
-
+        """Run a test case for _get_dyson_lmult_rule."""
+        lmult_rule = _get_dyson_lmult_rule(complete_index_multisets, perturbation_labels)
         self.assertMultRulesEqual(lmult_rule, expected)
 
     def assertMultRulesEqual(self, rule1, rule2):

test/dynamics/perturbation/test_multiset_utils.py

@@ -59,6 +59,15 @@ def test_case1(self):
         expected = [Multiset([1]), Multiset([0, 2]), Multiset([0, 0, 1]), Multiset([0, 1, 1])]
         self.assertTrue(output == expected)
 
+    def test_case2(self):
+        """Test case which would have caught a bug introduced by a previous implementation of the
+        ordering.
+        """
+        multisets = [Multiset([10]), Multiset([1])]
+        output = _sorted_multisets(multisets)
+        expected = [Multiset([1]), Multiset([10])]
+        self.assertTrue(output == expected)
+
 
 class TestToSortedList(QiskitDynamicsTestCase):
     """Test conversion to sorted list."""

test/dynamics/perturbation/test_solve_lmde_perturbation.py

@@ -501,6 +501,81 @@ def A1(t):
         self.assertAllClose(expected_D0001, results.perturbation_data.get_item([0, 0, 0, 1])[-1])
         self.assertAllClose(expected_D0011, results.perturbation_data.get_item([0, 0, 1, 1])[-1])
 
+    def test_dyson_analytic_case1_1d_relabeled(self):
+        """This is the same numerical test as test_dyson_analytic_case1_1d, however it relabels
+        the perturbation indices as 0 -> 1, and 1 -> 10. This is an integration test that would
+        have caught multiset sorting bug that broke solve_lmde_perturbation results computation.
+        """
+
+        def generator(t):
+            return Array([[1, 0], [0, 1]], dtype=complex).data
+
+        def A0(t):
+            return Array([[0, t], [t**2, 0]], dtype=complex).data
+
+        def A1(t):
+            return Array([[t, 0], [0, t**2]], dtype=complex).data
+
+        T = np.pi * 1.2341
+
+        results = solve_lmde_perturbation(
+            perturbations=[A0, A1],
+            t_span=[0, T],
+            generator=generator,
+            y0=np.array([0.0, 1.0], dtype=complex),
+            expansion_method="dyson",
+            expansion_order=2,
+            expansion_labels=[[1, 1, 10], [1, 1, 1, 10], [1, 1, 10, 10]],
+            perturbation_labels=[[1], [10]],
+            dyson_in_frame=False,
+            integration_method=self.integration_method,
+            atol=1e-13,
+            rtol=1e-13,
+        )
+
+        T2 = T**2
+        T3 = T * T2
+        T4 = T * T3
+        T5 = T * T4
+        T6 = T * T5
+        T7 = T * T6
+        T8 = T * T7
+        T9 = T * T8
+        T10 = T * T9
+        T11 = T * T10
+
+        U = expm(np.array(generator(0)) * T)
+
+        expected_D0 = U @ np.array([[T2 / 2], [0]], dtype=complex)
+        expected_D1 = U @ np.array([[0], [T3 / 3]], dtype=complex)
+        expected_D00 = U @ np.array([[0], [T5 / 10]], dtype=complex)
+        expected_D01 = U @ (
+            np.array([[T5 / 15], [0]], dtype=complex) + np.array([[T4 / 8], [0]], dtype=complex)
+        )
+        expected_D11 = U @ np.array([[0], [T6 / 18]], dtype=complex)
+        expected_D001 = U @ (
+            np.array([[0], [T8 / 120]], dtype=complex)
+            + np.array([[0], [T7 / 56]], dtype=complex)
+            + np.array([[0], [T8 / 80]], dtype=complex)
+        )
+        expected_D0001 = U @ np.array([[(T10 / 480) + (T9 / 280)], [0]], dtype=complex)
+        expected_D0011 = U @ np.array(
+            [
+                [0],
+                [(T11 / 396) + 23 * (T10 / 8400) + (T9 / 432)],
+            ],
+            dtype=complex,
+        )
+
+        self.assertAllClose(expected_D0, results.perturbation_data.get_item([1])[-1])
+        self.assertAllClose(expected_D1, results.perturbation_data.get_item([10])[-1])
+        self.assertAllClose(expected_D00, results.perturbation_data.get_item([1, 1])[-1])
+        self.assertAllClose(expected_D01, results.perturbation_data.get_item([1, 10])[-1])
+        self.assertAllClose(expected_D11, results.perturbation_data.get_item([10, 10])[-1])
+        self.assertAllClose(expected_D001, results.perturbation_data.get_item([1, 1, 10])[-1])
+        self.assertAllClose(expected_D0001, results.perturbation_data.get_item([1, 1, 1, 10])[-1])
+        self.assertAllClose(expected_D0011, results.perturbation_data.get_item([1, 1, 10, 10])[-1])
+
     def test_dyson_analytic_case1(self):
         """Analytic test of computing dyson terms.