Tutorial - Quantum support vector regression for disability insurance

B. Djehiche and B. Löfdahl – KTH, Department of Mathematics / SEB – Risks 2021, 9(12), 216; https://doi.org/10.3390/risks9120216

Hands-on assistance: https://enccs.se/anastasiia-andriievska

# install necessary packages via pip
# import the packages
from qiskit_iqm import IQMProvider
import networkx as nx
from qiskit import QuantumCircuit, QuantumRegister, execute, transpile
from qiskit.visualization import *
import pylatexenc
import pandas as pd
import numpy as np
from qiskit import Aer, execute
from qiskit.circuit import QuantumCircuit, QuantumRegister, Parameter, ParameterVector
from qiskit_machine_learning.kernels import QuantumKernel
from qiskit.utils import QuantumInstance
from qiskit.tools.monitor import job_monitor
from qiskit.tools.jupyter import *
from qiskit.visualization import *
#from ibm_quantum_widgets import *
import qiskit.quantum_info as qi

import matplotlib.pyplot as plt
from sklearn.svm import SVR

import qiskit
qiskit.__version__

'0.25.2'

First, we create some dummy data. The x-variable holds the available data. x[1] represents gender (0 or 1). x[2] represents age, normalized to the [0,1] interval. The y-variable is the value that we seek to model. It represents the logistic transition probability between two states. The w-variable represents the weight that we assign to each data point.

nsamp = 100
x = np.array( [[0.0, x/nsamp] for x in range(nsamp)] + [[1.0, x/nsamp] for x in range(nsamp)] )
y = np.array([-4.5 + 0.2*a[0] + a[1]*(1+a[0]*0.2) + 0.1*np.random.normal() for a in x])

p = 1/(1+np.exp(-y))

n = [int(200 + 600*xj[1]/0.7) if xj[1] <= 0.7 else int(800 - 300*(xj[1]-0.7)/0.3) for xj in x]

d = [nj*pj for nj,pj in zip(n,p)]

We can use the Pandas library to create a DataFrame and then save it as a CSV file, which support column headers. This will create a CSV file with each column labeled as ‘gender’, ‘age group (normalized)’, ‘population size’, and ‘number of transitions’ respectively. The index=False argument is used to prevent pandas from writing row indices into the file.

# Create a DataFrame
df = pd.DataFrame({
    'gender': [a[0] for a in x],
    'age group (normalized)': [a[1] for a in x],
    'population size': n,
    'number of transitions': d
})

# Write the DataFrame to a CSV file
df.to_csv('fdata.csv', index=False)  # index=False argument is used to prevent pandas from writing row indices into the CSV file

One can read the data back from the file like this:

data = pd.read_csv('fdata.csv')

Transform and inspect data

x = np.array( [[a,b] for a,b in zip(data['gender'], data['age group (normalized)'])] )
print(x)

