GSoC 2025 | Quantum Kolmogorov-Arnold Networks for High Energy Physics Analysis at the LHC (Part I)

I’m Ria Khatoniar, a contributor to ML4Sci through Google Summer of Code. My project focuses on Quantum Kolmogorov-Arnold Networks for High-Energy Physics Analysis at the Large Hadron Collider, and this is a mid-term update on my progress so far.

All development for this project is open and can be accessed at the following GitHub repository: https://github.com/Ria-K912/QMLHEP/tree/main/Quantum_KAN_for_HEP_Ria_Khatoniar

Part II — https://medium.com/@riakhatoniar1234/gsoc-2025-quantum-kolmogorov-arnold-networks-for-high-energy-physics-analysis-at-the-lhc-part-8b44f5616e6f

Why Quantum Kolmogorov–Arnold Networks for the LHC?

The Large Hadron Collider (LHC) is on the verge of a massive upgrade-the High-Luminosity LHC (HL-LHC) will soon start delivering petabytes of data at an unprecedented scale!!!

While this opens up a treasure trove of scientific possibilities, it also presents a huge challenge: how do we make sense of such complex, high-dimensional data to uncover new physics?

That’s the heart of my project!

As part of GSoC 2025, I’m exploring a new frontier in quantum machine learning: Quantum Kolmogorov–Arnold Networks (QKANs).

So why QKANs?

At its core, QKANs are based on a powerful mathematical result: the Kolmogorov–Arnold theorem, which states that any multivariate continuous function can be broken down into a sum and composition of univariate functions. In simpler terms, even the most complex patterns can be reconstructed from simpler, one-variable functions , kind of like building a beautiful mosaic out of tiny tiles!

This property is incredibly attractive for machine learning, especially in the quantum world. By bringing this idea into quantum circuits, QKANs offer a way to build compact, interpretable, and highly expressive quantum neural networks. And since quantum computers naturally excel at handling entangled, high-dimensional states, QKANs are a natural fit for analyzing collider data.

I’m developing the model using PennyLane, and training it on realistic datasets inspired by particle physics. The goal is not just to build a working model, but to benchmark QKAN against classical networks (like MLPs and spline-based KAN), evaluate its scalability, and understand its potential in actual HL-LHC data analysis pipelines.

Background Research

Before diving into building QKAN, I spent time understanding it’s roots and the landscape of models it evolved from. Here’s what I found fascinating:

Most neural networks rely on fixed nonlinear activation functions at each node (like ReLU or sigmoid). But KANs flip this idea. Proposed by **Liu et al. (2024) , KANs use learnable univariate functions **on the edges of the network, not the nodes. These functions are typically splines and can flexibly adapt to the data, making the model both interpretable and powerful.

MLP vs KAN (Adopted from [2])

KANs also align with the Kolmogorov–Arnold theorem, offering a principled way to break down complex functions into simple parts. They showed strong results on several science benchmarks while maintaining fewer parameters than large MLPs. This made me wonder: could a quantum version push things even further?

KANs are great, but spline-based activations can be computationally expensive. That’s where SineKAN by **Reinhardt et al. (2024) comes in. Instead of splines, it uses sinusoidal functions** with trainable frequency and phase. This change dramatically improves training time and memory efficiency, without losing the theoretical benefits of KANs.

SineKAN became one of the baselines I compared my QKAN model against, as we shall see later, especially because its sinusoidal activations are also inspired by Fourier-like encodings which is a common theme in quantum circuits.

This is where things get really exciting!!!

In 2024, Ivashkov et al. introduced QKAN: Quantum Kolmogorov–Arnold Networks : the first truly quantum version of KANs. Instead of using splines or sine functions, they suggested univariate functions using Chebyshev polynomials and quantum singular value transformation (QSVT). This method allows the model to process quantum states directly.

CHEBYSHEV-QKAN layer (Adopted from [4])

