Saltar al contenido principal

Mejore la clasificación de características utilizando núcleos cuánticos proyectados

Estimación de uso: 80 minutos en un procesador Heron r3 (NOTA: Esto es solo una estimación. Su tiempo de ejecución puede variar.)

En este tutorial, demostramos cómo ejecutar un núcleo cuántico proyectado (PQK, por sus siglas en inglés) con Qiskit en un conjunto de datos biológicos del mundo real, basado en el artículo Enhanced Prediction of CAR T-Cell Cytotoxicity with Quantum-Kernel Methods [1].

PQK es un método utilizado en el aprendizaje automático cuántico (QML, por sus siglas en inglés) para codificar datos clásicos en un espacio de características cuántico y proyectarlos de vuelta al dominio clásico, utilizando computadoras cuánticas para mejorar la selección de características. Implica codificar datos clásicos en estados cuánticos mediante un circuito cuántico, típicamente a través de un proceso llamado mapeo de características, donde los datos se transforman en un espacio de Hilbert de alta dimensión. El aspecto "proyectado" se refiere a la extracción de información clásica de los estados cuánticos, mediante la medición de observables específicos, para construir una matriz de núcleo que pueda ser utilizada en algoritmos clásicos basados en núcleos, como las máquinas de vectores de soporte. Este enfoque aprovecha las ventajas computacionales de los sistemas cuánticos para lograr potencialmente un mejor rendimiento en ciertas tareas en comparación con los métodos clásicos.

Este tutorial también asume familiaridad general con los métodos de QML. Para una exploración más profunda del QML, consulta el curso Quantum machine learning en IBM Quantum Learning.

Requisitos

Antes de comenzar este tutorial, asegúrate de tener instalado lo siguiente:

  • Qiskit SDK v2.0 o posterior, con soporte de visualización
  • Qiskit Runtime v0.40 o posterior (pip install qiskit-ibm-runtime)
  • Category encoders 2.8.1 (pip install category-encoders)
  • NumPy 2.3.2 (pip install numpy)
  • Pandas 2.3.2 (pip install pandas)
  • Scikit-learn 1.7.1 (pip install scikit-learn)
  • Tqdm 4.67.1 (pip install tqdm)

Configuración

# Added by doQumentation — required packages for this notebook
!pip install -q category-encoders numpy pandas qiskit qiskit-ibm-runtime scipy scikit-learn tqdm
import warnings

# Standard libraries
import os
import numpy as np
import pandas as pd

# Machine learning and data processing
import category_encoders as ce
from scipy.linalg import inv, sqrtm
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.model_selection import GridSearchCV, StratifiedKFold
from sklearn.svm import SVC

# Qiskit and Qiskit Runtime
from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.circuit.library import UnitaryGate, ZZFeatureMap
from qiskit.quantum_info import SparsePauliOp, random_unitary
from qiskit.transpiler import generate_preset_pass_manager
from qiskit_ibm_runtime import (
Batch,
EstimatorOptions,
EstimatorV2 as Estimator,
QiskitRuntimeService,
)

# Progress bar
import tqdm

warnings.filterwarnings("ignore")

Paso 1: Mapear las entradas clásicas a un problema cuántico

Preparación del conjunto de datos

En este tutorial utilizamos un conjunto de datos biológicos del mundo real para una tarea de clasificación binaria, generado por Daniels et al. (2022) y disponible para descarga desde el material complementario incluido con el artículo. Los datos consisten en células T CAR, que son células T modificadas genéticamente utilizadas en inmunoterapia para tratar ciertos tipos de cáncer. Las células T, un tipo de célula inmunitaria, se modifican en un laboratorio para expresar receptores de antígenos quiméricos (CAR) que se dirigen a proteínas específicas en las células cancerosas. Estas células T modificadas pueden reconocer y destruir las células cancerosas de manera más eficaz. Las características de los datos son los motivos de las células T CAR, que se refieren al componente estructural o funcional específico del CAR diseñado en las células T. Basándose en estos motivos, nuestra tarea es predecir la citotoxicidad de una célula T CAR dada, etiquetándola como tóxica o no tóxica. A continuación se muestran las funciones auxiliares para preprocesar este conjunto de datos.

def preprocess_data(dir_root, args):
"""
Preprocess the training and test data.
"""
# Read from the csv files
train_data = pd.read_csv(
os.path.join(dir_root, args["file_train_data"]),
encoding="unicode_escape",
sep=",",
)
test_data = pd.read_csv(
os.path.join(dir_root, args["file_test_data"]),
encoding="unicode_escape",
sep=",",
)

# Fix the last motif ID
train_data[train_data == 17] = 14
train_data.columns = [
"Cell Number",
"motif",
"motif.1",
"motif.2",
"motif.3",
"motif.4",
"Nalm 6 Cytotoxicity",
]
test_data[test_data == 17] = 14
test_data.columns = [
"Cell Number",
"motif",
"motif.1",
"motif.2",
"motif.3",
"motif.4",
"Nalm 6 Cytotoxicity",
]

