nate's suggestions + default built in for types + concatenate block i…

…nstead of copying np array

mitiq/lre/inference/multivariate_richardson.py

@@ -9,7 +9,7 @@
 
 import warnings
 from itertools import product
-from typing import Any, List, Optional, Tuple
+from typing import Any, Optional
 
 import numpy as np
 from cirq import Circuit
@@ -22,7 +22,7 @@
 
 def _full_monomial_basis_term_exponents(
     num_layers: int, degree: int
-) -> List[Tuple[int, ...]]:
+) -> list[tuple[int, ...]]:
     """Exponents of monomial terms required to create the sample matrix."""
     exponents = {
         exps
@@ -83,15 +83,15 @@ def sample_matrix(
     variable_exp = _full_monomial_basis_term_exponents(num_layers, degree)
     sample_matrix = np.empty((len(variable_exp), len(variable_exp)))
 
-    for rows, i in enumerate(scale_factor_vectors):
-        for cols, j in enumerate(variable_exp):
+    for i, scale_factors in enumerate(scale_factor_vectors):
+        for j, exponent in enumerate(variable_exp):
             evaluated_terms = []
-            for base, exp in zip(list(i), j):
+            for base, exp in zip(scale_factors, exponent):
                 # raise scale factor value by the exponent dict value
                 evaluated_terms.append(base**exp)
             # multiply both elements in the list to create an evaluated
             # monomial term
-            sample_matrix[rows, cols] = np.prod(evaluated_terms)
+            sample_matrix[i, j] = np.prod(evaluated_terms)
 
     return sample_matrix
 
@@ -101,10 +101,11 @@ def linear_combination_coefficients(
     degree: int,
     fold_multiplier: int,
     num_chunks: Optional[int] = None,
-) -> List[np.float64]:
+) -> list[float]:
     r"""
-    Defines the sample matrix required for multivariate extrapolation as
-    defined in :cite:`Russo_2024_LRE`.
+    Defines the function to find the linear combination coefficients from the
+    sample matrix as required for multivariate extrapolation (defined in
+    :cite:`Russo_2024_LRE`).
 
     Args:
         input_circuit: Circuit to be scaled.
@@ -118,11 +119,6 @@ def linear_combination_coefficients(
         List of the evaluated monomial basis terms using the scale factor
             vectors.
     """
-    num_layers = len(
-        _get_scale_factor_vectors(
-            input_circuit, degree, fold_multiplier, num_chunks
-        )
-    )
     input_sample_matrix = sample_matrix(
         input_circuit, degree, fold_multiplier, num_chunks
     )
@@ -138,7 +134,7 @@ def linear_combination_coefficients(
         sign, logdet = np.linalg.slogdet(  # pragma: no cover
             input_sample_matrix
         )
-        det = np.exp(logdet)  # pragma: no cover
+        det = sign * np.exp(logdet)  # pragma: no cover
 
     if np.isinf(det):
         raise ValueError(  # pragma: no cover
@@ -149,11 +145,26 @@ def linear_combination_coefficients(
     assert det != 0.0
 
     coeff_list = []
-    for i in range(num_layers):
-        # taken from https://stackoverflow.com/a/52866044
-        sample_matrix_minor = np.delete(
-            np.delete(input_sample_matrix, i, axis=0), i, axis=1
-        )
-        coeff_list.append(np.linalg.det(sample_matrix_minor) / det)
+    mat_row, mat_cols = input_sample_matrix.shape
+    assert mat_row == mat_cols
+    # replace a row of the sample matrix with [1, 0, 0, .., 0]
+    repl_row = np.array([[1] + [0] * (mat_cols - 1)])
+    for i in range(mat_row):
+        if i == 0:  # first row
+            new_mat = np.concatenate(
+                (repl_row, input_sample_matrix[1:]), axis=0
+            )
+        elif i == mat_row - 1:  # last row
+            new_mat = np.concatenate(
+                (input_sample_matrix[:i], repl_row), axis=0
+            )
+        else:
+            frst_sl = np.concatenate(
+                (input_sample_matrix[:i], repl_row), axis=0
+            )
+            sec_sl = input_sample_matrix[i + 1 :]
+            new_mat = np.concatenate((frst_sl, sec_sl), axis=0)
+
+        coeff_list.append(np.linalg.det(new_mat) / det)
 
     return coeff_list

mitiq/lre/tests/test_multivariate_richardson.py

@@ -68,7 +68,8 @@ def test_basis_exp(test_num_layers, test_degree, expected):
 @pytest.mark.parametrize(
     "test_num_layers, test_degree",
     [(1, 1), (2, 2), (3, 2), (10, 4)],
-    # Note need to add (100, 2) here. This makes the unit test very slow.
+    # TO DO: Note need to add (100, 2) here.
+    # This makes the unit test very slow.
 )
 def test_basis_exp_len(test_num_layers, test_degree):
     calc_dict = _full_monomial_basis_term_exponents(
@@ -116,9 +117,9 @@ def test_basis_exp_len(test_num_layers, test_degree):
     ],
 )
 def test_sample_matrix(test_circ, test_degree, expected_matrix):
-    assert (
-        expected_matrix - sample_matrix(test_circ, test_degree, 1)
-    ).all() <= 1e-3
+    assert np.allclose(
+        expected_matrix, sample_matrix(test_circ, test_degree, 1), atol=1e-3
+    )
 
 
 @pytest.mark.parametrize(
@@ -157,16 +158,21 @@ def test_sample_matrix(test_circ, test_degree, expected_matrix):
     ],
 )
 def test_coeffs(test_circ, test_degree, test_fold_multiplier, expected_matrix):
-    assert (
-        abs(
-            np.array(expected_matrix)
-            - np.array(
-                linear_combination_coefficients(
-                    test_circ, test_degree, test_fold_multiplier
-                )
+    assert np.allclose(
+        expected_matrix,
+        linear_combination_coefficients(
+            test_circ, test_degree, test_fold_multiplier
+        ),
+        atol=1e-3,
+    )
+
+    assert np.isclose(
+        sum(
+            linear_combination_coefficients(
+                test_circ, test_degree, test_fold_multiplier
             )
-        ).all()
-        <= 1e-3
+        ),
+        1.0,
     )
 
 
@@ -197,13 +203,17 @@ def test_invalid_degree_fold_multiplier_sample_matrix(
 
 
 def test_lre_inference_with_chunking():
+    """Verify the dimension of a chunked sample matrix for some input circuit
+    is smaller than the non-chunked sample matrix for the same input circuit.
+    """
     circ = test_circuit1 * 7
-    chunked_sample_matrix_dim = sample_matrix(circ, 2, 2, 4).shape
+    chunked_sample_matrix_dim = sample_matrix(circ, 2, 2, num_chunks=4).shape
     non_chunked_sample_matrix_dim = sample_matrix(circ, 2, 2).shape
     assert chunked_sample_matrix_dim[0] < non_chunked_sample_matrix_dim[0]
 
 
 def test_sample_matrix_numerical_stability():
+    """Verify sample matrix function works for very large circuits."""
     large_circuit = Circuit([ops.H.on(LineQubit(i)) for i in range(10000)])
     matrix = sample_matrix(large_circuit, 5, 10000)
     assert np.isfinite(matrix).all()
@@ -212,10 +222,13 @@ def test_sample_matrix_numerical_stability():
 
 @pytest.mark.parametrize("num_chunks", [2, 3])
 def test_eval(num_chunks):
+    """Verify the number of calculated linear combination coefficients matches
+    to the number of scaled chunked circuits."""
     coeffs = linear_combination_coefficients(
         7 * test_circuit2, 2, 2, num_chunks
     )
     multiple_scaled_circuits = multivariate_layer_scaling(
         7 * test_circuit2, 2, 2, num_chunks
     )
     assert len(coeffs) == len(multiple_scaled_circuits)
+    assert np.isclose(sum(coeffs), 1.0)