Tag: QHack 2023

QHack 2023: 4 week prep-challenge, day 23, Draper adder

I went for another QOSF task (Oct. 2022), there I tackled the seemingly easy task of implementing a Draper adder, as an intermediate step towards a multiplier. However, I got lost in index engineering. At this point, I am not satisfied with the solution, but the tests check out, so it is usable.

I added two quantum registers; instead of the much easier multiplication of a quantum with a classical register.

import pennylane as qml
from pennylane import numpy as np

reg_size = 4
n_wires = 2*reg_size

dev = qml.device("default.qubit", wires=n_wires)

def base_encode(a, wires_offset=0):
    a_bin = np.binary_repr(a, reg_size)
    a_bin = a_bin if wires_offset else reversed(a_bin)
    
    for i, ai in enumerate(a_bin):
        if ai == "1":
            qml.PauliX(i + wires_offset)

@qml.qnode(dev)
def circuit(a, b):
    
    # encode basis state
    base_encode(a)
    base_encode(b, wires_offset=reg_size)
    
    # apply QFT to one register
    qml.QFT(wires=range(reg_size, n_wires))
    
    # in this implementation, I pretend that I do not know the value of b, 
    # this could be the case, if this addition was amidst a longer calculation
    for i in range(reg_size): # ith significant bit (target)
        for j in range(reg_size-i): # control bits
            trgt_idx = n_wires - 1 - i
            ctrl_idx = trgt_idx - reg_size - j
            angle = np.pi/(2**j)
            qml.ControlledPhaseShift(angle, wires=[ctrl_idx, trgt_idx])
    
    # invert QFT
    qml.adjoint(qml.QFT)(wires=range(reg_size, n_wires))
    
    return qml.probs(wires=range(reg_size, n_wires))# [qml.expval(qml.PauliZ(i)) for i in range(reg_size, n_wires)]

As already mentioned, I wouldn’t trust the implementation, if the tests didn’t check out, but they do.

circ_draw = qml.draw(circuit)

res = np.array([])
for a in range(reg_size**2 - 1):
    for b in range(reg_size**2 - 1 -a):
        sample = circuit(a, b)
        expected_result = a+b
        if np.argmax(sample)!=expected_result:
            print(f"a: {a}, b: {b}, sample: {np.argmax(sample)}")
            print(circ_draw(a, b))
            break
        
        res = np.append(res, [np.argmax(sample)==a+b])
        
print(f"The algorithm is correct! {len(res)} tests successful!" if np.all(res) else "There is an error")

QHack 2023: 4 week prep-challenge, day 22, Constructing diagonal operators, pt. 2

After implementing the basic level task of creating a 2×2 diagonal unitary operator (i.e. only possible values in the diagonal are +1 and -1), today I am generalizing the algorithm to create arbitrary diagonal unitary operators, with complex eigenvalues. For this purpose I used the multi-controlled P-gate, which differs slightly from the RZ-gate:

\[ RZ(\phi) = \left( \begin{array}{rr} e^{-i \frac{\phi}{2}} & 0 \\ 0 & e^{i \frac{\phi}{2}} \end{array} \right) \text{, and} \\ P(\phi) = \left( \begin{array}{rr} 1 & 0 \\ 0 & e^{i \phi} \end{array} \right) \]

In code, most remained the same to last day’s solution, I just had to modify the marking function from a CZ to a multi-controlled P-gate.

import pennylane as qml
from pennylane import numpy as np

num_wires = 3
dev = qml.device("default.qubit", wires=num_wires)

def prepare_address(address):
    """Flip bits, to prepare addressation by oracle, or the like.
    """
    bin_address = np.binary_repr(address, width=num_wires)
    
    for i, addr_bit in enumerate(bin_address):
        if addr_bit == "0":
            qml.PauliX(i)


