Allow users to optionally specify Isometry tolerance (#6482)

* Allow users to optionally specify Isometry tolerance

Instead of hard-coding a fixed value for floating-point epsilon
tolerances, allow users to optionally specify what value of the
tolerance they would like to use instead.

Fixes #3789

* Added unit tests for custom isometry tolerance

* Added release notes for custom isometry tolerance

* Fixed lint errors

* Update test/python/circuit/test_isometry.py

Co-authored-by: Abby Mitchell <abby.mitchell@btinternet.com>
Co-authored-by: Luciano Bello <bel@zurich.ibm.com>

qiskit/extensions/quantum_initializer/isometry.py

@@ -54,6 +54,8 @@ class Isometry(Instruction):
             not accounted for here).
 
         num_ancillas_dirty (int): number of additional ancillas that start in an arbitrary state
+
+        epsilon (float) (optional): error tolerance of calculations
     """
 
     # Notation: In the following decomposition we label the qubit by
@@ -63,7 +65,7 @@ class Isometry(Instruction):
     # finally, we convert the labels back to the qubit numbering used in Qiskit
     # (using: _get_qubits_by_label)
 
-    def __init__(self, isometry, num_ancillas_zero, num_ancillas_dirty):
+    def __init__(self, isometry, num_ancillas_zero, num_ancillas_dirty, epsilon=_EPS):
         # Convert to numpy array in case not already an array
         isometry = np.array(isometry, dtype=complex)
 
@@ -74,6 +76,7 @@ def __init__(self, isometry, num_ancillas_zero, num_ancillas_dirty):
         self.num_ancillas_zero = num_ancillas_zero
         self.num_ancillas_dirty = num_ancillas_dirty
         self._inverse = None
+        self._epsilon = epsilon
 
         # Check if the isometry has the right dimension and if the columns are orthonormal
         n = np.log2(isometry.shape[0])
@@ -90,7 +93,7 @@ def __init__(self, isometry, num_ancillas_zero, num_ancillas_dirty):
             raise QiskitError(
                 "The input matrix has more columns than rows and hence " "it can't be an isometry."
             )
-        if not is_isometry(isometry, _EPS):
+        if not is_isometry(isometry, self._epsilon):
             raise QiskitError(
                 "The input matrix has non orthonormal columns and hence " "it is not an isometry."
             )
@@ -137,7 +140,7 @@ def _gates_to_uncompute(self):
             diag.append(remaining_isometry[column_index, 0])
             # remove first column (which is now stored in diag)
             remaining_isometry = remaining_isometry[:, 1:]
-        if len(diag) > 1 and not _diag_is_identity_up_to_global_phase(diag):
+        if len(diag) > 1 and not _diag_is_identity_up_to_global_phase(diag, self._epsilon):
             circuit.diagonal(np.conj(diag).tolist(), q_input)
         return circuit
 
@@ -167,9 +170,9 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
         index2 = (2 * _a(k, s + 1) + 1) * 2 ** s + _b(k, s + 1)
         target_label = n - s - 1
         # Check if a MCG is required
-        if _k_s(k, s) == 0 and _b(k, s + 1) != 0 and np.abs(v[index2, k_prime]) > _EPS:
+        if _k_s(k, s) == 0 and _b(k, s + 1) != 0 and np.abs(v[index2, k_prime]) > self._epsilon:
             # Find the MCG, decompose it and apply it to the remaining isometry
-            gate = _reverse_qubit_state([v[index1, k_prime], v[index2, k_prime]], 0)
+            gate = _reverse_qubit_state([v[index1, k_prime], v[index2, k_prime]], 0, self._epsilon)
             control_labels = [
                 i
                 for i, x in enumerate(_get_binary_rep_as_list(k, n))
@@ -189,7 +192,7 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
         # UCGate to disentangle a qubit:
         # Find the UCGate, decompose it and apply it to the remaining isometry
         single_qubit_gates = self._find_squs_for_disentangling(v, k, s)
-        if not _ucg_is_identity_up_to_global_phase(single_qubit_gates):
+        if not _ucg_is_identity_up_to_global_phase(single_qubit_gates, self._epsilon):
             control_labels = list(range(target_label))
             diagonal_ucg = self._append_ucg_up_to_diagonal(
                 circuit, q, single_qubit_gates, control_labels, target_label
@@ -225,6 +228,7 @@ def _find_squs_for_disentangling(self, v, k, s):
                     v[(2 * i + 1) * 2 ** s + _b(k, s), k_prime],
                 ],
                 _k_s(k, s),
+                self._epsilon,
             )
             for i in range(i_start, 2 ** (n - s - 1))
         ]
@@ -311,10 +315,10 @@ def inverse(self):
 
 # Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or
 # basis_state=1 respectively
-def _reverse_qubit_state(state, basis_state):
+def _reverse_qubit_state(state, basis_state, epsilon):
     state = np.array(state)
     r = np.linalg.norm(state)
-    if r < _EPS:
+    if r < epsilon:
         return np.eye(2, 2)
     if basis_state == 0:
         m = np.array([[np.conj(state[0]), np.conj(state[1])], [-state[1], state[0]]]) / r
@@ -523,8 +527,8 @@ def _k_s(k, s):
 # Check if a gate of a special form is equal to the identity gate up to global phase
 
 
-def _ucg_is_identity_up_to_global_phase(single_qubit_gates):
-    if not np.abs(single_qubit_gates[0][0, 0]) < _EPS:
+def _ucg_is_identity_up_to_global_phase(single_qubit_gates, epsilon):
+    if not np.abs(single_qubit_gates[0][0, 0]) < epsilon:
         global_phase = 1.0 / (single_qubit_gates[0][0, 0])
     else:
         return False
@@ -534,19 +538,25 @@ def _ucg_is_identity_up_to_global_phase(single_qubit_gates):
     return True
 
 
-def _diag_is_identity_up_to_global_phase(diag):
-    if not np.abs(diag[0]) < _EPS:
+def _diag_is_identity_up_to_global_phase(diag, epsilon):
+    if not np.abs(diag[0]) < epsilon:
         global_phase = 1.0 / (diag[0])
     else:
         return False
     for d in diag:
-        if not np.abs(global_phase * d - 1) < _EPS:
+        if not np.abs(global_phase * d - 1) < epsilon:
             return False
     return True
 
 
 def iso(
-    self, isometry, q_input, q_ancillas_for_output, q_ancillas_zero=None, q_ancillas_dirty=None
+    self,
+    isometry,
+    q_input,
+    q_ancillas_for_output,
+    q_ancillas_zero=None,
+    q_ancillas_dirty=None,
+    epsilon=_EPS,
 ):
     """
     Attach an arbitrary isometry from m to n qubits to a circuit. In particular,