# Adjust motif at the third position
if args["filter_for_spacer_motif_third_position"]:
train_data = train_data[
(train_data["motif.2"] == 14) | (train_data["motif.2"] == 0)
]
test_data = test_data[
(test_data["motif.2"] == 14) | (test_data["motif.2"] == 0)
]

train_data = train_data[
args["motifs_to_use"] + [args["label_name"], "Cell Number"]
]
test_data = test_data[
args["motifs_to_use"] + [args["label_name"], "Cell Number"]
]

# Adjust motif at the last position
if not args["allow_spacer_motif_last_position"]:
last_motif = args["motifs_to_use"][len(args["motifs_to_use"]) - 1]
train_data = train_data[
(train_data[last_motif] != 14) & (train_data[last_motif] != 0)
]
test_data = test_data[
(test_data[last_motif] != 14) & (test_data[last_motif] != 0)
]

# Get the labels
train_labels = np.array(train_data[args["label_name"]])
test_labels = np.array(test_data[args["label_name"]])

# For the classification task use the threshold to binarize labels
train_labels[train_labels > args["label_binarization_threshold"]] = 1
train_labels[train_labels < 1] = args["min_label_value"]
test_labels[test_labels > args["label_binarization_threshold"]] = 1
test_labels[test_labels < 1] = args["min_label_value"]

# Reduce data to just the motifs of interest
train_data = train_data[args["motifs_to_use"]]
test_data = test_data[args["motifs_to_use"]]

# Get the class and motif counts
min_class = np.min(np.unique(np.concatenate([train_data, test_data])))
max_class = np.max(np.unique(np.concatenate([train_data, test_data])))

num_class = max_class - min_class + 1
num_motifs = len(args["motifs_to_use"])
print(str(max_class) + ":" + str(min_class) + ":" + str(num_class))

train_data = train_data - min_class
test_data = test_data - min_class

return (
train_data,
test_data,
train_labels,
test_labels,
num_class,
num_motifs,
)

def data_encoder(args, train_data, test_data, num_class, num_motifs):
"""
Use one-hot or binary encoding for classical data representation.
"""
if args["encoder"] == "one-hot":
# Transform to one-hot encoding
train_data = np.eye(num_class)[train_data]
test_data = np.eye(num_class)[test_data]

train_data = train_data.reshape(
train_data.shape[0], train_data.shape[1] * train_data.shape[2]
)
test_data = test_data.reshape(
test_data.shape[0], test_data.shape[1] * test_data.shape[2]
)

elif args["encoder"] == "binary":
# Transform to binary encoding
encoder = ce.BinaryEncoder()

base_array = np.unique(np.concatenate([train_data, test_data]))
base = pd.DataFrame(base_array).astype("category")
base.columns = ["motif"]
for motif_name in args["motifs_to_use"][1:]:
base[motif_name] = base.loc[:, "motif"]
encoder.fit(base)

train_data = encoder.transform(train_data.astype("category"))
test_data = encoder.transform(test_data.astype("category"))

train_data = np.reshape(
train_data.values, (train_data.shape[0], num_motifs, -1)
)
test_data = np.reshape(
test_data.values, (test_data.shape[0], num_motifs, -1)
)

train_data = train_data.reshape(
train_data.shape[0], train_data.shape[1] * train_data.shape[2]
)
test_data = test_data.reshape(
test_data.shape[0], test_data.shape[1] * test_data.shape[2]
)

else:
raise ValueError("Invalid encoding type.")

return train_data, test_data

Puedes ejecutar este tutorial ejecutando la siguiente celda, que crea automáticamente la estructura de carpetas necesaria y descarga tanto los archivos de entrenamiento como los de prueba directamente en tu entorno. Si ya tienes estos archivos localmente, este paso los sobrescribirá de forma segura para garantizar la consistencia de la versión.

## Download dataset

# Create data directory if it doesn't exist
!mkdir -p data_tutorial/pqk

# Download the training and test sets from the official Qiskit documentation repo
!wget -q --show-progress -O data_tutorial/pqk/train_data.csv \
https://raw.githubusercontent.com/Qiskit/documentation/main/datasets/tutorials/pqk/train_data.csv

!wget -q --show-progress -O data_tutorial/pqk/test_data.csv \
https://raw.githubusercontent.com/Qiskit/documentation/main/datasets/tutorials/pqk/test_data.csv

!wget -q --show-progress -O data_tutorial/pqk/projections_train.csv \
https://raw.githubusercontent.com/Qiskit/documentation/main/datasets/tutorials/pqk/projections_train.csv

!wget -q --show-progress -O data_tutorial/pqk/projections_test.csv \
https://raw.githubusercontent.com/Qiskit/documentation/main/datasets/tutorials/pqk/projections_test.csv