@qml.qnode(dev)
def prepare_unitary(u_eigenvalues):
    """Prepare a circuit, that realizes a unitary transformation, with a diagonal matrix representation, with the provided eigenvalues.
    
    Args: 
        - u_eigenvalues (np.array(complex)): Array that holds a list of unitary eigenvalues of the form np.exp(1.j*alpha), i.e. complex, norm 1
    
    Return: 
        - qml.State
    """
    assert np.all(np.isclose(np.absolute(u_eigenvalues), 1.)), "Invalid eigenvalues provided, do not have norm of 1"
    
    angles = np.angle(u_eigenvalues)
    
    for i, angle in enumerate(angles):
        prepare_address(i)
        qml.ctrl(qml.PhaseShift(angle, num_wires-1), range(num_wires-1))
        prepare_address(i)
        
    return qml.state()


def adjust_global_phase_of_diag_unitary(trg_arr, ref_value=1.+0.j):
    """Rotate the target diagonal matrix such, that it's global aligns with that of a source matrix. 
    If no source matrix is provided, the target matrix is rotated such, that the first element is real.
    
    Args:
        - trg_arr (np.array): approx. diagonalized unitary matrix (i.e. the norm of the diagonal elements approximate to 1, all other elements approx. to 0)
        - ref_value (complex): the reference value, holding the target phase; defaults to 1.+0.j 
        
    Return:
    
    """
    assert np.all(np.isclose(np.absolute(trg_arr), np.identity(len(trg_arr)))), "Invalid diagonalized unitary provided."
    assert np.isclose(np.absolute(ref_value), 1.), "Invalid reference value provided; must be complex with norm 1"
    
    angle_trg = np.angle(trg_arr.flat[0])
    angle_src = np.angle(ref_value)
    return trg_arr * np.exp(1.j*(angle_src - angle_trg))

.. and a small adjustment, to the test data + checking procedure.

# prepare test data
rng = np.random.default_rng(seed=314159)
rand_angles = rng.random(2**num_wires) * 2. * np.pi
eigenvalues = np.exp(1.j*rand_angles)

get_matrix = qml.matrix(prepare_unitary)
matrix = get_matrix(eigenvalues)
phase_adjusted_matrix = adjust_global_phase_of_diag_unitary(matrix, eigenvalues[0])

result_str = "correct!" if np.all(np.isclose(phase_adjusted_matrix, np.diagflat(eigenvalues))) else "INCORRECT! "

print("The result is", result_str)

The tests check out, all fine. It remains to exchange the MCP-gate with a series of two qubit gates to solve the challenge in it’s entirety.

QHack 2023: 4 week prep-challenge, day 21, Constructing unitary diagonal operators, pt. 1

I dabbled into the first QOSF monthly challenge, first published in November 2020. The task is to create a circuit, that implements a unitary operator in it’s eigenbasis, i.e. a diagonal matrix, where the diagonal elements have norm 1.

The task is split into three difficulty levels, today I’m tackling the easiest level: Create an algorithm that constructs a circuit, that translates into a 4×4 diagonal matrix with real entries. Now, the only two allowed real values for the diagonal elements are +1 and -1.

I set out to solve the more difficult problem, to create a circuit that allows for complex values in the diagonal, therefore, the logic is formulated slightly more general (except for the crucial gates, there the limits are hard in place).

import pennylane as qml
from pennylane import numpy as np

num_wires = 2
dev = qml.device("default.qubit", wires=num_wires)

def prepare_address(address):
    """Flip bits, to prepare addressation by oracle, or the like.
    """
    bin_address = np.binary_repr(address, width=num_wires)
    
    for i, addr_bit in enumerate(bin_address):
        if addr_bit == "0":
            qml.PauliX(i)