[[0.   0.  ]
 [0.   0.01]
 [0.   0.02]
 [0.   0.03]
 [0.   0.04]
 [0.   0.05]
 [0.   0.06]
 [0.   0.07]
 [0.   0.08]
 [0.   0.09]
 [0.   0.1 ]
 [0.   0.11]
 [0.   0.12]
 [0.   0.13]
 [0.   0.14]
 [0.   0.15]
 [0.   0.16]
 [0.   0.17]
 [0.   0.18]
 [0.   0.19]
 [0.   0.2 ]
 [0.   0.21]
 [0.   0.22]
 [0.   0.23]
 [0.   0.24]
 [0.   0.25]
 [0.   0.26]
 [0.   0.27]
 [0.   0.28]
 [0.   0.29]
 [0.   0.3 ]
 [0.   0.31]
 [0.   0.32]
 [0.   0.33]
 [0.   0.34]
 [0.   0.35]
 [0.   0.36]
 [0.   0.37]
 [0.   0.38]
 [0.   0.39]
 [0.   0.4 ]
 [0.   0.41]
 [0.   0.42]
 [0.   0.43]
 [0.   0.44]
 [0.   0.45]
 [0.   0.46]
 [0.   0.47]
 [0.   0.48]
 [0.   0.49]
 [0.   0.5 ]
 [0.   0.51]
 [0.   0.52]
 [0.   0.53]
 [0.   0.54]
 [0.   0.55]
 [0.   0.56]
 [0.   0.57]
 [0.   0.58]
 [0.   0.59]
 [0.   0.6 ]
 [0.   0.61]
 [0.   0.62]
 [0.   0.63]
 [0.   0.64]
 [0.   0.65]
 [0.   0.66]
 [0.   0.67]
 [0.   0.68]
 [0.   0.69]
 [0.   0.7 ]
 [0.   0.71]
 [0.   0.72]
 [0.   0.73]
 [0.   0.74]
 [0.   0.75]
 [0.   0.76]
 [0.   0.77]
 [0.   0.78]
 [0.   0.79]
 [0.   0.8 ]
 [0.   0.81]
 [0.   0.82]
 [0.   0.83]
 [0.   0.84]
 [0.   0.85]
 [0.   0.86]
 [0.   0.87]
 [0.   0.88]
 [0.   0.89]
 [0.   0.9 ]
 [0.   0.91]
 [0.   0.92]
 [0.   0.93]
 [0.   0.94]
 [0.   0.95]
 [0.   0.96]
 [0.   0.97]
 [0.   0.98]
 [0.   0.99]
 [1.   0.  ]
 [1.   0.01]
 [1.   0.02]
 [1.   0.03]
 [1.   0.04]
 [1.   0.05]
 [1.   0.06]
 [1.   0.07]
 [1.   0.08]
 [1.   0.09]
 [1.   0.1 ]
 [1.   0.11]
 [1.   0.12]
 [1.   0.13]
 [1.   0.14]
 [1.   0.15]
 [1.   0.16]
 [1.   0.17]
 [1.   0.18]
 [1.   0.19]
 [1.   0.2 ]
 [1.   0.21]
 [1.   0.22]
 [1.   0.23]
 [1.   0.24]
 [1.   0.25]
 [1.   0.26]
 [1.   0.27]
 [1.   0.28]
 [1.   0.29]
 [1.   0.3 ]
 [1.   0.31]
 [1.   0.32]
 [1.   0.33]
 [1.   0.34]
 [1.   0.35]
 [1.   0.36]
 [1.   0.37]
 [1.   0.38]
 [1.   0.39]
 [1.   0.4 ]
 [1.   0.41]
 [1.   0.42]
 [1.   0.43]
 [1.   0.44]
 [1.   0.45]
 [1.   0.46]
 [1.   0.47]
 [1.   0.48]
 [1.   0.49]
 [1.   0.5 ]
 [1.   0.51]
 [1.   0.52]
 [1.   0.53]
 [1.   0.54]
 [1.   0.55]
 [1.   0.56]
 [1.   0.57]
 [1.   0.58]
 [1.   0.59]
 [1.   0.6 ]
 [1.   0.61]
 [1.   0.62]
 [1.   0.63]
 [1.   0.64]
 [1.   0.65]
 [1.   0.66]
 [1.   0.67]
 [1.   0.68]
 [1.   0.69]
 [1.   0.7 ]
 [1.   0.71]
 [1.   0.72]
 [1.   0.73]
 [1.   0.74]
 [1.   0.75]
 [1.   0.76]
 [1.   0.77]
 [1.   0.78]
 [1.   0.79]
 [1.   0.8 ]
 [1.   0.81]
 [1.   0.82]
 [1.   0.83]
 [1.   0.84]
 [1.   0.85]
 [1.   0.86]
 [1.   0.87]
 [1.   0.88]
 [1.   0.89]
 [1.   0.9 ]
 [1.   0.91]
 [1.   0.92]
 [1.   0.93]
 [1.   0.94]
 [1.   0.95]
 [1.   0.96]
 [1.   0.97]
 [1.   0.98]
 [1.   0.99]]

# raw transition frequency per group, used for plots
p_raw = [a/b for a,b in zip(d,n)]

# logistic probability, used as input to SVR optimization
y = np.array([np.log( p / (1-p) ) for p in p_raw])

# calculate weights for SVR based on population size
w = np.array([1./nj for nj in n])

# normalize weights (equivalent to changing the hyperparameters of the SVR)
w = w/sum(w)*len(x)

fig, axs = plt.subplots(1, 1, figsize=(10, 5))
axs.plot([xj[1] for xj,p in zip(x, p_raw) if xj[0] == 0], [p for xj,p in zip(x, p_raw) if xj[0] == 0], 'or')
axs.plot([xj[1] for xj,p in zip(x, p_raw) if xj[0] == 1], [p for xj,p in zip(x, p_raw) if xj[0] == 1], 'ob')
axs.set_title("raw transition frequencies")
plt.show()

