Skip to main content
IBM Quantum Platform

投影量子カーネルを用いた特徴分類の強化

使用時間の目安：Heron r3 プロセッサーで80分（注：これはあくまでも目安です。 ランタイムは異なるかもしれない)。

学習成果

  • 投影量子カーネル（PQK）の仕組みと、どのような場合に量子優位性が期待できるか。
  • 実際のデータセットを使用して、ハードウェア上でPQKを実行する方法。

前提条件

背景

このチュートリアルでは、論文 Enhanced Prediction of CAR T-Cell Cytotoxicity with Quantum-Kernel Methods [1] に基づき、Qiskitを用いて投影量子カーネル （PQK）を実際の生物学的データセット上で実行する方法を示します。

PQKは量子機械学習（QML）で使われる手法で、特徴選択を強化するために量子コンピュータを使うことで、古典的なデータを量子特徴空間に符号化し、古典的な領域に投影し直す。 これは、量子回路を使って古典的なデータを量子状態にエンコードするもので、一般的には特徴マッピングと呼ばれるプロセスを経て、データを高次元ヒルベルト空間に変換する。 投影」とは、特定の観測量を測定することで、量子状態から古典的な情報を抽出し、サポートベクターマシンのような古典的なカーネルベースのアルゴリズムで使用できるカーネル行列を構築することである。 このアプローチは、量子システムの計算上の利点を活用することで、古典的な方法と比較して特定のタスクでより優れた性能を達成できる可能性がある。

PQKの主要な構成要素は、量子特徴マップの射影測定を通じて得られる還元密度行列（RDM）である。 特に、通常は各量子ビットについて、単一量子ビットの縮約密度行列（1 RDM）を計算する。 これらの測定値は、その後、指数カーネルなどの古典的なカーネル関数の入力として用いられ、最終的なカーネル行列が構築される。

PQKは、特に短期的に実用化が見込まれる量子ハードウェアにおいて、標準的な量子カーネルに比べて潜在的な利点をもたらす。 標準的な量子カーネルは、通常、グローバルな状態の重なりを推定することに依存していますが、量子ビットの数が増えるにつれてこれを正確に測定することはますます困難になり、またノイズの影響を強く受けます。 対照的に、PQKでは単一量子ビットの縮約密度行列（1 RDM）などの局所的な観測量を用いるため、サンプリングのオーバーヘッドが低減され、ハードウェアノイズに対する堅牢性が向上し、スケーラビリティも高まる。 PQKは、古典的なカーネル関数を適用する前に量子状態を局所的な測定特徴に射影することで、有用な量子相関を維持しつつ、近未来のデバイスにおいてより実用的なものとなることができる。

要件

このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：

  • Qiskit SDK v2.0 またはそれ以降、 可視化サポート付き
  • Qiskit Runtime v0.40 またはそれ以降 (pip install qiskit-ibm-runtime)
  • カテゴリーエンコーダ 2.8.1 (pip install category-encoders)
  • NumPy 2.3.2 (pip install numpy)
  • パンダ 2.3.2 (pip install pandas)
  • Scikit-learn 1.7.1 (pip install scikit-learn)
  • Tqdm 4.67.1 (pip install tqdm)

セットアップ

import warnings

# Standard libraries
import os
import urllib.request
from pathlib import Path
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 IBM Quantum Compute Service
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")

小規模シミュレータの例

このチュートリアルでは、小規模なシミュレータの例は省略します。なぜなら、私たちの主な目的は、投影量子カーネルがより大規模なシステムや実際のハードウェアにどのように拡張できるかを実証することだからです。

大規模なハードウェアの例

ステップ1：古典的な入力を量子問題にマッピングする

データセットの準備

このチュートリアルでは、Danielsら(2022)によって作成され、論文に含まれる補足資料からダウンロードできる、バイナリ分類タスクのための実世界の生物学的データセットを使用します。 このデータはCAR T細胞から構成されており、CAR T細胞は特定の癌を治療する免疫療法に使用される遺伝子操作されたT細胞である。 免疫細胞の一種であるT細胞は、がん細胞上の特定のタンパク質を標的とするキメラ抗原受容体（CAR）を発現するように研究室で改変される。 これらの改変されたT細胞は、がん細胞をより効果的に認識し、破壊することができる。 データの特徴はCAR T細胞モチーフであり、これはT細胞に工学的に組み込まれたCARの特定の構造的または機能的構成要素を指す。 これらのモチーフに基づき、与えられたCAR T細胞の細胞毒性を予測し、有毒か無毒かのラベルを貼るのが我々の仕事である。