@qml.qnode(dev)
def prepare_unitary(u_eigenvalues):
    """Prepare a circuit, that realizes a unitary transformation, with a diagonal matrix representation, with the provided eigenvalues.
    
    Args: 
        - u_eigenvalues (np.array(complex)): Array that holds a list of unitary eigenvalues of the form np.exp(1.j*alpha), i.e. complex, norm 1
    
    Return: 
        - qml.State
    """
    assert np.all(np.isclose(np.absolute(u_eigenvalues), 1.)), "Invalid eigenvalues provided, do not have norm of 1"
    
    angles = np.angle(u_eigenvalues)
    
    for i, angle in enumerate(angles):
        prepare_address(i)
        if angle:
            qml.CZ(wires=range(num_wires))
        prepare_address(i)
        
    return qml.state()

My previous attempts lead me to write this little helper function, to adjust the global phase of a result matrix, for easier comparison….


def adjust_global_phase_of_diag_unitary(trg_arr, ref_value=1.+0.j):
    """Rotate the target diagonal matrix such, that it's global aligns with that of a source matrix. 
    If no source matrix is provided, the target matrix is rotated such, that the first element is real.
    
    Args:
        - trg_arr (np.array): approx. diagonalized unitary matrix (i.e. the norm of the diagonal elements approximate to 1, all other elements approx. to 0)
        - ref_value (complex): the reference value, holding the target phase; defaults to 1.+0.j 
        
    Return:
    
    """
    assert np.all(np.isclose(np.absolute(trg_arr), np.identity(len(trg_arr)))), "Invalid diagonalized unitary provided."
    assert np.isclose(np.absolute(ref_value), 1.), "Invalid reference value provided; must be complex with norm 1"
    
    angle_trg = np.angle(trg_arr.flat[0])
    angle_src = np.angle(ref_value)
    return trg_arr * np.exp(1.j*(angle_src - angle_trg))

Eventually, I tested the code, with a few configurations…

eigenvalues = np.array([1,-1,-1,1])

get_matrix = qml.matrix(prepare_unitary)
matrix = get_matrix(eigenvalues)
phase_adjusted_matrix = adjust_global_phase_of_diag_unitary(matrix, eigenvalues[0])

result_str = "correct!" if np.all(np.isclose(phase_adjusted_matrix, np.diagflat(eigenvalues))) else "INCORRECT! "

print("The result is", result_str)

… and I end up with correct results.

As a next step, I will generalize this algorithm to allow for arbitrary dimensions.

QHack 2023: 4 week prep challenge, day 20

Today, I tried to solve further exercises, but the remaining problem statements in the QHack 2022 Games category don’t play into my strengths. I will need to improve my theoretical understanding , to progress on them.

For the next few days, I will leave them aside, and try to apply myself to other problem sets. I believe, I can learn plenty from the QOSF monthly challenges.

QHack 2023: 4 week prep-challenge, day 19, A reading night

Recently, I was feeling quite comfortable with solving these challenges, to me this is a sign, that I can set the bar a little higher or mix-up my activities to increase my learning. The QHack comprises talks, the coding challenge and an open hackathon. The coding challenges will be of the same format, as my exercises of the last days, the hackathon invites for free-form contributions and offer much bigger rewards. Now, I deem it beneficial to broaden my view and get some inspiration reading up a little bit.

I read an insightful industry and market analysis by BCG, providing forecasts for market sizes and areas for near term applications.

Further, I found a consortium of big German companies, evaluating quantum computing: QUTAC. Some of the researched applications are outlined in this paper. Over the next few days, I want to study this overview more carefully.

QHack 2023: 4 week prep-challenge, day 18, QAOA with custom Hamiltonian

The final exercise on the QML track, included some bug-hunting for me, but this made it much more rewarding, to eventually solve it.

The basic task was conceivably easy: given some data, create a hamiltonian and use a QAOA optimization scheme to obtain a ground state approximation. The problem described some special graph covering problem. The mistakes that took some efforts to resolve where

  • Mistake 1: PennyLane is exceptionally greedy, when arrays appear along the optimization scheme. The way I constructed the hamiltonian included numpy arrays for indexing qubits, PennyLane tried to differentiate the problem by them and vary them subsequently – that lead to exceptions, and it was a real pain to find the source of the error.
  • Mistake 2: I tried to construct operators, instead of decomposing the occupation operator into Pauli matrices. This might work somehow, but I had an easier time decomposing the operators eventually.
  • Mistake 3: I omitted the constant terms, the optimal configuration is not affected by the constant term, but my tests relied on the final energy.