The first step of the modelling is to define the mapping from data space to our quantum feature space. This is done by setting up the qc circuit. The circuit below is a very simplified two-variable circuit. Feel free to experiment!

# set up 'kernel circuit'
p = ParameterVector('x', length=2)
qc = QuantumCircuit(2)

qc.ry(p[0]*np.pi, 0)
qc.ry(p[1]*np.pi, 1)

qc.draw('mpl')

Exercise 1

Experimenting with quantum circuits involves changing the structure of the circuit or the parameters used. Here are a few ways you can experiment with the given circuit:

1. Change the gates: Instead of using ry gates, you could use other gates like rx, rz, h (Hadamard), etc.

qc.rx(p[0]*np.pi, 0)
qc.rz(p[1]*np.pi, 1)

2. Add more gates: You can add more gates to your circuit to make it more complex.

qc.ry(p[0]*np.pi, 0)
qc.ry(p[1]*np.pi, 1)
qc.h(0)
qc.cx(0, 1) # CNOT gate

3. Change the parameters: Instead of multiplying the parameters by pi, you could multiply them by other constants or use more complex functions of the parameters.

qc.ry(p[0], 0)
qc.ry(np.sin(p[1]), 1)

4. Increase the number of qubits: You can increase the number of qubits in your circuit and see how that affects your results.

Remember, when experimenting with quantum circuits, it’s important to understand what each gate does and how changing the circuit might affect your results. Happy experimenting! 😊

A few general guidelines for experimenting with quantum circuits:

1. Understand the problem: Before you start experimenting, make sure you understand the problem you’re trying to solve and how a quantum circuit could help solve it.

2. Start simple: Start with a simple circuit and gradually add complexity. This makes it easier to understand what each part of your circuit is doing.

3. Use a variety of gates: Different gates can have different effects on your quantum state, so try using a variety of gates in your experiments.

4. Test on a simulator first: Before running your circuit on a real quantum computer, test it on a simulator to make sure it’s working as expected.

5. Analyze your results: After running your circuit, analyze the results to see if they make sense and if they’re helping you solve your problem.

6. Iterate: Experimenting with quantum circuits often involves a lot of trial and error. Don’t be afraid to iterate on your design and try new things.

In the context of the Quantum Support Vector Regression (QSVR) for disability insurance that you’re working on, you might want to experiment with different ways of encoding your data into the quantum state or different types of quantum kernels. Remember, the goal is to find a quantum circuit that can effectively model your data and help you make accurate predictions.

Hints

1. x[0] is a binary variable representing gender. Hence, flipping the first qubit from |0> to |1> based on gender will maximize the distance (wrt only the first qubit) between two subpopulations belong to different genders.

2. x[1] is a representing age. For two sub-populations of similar age, the rotations of the second qubit will be similar, causing a low distance (wrt only the second qubit).

3. One possible extension is to add a third qubit that would capture a first order interaction effect between x[0] and x[1], e.g. a cy or cx gate. We could also incorporate this by adding more gates to qubits 1 or 2.

# Test code here
# Then rerun the following cells 

The QuantumKernel function prepares the kernel estimation by setting up inner product circuits, based on each combination of two data points.

zz_kernel = QuantumKernel(feature_map=qc, quantum_instance=Aer.get_backend('statevector_simulator'))

Then we specify the quantum instance:

zz_circuit = zz_kernel.construct_circuit(ParameterVector('x_i', length=2),ParameterVector('x_j', length=2))
zz_circuit.decompose().draw(output='mpl')

We evaluate the kernel matrix based on our x data.

matrix_train = zz_kernel.evaluate(x_vec=x)

fig, axs = plt.subplots(1, 1, figsize=(10, 5))
pos = axs.imshow(np.asmatrix(matrix_train),
              interpolation='nearest', origin='upper', cmap='Blues')
axs.set_title("training kernel matrix")
fig.colorbar(pos, ax=axs)
plt.show()

Figure 1 displays the estimated kernel matrix. This matrix has an interesting structure: it is block-diagonal. This is due to the fact that the second quadrant of the matrix correspond to the inner products of the population groups with gender=0. These share the common characteristic gender=0, and each row is similar to its neighbours due to the encoding: similar ages are also similar in the quantum feature space. Analogously, the fourth quadrant of the matrix contains the gender=1 population groups. The first and third quadrants contain the inner products between gender=0 and gender=1 population groups and so are dissimilar in the quantum feature space.

