Deutsch Algorithm (Part 3): Implementation of oracle-based algorithms

Because the Deutsch algorithm is just a simple instance of the Deutsch-Jozsa algorithm where $n=2$, we’ll only discuss how to implement the general one. Let’s go.

Difficulty of reconstructing gates from matrices

This is the first time we encounter a so-called black-box function in an algorithm. It’s called a black box since we don’t how it calculate the output from the input, but we do know the output it returns. We can query the oracle if it is provided beforehand. Querying is easy but you shouldn’t take it for granted. In fact, creating an oracle from fundamental gate is, generally, extremely hard and a soft underbelly of quantum computing. This huge challenge tends to be overlooked because it’s ignored in the complexity evaluation of computation, which focuses on “classical” actions instead.

It’s known that oracles are indeed hidden operators, thus represented by hidden matrices. If we know the purpose of the operator we want to build, it’s not hard to write it out mathematically as a matrix. However, retaining a respective set of fundamental gates from the matrix is a monsterous challenge. There is yet a general process to perform reverse engineering to be discovered while verifying a matrix resulting from gates is a no difficult task. If needed, we have to deduce gates from matrices manually. For example, if we want to reverse engineer the matrix $$\frac{1}{2} \begin{bmatrix} 1 & 0 & 1 & 0 & 0 & 1 & 0 & -1 \\ 0 & 1 & 0 & 1 & 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 & 1 & 0 & 1 & 0 \\ 1 & 0 & -1 & 0 & 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 & 0 & -1 & 0 & 1 \\ 0 & 1 & 0 & 1 & -1 & 0 & 1 & 0 \\ 0 & 1 & 0 & -1 & -1 & 0 & -1 & 0 \\ 1 & 0 & -1 & 0 & 0 & -1 & 0 & -1 \end{bmatrix},$$ we can input it into Quirk, diagonalizing the matrix by adding some operations, then simplying them while hoping that the matrix of our interest isn’t created by any complex operations; “complex” means that it don’t contain rotation gates with odd angle parameters or any other kind of complicated transformations. In this case, we by luck retain the circuit for the mentioned matrix as follows

Deutsch-Jozsa oracle

As I want to implement and verify the Deutsch-Jozsa algorithm for various functions $f$, it’s not effective to construct gate-based oracles with each function reverse engineered. Now I want to illustrate a coding approach, more exploratory, that differs from my previous codes. I’ll use Python instead of Microsoft Q# to perform it.

Users are now allowed to input their own function $f$ as long as it’s either constant or balanced. The function maps each integer value from $0$ to $2^n$ to a binary value. For example, a function $f = [0,0,1,1]$ maps $00$ to $0$, $01$ to $0$, $10$ to $1$, and $11$ to $1$. This specific function is apparently balanced because the number of $0$ is equal the number of $1$. However, we’ll pretend not to know this and perform Deutsch-Jozsa algorithm to determine its constant-balanced property.

Manual coding

Since we cannot reconstruct operators equivalent to an arbitrary function, we have to implement the algorithm manually. First, I create the entry point for inputing and executing function.

1
2
3
4
5
6
7
8

def main():
    n = [2,3,3]  # Input the number of qubits
    f_map = [[0,0,1,1],
             [1,1,1,1,1,1,1,1],
             [1,0,0,1,1,0,1,0]]  # Input the mapping functions
    for index, value in enumerate(n):
        Deutsch_Jozsa(n[index], f_map[index])  # Algorithm executed here
 

1
2
3
4
5
6
7

def Deutsch_Jozsa(n, f_map): 
	num_qubits = n + 1  # Plus one qubit of the second register, can be called as ancilla qubit
	state_0 = np.array([[1],[0]])  # Standard state |0> as a column vector
	I_gate = np.array([[1,0], [0,1]])  # Identity gate
	X_gate = np.array([[0,1], [1,0]])  # NOT gate
	H_gate = np.array([[1,1], [1,-1]])/sqrt(2)  # Hadamard gate
 

1
2
3
4
5
6

def getTensor(matrices):
    product = matrices[0]
    for matrix in matrices[1:]:
        product = np.kron(product,matrix)  ## np.kron stands for Kronecker product, the official name for tensor product
    return product
 