このデータセットを前処理するヘルパー関数を以下に示す。

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"]),
        sep=",",
    )
    test_data = pd.read_csv(
        os.path.join(dir_root, args["file_test_data"]),
        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

以下のセルを実行することにより、このチュートリアルを実行することができます。これにより、必要なフォルダ構造が自動的に作成され、トレーニングファイルとテストファイルの両方があなたの環境に直接ダウンロードされます。 これらのファイルがすでにローカルにある場合は、この手順で安全に上書きし、バージョンの一貫性を確保する。

## Download dataset


def download_pqk_dataset(data_dir="data_tutorial/pqk"):
    """Download the four CSV files from the Qiskit documentation repo."""
    data_dir = Path(data_dir)
    data_dir.mkdir(parents=True, exist_ok=True)

    base_url = (
        "https://raw.githubusercontent.com/Qiskit/documentation/main/"
        "datasets/tutorials/pqk"
    )
    files = [
        "train_data.csv",
        "test_data.csv",
        "projections_train.csv",
        "projections_test.csv",
    ]

    for filename in files:
        url = f"{base_url}/{filename}"
        dest = data_dir / filename
        print(f"  {filename} ...", end=" ", flush=True)
        urllib.request.urlretrieve(url, dest)
        print(f"OK ({dest.stat().st_size:,} bytes)")

    print(f"\nAll files saved to {data_dir}/")
    return data_dir


DATA_DIR = download_pqk_dataset()

Output:

  train_data.csv ... OK (5,012 bytes)
  test_data.csv ... OK (2,194 bytes)
  projections_train.csv ... OK (779,730 bytes)
  projections_test.csv ... OK (335,529 bytes)

All files saved to data_tutorial/pqk/

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 points to the folder where the downloaded CSVs live
dir_root = str(DATA_DIR)

# 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
)

Output:

14:0:15

また、 11 が π/2\pi/2 として表現されるようにデータセットを変換し、スケーリングを行う。

# 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

トレーニングデータセットとテストデータセットのサイズと形状を検証する。

print(train_data.shape, train_labels.shape)
print(test_data.shape, test_labels.shape)

Output:

(172, 60) (172,)
(74, 60) (74,)

ステップ2：量子ハードウェア実行に向けた問題の最適化

量子回路

ここで、古典的なデータセットを高次元特徴空間に埋め込む特徴マップを構築する。 この埋め込みには ZZFeatureMap を使います。

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:

Output of the previous code cell

もう一つの量子埋め込みオプションは、 1D-Heisenberg ハミルトニアン進化アサッツである。 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:

Output of the previous code cell

ステップ3: Qiskit primitivesを使用して実行する

措置1-RDM

このステップでは、量子特徴マップの射影測定を通じて、すべての単一量子ビットの縮約密度行列（1-RDM）を取得します。これらは、後で古典的な指数カーネル関数に投入されます。

すべてのデータについて計算する前に、データセットから1つのデータ点を与えて1-RDMを計算する方法を見てみよう。 1-RDMは、すべての量子ビット上のパウリ X、 Y 、 Z 演算子の単一量子ビット測定のコレクションです。 なぜなら、1量子ビットのRDMは次のように完全に表現できるからである： ρ=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)

まず、使用するバックエンドを選択します。

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

それから量子回路を動かし、投影を測定する。 ゼロノイズ外挿（ZNE）を含むエラー緩和をオンにしていることに注意。

# 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])

次に結果を取り出す。

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)

Output:

[ 0.03530987 -0.06207794 -0.03529884 -0.1418671   0.00209782  0.0045834
  0.00407694  0.02528003  0.00233791  0.01800766  0.00718357  0.01927931
  0.0073651  -0.02009021  0.01144208  0.01333925  0.00521008  0.00535276
 -0.04354042 -0.0383848  -0.04472125  0.00641964 -0.03954627  0.03207479
  0.01823132  0.02546267 -0.          0.16288225  0.03246113  0.
  0.06107868  0.01082782  0.00240078  0.13147612  0.14033432  0.14925945
  0.11577918  0.00016128 -0.          0.00604693  0.02433089  0.02033885
  0.01492506  0.00494294  0.00926954  0.00569533  0.09867722  0.05662552
 -0.00001734  0.          0.          0.04625459 -0.02480763  0.01360688
  0.11511306  0.01260572 -0.01656313 -0.02510078 -0.03256272  0.00058607]