Next, we run a support vector regression (SVR) taking the kernel matrix, the weights, and the y-values as input.

f_matrix_train = matrix_train

zzpc_svr = SVR(kernel='precomputed')
zzpc_svr.fit(f_matrix_train, y, w)
zzpc_score = zzpc_svr.score(f_matrix_train, y, w)
print(f'score: {zzpc_score}')
pred = zzpc_svr.predict(f_matrix_train)

p_hat = 1/(1+np.exp(-pred))
p = 1/(1+np.exp(-y))
plt.figure()
plt.plot([a[1] for a in x if a[0] == 0], [p for a,p in zip(x,p) if a[0] == 0],'or')
plt.plot([a[1] for a in x if a[0] == 0], [p for a,p in zip(x,p_hat) if a[0] == 0],'r')
plt.plot([a[1] for a in x if a[0] == 1], [p for a,p in zip(x,p) if a[0] == 1],'ob')
plt.plot([a[1] for a in x if a[0] == 1], [p for a,p in zip(x,p_hat) if a[0] == 1],'b')
plt.legend(['Group 0', 'Group 0 model', 'Group 1', 'Group 1 model'])
plt.show()

score: 0.9315524163072193

The score you see is the R-squared score from the SVR model. This score represents the proportion of the variance in the dependent variable that is predictable from the independent variables1. In other words, it gives you a measure of how well your model is fitting the data. A score of 1.0 indicates a perfect fit.

We check the results against a classical SVR.

classical = SVR(kernel='rbf')
classical.fit(x, y, w)
classical_score = classical.score(x, y, w)
print(f'score: {classical_score}')
pred = classical.predict(x)

p_hat = 1/(1+np.exp(-pred))
p = 1/(1+np.exp(-y))
plt.figure()
plt.plot([a[1] for a in x if a[0] == 0], [p for a,p in zip(x,p) if a[0] == 0],'or')
plt.plot([a[1] for a in x if a[0] == 0], [p for a,p in zip(x,p_hat) if a[0] == 0],'r')
plt.plot([a[1] for a in x if a[0] == 1], [p for a,p in zip(x,p) if a[0] == 1],'ob')
plt.plot([a[1] for a in x if a[0] == 1], [p for a,p in zip(x,p_hat) if a[0] == 1],'b')
plt.legend(['Group 0', 'Group 0 model', 'Group 1', 'Group 1 model'])
plt.show()

score: 0.9390284403415132

The above analysis is done in-sample. Below, we use the leave-one-out cross-validation functionality from the sklearn module.

from sklearn.model_selection import LeaveOneOut
from sklearn.metrics import r2_score

zzpc_svr = SVR(kernel='precomputed')
classical = SVR(kernel='rbf')

cols = range(len(x))

loo = LeaveOneOut()
pred_loo = []
pred_loo_classic = []
used_K = matrix_train
for train_index, test_index in loo.split(cols):
    #print(x[train_index])
    zzpc_svr.fit(used_K[np.ix_(train_index,train_index)], y[train_index], w[train_index])
    pred = zzpc_svr.predict(used_K[np.ix_(test_index,train_index)])
    for p in pred:
        pred_loo.append(p)
    classical.fit(x[train_index], y[train_index], w[train_index])
    cpred = classical.predict(x[test_index])
    for p in cpred:
        pred_loo_classic.append(p)
pred_loo = np.array(pred_loo)
pred_loo_classic = np.array(pred_loo_classic)

from sklearn.model_selection import LeaveOneOut
from sklearn.metrics import r2_score

def predict_loo(K, y, w, zzpc_svr):
    cols = range(len(y))
    loo = LeaveOneOut()
    pred_loo = []
    for train_index, test_index in loo.split(cols):
        #print(x[train_index])
        zzpc_svr.fit(K[np.ix_(train_index,train_index)], y[train_index], w[train_index])
        pred = zzpc_svr.predict(K[np.ix_(test_index,train_index)])
        for p in pred:
            pred_loo.append(p)
    return np.array(pred_loo)

statevector_loo = predict_loo(matrix_train, y, w, SVR(kernel='precomputed'))