The code for the Hamiltonian

def hamiltonian_coeffs_and_obs(graph):
    """Creates an ordered list of coefficients and observables used to construct
    the UDMIS Hamiltonian.

    Args:
        - graph (list((float, float))): A list of x,y coordinates. e.g. graph = [(1.0, 1.1), (4.5, 3.1)]

    Returns:
        - coeffs (list): List of coefficients for elementary parts of the UDMIS Hamiltonian
        - obs (list(qml.ops)): List of qml.ops
    """

    num_vertices = len(graph)
    E, num_edges = edges(graph)
    u = 1.35
    obs = []
    coeffs = []

    # single terms
    for i in range(num_vertices):
        obs.append(qml.PauliZ(wires=i))
        coeffs.append(-.5)
    
    # constant term
    obs.append(qml.Identity(wires=0))
    coeffs.append(-.5*num_vertices)
    
    # interaction terms
    pairs = np.argwhere(E>0)
    for pair in pairs:
        i = int(pair[0])
        j = int(pair[1])
        obs.append(qml.PauliZ(wires=i) @ qml.PauliZ(wires=j))
        coeffs.append(.25*u)
        coeffs[i] += .25*u
        coeffs[j] += .25*u
        coeffs[num_vertices] += .25*u
    
    return coeffs, obs

The tests check out, success!

QHack 2023: 4 week prep challenge, day 17, QRAM

Today, I am working on the next exercise on the QML track. The goal is to create a QRAM, that entangles phase-information with address bits. A good description of QRAM that I know of, stems from an IBM Quantum hackathon from 2020. The main idea is to encode the addresses, and entangle the address register with the data register via a controlled rotation.

My code …

def prep_address(address):
    bin_address = np.binary_repr(address, width=3)
        
    for i, add_bit in enumerate(bin_address):
        if add_bit == "0":
            qml.PauliX(wires=i)

dev = qml.device("default.qubit", wires=range(4))

@qml.qnode(dev)
def circuit():
    qml.broadcast(qml.Hadamard, range(3), pattern="single")
        
    for i, theta in enumerate(thetas):
        prep_address(i)
        qml.ctrl(op=qml.RY(theta, wires=3), control=range(3))
        prep_address(i)
    return qml.state()

print(qml.draw(circuit)())

Instinctively, I would have reversed the address bit-order, but as long as you keep it consistent, it doesn’t matter. I also didn’t care about saving gates or anything, it’s just important to calculate the addresses correctly.

QHack 2023: 4 week prep challenge, day 16, Variational Quantum Classifier

I continued with the next exercise on the QML track, the task was to create a VQC (Variational Quantum Classifier) to classify Ising-spin systems into ordered and disordered states, with an accuracy of at least 90%.

The Ising model, is a simplified model of a ferromagnet, describing it’s magnetic properties. A ferromagnet displays spontaneous magnetization (i.e. parallel alignment of it’s atoms spins) below its critical temperature. In our dataset, the ferromagnet is represented by five spins arranged in a linear setup.

Variational Quantum Classifier

A Variational Quantum Classifier differs in two essential points from a conventional NN classifier. Firstly, the input data needs to be encoded into qubits. Secondly, the hypothesis function is a parameterized quantum circuit. These parameters are subsequently optimized using classical optimizers.

Implementation

My implementation heavily draws on the PennyLane VQC tutorial, I made a few adaption though.

First we need to create our hypothesis, similarly to the tutorial I chose strongly entangling layers. I only want to distinguish two classes, and I do not want to create a quantum-classical hybrid pipeline (apart from the classically added bias term). Therefore, I only need the expectation value of a single qubit, it doesn’t matter which one I choose.

