Controlled-Initialize instruction

Question

I'm trying to build a new instruction for my circuit. This instruction needs both a controller qubit qctl and an arbitrary register qreg. When qctl is set then the Qiskit's initialize function is applied to qreg.

The original Initialize gate (version 10.5) can be found in the official documentation or locally at path: /anaconda3/envs/<environment name>/lib/python3.7/site-packages/qiskit/extensions/initializer.py. It follows a particular procedure which consists of applying a sequence of RY and RZ gates in order to match the desired state.

The idea is:

• copy the Initialize instruction, naming it ControlledInitialize;
• pass an additional single qubit register qctl to ControlledInitialize;
• change all RZ, RY gates with CRZ, CRY gates (the first one is already available, the second one have to be made from scratch).

The problem

It seems that I have passed qctl register in the wrong way, in fact the below minimum example throws the error:

DAGCircuitError: '(qu)bit qctl[0] not found'

Minimum example

import numpy as np
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import BasicAer, execute

# copy here CRY and ControlledInitialize implementation

desired_vector = [ 1 / math.sqrt(2), 0, 0, 1 / math.sqrt(2) ]
qctl = QuantumRegister(1, "qctl")
qreg = QuantumRegister(2, "qreg")
creg = ClassicalRegister(2, "creg")
circuit = QuantumCircuit(qctl, qreg, creg)

circuit.x(qctl)
circuit.controlled_initialize(qctl, desired_vector, qreg)
circuit.measure(qreg, creg)

job = execute(circuit, BasicAer.get_backend('qasm_simulator'), shots=10000)
print('Counts: ', job.result().get_counts(circuit))

The implementation

The whole code can be seen here as .py and here as Jupyter notebook.

CRZ, CRY gates

CRZ is a standard gate, you can use it with

from qiskit.extensions.standard.crz import CrzGate

CRY is present in Aqua module as a function but not as a subclass of Gate class. You can easily derive the gate implementation:

from qiskit.circuit import CompositeGate
from qiskit.circuit import Gate
from qiskit.circuit import QuantumCircuit
from qiskit.circuit import QuantumRegister
from qiskit.circuit.decorators import _op_expand, _to_bits
from qiskit.extensions.standard.u3 import U3Gate
from qiskit.extensions.standard.cx import CnotGate

class CryGate(Gate):
    """controlled-rz gate."""

    def __init__(self, theta):
        """Create new cry gate."""
        super().__init__("cry", 2, [theta]) # 2 = number of qubits

    def _define(self):
        """
        self.u3(theta / 2, 0, 0, q_target)
        self.cx(q_control, q_target)
        self.u3(-theta / 2, 0, 0, q_target)
        self.cx(q_control, q_target)
        """
        definition = []
        q = QuantumRegister(2, "q")
        rule = [
            (U3Gate(self.params[0] / 2, 0, 0), [q[1]], []),
            (CnotGate(), [q[0], q[1]], []),
            (U3Gate(-self.params[0] / 2, 0, 0), [q[1]], []),
            (CnotGate(), [q[0], q[1]], [])
        ]
        for inst in rule:
            definition.append(inst)
        self.definition = definition

    def inverse(self):
        """Invert this gate."""
        return CrzGate(-self.params[0])


@_to_bits(2)
@_op_expand(2)
def cry(self, theta, ctl, tgt):
    """Apply crz from ctl to tgt with angle theta."""
    return self.append(CryGate(theta), [ctl, tgt], [])


QuantumCircuit.cry = cry
CompositeGate.cry = cry

ControlledInitialize instruction

Any modification of original Initialize instruction is denoted with WATCH ME comment. Here an overview:

• in __init__ I just save the single qubit control register;
• in _define, gates_to_uncompute, multiplexer the temporary circuit built will have also qctl register;
• in _define, gates_to_uncompute, multiplexer any append function call is enriched with qctl register in the list of qubits taken as second parameter;