It also opens the door for end-to-end quantum learning on high-dimensional data, using tools like block encodings and variational quantum layers. This became the technical backbone of my implementation.

Finally, I looked into how KANs can support quantum circuit design themselves. In KANQAS by **Kundu et al. (2024)** , researchers replaced MLPs in a reinforcement learning pipeline with KANs.

Schematic of the KANQAS algorithm (Adopted from [5])

The result?

The agent learned to generate shorter, more efficient quantum circuits which are crucial for today’s noisy quantum devices. This work showed me that KANs aren’t just great learners ; they can also learn how to build better quantum systems, making them doubly powerful in quantum machine learning.

Datasets Used

To thoroughly evaluate the performance of Quantum Kolmogorov–Arnold Networks (QKANs), I began by testing the model across a range of structured datasets. These included: Iris Dataset, Social Network Ads Dataset, Titanic Dataset and Digits Dataset.

As the project progressed, I moved toward more domain-relevant datasets from high energy physics:

Quark-Gluon Dataset-Derived from simulated jet images of quarks and gluons produced in high-energy proton-proton collisions. The dataset contains high-dimensional jet observables and calorimeter-based features that distinguish quark-initiated and gluon-initiated jets.

Jet Mass Dataset-A structured dataset focused on jet mass reconstruction, which is vital in tagging boosted objects like W/Z/H bosons or top quarks.

Higgs Boson Dataset-A benchmark dataset from a Kaggle competition hosted by CERN. It contains 28 features (some derived and some raw) and aims to classify whether a signal event corresponds to a Higgs boson decay or not. It reflects the real-world challenges of signal-to-noise separation in HEP.

Preprocessing

For all datasets, the input data was first transformed into a format compatible with quantum computation using diagonal block encoding.

What is Block-Encoding?

Block encoding allows a classical vector x ∈ R^N to be embedded as a subblock within a unitary matrix U, such that quantum operations can be performed on it efficiently. This is the foundation of quantum linear algebra methods like QSVT and LCU, explained later.

• For input dimensions N = 2^n, we typically require 2n qubits.
• The extra factor accounts for the data register and an ancilla/control register used to embed the data as a subblock of a larger unitary.

When working with quark-gluon and jet mass datasets, both of which are high-dimensional jet image representations, I experimented with several preprocessing strategies to reduce the number of input features:

Initial Hypothesis: Central Cropping

I first hypothesized that the center of the jet image contains the most meaningful information (i.e., high-energy jet substructure). I performed an 8×8×3 center crop on each original 125×125×3 image from the quark gluon dataset. However, the cropped images still contained many near-zero pixels, and accuracy remained poor even after normalization.

Empirical Analysis: Intensity Distributions

To investigate further, I analyzed the pixel intensity distributions per channel. What I found:

• Most pixels (across all channels) had very low values (< 10⁻¹).
• There were minor average intensity differences between quark and gluon classes, but many pixels were essentially noise.

Channel 0 pixel distribution for Quark and Gluon classes

Channel 1 pixel distribution for Quark and Gluon classes

Channel 2 pixel distribution for Quark and Gluon classes

This highlighted a key limitation of cropping alone: it didn’t necessarily preserve the most informative features.

Tried Tom Magorsch’s Discretization Approach

Inspired by an alternate method (Tom’s binning approach), I implemented:

• Cropping + resizing (to e.g., 8×8)
• Single channel selection (e.g., channel 1)
• Normalization and discretization (e.g., 2-bin quantization: 0 or 0.5)

with h5py.File("quark-gluon_train-set_n793900-001.hdf5", "r") as f:
    X_jets = f["X_jets"][:Nsamples]
    y = f["y"][:Nsamples]

def discretize_img(img, Nvals):
    bins = np.linspace(0, 1, Nvals + 1)
    vals = np.linspace(0, 1, Nvals)
    inds = np.digitize(img, bins)
    return vals[inds - 1]

