How do we code the matrix for a controlled operation knowing the control qubit, the target qubit and the $2\times 2$ unitary?

Question

Having n qubits, I want to have the unitary described a controlled operation. Say for example you get as input a unitary, an index for a controlled qubit and another for a target.

How would you code this unitary operation?

Answer 1

I would solve it like that: suppose you have CNOT gate; it's unitary matrix can be written as a sum of tensor products $$ \begin{pmatrix}1 & 0 \\ 0 & 0 \end{pmatrix}\otimes \begin{pmatrix}1 & 0 \\ 0 & 1 \end{pmatrix}+ \begin{pmatrix}0 & 0 \\ 0 & 1 \end{pmatrix}\otimes \begin{pmatrix}0 & 1 \\ 1 & 0 \end{pmatrix} $$ (this is another way to say: if first qubit is 0, the second qubit is unchanged; if first qubit is 1, the second qubit is swapped).

Now if we have more qubits, we need to insert additional identity matrices $$ \begin{pmatrix}1 & 0 \\ 0 & 1 \end{pmatrix} $$ into the tensor products from the left, from the right, or into the middle, depending on the indexes of control and target CNOT qubits.

PS: I also assumed that the index of control qubit is less than the index of target qubit; if this is not the case swap the terms in the tensor products in the first equation.

Answer 2

Here’s some pseudo code, where id(n) creates an $2^n\times 2^n$ identity matrix, and tensor(A,B,...) returns $A\otimes B\otimes\ldots$.

def cU(ctrl,targ,U,size):
    '''implement controlled-U with:
          control qubit ctrl,
          target qubit targ,
          within a set of size qubits'''
#check input ranges
    assert 1<=ctrl<=size
    assert 1<=targ<=size
    assert ctrl<>targ
    assert ctrl,targ,size ∊ ℤ
#ensure U is a 2x2 unitary
    assert U∊ℂ2x2
    assert U.U†=id(2)
#the actual code
    if ctrl<targ:
        return id(size)+tensor(id(ctrl-1),id(1)-Z,id(targ-1-ctrl),U-id(1),id(size-targ))/2
    else:
        return id(size)+tensor(id(targ-1),U-id(1),id(ctrl-1-targ),id(1)-Z,id(size-ctrl))/2

However, remember that usually you're trying to calculate the action of a unitary on some state vector. It will be far more memory efficient to calculate that application directly, rather than first calculating the unitary matrix and applying it to the state vector.

To understand where this formula came from, think about the two-qubit version, where the first qubit is the control qubit. You'd normally write the unitary as $$ |0\rangle\langle 0|\otimes\mathbb{I}+|1\rangle\langle 1|\otimes U. $$ Let's rewrite this as $$ (\mathbb{I}-|1\rangle\langle 1|)\otimes\mathbb{I}+|1\rangle\langle 1|\otimes U=\mathbb{I}\otimes\mathbb{I}+|1\rangle\langle 1|\otimes (U-\mathbb{I}). $$ It can be easier to write things in terms of Pauli matrices, so $$ |1\rangle\langle 1|=(\mathbb{I}-Z)/2. $$ To get the same unitary on a different number of qubits, you just need to pad with identity matrices everywhere.