• in gates_to_uncompute just substitute RYGate/RZGate with CryGate/CrzGate.

  class ControlledInitialize(Instruction):
  
  """Complex amplitude initialization.
  
  Class that implements the (complex amplitude) initialization of some
  flexible collection of qubit registers (assuming the qubits are in the
  zero state).
  """
  
  def __init__(self, controlled_qubit, params):
      """Create new initialize composite.
  
      params (list): vector of complex amplitudes to initialize to
      """
  
      # WATCH ME: save controlled qubit register
      self.controlled_qubit = controlled_qubit
  
      num_qubits = math.log2(len(params))
  
      # Check if param is a power of 2
      if num_qubits == 0 or not num_qubits.is_integer():
          raise QiskitError("Desired statevector length not a positive power of 2.")
  
      # Check if probabilities (amplitudes squared) sum to 1
      if not math.isclose(sum(np.absolute(params) ** 2), 1.0,
                          abs_tol=_EPS):
          raise QiskitError("Sum of amplitudes-squared does not equal one.")
  
      num_qubits = int(num_qubits)
  
      super().__init__("controlledinitialize", num_qubits, 0, params) # +1 per il controllo
  
  def _define(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
  
      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()
  
      # invert the circuit to create the desired vector from zero (assuming
      # the qubits are in the zero state)
      initialize_instr = disentangling_circuit.to_instruction().inverse()
  
      q = QuantumRegister(self.num_qubits, 'q')
      initialize_circuit = QuantumCircuit(self.controlled_qubit, q, name='init_def')
      for qubit in q:
          initialize_circuit.append(Reset(), [qubit])
  
      # WATCH ME: cambiati registri
      temp_qubitsreg = [ self.controlled_qubit[0] ] + q[:]
      # initialize_circuit.append(initialize_instr, q[:])
      initialize_circuit.append(initialize_instr, temp_qubitsreg)
  
      self.definition = initialize_circuit.data
  
  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 |00..0>
      """
      q = QuantumRegister(self.num_qubits)
      # WATCH ME: aggiunto registro controlled_qubit
      circuit = QuantumCircuit(self.controlled_qubit, 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) = ControlledInitialize._rotations_to_disentangle(remaining_param)
          # WATCH ME: Initialize._rotations_to_disentangle diventa ControlledInitialize._rotations_to_disentangle
  
          # perform the required rotations to decouple the LSB qubit (so that
          # it can be "factored" out, leaving a shorter amplitude vector to peel away)
  
          # WATCH ME: substitute RZ with CRZ
          # rz_mult = self._multiplex(RZGate, phis)
          rz_mult = self._multiplex(CrzGate, phis)
  
          # WATCH ME: substitute RY with CRY
          # ry_mult = self._multiplex(RYGate, thetas)
          ry_mult = self._multiplex(CryGate, thetas)
  
          # WATCH ME: cambiati registri
          temp_qubitsreg = [ self.controlled_qubit[0] ] + q[i:self.num_qubits]
          # circuit.append(rz_mult.to_instruction(), q[i:self.num_qubits])
          # circuit.append(ry_mult.to_instruction(), q[i:self.num_qubits])
          circuit.append(rz_mult.to_instruction(), temp_qubitsreg)
          circuit.append(ry_mult.to_instruction(), temp_qubitsreg)
  
          print("Z: ", phis, " | Y: ", thetas)
  
      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) = ControlledInitialize._bloch_angles(local_param[2 * i: 2 * (i + 1)])
          # WATCH ME: Initialize._bloch_angles diventa ControlledInitialize._bloch_angles
  
          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 = np.absolute(a_complex)
      final_r = float(np.sqrt(mag_a ** 2 + np.absolute(b_complex) ** 2))
      if final_r < _EPS:
          theta = 0
          phi = 0
          final_r = 0
          final_t = 0
      else:
          theta = float(2 * np.arccos(mag_a / final_r))
          a_arg = np.angle(a_complex)
          b_arg = np.angle(b_complex)
          final_t = a_arg + b_arg
          phi = b_arg - a_arg
  
      return final_r * np.exp(1.J * final_t / 2), theta, phi
  
  def _multiplex(self, target_gate, list_of_angles):
      """
      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
  
      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)
      # WATCH ME: aggiunto registro controlled_qubit
      circuit = QuantumCircuit(self.controlled_qubit, q, name="multiplex" + local_num_qubits.__str__())
  
      lsb = q[0]
      msb = q[local_num_qubits - 1]
  
      # case of no multiplexing: base case for recursion
      if local_num_qubits == 1:
          temp_qubitsreg = [ self.controlled_qubit[0], q[0] ]
          circuit.append(target_gate(list_of_angles[0]), temp_qubitsreg)
          return circuit
  
  
      # calc angle weights, assuming recursion (that is the lower-level
      # requested angles have been correctly implemented by recursion
      angle_weight = scipy.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)])
      temp_qubitsreg =  [ self.controlled_qubit[0] ] + q[0:-1]
      circuit.append(multiplex_1.to_instruction(), temp_qubitsreg)
  
      # attach CNOT as follows, thereby flipping the LSB qubit
      circuit.append(CnotGate(), [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):])
      temp_qubitsreg = [ self.controlled_qubit[0] ] + q[0:-1]
      if list_len > 1:
          circuit.append(multiplex_2.to_instruction().mirror(), temp_qubitsreg)
      else:
          circuit.append(multiplex_2.to_instruction(), temp_qubitsreg)
  
      # attach a final CNOT
      circuit.append(CnotGate(), [msb, lsb])
  
      return circuit
  

Qiskit version

Latest version is used:

import qiskit
qiskit.__qiskit_version__

{'qiskit': '0.10.5',
 'qiskit-terra': '0.8.2',
 'qiskit-ignis': '0.1.1',
 'qiskit-aer': '0.2.1',
 'qiskit-ibmq-provider': '0.2.2',
 'qiskit-aqua': '0.5.2'}

Answer 1

This is occurring because you declare your definition rule on two registers, but the way nodes are added to the DAG, only one register will be added. It is defined over both the QuantumRegister "q" you define in the method, and also the register passed in to self.params.

To fix this therefore you need to update your definition to work on only one register. For example, the first qubit could be the control and the other qubits can be initialised.