def crop_and_resize_opencv(x, channel, crop_fraction, scale):
    cropped_imgs = []
    for i in range(x.shape[0]):
        img = x[i, :, :, channel]
        h, w = img.shape
        crop_size = int(h * crop_fraction)
        start = (h - crop_size) // 2
        cropped = img[start:start+crop_size, start:start+crop_size]
        resized = cv2.resize(cropped, (scale, scale), interpolation=cv2.INTER_LINEAR)
        cropped_imgs.append(resized)
    return np.array(cropped_imgs)

88 (channel 2) images plotted after using discretization approach*

This helped slightly, but results plateaued quickly. As can be seen in the above images, the cropping around the center approach has produced sparse representations, while reducing continuous pixel values to coarse bins (discretization) oversimplifies the data, resulting in the loss of critical substructure.

Final Approach: T-Test–Based Feature Selection

I finally moved to a statistically grounded method:

• Flattened each image (125×125×3 → 46,875 features)
• To identify which parts of the image were most useful in distinguishing quark and gluon jets, I conducted a statistical analysis on each pixel location
• Specifically, I calculated a t-statistic comparing the pixel intensity distributions between the two classes. This helped test whether the average pixel values at each location were statistically different across quark and gluon samples.

162 pixels found with p value less than 0.01

Out of all pixels, 162 were found to be significantly different with p-value < 0.01.Interestingly, many of these pixels were located outside the central region, contradicting the common assumption that central pixels carry the most signal. Furthermore, as seen from the pixel (i, j, channel) locations in the image, Channel 2 (HCAL subdetector) was the most informative, followed by Channel 1(ECAL), and then Channel 0 (tracks).

from scipy.stats import ttest_ind
import numpy as np

X_flat = X.reshape(X.shape[0], -1)  
quark_flat = X_flat[y == 0]
gluon_flat = X_flat[y == 1]
t_vals, p_vals = ttest_ind(quark_flat, gluon_flat, axis=0, equal_var=False)

top_indices = np.argsort(p_vals)[:192]  # top 8x8x3 = 192 pixel indices

X_selected = X_flat[:, top_indices]  # shape: (1000, 192)

from sklearn.model_selection import train_test_split
import torch

X_train, X_test, y_train, y_test = train_test_split(
    X_selected, y, test_size=0.2, random_state=42, stratify=y
)
X_train = torch.tensor(X_train, dtype=torch.float32)
X_test = torch.tensor(X_test, dtype=torch.float32)
y_train = torch.tensor(y_train, dtype=torch.long)
y_test = torch.tensor(y_test, dtype=torch.long)

By selecting pixels with low p-values, I retained only those features that were statistically discriminative between classes. This approach proved far more effective than assumptions about centrality or channel dominance and laid the foundation for better generalization in the QKAN architecture.

Model Architecture

Hybrid Quantum KAN

The architectural foundation of this hybrid model draws direct inspiration from Ivashkov et al. (2024), who proposed QKAN: Quantum Kolmogorov–Arnold Networks.

In particular, I adopted and extended the core idea of:

• Quantum polynomial function encoding via Quantum Singular Value Transformation (QSVT)
• Quantum-native aggregation mechanisms using controlled unitaries and Hadamard tests
• End-to-end training via variational parameters, integrated into a PyTorch-compatible pipeline.

While Ivashkov et al. presented the theoretical framework, I focused on scaling and applying this architecture to real-world HEP datasets, and additionally, I designed modular blocks for QSVT, LCU(Linear Combination of Unitaries), and quantum summation that are reusable and task-specific.

Step 1: QSVT-Based Quantum Encoding Layer

For each input feature, I apply a Quantum Signal Processing circuit that maps a scalar input into a vector of polynomial-transformed outputs. These quantum circuits are initialized to approximate composite nonlinear functions, and their parameters (angles) are learnable via gradient descent.