p_hat = 1/(1+np.exp(-statevector_loo))
p = 1/(1+np.exp(-y))
score = r2_score(p, p_hat, sample_weight=w)
print(f'score: {score}')
fig, ax = plt.subplots(1, 1)
plt.plot(np.array([a[1] for a in x if a[0] == 0])*100, [p for a,p in zip(x,p) if a[0] == 0],'or')
plt.plot(np.array([a[1] for a in x if a[0] == 1])*100, [p for a,p in zip(x,p) if a[0] == 1],'ob')
plt.plot(np.array([a[1] for a in x if a[0] == 0])*100, [p for a,p in zip(x,p_hat) if a[0] == 0],'r')
plt.plot(np.array([a[1] for a in x if a[0] == 1])*100, [p for a,p in zip(x,p_hat) if a[0] == 1],'b')
plt.legend(['Group 0', 'Group 1', 'Group 0, statevector', 'Group 1, statevector'])
ax.axes.yaxis.set_ticklabels([])
plt.show()

score: 0.917110541114462

Quantum Support Vector Regression (QSVR) Score (0.9155345528008413): This score suggests that your QSVR model is doing a good job of capturing the relationship between your features and your target variable.

Classical Support Vector Regression (SVR) Score (0.9295558682302361): This score is slightly higher than your QSVR score, suggesting that your classical SVR model might be fitting the data slightly better than your QSVR model.

Leave-One-Out Cross-Validation Score (0.8977158412990511): This score is lower than both your QSVR and classical SVR scores. This isn’t surprising, as leave-one-out cross-validation can often result in lower scores because it’s a more rigorous form of validation. It’s a good way to get an unbiased estimate of how well your model will perform on unseen data.

To compare these results, you can look at how much each score deviates from 1.0. The closer the score is to 1.0, the better the model is fitting your data. However, keep in mind that a high R-squared score doesn’t always mean your model will perform well on unseen data. It’s important to also look at other metrics and plots (like residuals) to assess your model’s performance.

Exercise 2: Residuals

In the context of regression analysis, residuals are the difference between the observed (actual) and predicted values for a data point. In other words, a residual is an error in prediction.

Mathematically, if you have an observed value $\(y_{\text{obs}}\)\( and a predicted value \)\(y_{\text{pred}}\)\( for a data point, the residual \)\(e\)$ for that data point is calculated as:

\[e = y_{\text{obs}} - y_{\text{pred}}\]

Residuals are used to understand the accuracy of a regression model. If a model fits the data well, the residuals should be small and randomly distributed around zero. Patterns or trends in the residuals can indicate that the model is not adequately capturing some aspect of the data.

In our case, we will plot a histogram of residuals from our QSVR model. This visualization can help us understand the distribution of errors made by the model. If the model is performing well, we would expect to see most residuals close to zero, indicating that the predicted and actual values are close. Large residuals indicate that the model’s predictions are far off from the actual values.

Hints

1. use residuals = y - pred

2. in plt.histuse bins=20 # bins parameter determines the number of bins (or bars) to divide your data into

3. use appropriate lables

4. Advanced – use p_hat-p for residuals

You’ve done a great job implementing the Quantum Support Vector Regression (QSVR) and comparing it with a classical SVR. You’ve also used leave-one-out cross-validation to evaluate the performance of your models, which is a robust method for estimating the generalization error.

In terms of the residuals, you’ve plotted two histograms:

1. The first histogram represents the residuals calculated as y - pred, where y is the observed value and pred is the predicted value from your QSVR model. This gives you an idea of how well your model is predicting the actual y values. If your model is accurate, most residuals should be close to zero, indicating that the predicted and actual values are close. Large residuals indicate that the model’s predictions are far off from the actual values.

2. The second histogram represents the residuals calculated as p_hat - p, where p_hat is the predicted probability and p is the actual probability. This gives you an idea of how well your model is predicting the actual probabilities. Again, if your model is accurate, most residuals should be close to zero.

To describe the results, you would typically look at the distribution of residuals in your histograms. If your model is performing well, you would expect to see a distribution centered around zero, indicating that most predictions are close to the actual values or probabilities. If you see a skewed distribution or a distribution with a long tail, it might indicate that your model is systematically overestimating or underestimating the values or probabilities.