@qml.qnode(dev)
def circuit(weights, x):
    """VQC architecture, encodes input data and performs a series of strongly entangling layers, to perform the classification.
    
    Args:
        - weights (np.ndarray): weights for the entangling layers
        - x (np.ndarray): spin configurations

    Returns:
        - prediction (float): expectation value of one qubit, after the entangling layers
    """
    qml.BasisState(x, wires=range(num_wires))
    qml.StronglyEntanglingLayers(weights=weights, wires=range(4))

    return qml.expval(qml.PauliZ(0))

def variational_classifier(weights, bias, x):
    return circuit(weights, x) + bias

Next, I initialize my parameters and my optimizer. I opted for three layers, and a smaller step size.

np.random.seed(0)
n_layers = 3
shape = qml.StronglyEntanglingLayers.shape(n_layers=n_layers, n_wires=num_wires)
weights = np.random.random(size=shape)
bias = np.array(0.0, requires_grad=True)
    
opt = NesterovMomentumOptimizer(0.1)
batch_size = 10

Finally, the optimization loop …

for it in range(25):

    # Update the weights by one optimizer step
    batch_index = np.random.randint(0, len(ising_configs), (batch_size,))
    X_batch = ising_configs[batch_index]
    Y_batch = labels[batch_index]
    weights, bias, _, _ = opt.step(cost, weights, bias, X_batch, Y_batch)

    # Compute accuracy
    predictions = [np.sign(variational_classifier(weights, bias, x)) for x in ising_configs]
    acc = accuracy(labels, predictions)

    print(
        "Iter: {:5d} | Cost: {:0.7f} | Accuracy: {:0.7f} ".format(
            it + 1, cost(weights, bias, ising_configs, labels), acc
        )
    )

I reach an accuracy of 93.6%, a good result for this toy problem.

QHack 2023: 4 week prep challenge, day 15, SWAP Test

I did another basic exercise on the QML track, but honestly, I was struggling a lot with this one. The hard-earned victories are both satisfying and instructive. But first a quick outline of the exercise.

The challenge was to estimate the inner product of two feature vectors. A distance measure is derived from this quantity and is used to perform a k-NN (k-nearest neighbor) algorithm – in this configuration an artificial use case.

The inner product can be estimated with the SWAP-test.

SWAP test circuit, taken from Wikipedia.

The inner product estimate is calculated from

\[ P(\text{control qubit} = 0) = \frac{1}{2} + \frac{1}{2}\lvert\langle \psi | \phi \rangle\rvert^2 \]

Implementation

In this implementation, each feature vector is normalized, and since each vector has dimension two, a single qubit suffices to encode the data.

from pennylane import numpy as np
import pennylane as qml


def distance(A, B):
    """Function that returns the distance between two vectors.

    Args:
        - A (list[int]): person's information: [age, minutes spent watching TV].
        - B (list[int]): person's information: [age, minutes spent watching TV].

    Returns:
        - (float): distance between the two feature vectors.
    """

    def prep_state():
        """Function that encodes the two 2D feature vectors into the amplitudes of one qubit each."""
        
        tan_thetaA = A[0]/A[1]
        thetaA = np.arctan(tan_thetaA)
        tan_thetaB = B[0]/B[1]
        thetaB = np.arctan(tan_thetaB)
        
        qml.RY(2.*thetaA, wires=0)
        qml.RY(2.*thetaB, wires=1)
    
    
    dev = qml.device("default.qubit", wires=3)
    @qml.qnode(dev)
    def circuit():
        """Quantum circuit that encodes two 2D-feature vectors and performs a SWAP test.
        
        Returns:
            - ([float, float]): probabilities of states measured on the control qubit.
        """
        prep_state()
        
        qml.Hadamard(wires=2)
        qml.CSWAP(wires=range(2,-1,-1))
        qml.Hadamard(wires=2)
        
        return qml.probs(2)

    # print(qml.draw(circuit)())
    
    tmp = circuit()[0] # .5 + .5<A,B>^2
    A_times_B = np.sqrt(2.*tmp-1.)
    
    return np.sqrt(2*(1.-A_times_B))