class QSVT(nn.Module):
    def __init__(self, wires=4, degree=5):
        super().__init__()
        self.wires = wires
        self.dev = qml.device("default.qubit", wires=wires)

        xs = np.linspace(-1, 1, 200)
        fx = sum(np.sin(k * np.pi * xs) for k in [1, 3, 5])
        poly = np.polyfit(xs, fx, deg=degree)
        poly /= np.max(np.abs(np.polyval(poly, xs))) + 1e-6
        odd_poly = [c if i % 2 == 1 else 0 for i, c in enumerate(poly[::-1])]
        self.qsvt_phis = nn.Parameter(torch.tensor(odd_poly, dtype=torch.float32), requires_grad=True)
        self.theta = nn.Parameter(torch.randn(wires))

        @qml.qnode(self.dev, interface="torch", diff_method="backprop")
        def circuit(x, phis, theta):
            for i in range(wires):
                qml.RY(x[i], wires=i)
            A = torch.diag(x.detach())
            ph_np = phis.detach().cpu().numpy()
            ph_np /= np.max(np.abs(ph_np)) + 1e-6
            qml.qsvt(A, ph_np, encoding_wires=list(range(wires)), block_encoding="embedding")
            for i in range(wires):
                qml.RZ(theta[i], wires=i)
            return [qml.expval(qml.PauliZ(i)) for i in range(wires)]

        self.qnode = circuit

    def forward(self, x):
        return self.qnode(x, self.qsvt_phis, self.theta)

Step 2: Quantum LCU Block

Once I had polynomial-transformed outputs per feature, the next challenge was: how do we weight these outputs meaningfully?

I used a Quantum Linear Combination of Unitaries (LCU) approach. Here:

• Each polynomial component is weighted quantumly via state preparation.
• A trainable vector of weights is encoded as a quantum state, which is then used to control parameterized rotations over a target wire.
• This block performs an elementwise multiplication between the QSVT polynomial outputs and their corresponding quantum weights for each feature.

def quantum_lcu_block(qsvt_vals, weight_vals):
    """
    qsvt_vals: tensor of shape (P,)
    weight_vals: tensor of shape (P,)
    Output: tensor scalar (Z expectation after quantum weighting)
    """
    P = len(qsvt_vals)
    n_ctrl = math.ceil(math.log2(P))
    wires = list(range(n_ctrl + 1))  # control + 1 target
    dev = qml.device("default.qubit", wires=len(wires))

    alpha = torch.sqrt(torch.abs(weight_vals) + 1e-8)
    alpha = alpha / torch.norm(alpha + 1e-8)

    @qml.qnode(dev, interface="torch", diff_method="backprop")
    def circuit():
        qml.StatePrep(alpha, wires=wires[:-1])

        for i in range(P):
            ctrl_bin = [int(b) for b in f"{i:0{n_ctrl}b}"]
            qml.ctrl(qml.RY, control=wires[:-1], control_values=ctrl_bin)(2 * qsvt_vals[i], wires=wires[-1])

        qml.adjoint(qml.StatePrep(alpha, wires=wires[:-1]))

        return qml.expval(qml.PauliZ(wires[-1]))

    return circuit()

Step 3: Quantum Summation using Hadamard Test

Each feature now has a quantum-weighted set of polynomial outputs. But we still need to aggregate them into a single scalar per feature to proceed further.

For this, I implemented a Hadamard test–based summation block:

• It uses an ancilla qubit and applies controlled phase shifts based on the weights and polynomial values.
• The result is a scalar expectation value that reflects the sum of weighted components via quantum interference.
• Each feature gets its instance of this summation block with separate trainable weights.