In terms of what you’ve done here, you’ve implemented a QSVR for modeling transition probabilities in disability insurance, compared it with a classical SVR, and evaluated their performance using leave-one-out cross-validation and residual analysis. This is a comprehensive approach to developing and evaluating predictive models. Great work! 😊

We used the statevector simulator and now we’ll run it on Helmi.

provider = IQMProvider("https://qc.vtt.fi/cocos")
backend = provider.get_backend()

print(f'Native operations: {backend.operation_names}')
print(f'Number of qubits: {backend.num_qubits}')
print(f'Coupling map: {backend.coupling_map}')

Native operations: ['r', 'id', 'cz', 'measure']
Number of qubits: 5
Coupling map: [[0, 2], [2, 0], [1, 2], [2, 1], [2, 3], [3, 2], [2, 4], [4, 2]]

The feature_map function is a way to encode classical data into a quantum state. In this case, it’s creating a quantum circuit qc with 2 qubits. It then applies a rotation around the y-axis (ry) on each qubit. The rotation angle is determined by the input data x multiplied by pi. The ParameterVector is a way to create a list of parameters.

First, we create some dummy data, but smaller (fewer data points) than it was before.

nsamp = 10
x = np.array( [[0.0, x/nsamp] for x in range(nsamp)] + [[1.0, x/nsamp] for x in range(nsamp)] )
y = np.array([-4.5 + 0.2*a[0] + a[1]*(1+a[0]*0.2) + 0.1*np.random.normal() for a in x])

p = 1/(1+np.exp(-y))

n = [int(200 + 600*xj[1]/0.7) if xj[1] <= 0.7 else int(800 - 300*(xj[1]-0.7)/0.3) for xj in x]

d = [nj*pj for nj,pj in zip(n,p)]

# Create a DataFrame
df1 = pd.DataFrame({
    'gender': [a[0] for a in x],
    'age group (normalized)': [a[1] for a in x],
    'population size': n,
    'number of transitions': d
})

# Write the DataFrame to a CSV file
df1.to_csv('fdata1.csv', index=False)  # index=False argument is used to prevent pandas from writing row indices into the CSV file

data1 = df1
# data1 = pd.read_csv('fdata1.csv')

# Transform and inspect data
x1 = np.array( [[a,b] for a,b in zip(data1['gender'], data1['age group (normalized)'])] )
print(x1)

[[0.  0. ]
 [0.  0.1]
 [0.  0.2]
 [0.  0.3]
 [0.  0.4]
 [0.  0.5]
 [0.  0.6]
 [0.  0.7]
 [0.  0.8]
 [0.  0.9]
 [1.  0. ]
 [1.  0.1]
 [1.  0.2]
 [1.  0.3]
 [1.  0.4]
 [1.  0.5]
 [1.  0.6]
 [1.  0.7]
 [1.  0.8]
 [1.  0.9]]

# raw transition frequency per group, used for plots
p_raw = [a/b for a,b in zip(d,n)]

# logistic probability, used as input to SVR optimization
y = np.array([np.log( p / (1-p) ) for p in p_raw])

# calculate weights for SVR based on population size
w = np.array([1./nj for nj in n])

# normalize weights (equivalent to changing the hyperparameters of the SVR)
w = w/sum(w)*len(x)

fig, axs = plt.subplots(1, 1, figsize=(10, 5))
axs.plot([x1j[1] for x1j,p in zip(x1, p_raw) if x1j[0] == 0], [p for x1j,p in zip(x1, p_raw) if x1j[0] == 0], 'or')
axs.plot([x1j[1] for x1j,p in zip(x1, p_raw) if x1j[0] == 1], [p for x1j,p in zip(x1, p_raw) if x1j[0] == 1], 'ob')
axs.set_title("raw transition frequencies")
plt.show()

Exercise 3: Implement and Visualize a Quantum Feature Map

In this exercise, you will implement a quantum feature map using Qiskit and visualize the resulting quantum circuit. The feature map is a way of encoding classical data into quantum states. In this case, we will use a simple feature map that applies a rotation around the y-axis on each qubit.

Task 1: Implement the Feature Map

First, let’s implement the feature map. Here is the function you need to complete:

def feature_map(x1):
    p = ParameterVector('x1', length=2)
    qc = QuantumCircuit(2)

    # TODO: Apply a rotation around the y-axis on each qubit
    # Hint: Use the ry gate and the input data x1

    return qc

Task 2: Visualize the Circuit