The provided tests check out.

Lessons learned

I will seize this opportunity to reflect upon my solution strategies. The problem statement was a little less concrete than the ones of previous exercises. In the starter code two hints were pointing at the SWAP Test and the PennyLane amplitude encoding functionality. After an initial implementation failed, I did some more research and found a SWAP-test implementation on a PennyLane forum. This was more confusing then helpful, I believe the core confusion being how the normalization is done in this encoding strategy.

Anyways, the clean bottom up solution, along with some debugging print-outs eventually worked out.

QHack 2023: 4 week prep challenge, day 14, First steps in QAOA

Today, I’m not doing an exercise from the repository, but I looking into one of the basic goto algorithms in the world of QC, the Quantum Approximate Optimization Algorithm (QAOA). The code snippets are taken from the very instructive tutorial published here.

A common application of this algorithm is to solve so called Quadratic Unbound Binary Optimization (QUBO) problems. The goal of these problems is to find the binary feature vector x, that minimizes an at most quadratic function f:

\[ f_Q(x) = x^TQx = \sum^n_{i=1}\sum^n_{j=i}Q_{ij}x_ix_j \]

Therefore, the solution to the problem is a series of zeros and ones. As a matter of fact, a large class of problems can be formulated as such a QUBO problem.

QAOA algorithm

The QAOA algorithm employs a very simple idea. A quantum system is first initialized in a promising state (ideally close to the expected global minimum). Then a series of time evolutions is performed, the time evolutions apply alternating the Hamiltonian encoding the problem and a mixer Hamiltonian, to leave meta-stable local optima states.

\[ U(\gamma, \alpha) = \exp(-\alpha_nH_M)\exp(-i\gamma_nH_C)\cdots\exp(-\alpha_1H_M)\exp(-i\gamma_1H_C) \]

In the tutorial, the parameters Alpha and Gamma were optimized using a VQE scheme.

Implementation in PennyLane

First, you’ve got to define a Hamiltonian that characterizes your problem, here’s a toy Hamiltonian

import pennylane as qml

H_P = qml.Hamiltonian(
    [1, 1, 0.5],
    [qml.PauliX(0), qml.PauliZ(1), qml.PauliX(0) @ qml.PauliX(1)]
)
\[ H_P = \sigma_x^{(0)} + \sigma_z^{(1)} + \frac{1}{2} \sigma_x^{(0)} \otimes \sigma_x^{(1)} \]

and your mixer Hamiltonian

\[ H_M = \frac{1}{2} \left( \sigma_x^{(0)} + \sigma_x^{(1)} \right) \]

Then define a time evolution step

def qaoa_layer(params):
    qaoa.cost_layer(params[0], H_P)
    qaoa.mixer_layer(params[1], H_C)

Now, you can combine the initialization and the time evolution sequence to a single circuit

wires = range(4)
depth = 2

def circuit(params, **kwargs):
    for w in wires:
        qml.Hadamard(wires=w)
    qml.layer(qaoa_layer, depth, params[0], params[1])

And having a cost function in place, that evaluates the problem function …

dev = qml.device("qulacs.simulator", wires=wires)

@qml.qnode(dev)
def cost_function(params):
    circuit(params)
    return qml.expval(cost_h)

you can run an optimization procedure, to find a good parameterization for the time evolutions.

optimizer = qml.GradientDescentOptimizer()
steps = 70
params = np.array([[0.5, 0.5], [0.5, 0.5]], requires_grad=True)

for i in range(steps):
    params = optimizer.step(cost_function, params)

print("Optimal Parameters")
print(params)

Interesting observation

The procedure outlined in the tutorial yield good results, providing the highest probabilities for sensible configurations. However, after I went on to modify the algorithm, and increased the number of evolutions from 2 to 5, a less favorable configuration yielded the highest probability. I’ll have to do some research, into why this happens.