@@ -568,6 +578,8 @@ def iso(
             which are assumed to start in the zero state. Default is q_ancillas_zero = None.
         q_ancillas_dirty (QuantumRegister|list[Qubit]): list of ancilla qubits
             which can start in an arbitrary state. Default is q_ancillas_dirty = None.
+        epsilon (float): error tolerance of calculations.
+            Default is epsilon = _EPS.
 
     Returns:
         QuantumCircuit: the isometry is attached to the quantum circuit.
@@ -595,7 +607,7 @@ def iso(
         q_ancillas_dirty = q_ancillas_dirty[:]
 
     return self.append(
-        Isometry(isometry, len(q_ancillas_zero), len(q_ancillas_dirty)),
+        Isometry(isometry, len(q_ancillas_zero), len(q_ancillas_dirty), epsilon=epsilon),
         q_input + q_ancillas_for_output + q_ancillas_zero + q_ancillas_dirty,
     )
 

releasenotes/notes/custom-isometry-tolerance-ee6dcb4cafce9ad4.yaml

@@ -0,0 +1,21 @@
+---
+features:
+  - |
+    Allowed users to set a custom isometry tolerance. Instead of hard-coding a 
+    fixed value for floating-point epsilon tolerances, allow users to 
+    optionally specify what value of the tolerance they would like to use 
+    instead. For example::
+
+        import numpy as np
+        from qiskit import QuantumRegister, QuantumCircuit
+
+        tolerance = 1e-8
+        iso = np.eye(2,2)
+        num_q_output = int(np.log2(iso.shape[0]))
+        num_q_input = int(np.log2(iso.shape[1]))
+        q = QuantumRegister(num_q_output)
+        qc = QuantumCircuit(q)
+
+        qc.isometry(iso, q[:num_q_input], q[num_q_input:], epsilon=tolerance)
+
+

test/python/circuit/test_isometry.py

@@ -71,6 +71,45 @@ def test_isometry(self, iso):
         iso_desired = iso
         self.assertTrue(matrix_equal(iso_from_circuit, iso_desired, ignore_phase=True))
 
+    @data(
+        np.eye(2, 2),
+        random_unitary(2).data,
+        np.eye(4, 4),
+        random_unitary(4).data[:, 0],
+        np.eye(4, 4)[:, 0:2],
+        random_unitary(4).data,
+        np.eye(4, 4)[:, np.random.permutation(np.eye(4, 4).shape[1])][:, 0:2],
+        np.eye(8, 8)[:, np.random.permutation(np.eye(8, 8).shape[1])],
+        random_unitary(8).data[:, 0:4],
+        random_unitary(8).data,
+        random_unitary(16).data,
+        random_unitary(16).data[:, 0:8],
+    )
+    def test_isometry_tolerance(self, iso):
+        """Tests for the decomposition of isometries from m to n qubits with a custom tolerance"""
+        if len(iso.shape) == 1:
+            iso = iso.reshape((len(iso), 1))
+        num_q_output = int(np.log2(iso.shape[0]))
+        num_q_input = int(np.log2(iso.shape[1]))
+        q = QuantumRegister(num_q_output)
+        qc = QuantumCircuit(q)
+
+        # Compute isometry with custom tolerance
+        qc.isometry(iso, q[:num_q_input], q[num_q_input:], epsilon=1e-3)
+
+        # Verify the circuit can be decomposed
+        self.assertIsInstance(qc.decompose(), QuantumCircuit)
+
+        # Decompose the gate
+        qc = transpile(qc, basis_gates=["u1", "u3", "u2", "cx", "id"])
+
+        # Simulate the decomposed gate
+        simulator = BasicAer.get_backend("unitary_simulator")
+        result = execute(qc, simulator).result()
+        unitary = result.get_unitary(qc)
+        iso_from_circuit = unitary[::, 0 : 2 ** num_q_input]
+        self.assertTrue(matrix_equal(iso_from_circuit, iso, ignore_phase=True))
+
     @data(
         np.eye(2, 2),
         random_unitary(2).data,