# Check the files have been downloaded
!echo "Dataset files downloaded:"
!ls -lh data_tutorial/pqk/*.csv
args = {
"file_train_data": "train_data.csv",
"file_test_data": "test_data.csv",
"motifs_to_use": ["motif", "motif.1", "motif.2", "motif.3"],
"label_name": "Nalm 6 Cytotoxicity",
"label_binarization_threshold": 0.62,
"filter_for_spacer_motif_third_position": False,
"allow_spacer_motif_last_position": True,
"min_label_value": -1,
"encoder": "one-hot",
}
dir_root = "./"

# Preprocess data
train_data, test_data, train_labels, test_labels, num_class, num_motifs = (
preprocess_data(dir_root=dir_root, args=args)
)

# Encode the data
train_data, test_data = data_encoder(
args, train_data, test_data, num_class, num_motifs
)
14:0:15

También transformamos el conjunto de datos de manera que 11 se represente como π/2\pi/2 con fines de escalado.

# Change 1 to pi/2
angle = np.pi / 2

tmp = pd.DataFrame(train_data).astype("float64")
tmp[tmp == 1] = angle
train_data = tmp.values

tmp = pd.DataFrame(test_data).astype("float64")
tmp[tmp == 1] = angle
test_data = tmp.values

Verificamos los tamaños y las formas de los conjuntos de datos de entrenamiento y prueba.

print(train_data.shape, train_labels.shape)
print(test_data.shape, test_labels.shape)
(172, 60) (172,)
(74, 60) (74,)

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

Circuito cuántico

Ahora construimos el mapa de características que incorpora nuestro conjunto de datos clásico en un espacio de características de mayor dimensión. Para esta incorporación, utilizamos el ZZFeatureMap de Qiskit.

feature_dimension = train_data.shape[1]
reps = 24
insert_barriers = True
entanglement = "pairwise"

# ZZFeatureMap with linear entanglement and a repetition of 2
embed = ZZFeatureMap(
feature_dimension=feature_dimension,
reps=reps,
entanglement=entanglement,
insert_barriers=insert_barriers,
name="ZZFeatureMap",
)
embed.decompose().draw(output="mpl", style="iqp", fold=-1)

Output of the previous code cell

Otra opción de incorporación cuántica es el ansatz de evolución del Hamiltoniano de Heisenberg unidimensional. Puede omitir la ejecución de esta sección si desea continuar con el ZZFeatureMap.

feature_dimension = train_data.shape[1]
num_qubits = feature_dimension + 1
embed2 = QuantumCircuit(num_qubits)
num_trotter_steps = 6
pv_length = feature_dimension * num_trotter_steps
pv = ParameterVector("theta", pv_length)

# Add Haar random single qubit unitary to each qubit as initial state
np.random.seed(42)
seeds_unitary = np.random.randint(0, 100, num_qubits)
for i in range(num_qubits):
rand_gate = UnitaryGate(random_unitary(2, seed=seeds_unitary[i]))
embed2.append(rand_gate, [i])

def trotter_circ(feature_dimension, num_trotter_steps):
num_qubits = feature_dimension + 1
circ = QuantumCircuit(num_qubits)
# Even
for i in range(0, feature_dimension, 2):
circ.rzz(2 * pv[i] / num_trotter_steps, i, i + 1)
for i in range(0, feature_dimension, 2):
circ.rxx(2 * pv[i] / num_trotter_steps, i, i + 1)
for i in range(0, feature_dimension, 2):
circ.ryy(2 * pv[i] / num_trotter_steps, i, i + 1)
# Odd
for i in range(1, feature_dimension, 2):
circ.rzz(2 * pv[i] / num_trotter_steps, i, i + 1)
for i in range(1, feature_dimension, 2):
circ.rxx(2 * pv[i] / num_trotter_steps, i, i + 1)
for i in range(1, feature_dimension, 2):
circ.ryy(2 * pv[i] / num_trotter_steps, i, i + 1)
return circ

# Hamiltonian evolution ansatz
for step in range(num_trotter_steps):
circ = trotter_circ(feature_dimension, num_trotter_steps)
if step % 2 == 0:
embed2 = embed2.compose(circ)
else:
reverse_circ = circ.reverse_ops()
embed2 = embed2.compose(reverse_circ)

embed2.draw(output="mpl", style="iqp", fold=-1)

Output of the previous code cell

Paso 3: Ejecutar utilizando las primitivas de Qiskit

Medir las 1-RDM

Los principales bloques de construcción de los núcleos cuánticos proyectados son las matrices de densidad reducidas (RDM), que se obtienen mediante mediciones proyectivas del mapa de características cuántico. En este paso, obtenemos todas las matrices de densidad reducidas de un solo cúbit (1-RDM), que posteriormente se proporcionarán a la función de núcleo exponencial clásica. Veamos cómo calcular la 1-RDM dado un único punto de datos del conjunto de datos antes de ejecutar sobre todos los datos. Las 1-RDM son una colección de mediciones de un solo cúbit de los operadores de Pauli X, Y y Z en todos los cúbits. Esto se debe a que una RDM de un solo cúbit puede expresarse completamente como: ρ=12(I+σxσx+σyσy+σzσz)\rho = \frac{1}{2} \big( I + \braket \sigma_x \sigma_x + \braket \sigma_y \sigma_y + \braket \sigma_z \sigma_z \big) Primero seleccionamos el backend a utilizar.

service = QiskitRuntimeService()
backend = service.least_busy(
operational=True, simulator=False, min_num_qubits=133
)
target = backend.target

Luego ejecutamos el circuito cuántico y medimos las proyecciones. Ten en cuenta que activamos la mitigación de errores, incluyendo la Extrapolación de Ruido Cero (ZNE, por sus siglas en inglés).

# Let's select the ZZFeatureMap embedding for this example
qc = embed
num_qubits = feature_dimension

# Identity operator on all qubits
id = "I" * num_qubits

# Let's select the first training datapoint as an example
parameters = train_data[0]

# Bind parameter to the circuit and simplify it
qc_bound = qc.assign_parameters(parameters)
transpiler = generate_preset_pass_manager(
optimization_level=3, basis_gates=["u3", "cz"]
)
transpiled_circuit = transpiler.run(qc_bound)

# Transpile for hardware
transpiler = generate_preset_pass_manager(optimization_level=3, target=target)
transpiled_circuit = transpiler.run(transpiled_circuit)

# We group all commuting observables
# These groups are the Pauli X, Y and Z operators on individual qubits
observables_x = [
SparsePauliOp(id[:i] + "X" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]
observables_y = [
SparsePauliOp(id[:i] + "Y" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]
observables_z = [
SparsePauliOp(id[:i] + "Z" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]

# We define the primitive unified blocs (PUBs) consisting of the embedding circuit,
# set of observables and the circuit parameters
pub_x = (transpiled_circuit, observables_x)
pub_y = (transpiled_circuit, observables_y)
pub_z = (transpiled_circuit, observables_z)

# Experiment options for error mitigation
num_randomizations = 300
shots_per_randomization = 100
noise_factors = [1, 3, 5]

experimental_opts = {}
experimental_opts["resilience"] = {
"measure_mitigation": True,
"zne_mitigation": True,
"zne": {
"noise_factors": noise_factors,
"amplifier": "gate_folding",
"extrapolated_noise_factors": [0] + noise_factors,
},
}
experimental_opts["twirling"] = {
"num_randomizations": num_randomizations,
"shots_per_randomization": shots_per_randomization,
"strategy": "active-accum",
}

# We define and run the estimator to obtain <X>, <Y> and <Z> on all qubits
estimator = Estimator(mode=backend, options=experimental_opts)

job = estimator.run([pub_x, pub_y, pub_z])

A continuación recuperamos los resultados.

job_result_x = job.result()[0].data.evs
job_result_y = job.result()[1].data.evs
job_result_z = job.result()[2].data.evs
print(job_result_x)
print(job_result_y)
print(job_result_z)
[ 3.67865951e-03  1.01158571e-02 -3.95790878e-02  6.33984326e-03
1.86035759e-02 -2.91533268e-02 -1.06374793e-01 4.48873518e-18
4.70201764e-02 3.53997968e-02 2.53130819e-02 3.23903401e-02
6.06327843e-03 1.16313667e-02 -1.12387504e-02 -3.18457725e-02
-4.16445718e-04 -1.45609602e-03 -4.21737114e-01 2.83705669e-02
6.91332890e-03 -7.45363001e-02 -1.20139326e-02 -8.85566135e-02
-3.22648394e-02 -3.24228074e-02 6.20431299e-04 3.04225434e-03
5.72795792e-03 1.11288428e-02 1.50395861e-01 9.18380197e-02
1.02553163e-01 2.98312847e-02 -3.30298912e-01 -1.13979648e-01
4.49159340e-03 8.63861493e-02 3.05666566e-02 2.21463145e-04
1.45946735e-02 8.54537275e-03 -8.09805979e-02 -2.92608104e-02
-3.91243644e-02 -3.96632760e-02 -1.41187613e-01 -1.07363243e-01
1.81089440e-02 2.70778895e-02 1.45139414e-02 2.99480458e-02
4.99137134e-02 7.08789852e-02 4.30565759e-02 8.71287156e-02
1.04334798e-01 7.72191962e-02 7.10059720e-02 1.04650403e-01]
[-7.31765102e-05 7.42669174e-03 9.82277344e-03 5.92638249e-02
4.24120486e-02 -9.06473416e-03 4.55057675e-03 8.43494094e-03
6.92097339e-02 -6.82234424e-02 6.13509008e-02 3.94200491e-02
-1.24037979e-02 1.01976642e-01 7.90538600e-03 -7.19726160e-02
-1.19501703e-16 -1.03796614e-02 7.37382463e-02 1.97238568e-01
-3.59250635e-02 -2.67554009e-02 3.55010633e-02 7.68877990e-02
6.50677589e-05 -6.59298767e-03 -1.23719487e-02 -6.41938151e-02
1.95603072e-02 -2.48448551e-02 5.17784810e-02 -5.93767100e-02
3.11897681e-02 -3.91959720e-18 -4.47769148e-03 1.39202197e-01
-6.56387523e-02 -5.85665483e-02 9.52905894e-03 -8.61460731e-02
3.91790656e-02 -1.27544375e-01 1.63712244e-01 3.36816934e-04
2.26230028e-02 -2.45023393e-05 4.95635588e-03 1.44779564e-01
3.71625177e-02 3.65675948e-03 2.83694017e-02 -7.10500602e-02
-1.15467702e-01 6.21712129e-03 -4.80958959e-02 2.21021066e-02
7.99062499e-02 -1.87164076e-02 -3.67100369e-02 -2.38923731e-02]
[ 6.85871605e-01 5.07725024e-01 8.71024642e-03 3.34823455e-02
4.58684961e-02 9.44384189e-17 -4.46829296e-02 -2.91296778e-02
4.15466461e-02 2.89628330e-02 1.88624017e-03 5.37110446e-02
2.59579053e-03 1.39327071e-02 -2.90781778e-02 5.07209866e-03
5.83403000e-02 2.60764440e-02 4.45999706e-17 -6.66701417e-03
3.03215873e-01 2.26172533e-02 2.43105960e-02 4.98861041e-18
-2.45530791e-02 6.26940708e-02 1.21058073e-02 2.76675948e-04
2.63980996e-02 2.58302364e-02 7.47856723e-02 8.42728943e-02
5.70989097e-02 6.92955086e-02 -5.68313712e-03 1.32199452e-01
8.90511238e-02 -3.45204621e-02 -1.05445836e-01 6.03864150e-03
2.16291384e-02 8.22303162e-03 1.00856715e-02 6.28973151e-02
6.26727169e-02 6.15399206e-02 9.67320897e-02 1.03045269e-16
1.79688783e-01 -1.59960520e-02 -1.15422952e-02 9.60200470e-03
6.58396672e-02 7.78329830e-03 6.53226955e-02 2.45778685e-03
4.36694753e-03 5.75098762e-03 -2.48896201e-02 8.33740755e-05]

Imprimimos el tamaño del circuito y la profundidad de las compuertas de dos cúbits.

print(f"qubits: {qc.num_qubits}")
print(
f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}"
)
print(
f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}"
)
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw("mpl", fold=-1, style="clifford", idle_wires=False)
qubits: 60
2q-depth: 64
2q-size: 1888
Operator counts: OrderedDict({'rz': 6016, 'sx': 4576, 'cz': 1888, 'x': 896, 'barrier': 31})

Output of the previous code cell

Ahora podemos iterar sobre todo el conjunto de datos de entrenamiento para obtener todas las 1-RDM. También proporcionamos los resultados de un experimento que ejecutamos en hardware cuántico. Puede ejecutar el entrenamiento tú mismo estableciendo la bandera a continuación en True, o utilizar los resultados de proyección que proporcionamos.

# Set this to True if you want to run the training on hardware
run_experiment = False
# Identity operator on all qubits
id = "I" * num_qubits

# projections_train[i][j][k] will be the expectation value of the j-th Pauli operator (0: X, 1: Y, 2: Z)
# of datapoint i on qubit k
projections_train = []
jobs_train = []

# Experiment options for error mitigation
num_randomizations = 300
shots_per_randomization = 100
noise_factors = [1, 3, 5]

experimental_opts = {}
experimental_opts["resilience"] = {
"measure_mitigation": True,
"zne_mitigation": True,
"zne": {
"noise_factors": noise_factors,
"amplifier": "gate_folding",
"return_all_extrapolated": True,
"return_unextrapolated": True,
"extrapolated_noise_factors": [0] + noise_factors,
},
}
experimental_opts["twirling"] = {
"num_randomizations": num_randomizations,
"shots_per_randomization": shots_per_randomization,
"strategy": "active-accum",
}
options = EstimatorOptions(experimental=experimental_opts)

if run_experiment:
with Batch(backend=backend):
for i in tqdm.tqdm(
range(len(train_data)), desc="Training data progress"
):
# Get training sample
parameters = train_data[i]

# Bind parameter to the circuit and simplify it
qc_bound = qc.assign_parameters(parameters)
transpiler = generate_preset_pass_manager(
optimization_level=3, basis_gates=["u3", "cz"]
)
transpiled_circuit = transpiler.run(qc_bound)

# Transpile for hardware
transpiler = generate_preset_pass_manager(
optimization_level=3, target=target
)
transpiled_circuit = transpiler.run(transpiled_circuit)

# We group all commuting observables
# These groups are the Pauli X, Y and Z operators on individual qubits
observables_x = [
SparsePauliOp(id[:i] + "X" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]
observables_y = [
SparsePauliOp(id[:i] + "Y" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]
observables_z = [
SparsePauliOp(id[:i] + "Z" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]

# We define the primitive unified blocs (PUBs) consisting of the embedding circuit,
# set of observables and the circuit parameters
pub_x = (transpiled_circuit, observables_x)
pub_y = (transpiled_circuit, observables_y)
pub_z = (transpiled_circuit, observables_z)

# We define and run the estimator to obtain <X>, <Y> and <Z> on all qubits
estimator = Estimator(options=options)

job = estimator.run([pub_x, pub_y, pub_z])
jobs_train.append(job)
Training data progress: 100%|██████████| 172/172 [13:03<00:00,  4.55s/it]

Una vez que los trabajos están completos, podemos recuperar los resultados.

if run_experiment:
for i in tqdm.tqdm(
range(len(train_data)), desc="Retrieving training data results"
):
# Completed job
job = jobs_train[i]

# Job results
job_result_x = job.result()[0].data.evs
job_result_y = job.result()[1].data.evs
job_result_z = job.result()[2].data.evs

# Record <X>, <Y> and <Z> on all qubits for the current datapoint
projections_train.append([job_result_x, job_result_y, job_result_z])

Repetimos esto para el conjunto de prueba.

# Identity operator on all qubits
id = "I" * num_qubits

# projections_test[i][j][k] will be the expectation value of the j-th Pauli operator (0: X, 1: Y, 2: Z)
# of datapoint i on qubit k
projections_test = []
jobs_test = []

# Experiment options for error mitigation
num_randomizations = 300
shots_per_randomization = 100
noise_factors = [1, 3, 5]

experimental_opts = {}
experimental_opts["resilience"] = {
"measure_mitigation": True,
"zne_mitigation": True,
"zne": {
"noise_factors": noise_factors,
"amplifier": "gate_folding",
"return_all_extrapolated": True,
"return_unextrapolated": True,
"extrapolated_noise_factors": [0] + noise_factors,
},
}
experimental_opts["twirling"] = {
"num_randomizations": num_randomizations,
"shots_per_randomization": shots_per_randomization,
"strategy": "active-accum",
}
options = EstimatorOptions(experimental=experimental_opts)

if run_experiment:
with Batch(backend=backend):
for i in tqdm.tqdm(range(len(test_data)), desc="Test data progress"):
# Get test sample
parameters = test_data[i]

# Bind parameter to the circuit and simplify it
qc_bound = qc.assign_parameters(parameters)
transpiler = generate_preset_pass_manager(
optimization_level=3, basis_gates=["u3", "cz"]
)
transpiled_circuit = transpiler.run(qc_bound)

# Transpile for hardware
transpiler = generate_preset_pass_manager(
optimization_level=3, target=target
)
transpiled_circuit = transpiler.run(transpiled_circuit)

# We group all commuting observables
# These groups are the Pauli X, Y and Z operators on individual qubits
observables_x = [
SparsePauliOp(id[:i] + "X" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]
observables_y = [
SparsePauliOp(id[:i] + "Y" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]
observables_z = [
SparsePauliOp(id[:i] + "Z" + id[(i + 1) :]).apply_layout(
transpiled_circuit.layout
)
for i in range(num_qubits)
]

# We define the primitive unified blocs (PUBs) consisting of the embedding circuit,
# set of observables and the circuit parameters
pub_x = (transpiled_circuit, observables_x)
pub_y = (transpiled_circuit, observables_y)
pub_z = (transpiled_circuit, observables_z)

# We define and run the estimator to obtain <X>, <Y> and <Z> on all qubits
estimator = Estimator(options=options)

job = estimator.run([pub_x, pub_y, pub_z])
jobs_test.append(job)
Test data progress: 100%|██████████| 74/74 [00:13<00:00,  5.56it/s]

Podemos recuperar los resultados como antes.

if run_experiment:
for i in tqdm.tqdm(
range(len(test_data)), desc="Retrieving test data results"
):
# Completed job
job = jobs_test[i]

# Job results
job_result_x = job.result()[0].data.evs
job_result_y = job.result()[1].data.evs
job_result_z = job.result()[2].data.evs

# Record <X>, <Y> and <Z> on all qubits for the current datapoint
projections_test.append([job_result_x, job_result_y, job_result_z])

Paso 4: Posprocesar y devolver el resultado en el formato clásico deseado

Definir el núcleo cuántico proyectado

El núcleo cuántico proyectado se define con la siguiente función de núcleo: kPQ(xi,xj)=exp(γkP{X,Y,Z}(Tr[Pρk(xi)]Tr[Pρk(xj)])2)k^{\textrm{PQ}}(x_i, x_j) = \textrm{exp} \Big(-\gamma \sum_k \sum_{P \in \{ X,Y,Z \}} (\textrm{Tr}[P \rho_k(x_i)] - \textrm{Tr}[P \rho_k(x_j)])^2 \Big) En la ecuación anterior, γ>0\gamma>0 es un hiperparámetro ajustable. Los KijPQ=kPQ(xi,xj)K^{\textrm{PQ}}_{ij} = k^{\textrm{PQ}}(x_i, x_j) son las entradas de la matriz de núcleo KPQK^{\textrm{PQ}}. Utilizando la definición de las 1-RDM, podemos ver que los términos individuales dentro de la función de núcleo pueden evaluarse como Tr[Pρk(xi)]=P\textrm{Tr}[P \rho_k (x_i)] = \braket P, donde P{X,Y,Z}P \in \{ X,Y,Z \}. Estos valores esperados son precisamente lo que medimos anteriormente. Utilizando scikit-learn, podemos de hecho calcular el núcleo de manera aún más sencilla. Esto se debe a la función de base radial ('rbf') disponible: exp(γxx2) \textrm{exp} (-\gamma \lVert x - x' \rVert^2). Primero, simplemente necesitamos reorganizar los nuevos conjuntos de datos de entrenamiento y prueba proyectados en arreglos bidimensionales. Ten en cuenta que recorrer todo el conjunto de datos puede llevar aproximadamente 80 minutos en la QPU. Para asegurarnos de que el resto del tutorial sea fácilmente ejecutable, además proporcionamos proyecciones de un experimento previamente ejecutado (que se incluyen en los archivos que descargó en el bloque de código Download dataset). Si realizó el entrenamiento tú mismo, puede continuar el tutorial con sus propios resultados.

if run_experiment:
projections_train = np.array(projections_train).reshape(
len(projections_train), -1
)
projections_test = np.array(projections_test).reshape(
len(projections_test), -1
)
else:
projections_train = np.loadtxt("projections_train.txt")
projections_test = np.loadtxt("projections_test.txt")

Máquina de Vectores de Soporte (SVM)

Ahora podemos ejecutar una SVM clásica sobre este núcleo precalculado, y utilizar el núcleo entre los conjuntos de prueba y entrenamiento para la predicción.

# Range of 'C' and 'gamma' values as SVC hyperparameters
C_range = [0.001, 0.005, 0.007]
C_range.extend([x * 0.01 for x in range(1, 11)])
C_range.extend([x * 0.25 for x in range(1, 60)])
C_range.extend(
[
20,
50,
100,
200,
500,
700,
1000,
1100,
1200,
1300,
1400,
1500,
1700,
2000,
]
)

gamma_range = ["auto", "scale", 0.001, 0.005, 0.007]
gamma_range.extend([x * 0.01 for x in range(1, 11)])
gamma_range.extend([x * 0.25 for x in range(1, 60)])
gamma_range.extend([20, 50, 100])

param_grid = dict(C=C_range, gamma=gamma_range)

# Support vector classifier
svc = SVC(kernel="rbf")

# Define the cross validation
cv = StratifiedKFold(n_splits=10)

# Grid search for hyperparameter tuning (q: quantum)
grid_search_q = GridSearchCV(
svc, param_grid, cv=cv, verbose=1, n_jobs=-1, scoring="f1_weighted"
)
grid_search_q.fit(projections_train, train_labels)

# Best model with best parameters
best_svc_q = grid_search_q.best_estimator_
print(
f"The best parameters are {grid_search_q.best_params_} with a score of {grid_search_q.best_score_:.4f}"
)

# Test accuracy
accuracy_q = best_svc_q.score(projections_test, test_labels)
print(f"Test accuracy with best model: {accuracy_q:.4f}")
Fitting 10 folds for each of 6622 candidates, totalling 66220 fits
The best parameters are {'C': 8.5, 'gamma': 0.01} with a score of 0.6980
Test accuracy with best model: 0.8108

Evaluación comparativa clásica

Podemos ejecutar una SVM clásica con la función de base radial como núcleo sin realizar una proyección cuántica. Este resultado es nuestra referencia clásica.

# Support vector classifier
svc = SVC(kernel="rbf")

# Grid search for hyperparameter tuning (c: classical)
grid_search_c = GridSearchCV(
svc, param_grid, cv=cv, verbose=1, n_jobs=-1, scoring="f1_weighted"
)
grid_search_c.fit(train_data, train_labels)

# Best model with best parameters
best_svc_c = grid_search_c.best_estimator_
print(
f"The best parameters are {grid_search_c.best_params_} with a score of {grid_search_c.best_score_:.4f}"
)

# Test accuracy
accuracy_c = best_svc_c.score(test_data, test_labels)
print(f"Test accuracy with best model: {accuracy_c:.4f}")
Fitting 10 folds for each of 6622 candidates, totalling 66220 fits
The best parameters are {'C': 10.75, 'gamma': 0.04} with a score of 0.7830
Test accuracy with best model: 0.7432

Apéndice: Verificar el potencial de ventaja cuántica del conjunto de datos en tareas de aprendizaje

No todos los conjuntos de datos ofrecen una ventaja potencial con el uso de PQK. Existen algunos límites teóricos que se pueden utilizar como prueba preliminar para determinar si un conjunto de datos particular puede beneficiarse de los PQK. Para cuantificar esto, los autores de Power of data in quantum machine learning [2] definen cantidades denominadas complejidades de modelo clásico y cuántico, y separación geométrica de los modelos clásico y cuántico. Para esperar una ventaja cuántica potencial de los PQK, la separación geométrica entre los núcleos clásicos y cuánticos proyectados debe ser aproximadamente del orden de N\sqrt{N}, donde NN es el número de muestras de entrenamiento. Si se satisface esta condición, procedemos a verificar las complejidades de los modelos. Si la complejidad del modelo clásico es del orden de NN mientras que la complejidad del modelo cuántico proyectado es sustancialmente menor que NN, podemos esperar una ventaja potencial del PQK. La separación geométrica se define de la siguiente manera (F19 en [2]): gcq=g(KcKq)=KqKc(Kc+λI)2KcKqg_{cq} = g(K^c \Vert K^q) = \sqrt{\Vert \sqrt{K^q} \sqrt{K^c} (K^c + \lambda I)^{-2} \sqrt{K^c} \sqrt{K^q}\Vert_{\infty}}

# Gamma values used in best models above
gamma_c = grid_search_c.best_params_["gamma"]
gamma_q = grid_search_q.best_params_["gamma"]

# Regularization parameter used in the best classical model above
C_c = grid_search_c.best_params_["C"]
l_c = 1 / C_c

# Classical and quantum kernels used above
K_c = rbf_kernel(train_data, train_data, gamma=gamma_c)
K_q = rbf_kernel(projections_train, projections_train, gamma=gamma_q)

# Intermediate matrices in the equation
K_c_sqrt = sqrtm(K_c)
K_q_sqrt = sqrtm(K_q)
K_c_inv = inv(K_c + l_c * np.eye(K_c.shape[0]))
K_multiplication = (
K_q_sqrt @ K_c_sqrt @ K_c_inv @ K_c_inv @ K_c_sqrt @ K_q_sqrt
)

# Geometric separation
norm = np.linalg.norm(K_multiplication, ord=np.inf)
g_cq = np.sqrt(norm)
print(
f"Geometric separation between classical and quantum kernels is {g_cq:.4f}"
)

print(np.sqrt(len(train_data)))
Geometric separation between classical and quantum kernels is 1.5440
13.114877048604

La complejidad del modelo se define de la siguiente manera (M1 en [2]): sK,λ(N)=λ2i=1Nj=1N(K+λI)ij2yiyjN+i=1Nj=1N((K+λI)1K(K+λI)1)ijyiyjN s_{K, \lambda}(N) = \sqrt{\frac{\lambda^2 \sum_{i=1}^N \sum_{j=1}^N (K+\lambda I)^{-2}_{ij} y_i y_j}{N}} + \sqrt{\frac{\sum_{i=1}^N \sum_{j=1}^N ((K+\lambda I)^{-1}K(K+\lambda I)^{-1})_{ij} y_i y_j}{N}}

# Model complexity of the classical kernel

# Number of training data
N = len(train_data)

# Predicted labels
pred_labels = best_svc_c.predict(train_data)
pred_matrix = np.outer(pred_labels, pred_labels)

# Intermediate terms
K_c_inv = inv(K_c + l_c * np.eye(K_c.shape[0]))

# First term
first_sum = np.sum((K_c_inv @ K_c_inv) * pred_matrix)
first_term = l_c * np.sqrt(first_sum / N)

# Second term
second_sum = np.sum((K_c_inv @ K_c @ K_c_inv) * pred_matrix)
second_term = np.sqrt(second_sum / N)

# Model complexity
s_c = first_term + second_term
print(f"Classical model complexity is {s_c:.4f}")
Classical model complexity is 1.3578
# Model complexity of the projected quantum kernel

# Number of training data
N = len(projections_train)

# Predicted labels
pred_labels = best_svc_q.predict(projections_train)
pred_matrix = np.outer(pred_labels, pred_labels)

# Regularization parameter used in the best classical model above
C_q = grid_search_q.best_params_["C"]
l_q = 1 / C_q

# Intermediate terms
K_q_inv = inv(K_q + l_q * np.eye(K_q.shape[0]))

# First term
first_sum = np.sum((K_q_inv @ K_q_inv) * pred_matrix)
first_term = l_q * np.sqrt(first_sum / N)

# Second term
second_sum = np.sum((K_q_inv @ K_q @ K_q_inv) * pred_matrix)
second_term = np.sqrt(second_sum / N)

# Model complexity
s_q = first_term + second_term
print(f"Quantum model complexity is {s_q:.4f}")
Quantum model complexity is 1.5806

Referencias

  1. Utro, Filippo, et al. "Enhanced Prediction of CAR T-Cell Cytotoxicity with Quantum-Kernel Methods." arXiv preprint arXiv:2507.22710 (2025).
  2. Huang, Hsin-Yuan, et al. "Power of data in quantum machine learning." Nature communications 12.1 (2021): 2631.
  3. Daniels, Kyle G., et al. "Decoding CAR T cell phenotype using combinatorial signaling motif libraries and machine learning." Science 378.6625 (2022): 1194-1200.