In [1]:
%matplotlib inline
import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt, mpld3
import operator
plt.rcParams['figure.figsize'] = 12, 8
mpld3.enable_notebook()

Quantum Computer

In [2]:
def kron(*Xs):
    Y = np.eye(1)
    for X in Xs:
        Y = np.kron(Y, X)
    return Y

def prod(iterable):
    return reduce(operator.mul, iterable, 1)

def randpsi(d):
    psi = np.random.normal(size=d) + 1j * np.random.normal(size=d)
    return psi / np.linalg.norm(psi)
In [3]:
class QC(object):
    def __init__(self, num_qubits):
        """starts out in state |0....0>"""
        self.num_qubits = num_qubits
        self.psi = np.zeros(2**num_qubits)
        self.psi[0] = 1
        
    def _tensor(self, indices, Xs):
        """compute tensor product of operators Xs at indices"""
        assert len(indices) == len(set(indices)), 'no repetitions'
        assert all(i < self.num_qubits for i in indices), 'index out of bounds'
        X = np.eye(1)
        for i in range(self.num_qubits):
            if i in indices:
                X = np.kron(X, Xs[indices.index(i)])
            else:
                X = np.kron(X, np.eye(2))
        return X
    
    def rdm(self, *indices):
        """reduced density matrix"""
        indices = list(indices)
        compl = [i for i in range(self.num_qubits) if i not in indices]
        A = self.psi.reshape([2] * self.num_qubits).transpose(indices + compl).reshape([2**len(indices), 2**len(compl)])
        return np.dot(A, A.T.conj())
    
    def H(self, idx):
        """Hadamard gate"""
        H = [[1, 1], [1, -1]] / np.sqrt(2)
        U = self._tensor([idx], [H])
        self.psi = np.dot(U, self.psi)
    
    def NOT(self, idx):
        """NOT gate"""
        NOT = [[0, 1], [1, 0]]
        U = self._tensor([idx], [NOT])
        self.psi = np.dot(U, self.psi)

    def Z(self, idx):
        """Pauli Z-operator"""
        Z = [[1, 0], [0, -1]]
        U = self._tensor([idx], [Z])
        self.psi = np.dot(U, self.psi)

    def C(self, U, idx, control_idx):
        """controlled unitary"""
        proj0 = [[1, 0], [0, 0]]
        proj1 = [[0, 0], [0, 1]]
        U = self._tensor([idx, control_idx], [np.eye(2), proj0]) + self._tensor([idx, control_idx], [U, proj1])
        self.psi = np.dot(U, self.psi)        
    
    def CNOT(self, idx, control_idx):
        """controller NOT"""
        NOT = [[0, 1], [1, 0]]
        self.C(NOT, idx, control_idx)
    
    def measure(self, idx):
        """measure a qubit"""
        proj0 = [[1, 0], [0, 0]]
        psi0 = np.dot(self._tensor([idx], [proj0]), self.psi)
        prob0 = np.linalg.norm(psi0)**2
        if np.random.random() < prob0:
            self.psi = psi0 / np.sqrt(prob0)
            return 0
        else:
            self.psi = (self.psi - psi0) / np.sqrt(1 - prob0)
            return 1

Teleportation

In [4]:
KET_EPR = np.array([1, 0, 0, 1]) / np.sqrt(2)
    
def teleport(qc, ALICE_1, ALICE_2, BOB):
    # teleportation protocol
    assert np.allclose(qc.rdm(ALICE_2, BOB), np.outer(KET_EPR, KET_EPR))
    qc.CNOT(ALICE_2, control_idx=ALICE_1)
    qc.H(ALICE_1)
    bits = [qc.measure(ALICE_1), qc.measure(ALICE_2)]   # ...these are sent to bob...
    if bits[1]:
        qc.NOT(BOB)
    if bits[0]:
        qc.Z(BOB)
    return qc.rdm(BOB)

# simple test
qc = QC(num_qubits=3)

# prepare initial state of Alice' first qubit
qc.H(0)
expected = qc.rdm(0)

# prepare maximally entangled state between Alice' other qubit and Bob's qubit
qc.H(1)
qc.CNOT(2, control_idx=1)
    
got = teleport(qc, 0, 1, 2)
assert np.allclose(expected, got)

Teleport random states:

In [5]:
%%time
qc = QC(num_qubits=3)

for _ in range(100):
    # initialize QC with random state on first qubit and EPR pair between the other two
    qc.psi = np.kron(randpsi(2), KET_EPR)
    expected = qc.rdm(0)
    
    # teleport
    teleport(qc, 0, 1, 2)
    
    # verify
    got = qc.rdm(2)
    assert np.allclose(expected, got)
Wall time: 119 ms

Entanglement Swapping

Note that the "entanglement" (in fact, the entire reduced state) is preserved by the teleportation procedure.

In [6]:
%%time
qc = QC(num_qubits=4)

for _ in range(100):
    # initialize QC with random state on first two qubits and EPR pair on the other two
    qc.psi = np.kron(randpsi(4), KET_EPR)
    expected = qc.rdm(0, 1)
    
    # teleport
    teleport(qc, 1, 2, 3)
    
    # verify
    got = qc.rdm(0, 3)
    assert np.allclose(expected, got)
Wall time: 153 ms

Thus given an EPR pair between Alice and Bob and another one between Bob and Charly we can get an EPR pair between Alice and Charly!

In [7]:
qc = QC(num_qubits=4)
ALICE, BOB_1, BOB_2, CHARLY = range(4)

# initialize with two EPR pairs between Alice and Bob1 & between Bob2 and Charly, respectively
qc.psi = np.kron(KET_EPR, KET_EPR)

# teleport Bob1 via Bob2--Charly
teleport(qc, BOB_1, BOB_2, CHARLY)

# Alice--Charly are now entangled (and can do stuff like QKD, etc.)
assert np.allclose(qc.rdm(0, 3), np.outer(KET_EPR, KET_EPR))

Deutsch-Josza Algorithm

Truth tables of all constant and balanced functions of $n$ bits:

In [8]:
from itertools import combinations

def constant_fns(n):
    yield [0] * (2**n)
    yield [1] * (2**n)

def balanced_fns(n):
    inputs = range(2**n)
    for s in combinations(inputs, 2**(n-1)):
        yield [int(i in s) for i in inputs]

def rand_balanced_fn(n):
    inputs = range(2**n)
    s = np.random.choice(inputs, 2**(n-1), replace=False)
    return [int(i in s) for i in inputs]

print list(constant_fns(2))
print list(balanced_fns(2))
print rand_balanced_fn(2)
[[0, 0, 0, 0], [1, 1, 1, 1]]
[[1, 1, 0, 0], [1, 0, 1, 0], [1, 0, 0, 1], [0, 1, 1, 0], [0, 1, 0, 1], [0, 0, 1, 1]]
[0, 1, 0, 1]

Unitary oracle for a boolean function $f$:

In [9]:
def U_oracle(n, f):
    U = np.zeros(shape=(2**(n+1), 2**(n+1)))
    for x in range(2**n):
        for y in range(2):
            z = (y + f[x]) % 2
            U[x*2 + z, x*2 + y] = 1
    return U

U_oracle(1, [1, 0])
Out[9]:
array([[ 0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

Deutsch--Josza algorithm:

In [10]:
def deutsch_josza(n, f):
    # construct unitary oracle
    U = U_oracle(n, f)
    
    # initialize in state |0>^{\otimes n}  \otimes  |1>
    qc = QC(n+1)
    qc.NOT(n)
    
    # hadamards on all qubits
    for i in range(n+1):
        qc.H(i)
    
    # oracle
    qc.psi = np.dot(U, qc.psi)
    
    # hadamards on all qubits (except for last one, which does not matter anymore)
    for i in range(n):
        qc.H(i)
    
    # measure...
    bits = [qc.measure(i) for i in range(n)]
    is_const = all(b == 0 for b in bits)
    return is_const

assert deutsch_josza(2, [0, 0, 0, 0]) == True
assert deutsch_josza(2, [0, 1, 0, 1]) == False

Systematic Test:

In [12]:
%%time
n = 3
for f in constant_fns(n):
    assert deutsch_josza(n, f) == True
for f in balanced_fns(n):
    assert deutsch_josza(n, f) == False
    
n = 5
for _ in range(100):
    f = rand_balanced_fn(n)
    assert deutsch_josza(n, f) == False
Wall time: 994 ms