[-0.0756078  -0.05445208 -0.0228333  -0.00015029  0.00006226  0.02925132
 -0.00325556 -0.00889965  0.0177611  -0.00437065  0.01682502 -0.00229805
 -0.01041899 -0.03208967 -0.03515749  0.17477371  0.03783633  0.2126005
  0.          0.          0.00754466 -0.08242599  0.          0.03263675
  0.00399151 -0.01984418 -0.02106749 -0.02580491  0.03973411 -0.02037816
 -0.01769352 -0.09720746  0.00098896 -0.11840454  0.14392615  0.13647983
  0.08683845  0.04492138  0.0046172   0.04171398 -0.0000869  -0.00270916
 -0.0019876  -0.00440696  0.0307905  -0.0284622   0.11237189  0.15042867
  0.1020601  -0.03812461  0.00302523 -0.05240398 -0.01304566 -0.00403933
 -0.01324601 -0.03658085  0.00934269 -0.00105112 -0.          0.01761827]
[ 0.57921657  0.2865493   0.          0.00028356  0.03177571  0.01152152
  0.00843001  0.02320127  0.00273558  0.00976802  0.00060077  0.00942531
  0.00096361 -0.03950026  0.00560635  0.00591487  0.00788236  0.01346192
  0.60752971  0.80203507  0.65649176  0.00069473  0.06010304  0.05922109
  0.01670672  0.02900743  0.0162253   0.0668811   0.01573204 -0.00288162
  0.04216451  0.00848301  0.00052577 -0.33798808  0.68075471  0.89471233
  0.72272544  0.08096828  0.02387351  0.01723619  0.00774532  0.05513527
  0.08285531  0.08102448  0.10677406  0.27778995  0.28883482  0.21497224
  0.17569826  0.00063149  0.0320076   0.06735008 -0.00053637 -0.0006907
  0.00991596  0.00414575 -0.08425133 -0.09569482  0.00219474  0.00241873]

回路サイズと2量子ビットのゲート深さをプリントアウトする。

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)

Output:

qubits: 60
2q-depth: 96
2q-size: 2832
Operator counts: OrderedDict([('rz', 8640), ('sx', 7104), ('cz', 2832), ('x', 720), ('barrier', 47)])
Output of the previous code cell

これで訓練データセット全体をループし、すべての1-RDMを得ることができる。

また、量子ハードウェア上で行った実験の結果も提供する。 以下のフラグを True に設定することで、ご自身でトレーニングを実行することもできますし、弊社が提供する投影結果を使用することもできます。

# 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)

ジョブが完了したら、結果を取り出すことができる。

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])

これをテストセットでも繰り返す。

# 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)

以前と同じように結果を取り出すことができる。

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])

ステップ4：後処理を行い、結果を希望の古典形式で返す

予測量子カーネルを定義する