class QuantumSumBlock(nn.Module):
    def __init__(self, num_polynomials):
        """
        Quantum summation over N polynomial outputs using Hadamard test.
        Supports arbitrary N (must be <= 2^num_index_qubits).

        Args:
            phi_vals (tensor): shape (N,), tensor of real-valued polynomial outputs.

        Returns:
            torch.Tensor: scalar value (quantum sum approximation via Hadamard test)
        """
        super().__init__()
        self.num_polynomials = num_polynomials
        self.weights = nn.Parameter(torch.randn(num_polynomials))  

        self.num_index_qubits = math.ceil(math.log2(num_polynomials))
        self.total_wires = self.num_index_qubits + 2  # index + hadamard + target
        self.wires = list(range(self.total_wires))

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

        @qml.qnode(self.dev, interface="torch", diff_method="backprop")
        def circuit(phi_vals, wts):
            index_wires = self.wires[:self.num_index_qubits]
            hadamard_wire = self.wires[self.num_index_qubits]
            target_wire = self.wires[-1]

            for i in index_wires:
                qml.Hadamard(wires=i)
            qml.Hadamard(wires=hadamard_wire)

            for p in range(self.num_polynomials):
                ctrl_bin = [int(b) for b in f"{p:0{self.num_index_qubits}b}"]
                ctrl_values = [1] + ctrl_bin
                qml.ctrl(qml.RZ, control=[hadamard_wire] + index_wires, control_values=ctrl_values)(
                    2 * wts[p] * phi_vals[p], wires=target_wire
                )

            qml.Hadamard(wires=hadamard_wire)
            return qml.expval(qml.PauliZ(hadamard_wire))

        self.circuit = circuit

    def forward(self, phi_vals):
        return self.circuit(phi_vals, self.weights)

Step 4: Final Prediction via Classical KAN Layer

Once all feature-level quantum blocks are complete, I’m left with a final vector of size equal to the number of input features. This acts as a quantum-processed representation of the sample.

To extract the final prediction, I pass this vector through a classical KAN (SineKAN/SplineKAN) layer.

class QuantumKANClassifier(nn.Module):
    def __init__(self, num_features, degree=3, num_classes=2):
        super().__init__()
        self.num_features = num_features
        self.degree = degree

        self.qsvt = QSVT(wires=1, degree=degree, depth=2)
        self.lcu_weights = nn.Parameter(torch.rand(num_features, degree))  # (F, P)
        self.sum_blocks = nn.ModuleList([QuantumSumBlock(degree) for _ in range(num_features)])
        self.kan = KANLayer(in_features=num_features, out_features=num_classes) 

    def forward(self, X):
        """
        Input: X of shape (B, F)
        Output: logits of shape (B, num_classes)
        """
        B = X.size(0)
        all_features = []

        for i in range(B):
            xi = X[i]  # shape: (F,)
            qsvt_vecs = [self.qsvt(xi[f]) for f in range(self.num_features)]  # each: (P,)
            lcu_vals = [quantum_lcu_block(qsvt_vecs[f], self.lcu_weights[f]) for f in range(self.num_features)]
            summed = [self.sum_blocks[f](lcu_vals[f]) for f in range(self.num_features)]
            all_features.append(torch.stack(summed))  # shape: (F,)

        features = torch.stack(all_features)  # shape: (B, F)
        return self.kan(features)  # shape: (B, num_classes)

Experimental Results

Regression Tasks

To evaluate the expressive power of the QuantumKAN-inspired model, I conducted experiments on a variety of nonlinear one-dimensional regression tasks. These functions were chosen to represent diverse complexity-sharp transitions, oscillations, exponential decay, and compound nonlinearity.

The tested functions include:

• Composite nonlinear functions f(x)=tanh(10x+0.5+ReLU(x2)⋅10) h(x)=σ(3x)+ReLU(sin⁡(2x)+x3) m(x)=Softplus(x2−1)+tanh⁡(4x+0.1) n(x)=e^(−x^2+0.3x)+ReLU(tanh⁡(2x−1))
• Canonical benchmark functions bump(x) = exp(−10(x − 0.2)²) + exp(−50(x + 0.5)²) runge(x) = 1 / (1 + 25x²) exp_sin(x) = exp(−x²) · sin(5x) noisy_step(x) = Heaviside(x) + 0.1 · sin(20x) sigmoid_bumps(x) = σ(8(x − 0.5)) + σ(−10(x + 0.3)) − 1 sawtooth(x) = 2(x − floor(x + 0.5))

