Replace initialization method by Isometry in StatePreparation (#12178)

* replace initializetion method by Isometry in StatePreparation

* add names to QuantumRegister and QuantumCircuit

* update docs to the new reference

* remove old initialization code based on Shende et al

* fix reference in docs

* add release notes

* fix lint errors

* fix references in docs

* add a benchmark for state preparation

* update circuit name following review

qiskit/circuit/library/data_preparation/initializer.py

@@ -36,6 +36,14 @@ class Initialize(Instruction):
     the :class:`~.library.StatePreparation` class.
     Note that ``Initialize`` is an :class:`~.circuit.Instruction` and not a :class:`.Gate` since it
     contains a reset instruction, which is not unitary.
+
+    The initial state is prepared based on the :class:`~.library.Isometry` synthesis described in [1].
+
+    References:
+        1. Iten et al., Quantum circuits for isometries (2016).
+           `Phys. Rev. A 93, 032318
+           <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.93.032318>`__.
+
     """
 
     def __init__(

qiskit/circuit/library/data_preparation/state_preparation.py

@@ -11,7 +11,6 @@
 # that they have been altered from the originals.
 """Prepare a quantum state from the state where all qubits are 0."""
 
-import cmath
 from typing import Union, Optional
 
 import math
@@ -21,11 +20,10 @@
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.quantumregister import QuantumRegister
 from qiskit.circuit.gate import Gate
-from qiskit.circuit.library.standard_gates.x import CXGate, XGate
+from qiskit.circuit.library.standard_gates.x import XGate
 from qiskit.circuit.library.standard_gates.h import HGate
 from qiskit.circuit.library.standard_gates.s import SGate, SdgGate
-from qiskit.circuit.library.standard_gates.ry import RYGate
-from qiskit.circuit.library.standard_gates.rz import RZGate
+from qiskit.circuit.library.generalized_gates import Isometry
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.quantum_info.states.statevector import Statevector  # pylint: disable=cyclic-import
 
@@ -71,13 +69,13 @@ def __init__(
         Raises:
             QiskitError: ``num_qubits`` parameter used when ``params`` is not an integer
 
-        When a Statevector argument is passed the state is prepared using a recursive
-        initialization algorithm, including optimizations, from [1], as well
-        as some additional optimizations including removing zero rotations and double cnots.
+        When a Statevector argument is passed the state is prepared based on the
+        :class:`~.library.Isometry` synthesis described in [1].
 
-        **References:**
-        [1] Shende, Bullock, Markov. Synthesis of Quantum Logic Circuits (2004)
-        [`https://arxiv.org/abs/quant-ph/0406176v5`]
+        References:
+            1. Iten et al., Quantum circuits for isometries (2016).
+               `Phys. Rev. A 93, 032318
+               <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.93.032318>`__.
 
         """
         self._params_arg = params
@@ -119,7 +117,7 @@ def _define(self):
         elif self._from_int:
             self.definition = self._define_from_int()
         else:
-            self.definition = self._define_synthesis()
+            self.definition = self._define_synthesis_isom()
 
     def _define_from_label(self):
         q = QuantumRegister(self.num_qubits, "q")
@@ -168,29 +166,18 @@ def _define_from_int(self):
         # we don't need to invert anything
         return initialize_circuit
 
-    def _define_synthesis(self):
-        """Calculate a subcircuit that implements this initialization
-
-        Implements a recursive initialization algorithm, including optimizations,
-        from "Synthesis of Quantum Logic Circuits" Shende, Bullock, Markov
-        https://arxiv.org/abs/quant-ph/0406176v5
+    def _define_synthesis_isom(self):
+        """Calculate a subcircuit that implements this initialization via isometry"""
+        q = QuantumRegister(self.num_qubits, "q")
+        initialize_circuit = QuantumCircuit(q, name="init_def")
 
-        Additionally implements some extra optimizations: remove zero rotations and
-        double cnots.
-        """
-        # call to generate the circuit that takes the desired vector to zero
-        disentangling_circuit = self._gates_to_uncompute()
+        isom = Isometry(self._params_arg, 0, 0)
+        initialize_circuit.append(isom, q[:])
 
         # invert the circuit to create the desired vector from zero (assuming
         # the qubits are in the zero state)
-        if self._inverse is False:
-            initialize_instr = disentangling_circuit.to_instruction().inverse()
-        else:
-            initialize_instr = disentangling_circuit.to_instruction()
-
-        q = QuantumRegister(self.num_qubits, "q")
-        initialize_circuit = QuantumCircuit(q, name="init_def")
-        initialize_circuit.append(initialize_instr, q[:])
+        if self._inverse is True:
+            return initialize_circuit.inverse()
 
         return initialize_circuit
 
@@ -253,164 +240,3 @@ def validate_parameter(self, parameter):
 
     def _return_repeat(self, exponent: float) -> "Gate":
         return Gate(name=f"{self.name}*{exponent}", num_qubits=self.num_qubits, params=[])
-
-    def _gates_to_uncompute(self):
-        """Call to create a circuit with gates that take the desired vector to zero.
-
-        Returns:
-            QuantumCircuit: circuit to take self.params vector to :math:`|{00\\ldots0}\\rangle`
-        """
-        q = QuantumRegister(self.num_qubits)
-        circuit = QuantumCircuit(q, name="disentangler")
-
-        # kick start the peeling loop, and disentangle one-by-one from LSB to MSB
-        remaining_param = self.params
-
-        for i in range(self.num_qubits):
-            # work out which rotations must be done to disentangle the LSB
-            # qubit (we peel away one qubit at a time)
-            (remaining_param, thetas, phis) = StatePreparation._rotations_to_disentangle(
-                remaining_param
-            )
-
-            # perform the required rotations to decouple the LSB qubit (so that
-            # it can be "factored" out, leaving a shorter amplitude vector to peel away)
-
-            add_last_cnot = True
-            if np.linalg.norm(phis) != 0 and np.linalg.norm(thetas) != 0:
-                add_last_cnot = False
-
-            if np.linalg.norm(phis) != 0:
-                rz_mult = self._multiplex(RZGate, phis, last_cnot=add_last_cnot)
-                circuit.append(rz_mult.to_instruction(), q[i : self.num_qubits])
-
-            if np.linalg.norm(thetas) != 0:
-                ry_mult = self._multiplex(RYGate, thetas, last_cnot=add_last_cnot)
-                circuit.append(ry_mult.to_instruction().reverse_ops(), q[i : self.num_qubits])
-        circuit.global_phase -= np.angle(sum(remaining_param))
-        return circuit
-
-    @staticmethod
-    def _rotations_to_disentangle(local_param):
-        """
-        Static internal method to work out Ry and Rz rotation angles used
-        to disentangle the LSB qubit.
-        These rotations make up the block diagonal matrix U (i.e. multiplexor)
-        that disentangles the LSB.
-
-        [[Ry(theta_1).Rz(phi_1)  0   .   .   0],
-        [0         Ry(theta_2).Rz(phi_2) .  0],
-                                    .
-                                        .
-        0         0           Ry(theta_2^n).Rz(phi_2^n)]]
-        """
-        remaining_vector = []
-        thetas = []
-        phis = []
-
-        param_len = len(local_param)
-
-        for i in range(param_len // 2):
-            # Ry and Rz rotations to move bloch vector from 0 to "imaginary"
-            # qubit
-            # (imagine a qubit state signified by the amplitudes at index 2*i
-            # and 2*(i+1), corresponding to the select qubits of the
-            # multiplexor being in state |i>)
-            (remains, add_theta, add_phi) = StatePreparation._bloch_angles(
-                local_param[2 * i : 2 * (i + 1)]
-            )
-
-            remaining_vector.append(remains)
-
-            # rotations for all imaginary qubits of the full vector
-            # to move from where it is to zero, hence the negative sign
-            thetas.append(-add_theta)
-            phis.append(-add_phi)
-
-        return remaining_vector, thetas, phis
-
-    @staticmethod
-    def _bloch_angles(pair_of_complex):
-        """
-        Static internal method to work out rotation to create the passed-in
-        qubit from the zero vector.
-        """
-        [a_complex, b_complex] = pair_of_complex
-        # Force a and b to be complex, as otherwise numpy.angle might fail.
-        a_complex = complex(a_complex)
-        b_complex = complex(b_complex)
-        mag_a = abs(a_complex)
-        final_r = math.sqrt(mag_a**2 + abs(b_complex) ** 2)
-        if final_r < _EPS:
-            theta = 0
-            phi = 0
-            final_r = 0
-            final_t = 0
-        else:
-            theta = 2 * math.acos(mag_a / final_r)
-            a_arg = cmath.phase(a_complex)
-            b_arg = cmath.phase(b_complex)
-            final_t = a_arg + b_arg
-            phi = b_arg - a_arg
-
-        return final_r * cmath.exp(1.0j * final_t / 2), theta, phi
-
-    def _multiplex(self, target_gate, list_of_angles, last_cnot=True):
-        """
-        Return a recursive implementation of a multiplexor circuit,
-        where each instruction itself has a decomposition based on
-        smaller multiplexors.
-
-        The LSB is the multiplexor "data" and the other bits are multiplexor "select".
-
-        Args:
-            target_gate (Gate): Ry or Rz gate to apply to target qubit, multiplexed
-                over all other "select" qubits
-            list_of_angles (list[float]): list of rotation angles to apply Ry and Rz
-            last_cnot (bool): add the last cnot if last_cnot = True
-
-        Returns:
-            DAGCircuit: the circuit implementing the multiplexor's action
-        """
-        list_len = len(list_of_angles)
-        local_num_qubits = int(math.log2(list_len)) + 1
-
-        q = QuantumRegister(local_num_qubits)
-        circuit = QuantumCircuit(q, name="multiplex" + str(local_num_qubits))
-
-        lsb = q[0]
-        msb = q[local_num_qubits - 1]
-
-        # case of no multiplexing: base case for recursion
-        if local_num_qubits == 1:
-            circuit.append(target_gate(list_of_angles[0]), [q[0]])
-            return circuit
-
-        # calc angle weights, assuming recursion (that is the lower-level
-        # requested angles have been correctly implemented by recursion
-        angle_weight = np.kron([[0.5, 0.5], [0.5, -0.5]], np.identity(2 ** (local_num_qubits - 2)))
-
-        # calc the combo angles
-        list_of_angles = angle_weight.dot(np.array(list_of_angles)).tolist()
-
-        # recursive step on half the angles fulfilling the above assumption
-        multiplex_1 = self._multiplex(target_gate, list_of_angles[0 : (list_len // 2)], False)
-        circuit.append(multiplex_1.to_instruction(), q[0:-1])
-
-        # attach CNOT as follows, thereby flipping the LSB qubit
-        circuit.append(CXGate(), [msb, lsb])
-
-        # implement extra efficiency from the paper of cancelling adjacent
-        # CNOTs (by leaving out last CNOT and reversing (NOT inverting) the
-        # second lower-level multiplex)
-        multiplex_2 = self._multiplex(target_gate, list_of_angles[(list_len // 2) :], False)
-        if list_len > 1:
-            circuit.append(multiplex_2.to_instruction().reverse_ops(), q[0:-1])
-        else:
-            circuit.append(multiplex_2.to_instruction(), q[0:-1])
-
-        # attach a final CNOT
-        if last_cnot:
-            circuit.append(CXGate(), [msb, lsb])
-
-        return circuit

qiskit/circuit/library/generalized_gates/isometry.py

@@ -45,10 +45,10 @@ class Isometry(Instruction):
 
     The decomposition is based on [1].
 
-    **References:**
-
-    [1] Iten et al., Quantum circuits for isometries (2016).
-        `Phys. Rev. A 93, 032318 <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.93.032318>`__.
+    References:
+        1. Iten et al., Quantum circuits for isometries (2016).
+           `Phys. Rev. A 93, 032318
+           <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.93.032318>`__.
 
     """
 
@@ -123,8 +123,8 @@ def _define(self):
         #  later here instead.
         gate = self.inv_gate()
         gate = gate.inverse()
-        q = QuantumRegister(self.num_qubits)
-        iso_circuit = QuantumCircuit(q)
+        q = QuantumRegister(self.num_qubits, "q")
+        iso_circuit = QuantumCircuit(q, name="isometry")
         iso_circuit.append(gate, q[:])
         self.definition = iso_circuit
 
@@ -139,8 +139,8 @@ def _gates_to_uncompute(self):
         Call to create a circuit with gates that take the desired isometry to the first 2^m columns
          of the 2^n*2^n identity matrix (see https://arxiv.org/abs/1501.06911)
         """
-        q = QuantumRegister(self.num_qubits)
-        circuit = QuantumCircuit(q)
+        q = QuantumRegister(self.num_qubits, "q")
+        circuit = QuantumCircuit(q, name="isometry_to_uncompute")
         (
             q_input,
             q_ancillas_for_output,

qiskit/circuit/library/generalized_gates/mcg_up_to_diagonal.py

@@ -68,8 +68,8 @@ def __init__(
     def _define(self):
         mcg_up_diag_circuit, _ = self._dec_mcg_up_diag()
         gate = mcg_up_diag_circuit.to_instruction()
-        q = QuantumRegister(self.num_qubits)
-        mcg_up_diag_circuit = QuantumCircuit(q)
+        q = QuantumRegister(self.num_qubits, "q")
+        mcg_up_diag_circuit = QuantumCircuit(q, name="mcg_up_to_diagonal")
         mcg_up_diag_circuit.append(gate, q[:])
         self.definition = mcg_up_diag_circuit
 
@@ -108,8 +108,8 @@ def _dec_mcg_up_diag(self):
             q=[q_target,q_controls,q_ancilla_zero,q_ancilla_dirty]
         """
         diag = np.ones(2 ** (self.num_controls + 1)).tolist()
-        q = QuantumRegister(self.num_qubits)
-        circuit = QuantumCircuit(q)
+        q = QuantumRegister(self.num_qubits, "q")
+        circuit = QuantumCircuit(q, name="mcg_up_to_diagonal")
         (q_target, q_controls, q_ancillas_zero, q_ancillas_dirty) = self._define_qubit_role(q)
         # ToDo: Keep this threshold updated such that the lowest gate count is achieved:
         # ToDo: we implement the MCG with a UCGate up to diagonal if the number of controls is

qiskit/circuit/library/generalized_gates/uc.py

@@ -148,10 +148,10 @@ def _dec_ucg(self):
         the diagonal gate is also returned.
         """
         diag = np.ones(2**self.num_qubits).tolist()
-        q = QuantumRegister(self.num_qubits)
+        q = QuantumRegister(self.num_qubits, "q")
         q_controls = q[1:]
         q_target = q[0]
-        circuit = QuantumCircuit(q)
+        circuit = QuantumCircuit(q, name="uc")
         # If there is no control, we use the ZYZ decomposition
         if not q_controls:
             circuit.unitary(self.params[0], [q])

qiskit/circuit/library/generalized_gates/uc_pauli_rot.py

@@ -69,7 +69,7 @@ def __init__(self, angle_list: list[float], rot_axis: str) -> None:
     def _define(self):
         ucr_circuit = self._dec_ucrot()
         gate = ucr_circuit.to_instruction()
-        q = QuantumRegister(self.num_qubits)
+        q = QuantumRegister(self.num_qubits, "q")
         ucr_circuit = QuantumCircuit(q)
         ucr_circuit.append(gate, q[:])
         self.definition = ucr_circuit
@@ -79,7 +79,7 @@ def _dec_ucrot(self):
         Finds a decomposition of a UC rotation gate into elementary gates
         (C-NOTs and single-qubit rotations).
         """
-        q = QuantumRegister(self.num_qubits)
+        q = QuantumRegister(self.num_qubits, "q")
         circuit = QuantumCircuit(q)
         q_target = q[0]
         q_controls = q[1:]

releasenotes/notes/replace-initialization-algorithm-by-isometry-41f9ffa58f72ece5.yaml

@@ -0,0 +1,7 @@
+---
+features_circuits:
+  - |
+    Replacing the internal synthesis algorithm of :class:`~.library.StatePreparation`
+    and :class:`~.library.Initialize` of Shende et al. by the algorithm given in
+    :class:`~.library.Isometry` of Iten et al.
+    The new algorithm reduces the number of CX gates and the circuit depth by a factor of 2.