投影量子カーネルは以下のカーネル関数で定義される： kPQ(xi,xj)=exp(−γ∑k∑P∈{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) 上式において、 γ>0\gamma>0 は調整可能なハイパーパラメータである。 KijPQ=kPQ(xi,xj)K^{\textrm{PQ}}_{ij} = k^{\textrm{PQ}}(x_i, x_j) はカーネル行列 KPQK^{\textrm{PQ}} のエントリである。

1-RDMの定義を使えば、カーネル関数内の個々の項は Tr[Pρk(xi)]=⟨P⟩\textrm{Tr}[P \rho_k (x_i)] = \braket P、 P∈{X,Y,Z}P \in \{ X,Y,Z \} として評価できることがわかる。これらの期待値は、まさに我々が上記で測定したものである。

scikit-learn を使えば、カーネルをより簡単に計算することができる。 これは、すぐに利用できる放射基底関数('rbf')カーネルによるものである： exp(−γ∥x−x′∥2) \textrm{exp} (-\gamma \lVert x - x' \rVert^2)。まず、新しい投影された訓練データセットとテストデータセットを2次元配列に再形成する必要がある。

QPUでは、全データセットに約80分かかる。 チュートリアルの残りの部分を簡単に実行できるようにするために、以前に実行した実験からの投影を追加で提供します（これは、 Download dataset コードブロックでダウンロードしたファイルに含まれています）。 自分でトレーニングを行った場合は、自分の結果を使ってチュートリアルを続けることができる。

# ---------------------------------------------------------------------------
# Load projections — either from the hardware run above, or from the
# pre-computed CSVs that were downloaded alongside the motif data.
# ---------------------------------------------------------------------------

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(DATA_DIR / "projections_train.csv")
    projections_test = np.loadtxt(DATA_DIR / "projections_test.csv")

Support Vector Machine (SVM)

この事前計算されたカーネルで古典的なSVMを実行し、テストセットとトレーニングセットの間のカーネルを予測に使うことができる。

# Range of 'C' and 'gamma' values as SVC hyperparameters.
#
# This is a reduced grid so the tutorial runs quickly (154 candidates).
# The optimal (C, gamma) reported below lie within it, so the results
# are unchanged. The full grid used originally had 6622 candidates and
# took roughly one hours to search:
#
#   C_range = [0.001, 0.005, 0.007]
#   C_range.extend([x * 0.01 for x in range(1, 11)])   # 0.01 .. 0.10
#   C_range.extend([x * 0.25 for x in range(1, 60)])   # 0.25 .. 14.75
#   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])

C_range = [
    0.001,
    0.01,
    0.1,
    0.5,
    1.0,
    2.0,
    4.0,
    6.0,
    8.5,
    10.75,
    14.0,
    20,
    50,
    100,
]
gamma_range = [
    0.001,
    0.005,
    0.007,
    0.01,
    0.02,
    0.03,
    0.04,
    0.05,
    0.1,
    0.5,
    1.0,
]

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}")

Output:

Fitting 10 folds for each of 154 candidates, totalling 1540 fits
The best parameters are {'C': 8.5, 'gamma': 0.01} with a score of 0.6980
Test accuracy with best model: 0.8108

古典的ベンチマーキング

量子射影を行わなくても、放射基底関数をカーネルとする古典的なSVMを実行することができる。 この結果は、我々の古典的なベンチマークである。

# 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}")

Output:

Fitting 10 folds for each of 154 candidates, totalling 1540 fits
The best parameters are {'C': 10.75, 'gamma': 0.04} with a score of 0.7830
Test accuracy with best model: 0.7432

付録：学習タスクにおけるデータセットの潜在的な量子優位性の検証

すべてのデータセットがPQKの使用から潜在的な利点を得られるわけではない。 特定のデータセットがPQKの恩恵を受けられるかどうかの予備テストとして使える理論的な境界がいくつかある。 これを定量化するために、『 Power of data in quantum machine learning [2] 』の著者は、古典モデルと量子モデルの複雑さ、および古典モデルと量子モデルの幾何学的分離と呼ばれる量を定義している。 PQKから量子的な利点を期待するためには、古典カーネルと量子射影カーネルの幾何学的な分離はおよそ N\sqrt{N} のオーダーになるはずである。 NN はトレーニングサンプルの数である。 この条件が満たされれば、モデルの複雑さのチェックに移る。 もし古典的なモデルの複雑さが NN のオーダーであるのに対し、量子予測モデルの複雑さが NN よりも大幅に小さい場合、PQKの潜在的な利点が期待できる。

幾何学的分離は以下のように定義される（ [2] の F19 ）： gcq=g(Kc∥Kq)=∥KqKc(Kc+λI)−2KcKq∥∞g_{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)))

Output:

Geometric separation between classical and quantum kernels is 1.5440
13.114877048604

モデルの複雑さは以下のように定義される（ [2] の M1 ）： sK,λ(N)=λ2∑i=1N∑j=1N(K+λI)ij−2yiyjN+∑i=1N∑j=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}")

Output:

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}")

Output:

Quantum model complexity is 1.5806

次のステップ

推奨事項

この作品に興味を持たれた方は、以下の資料もご参照ください：

参照

  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. ダニエルズ カイル・G et al. "Decoding CAR T cell phenotype using combinatorial signaling motif libraries and machine learning." Science 378.6625 (2022): 1194-1200.
このページは役に立ちましたか？
バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。