The Hybrid QKAN model can capture complex shapes, including sharp non-differentiabilities with stable convergence on both training and test sets.

Classification Tasks

To further assess the generalization ability of our hybrid model, I conducted experiments on multiple tabular classification datasets of varying dimensionality and difficulty. The architecture utilized spline-based Kolmogorov–Arnold Networks (KANs).

Accuracy Comparison

The slight dip in accuracy with QKAN can be attributed to quantum noise and hardware constraints, but results remain competitive and demonstrate the feasibility of hybrid learning on tabular data.

I further evaluate the models on particle physics datasets, comparing MLP, SineKAN, and Hybrid QKAN (utilising SineKAN layer ) architectures.

Quark Gluon Dataset Accuracy Comparision

However, due to quantum circuit overhead and memory limitations, I was unable to scale the Hybrid QKAN architecture beyond this point for the Quark Gluon Dataset, as larger configurations resulted in kernel crashes during training.

Higgs Boson Dataset Accuracy Comparision

Jet Mass Dataset Accuracy Comparision

Across all three high-dimensional physics datasets, the SineKAN model consistently performed competitively with a standard MLP baseline, demonstrating its expressive power even in complex domains. While the Hybrid QKAN model shows slightly lower accuracy, it maintains reasonable performance, validating the feasibility of integrating quantum components in real-world scenarios.

Conclusion & Future Directions

As the project evolves, HybridQKAN has emerged as a promising architecture, effectively combining QSVT-based quantum encodings with trainable quantum aggregation blocks and a classical KAN layer. While it currently incurs quantum overhead, the model consistently captures complex nonlinear patterns across diverse datasets, highlighting its potential for high-energy physics analysis.

Moving forward in the remainder of GSoC, I aim to make HybridQKAN more efficient and autonomous by integrating the t-test-based pixel selection directly into the model pipeline. This will allow the network to prioritize informative features dynamically. Additionally, I plan to explore the design of a fully quantum Kolmogorov–Arnold Network, paving the way for fully end-to-end quantum-native learning.

Acknowledgements

I would like to express my sincere gratitude to Sergei Gleyzer for his exceptional leadership and coordination of the ML4Sci program. I am also deeply thankful to everyone who shared valuable feedback and guidance throughout this journey, with special appreciation for my mentors-Eric Reinhardt, Roy T. Forestano, KC Kong, Konstantin Matchev, Katia Matcheva, Jogi Suda, and Amey Bhatuse for their constant support and insights.

References

[1] CERN. High-Luminosity LHC (HL-LHC): The Future of the LHC. Retrieved from https://hilumilhc.web.cern.ch/

[2] Liu, Z., Wang, Y., Vaidya, S., Ruehle, F., Halverson, J., Soljaˇci ́c, M., Hou, T.Y.,& Tegmark, M. (2025). KAN: Kolmogorov–Arnold Networks. Retrieved from https://arxiv.org/abs/2404.19756

[3] Reinhardt, E., Ramakrishnan, D., & Gleyzer, S. (2024). SineKAN: Kolmogorov–Arnold Networks using Sinusoidal Activation Functions. Retrieved from https://www.frontiersin.org/journals/artificial-intelligence/articles/10.3389/frai.2024.1462952/full

[4] Ivashkov, P., Huang, P.-W., Koor, K., Pira, L., & Rebentrost, P. (2024). QKAN: Quantum Kolmogorov-Arnold Networks. Retrieved from https://arxiv.org/html/2410.04435v1

[5] Kundu, A., Sarkar, A., & Sadhu, A. (2024). KANQAS: Kolmogorov-Arnold Networks for Quantum Architecture Search. Retrieved from https://arxiv.org/abs/2406.17630