After implementing the feature map, you visualize the resulting quantum circuit, as we did at the beginning of this lesson.

Good luck with your exercise! If you have any questions or need further clarification, feel free to ask. 😊

The following part of the code implements a list of all possible quantum circuits for the given data x1. For each pair of data points (x1[i], x1[j]), it creates two circuits using the feature_map function. The feature map is a way of encoding classical data into quantum states. It then appends the inverse of qc_j to qc_i using the compose function. This function combines two circuits, and with .inverse(), it applies the inverse of the second circuit. After composing the circuits, it adds a measurement on all qubits using measure_all(). The circuit is then transpiled for a specific backend using the transpile function.

# from qiskit import transpile

# Initialize an empty list to store all circuits
all_circuits = []

for i in range(len(x1)):
    for j in range(i):
        # Substitute parameters with values from your data
        qc_i = feature_map(x1[i])
        qc_j = feature_map(x1[j])
        
        # Append the inverse of qc_j to qc_i and measure
        qc_ij = qc_i.compose(qc_j.inverse())
        qc_ij.measure_all()
        
        # Transpile the circuit for the backend
        qc_ij_transpiled = transpile(qc_ij, backend=backend)
        
        # Add the transpiled circuit to the list
        all_circuits.append(qc_ij_transpiled)

print(len(all_circuits))

190

This block of code defines a batch size and splits the transpiled circuits into batches. It then initializes a list to store job IDs. Each batch of circuits is executed on the backend with 1024 shots per circuit execution. Also, it can take some time to execute.

# Define the batch size
batch_size = 20

# Split the transpiled circuits into batches
batches = [all_circuits[i:i + batch_size] for i in range(0, len(all_circuits), batch_size)]

# Initialize a list to store job IDs
job_ids = []

# Execute each batch of circuits
for batch in batches:
    job = backend.run(batch, shots=1024)
    job_ids.append(job.job_id())

# This block of code executes each batch until it's successful or an exception occurs
# If an exception occurs, it waits for a random amount of time between 1 and 10 seconds before retrying

import time
import random

# Initialize a list to store results
results = []
n_shots = 1024

# Initialize a counter for the job number
job_number = 1

for batch in batches:
    success = False
    while not success:
        try:
            # Execute the batch of circuits
            job = backend.run(batch, shots=n_shots)
            
            # Print job execution status
            print(f"Tracking execution of job{job_number}:")
            job_monitor(job)
            
            # Wait for the job to finish and get the result
            result = job.result()
            
            # If successful, append the result and break the loop
            results.append(result)
            success = True
            
        except Exception as e:
            print(f"Execution failed with error: {e}")
            
            # Sleep for a bit (random backoff) before retrying
            time.sleep(random.randint(1, 10))
    
    # Increment the job number for the next iteration
    job_number += 1

Tracking execution of job1:
Job Status: job has successfully run
Tracking execution of job2:
Job Status: job has successfully run
Tracking execution of job3:
Job Status: job has successfully run
Tracking execution of job4:
Job Status: job has successfully run
Tracking execution of job5:
Job Status: job has successfully run
Tracking execution of job6:
Job Status: job has successfully run
Tracking execution of job7:
Job Status: job has successfully run
Tracking execution of job8:
Job Status: job has successfully run
Tracking execution of job9:
Job Status: job has successfully run
Tracking execution of job10:
Job Status: job has successfully run

Job Retrieval and Frequency Calculation: The script starts by iterating over a list of job_ids. For each job_id, it retrieves the job result from the backend. It then calculates the frequency of the ‘00’ state in the results and appends this frequency to the estimates list.

# Initialize a counter for the job number
job_number = 1

for job_id in job_ids:
    print(f"Job{job_number}: {job_id}")
    job_number += 1

Job1: ea63caca-9859-4988-afcc-4617a4422261
Job2: 4ea8454c-0146-47e3-a4f3-22e7098bf8b5
Job3: 29c17939-328b-46bd-bcf6-e2df758b4efd
Job4: 80e39182-931b-4bc9-83d7-f100a877d81b
Job5: bae16e48-b12a-4cea-bdf4-287bce4e11ab
Job6: ed7f15be-36d1-41c4-b0e1-794746789a48
Job7: b07bd0df-a786-4389-8def-ea0926f2fd17
Job8: d346a014-4afc-43ee-8ee4-d16e229e7c1e
Job9: baf11441-c6d0-4b0b-8641-068a929c72c4
Job10: 1895b727-95b8-4925-bc95-e7e2347493f7