Most important is creating a matrix $U_f$ for the entered $f$.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

def U(f_map):
    """Generate an oracle matrix based on the given function mapping."""
    # INSPIRED BY https://github.com/meownoid/quantum-python/blob/master/quantum.py

    U = np.zeros((2**num_qubits, 2**num_qubits)) # Start with a matrix of zeroes.
    
    # Quantum state looks like IN-IN-...-IN-ANCILLA
    for input_state in range(2**num_qubits): # For each possible input
        input_string = input_state >> 1 # remove ANCILLA
        output_qubit = (input_state & 1) ^ (f_map[input_string]) # remove IN, XOR with f(IN)
        output_state = (input_string << 1) + output_qubit # the full state, with new ANCILLA
        U[input_state, output_state] = 1 # set that part of U to 1
    return U
 

Let me elaborate on the above definition. Because of the ancilla qubit, the oracle is a $2^{n+1} \times 2^{n+1}$. We initial it as a zero matrix. Because all quantum states now is in the form of $|x_1x_1\dots x_n y\rangle$, we can get the input string $x_1x_1\dots x_n$ by shifting the state to the right by $1$ place. Next we can obtain $y$ by bitwise ANDing $x_1x_1\dots x_n y$ with $00\dots1$, then XOR $y$ with $f(x_1x_1\dots x_n)$ to get $y \oplus f(\mathbf{x})$. The result is pasted to the ancilla’s position by adding it with the input string shifted to the left $1$ place: $x_1x_1\dots x_n 0 + \left(y \oplus f(\mathbf{x})\right) = x_1x_1\dots x_n \left(y \oplus f(\mathbf{x}) \right)$. This is the state $|\mathbf{x}\rangle|y\oplus f(\mathbf{x})\rangle$ we want to create using the oracle. To finish, we assign U[input_state, output_state] = 1 to indicate that the “input state” is transformed into the “output state” with absolute certainty. The matrix $U$ is completed by the time we finish repeating the procedure traversing through all quantum states possibly created by $n+1$ qubits.

The final subroutine in our program is used for measuring and displaying result.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

def measure(n, state):
    measurement = np.zeros(2**n)  # Initialize measurement result for n qubits in the first register
    for index, value in enumerate(state):
        measurement[index >> 1] += value * value  ## As the ancilla qubit is discarded, probabilities of the same kind, ie 100 and 101 will be combined

    # Last step: Determine the type of function f
    # f is constant if the probability of measuring |0> is positive
    if (abs(measurement[0]) > 1e-10): 
        print("The function is constant.")
    else:
        print("The function is balanced.")
 

The definitions listed above are sufficient for us to perform the algorithm. You may need to see the circuit for the algorithm before getting your hand in. Please note that there are some trivial differences between the states on the circuit and in the implementation code.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

def Deutsch_Jozsa(n, f_map):
    num_qubits = n + 1  # Plus one qubit and the second register, can be called as ancilla qubit
    state_0 = np.array([[1],[0]])  # Standard state |0> as a column vector
    I_gate = np.array([[1,0], [0,1]])  # Identity gate
    X_gate = np.array([[0,1], [1,0]])  # NOT gate
    H_gate = np.array([[1,1], [1,-1]])/sqrt(2)  # Hadamard gate
    
    ancilla = np.dot(X_gate, state_0)  # Create state |1> assigned to the ancilla
    
    # Create the a Hadamard transformation for all qubits and the state |ψ_0> 
    listStates = []
    listGates_H = []
    for i in range(n):
        listStates.append(state_0)
        listGates_H.append(H_gate)
    listStates.append(ancilla)
    listGates_H.append(H_gate)
    psi_0 = getTensor(listStates)
    composite_H = getTensor(listGates_H)
    
    # |ψ_1> is the dot product of the Hadamard transformation and |ψ_0>  
    psi_1 = np.dot(composite_H, psi_0)
    
    # Apply the oracle to |ψ_1>
    psi_2 = np.dot(U(n, f_map), psi_1)
    
    # H on all again
    psi_3 = np.dot(composite_H, psi_2)
    
    measure(n, psi_3)
 

It’s all done. Now let’s check if everything is okay.

It works as expected. Just free to experiment with other functions $f$ on yourself. Check out the source code archived in my GitHub directory.