estimates = []

for job_id in job_ids:
    job_res = backend.retrieve_job(job_id).result()
    for res in job_res.results:
        if '00' in res.data.counts:
            freq = res.data.counts['00']/n_shots
        else:
            freq = 0.0
        estimates.append(freq)

Kernel Matrix Construction: The script constructs a kernel matrix K_helmi using the estimates. The kernel matrix is a symmetric matrix where each element represents a measure of similarity between two data points.

K_helmi = np.identity(len(x1))

cnt = 0
for i in range(len(x1)):
    for j in range(i):
        #print(i, ' ', j)
        K_helmi[i,j] = estimates[cnt]
        K_helmi[j,i] = estimates[cnt]
        cnt += 1

Kernel Matrix Visualization: The kernel matrix is then visualized using matplotlib’s imshow function. The color map ‘Blues’ is used, and a color bar is added for reference.

fig, axs = plt.subplots(1, 1, figsize=(10, 5))
pos = axs.imshow(np.asmatrix(K_helmi),
              interpolation='nearest', origin='upper', cmap='Blues')
axs.set_title("kernel matrix from Helmi")
fig.colorbar(pos, ax=axs)
plt.show()

Support Vector Regression (SVR): The script uses the Support Vector Regression (SVR) model from scikit-learn with a precomputed kernel. It fits the model to the data (f_matrix_train and y) and calculates the score of the model using the score method.

Prediction and Plotting: The fitted model is used to predict outcomes, which are then plotted against the actual outcomes for two groups (Group 0 and Group 1).

f_matrix_train = K_helmi

zzpc_svr = SVR(kernel='precomputed')
zzpc_svr.fit(f_matrix_train, y, w)
zzpc_score = zzpc_svr.score(f_matrix_train, y, w)
print(f'score: {zzpc_score}')
pred = zzpc_svr.predict(f_matrix_train)

p_hat = 1/(1+np.exp(-pred))
p = 1/(1+np.exp(-y))
plt.figure()
plt.plot([a[1] for a in x if a[0] == 0], [p for a,p in zip(x,p) if a[0] == 0],'or')
plt.plot([a[1] for a in x if a[0] == 0], [p for a,p in zip(x,p_hat) if a[0] == 0],'r')
plt.plot([a[1] for a in x if a[0] == 1], [p for a,p in zip(x,p) if a[0] == 1],'ob')
plt.plot([a[1] for a in x if a[0] == 1], [p for a,p in zip(x,p_hat) if a[0] == 1],'b')
plt.legend(['Group 0', 'Group 0 model', 'Group 1', 'Group 1 model'])
plt.show()

score: 0.9001754733534563

Error Calculation and Visualization: Finally, the script calculates the absolute errors between K_helmi and matrix_train, visualizes these errors using imshow, and displays this with a color bar for reference.

# redo the steps for statevector kernel to compare with the Helmi results

zz_kernel = QuantumKernel(feature_map=qc, quantum_instance=Aer.get_backend('statevector_simulator'))
zz_circuit = zz_kernel.construct_circuit(ParameterVector('x_i', length=2),ParameterVector('x_j', length=2))
zz_circuit.decompose().draw(output='mpl')
matrix_train = zz_kernel.evaluate(x_vec=x)

errors = np.abs(K_helmi - matrix_train)

fig, axs = plt.subplots(1, 1, figsize=(10, 5))
pos = axs.imshow(np.asmatrix(errors),
              interpolation='nearest', origin='upper', cmap='Reds')
axs.set_title("erorrs")
fig.colorbar(pos, ax=axs)
plt.show()

# solution ex2
residuals = y - pred
plt.hist(residuals, bins=20)
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.title('Histogram of Residuals')
plt.show()

# solution1 ex2
residuals1 = p_hat-p
plt.hist(residuals1, bins=20)
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.title('Histogram of Residuals')
plt.show()

# solution ex3

def feature_map(x1):
    p = ParameterVector('x1', length=2)
    qc = QuantumCircuit(2)

    qc.ry(x1[0]*np.pi, 0)
    qc.ry(x1[1]*np.pi, 1)

    return qc

# solution1 ex3
# Draw the circuit
qc.draw('mpl')

TEST CODE