# Limited Memory Online Gradient Descent for Kernelized Pairwise Learning
with Dynamic Averaging

Hilal AlQuabeh1, William de Vazelhes1, Bin Gu1,2 Corresponding author.

###### Abstract

Pairwise learning, an important domain within machine learning, addresses loss
functions defined on pairs of training examples, including those in metric
learning and AUC maximization. Acknowledging the quadratic growth in
computation complexity accompanying pairwise loss as the sample size grows,
researchers have turned to online gradient descent (OGD) methods for enhanced
scalability. Recently, an OGD algorithm emerged, employing gradient
computation involving prior and most recent examples, a step that effectively
reduces algorithmic complexity to $O(T)$, with $T$ being the number of
received examples. This approach, however, confines itself to linear models
while assuming the independence of example arrivals. We introduce a
lightweight OGD algorithm that does not require the independence of examples
and generalizes to kernel pairwise learning. Our algorithm builds the gradient
based on a random example and a moving average representing the past data,
which results in a sub-linear regret bound with a complexity of $O(T)$.
Furthermore, through the integration of $O(\sqrt{T}{\log{T}})$ random Fourier
features, the complexity of kernel calculations is effectively minimized.
Several experiments with real-world datasets show that the proposed technique
outperforms kernel and linear algorithms in offline and online scenarios.

## Introduction

The concept of online pairwise learning has attracted significant attention
due to its unique characteristics and potential applications. Unlike
traditional machine learning, online pairwise learning delves into
understanding connections between pairs of entities or data points, rather
than isolating individual instances. This methodology holds particular
relevance in domains where relative evaluations between entities hold
substantial import, as seen in Area Under the ROC Curve (AUC) maximization
(e.g. in unbalanced-data classification) (Zhang et al. 2016), metric learning
(e.g. in image retrieval) (Kulis et al. 2012), and bipartite ranking (Du et
al. 2016). Given that the loss functions in these problems require pairs of
examples, a distinctive analysis is needed when compared to other machine
learning problems.

At the core of online pairwise learning lies the challenge of efficiently
updating models in real-time. This challenge emerges from the cumulative
influence of previous instances on the comprehension of new pairs. In this
landscape, Online Gradient Descent (OGD) stands as a significant tool within
the realm of online learning. OGD takes on a distinctive role within the
sphere of online pairwise learning, entailing a nuanced examination of
evolving loss gradients and the impact of prior examples with each incoming
instance. This nuanced viewpoint results in a computational complexity
translating to $O(T^{2})$ of gradient computations.

OGD has garnered considerable attention in the field of online pairwise
learning, leading to the development of various methodologies. These
approaches include online buffer learning (Zhao et al. 2011; Kar et al. 2013),
second-order statistic online learning (Gao et al. 2013), and methods rooted
in solving the saddle point problem (Ying, Wen, and Lyu 2016; Reddi et al.
2016), all of which have predominantly harnessed linear models. However,
exploration into non-linear models within this domain has been somewhat
limited, with kernelized learning receiving minimal attention (Ying and Zhou
2015; Du et al. 2016).

A notable avenue in the realm of pairwise learning is online buffer learning,
as introduced by Zhao et al. (Zhao et al. 2011). This technique employs a
finite buffer alongside reservoir sampling to mitigate time complexity,
reducing it to $O(sT)$, where $s$ signifies the buffer size. By retaining a
subset of data and ensuring uniform sampling within the buffer, this strategy
effectively mitigates the computational load. Moreover, Yang et al. (Yang et
al. 2021) achieved optimal generalization with a buffer size of $s=1$, all
while assuming the sequential data to be iid instances (e.g., stochastic
environments). This accomplishment marks a significant milestone in the study
of pairwise learning.

Algorithm | Problem | Model | Scheme | iid | Time Complexity | Space Complexity
---|---|---|---|---|---|---
(Natole, Ying, and Lyu 2018) | AUC | Linear | Online | Yes | $O(T)$ | $O(d)$
(Ying, Wen, and Lyu 2016) | AUC | Linear | Online | Yes | $O(T)$ | $O(d)$
(Zhao et al. 2011) | AUC | Linear | Online | Yes | $O(sT)$ | $O(ds)$
(Gao et al. 2013) | AUC | Linear | Online | Yes | $O(T)$ | $O(d^{2})$
(Gu, Huo, and Huang 2019) | AUC | Linear | Offline | Yes | $O(T)$ | $O(d)$
(Yang et al. 2021) | General | Linear | Online | Yes | $O(T)$ | $O(d)$
(Lin et al. 2017) | General | kernel | Online | No | $O(T^{2})$ | $O(dT)$
(Kakkar et al. 2017) | AUC | Kernel | Offline | Yes | $O(T\log{T})$ | $O(dT^{2})$
(Kar et al. 2013) | General | Linear | Online | Yes | $O(sT)$ | $O(ds)$
AOGD | General | Kernel | Online | No | $O(\frac{D}{d}T)$ | $O(D)$
Table 1: Recent pairwise learning algorithms ($T$ is the iteration number,
$d$: the dimension, $D$ is random features, and $s$ is a buffer size), note
the time complexity is w.r.t. gradients computations.

The prevailing frameworks in existing literature have predominantly
concentrated on data that is linearly separable, disregarding the complexities
linked with non-linear pairwise learning. Furthermore, the online buffer
techniques introduced thus far have not adequately tackled the sensitivity of
generalization to the non-iid data, particularly in adversarial environments,
constraining their capacity to grasp the intricacies inherent in real-world
scenarios.

Furthermore, the exploration of non-linear pairwise learning, particularly in
the realm of kernel approximation, has been relatively limited. While non-
linear methods offer enhanced expressive capability, the computational expense
attributed to kernel calculation – scaling as $O(T^{2})$ (Lin et al. 2017;
Kakkar et al. 2017) – poses hurdles to their practical applicability in terms
of scalability and efficiency.

Transitioning to the examination of generalization bounds, prior research
extensively investigated online pairwise gradient descent with buffers and
linear models (Wang et al. 2012; Kar et al. 2013). These studies establish a
bound of $O(1/s+1/\sqrt{T})$ for this approach. However, it’s imperative to
acknowledge that this bound attains optimality primarily when the buffer size
$s$ approximates $O(\sqrt{T})$. This aspect poses challenges in scenarios
where a smaller buffer size is preferred.

Collectively, these vulnerabilities underscore the necessity for further
research and advancement in the field of online pairwise learning. This effort
should encompass addressing the limitations of iid frameworks, conducting a
more comprehensive exploration of non-linear methods, and surmounting the
computational obstacles tied to kernel computation. Our method expands the
scope of online pairwise learning to accommodate nonlinear data by
incorporating kernelization into the input space. We tackle the ramifications
of non-iid data on regret by evaluating the gradient of both the moving
average and a random sample drawn from the past history. Leveraging random
Fourier features, we efficiently estimate the kernel while maintaining a
sublinear error bound, thus achieving computational efficiency without
compromising performance. Through the integration of kernelization, kernel
approximation, and dynamic averaging, our approach transcends linear
constraints, effectively addresses non-iid data challenges, and mitigates
kernel complexity. As a result, we establish a robust online pairwise learning
paradigm (as illustrated in Table 1). Our principal contributions can be
summarized as follows:

  * •

We introduce an online pairwise algorithm for non-linear models, characterized
by the ability to handle non-iid data. In scenarios featuring a wide margin in
the kernel space, our algorithm attains sublinear regret with a buffer size of
$O(1)$.

  * •

We confront the influence of non-iid data on regret and introduce a random
Bernoulli sampling strategy to manage and enhance the regret rate.

  * •

In the case of Gaussian kernel, we approximate the pairwise kernel function
utilizing a mere $O(\sqrt{T}\log{T})$ features, a significant reduction
compared to $O(T)$ features in literature, while keeping a sublinear regret
bound.

  * •

We validate the effectiveness of our proposed methodology across a range of
real-world datasets and compare its performance against state-of-the-art
methods for AUC maximization.

The subsequent sections are structured as follows: Section 2 delineates the
problem setting, Section 3 expounds on our proposed method, Section 4
furnishes the analysis of regret, Section 5 delves into related work, followed
by Section 6 presents experimental results, and lastly, Section 7 concludes
the paper.

## Problem Setting

The notion of pairwise learning pertains to a specific subset
$\mathcal{X}\subset\mathbb{R}^{d}$ and a label space
$\mathcal{Y}\subset\mathbb{R}$. This concept branches into two categories, as
discussed in more detail by (AlQuabeh and Abdurahimov 2022). Our focus lies in
analyzing pairwise loss functions from both branches, connecting the pairwise
kernel associated with pairwise hypotheses to regular kernels. This link
allows us to explore the characteristics of pairwise loss functions within the
regular kernel framework. In our study, we examine an algorithmic process that
learns from a set of instances denoted as
${z_{i}:=(x_{i},y_{i})\in\mathcal{Z}:=\mathcal{X}\times\mathcal{Y}}$, where
the index $i$ ranges from $1$ to $T$, encompassing the total number of
received examples. Assuming that the function $f$ belongs to a specific space
denoted as $\mathcal{H}$, we observe that the pairwise loss function plays a
crucial role in evaluating performance. This function is denoted as
$\ell:\mathcal{H}\times\mathcal{Z}^{2}\rightarrow\mathbb{R}+$, signifying that
it takes pairs of examples from the combined domain $\mathcal{Z}^{2}$, and
yields values within the positive real numbers $\mathbb{R}+$.

In online learning with pairwise losses, upon receiving a new data point
$z_{t}$, a localized error is used as a performance measure (Zhang et al.
2016; Boissier et al. 2016). This error represents the pairing of the new data
point with all preceding $t-1$ points. The properties of this localized error
depend on the selected pairwise loss function, illustrated as:

${L}_{t}(w_{t-1})=\frac{1}{t-1}\sum_{i=1}^{t-1}\ell(w_{t-1},z_{t},z_{i})$ (1)

where $w\in\mathcal{H}$ represents a linear function mapping, i.e.
$f_{w}(x)=\langle w,x\rangle$. To tackle the memory constraints, (Zhao et al.
2011) introduced a buffer-based mechanism for local error computation denoted
as $\bar{L}_{t}(w_{t-1})$, as explicated in Equation 2. At each iteration $t$,
the buffer, designated as $B_{t}$, accommodates a finite set of historical-
examples indices, and the size of this buffer is $|B_{t}|$ (or $s$ in
literature).

$\bar{L}_{t}(w_{t-1})=\frac{1}{|B_{t}|}\sum_{i\in
B_{t}}\ell(w_{t-1},z_{t},z_{i})$ (2)

This buffer has an important role in the learning journey, its contents
updated at each step through diverse strategies that span from randomized
methodologies such as reservoir sampling (Zhao et al. 2011; Kar et al. 2013)
to non-randomized approaches exemplified by FIFO (Yang et al. 2021). However,
despite the wide adoption of these sampling techniques, there exists a
discernible research void in understanding their implications within the
context of non-iid online examples. It’s imperative to highlight that the
issue of buffer size and its update methods becomes pronounced when dealing
with adversarial environment scenarios. The choice of buffer size and sampling
strategy can significantly impact the model’s performance and ability to
generalize in the case of non-iid data streams.

In our pursuit of addressing the intricacies inherent in complex real-world
data, our online pairwise approach takes into consideration the mapping of
both the hypothesis and the data to a Reproducing Kernel Hilbert Space (RKHS),
denoted as $\mathcal{H}$. The associated Mercer pairwise kernel function,
denoted as $k:\mathcal{X}^{4}\mapsto\mathbb{R}$, adheres to the reproducing
property $\langle k_{(x,x^{\prime})},g\rangle=g(x,x^{\prime})$, where
$x,x^{\prime}\in\mathcal{X}^{2}$ and $g\in\mathcal{H}$. Notably, for scenarios
involving pointwise hypotheses with pairwise losses, such as AUC loss, the
kernel function is streamlined to $k:\mathcal{X}^{2}\mapsto\mathbb{R}$. The
expanse of $\mathcal{H}$ encompasses a spectrum of linear combinations
encompassing functional mappings
${k_{(x,x^{\prime})}|(x,x^{\prime})\in\mathcal{X}^{2}}$ and their limit
points.

In our quest to circumvent the computational intricacies associated with
kernelization in the online realm, we employ the resourceful apparatus of
random Fourier features (RFF) as a pragmatic approximation for the Mercer
kernel function. RFF introduces a lower-dimensional mapping $r(\cdot)$,
serving as an approximation of the kernel function, denoted as
$\bar{k}(\cdot)$. This approximation empowers us to execute computations
through linear operations, culminating in a significant reduction in
computational overhead. The domain defined by the novel kernel functions finds
its place in the realm of $\bar{\mathcal{H}}$. It’s worth mentioning that
earlier investigations have explored the precision of random Fourier
approximation within pointwise and offline settings. As we navigate the online
landscape, the minimum number of random features necessary to ensure sublinear
regret has been ascertained to be $O(T)$. Our method introduces a novel error
bound for pairwise predicaments, hinging on the utilization of merely
$O(\sqrt{T}\log{T})$ random features (elucidated in Section 5 for
comprehensive details).

### Assumptions

Prior to unveiling our principal theorems, we lay down a foundation comprising
a collection of widely embraced assumptions pertaining to the attributes of
both the loss function and kernels.

######  Assumption 1 (Lipschitz Continuity).

Assume for any $z,z^{\prime}\in\mathcal{Z}$, the loss function
$\ell(\cdot,z,z^{\prime})$ is G-Lipschitz continuous, i.e. $\forall
w,w^{\prime}\in\mathcal{H}$,

$\displaystyle|\ell(w,z,z^{\prime})-\ell(w^{\prime},z,z^{\prime})|\leq
G\left\|w-w^{\prime}\right\|_{2}.$

######  Assumption 2 (M-smoothness ).

Assume the gradient of the loss function $\ell(\cdot,z,z^{\prime})$ is
M-Lipschitz continuous, i.e. $\forall w,w^{\prime}\in\mathcal{H}$ and for any
$z,z^{\prime}\in\mathcal{Z}$,

$\displaystyle\ell(w)\geq\ell(w^{\prime})$
$\displaystyle+\nabla\ell(w^{\prime})^{T}(w-w^{\prime})$
$\displaystyle+(2M)^{-1}\|\nabla\ell(w)-\nabla\ell(w^{\prime})\|^{2}$

######  Assumption 3 (Convexity).

Assume for any $z,z^{\prime}\in\mathcal{Z}$, the loss function
$\ell(\cdot,z,z^{\prime})$ is convex function, i.e. $\forall
w,w^{\prime}\in\mathcal{H}$,

$\displaystyle\ell(w,z,z^{\prime})\geq\ell(w^{\prime},z,z^{\prime})+\nabla\ell(w^{\prime},z,z^{\prime})^{T}(w-w^{\prime}).$

######  Assumption 4 (Finite Kernel).

Assume for any $\rho$-probability measure on $\mathcal{X}^{2}$ the positive
kernel function $k:\mathcal{X}^{2}\times\mathcal{X}^{2}\rightarrow\mathbb{R}$
is $\rho$-integrable, i.e. for any $(x,x^{\prime})\in\mathcal{X}^{2}$,

$\displaystyle\int\int_{\mathcal{X}^{2}}k((x,x^{\prime}),(\hat{x},\hat{x}^{\prime}))d\rho(\hat{x})d\rho(\hat{x}^{\prime})<\infty.$

## Methodology

In light of the core problem we aim to address, namely the computational
complexity arising from the necessity to accommodate all preceding examples in
the local error computation, our proposed methodology centers around
mitigating this challenge. We begin by suggesting the utilization of a moving
average of previous examples, thus alleviating the burden of handling the
complete historical dataset.

However, it’s imperative to note that the simple average approach may lead to
misleading outcomes, especially when the data exhibits high variance
distributions. Recognizing this potential limitation, we introduce a
corrective measure to enhance the reliability of the average-based gradient.
This correction mechanism involves introducing a random point from the
historical loss, enhancing the robustness of the gradient estimation, and
ameliorating the impact of potential average imbalances. This comprehensive
strategy effectively tackles the issues associated with handling the entire
historical dataset, while also accounting for data distribution complexities
that can skew outcomes, and most importantly not requiring the iid assumption.
Another parallel research addressing the high variance in the buffer has been
proposed recently (AlQuabeh, Mukhoty, and Gu 2023), however with $O(s)$
complexity.

Indeed, the loss function we seek to minimize can be formulated as a linear
combination of two distinct losses: one based on the moving average and the
other on a randomly chosen example. However, it’s important to recognize that
the minimization of this combined loss doesn’t inherently guarantee the
minimization of the true local loss. This discrepancy arises because, at the
optimal solution where the gradient is zero, the gradients of both components
remain nonzero, depending upon the linear coefficient constant.

To navigate this challenge, an alternating strategy to minimize both losses is
studied. In essence, we proceed by iteratively taking steps using the average-
based gradient and then employing the random gradient on the transition model
derived from the previous step. This alternating optimization strategy
accounts for the intricate balance between the two gradient components and
fosters a dynamic convergence process. Moreover, the iterative method ensures
that both the average-based and random-based gradients are appropriately
accounted for in the optimization process. The loss minimization upon
receiving new examples $z_{t}$ can be formulated using the average
$\bar{z}=(\sum_{i=1}^{t-1}\frac{x_{i}}{t-1},y_{t-1}|{y_{t-1}\neq y_{t}})$ and
the random sample $\hat{z}$ from $\\{z_{1},\dots,z_{t-1}\\}$, as follow;

$\min_{w}\ell\left(w-\eta_{t}\nabla\ell(w,z_{t},\bar{z}),z_{t},\hat{z}\right)$
(3)

where $\eta_{t}$ is the step size. As detailed in the subsequent section, the
magnitude of the step taken in the second update is contingent upon both the
variance of the data and the smoothness of the loss function. Interestingly,
when the variance is negligible, the step taken based on the average gradient
is deemed adequate, rendering the random step unnecessary. The objective
function can give more sense if we use second-order Taylor expansion after
denoting $\hat{\ell}(\cdot)=\ell(\cdot,z_{t},\hat{z})$ and
$\bar{\ell}(\cdot)=\ell(\cdot,z_{t},\bar{z})$ as follow,

$\displaystyle\ell(w-\eta_{t}\nabla\ell(w,z_{t},\bar{z}),z_{t},\hat{z})$
$\displaystyle\approx\hat{\ell}(w)-\eta_{t}\nabla\hat{\ell}(w)^{T}\nabla\bar{\ell}(w)$
$\displaystyle+\eta_{t}^{2}$
$\displaystyle\nabla\bar{\ell}(w)^{T}\nabla^{2}\hat{\ell}(w)\nabla\bar{\ell}(w)$
(4)

Hence, the process of minimizing the objective function involves addressing
three integral components. First, there is the endeavor to minimize the loss
associated with randomly selected examples. Then, leveraging the convex nature
of the loss function and the positive definiteness of its Hessian, the third
term—strictly non-negative in nature—ensures that the gradient of the average-
based loss is minimized. Of paramount significance is the role of the second
term, which guarantees that the gradients of losses, both average-based and
random-based, exhibit a maximum inner product. This intricate interplay is
pivotal in correcting the moving average gradient, particularly in scenarios
involving highly skewed distributions.

In order to handle the intricate patterns of real-world data, our approach to
online pairwise learning involves a fundamental assumption: both the
hypotheses and the data are transformed into a specialized space known as a
Hilbert space, denoted as $\mathcal{H}$. This space is equipped with a
specific type of kernel called a Mercer kernel. Essentially, this kernel is a
mathematical function that captures the relationships between data points.
This Hilbert space has two key attributes: first, it can reproduce the values
of the kernel function, and second, it encompasses various combinations of
these kernel functions along with their complete structures.

When we apply traditional buffer-based pairwise algorithms in an online
setting, the process of kernelization comes with a complexity that grows as
$O(T^{2}s)$, where $s$ represents the size of the buffer—a rather unwieldy
task for larger problems. To address this challenge, we introduce a method
that approximates the kernel function using a simpler approach. We assume the
existence of a lower-dimensional mapping, represented as
$r:\mathcal{X}^{2}\mapsto R^{D}$. This mapping enables us to simplify the
kernel function’s calculations, significantly reducing the complexity of the
process. This simplified approach is achieved using what are known as random
Fourier features, which are derived from the mathematical Fourier transform.
Although these techniques have been studied in previous works, our focus is on
adapting them to the online setting, and we establish that as few as
$O(\sqrt{T}\log{T})$ random features are enough to achieve effective results.

Our method is introduced in algorithm 1 for general pairwise learning, and it
includes the Gaussian kernel approximation using random Fourier features. In
the upcoming section, we delve into the performance of this online algorithm
and dissect its regret into two distinct components: the first is a direct
result of learning a model from an online stream of examples, and the second
arises from the approximation of the kernel function.

Algorithm 1 Average Online Gradient Descent (AOGD)

0: Random feature size $D$, initial solution $w_{1}\in\bar{\mathcal{H}}$,
initial average and random example $\bar{z}=\hat{z}\in R^{D}$, kernel function
$k$ and its Fourier transform $P(u)$, probability $p$, step sizes
$\eta_{t},\gamma_{t}$ and $\lambda$.

1: Sample Fourier feature $\\{u_{i}\\}_{i=1}^{D/2}\sim P(u_{i})$

2: for $t=2,\dots,T$ do

3: Receive a new example $z_{t}=(x_{t},y_{t})\in\mathcal{Z}$

4: Map to the new space:
$z_{t}=(\frac{2}{\sqrt{D}}[\sin(u_{i}^{T}x_{t}),\cos(u_{i}^{T}x_{t})]_{i=1}^{D/2},y_{t})$

5: Suffer loss $\bar{\ell}(w_{t-1},z_{t},\bar{z})$

6: Update model
$w^{\prime}_{t}=w_{t-1}-\eta_{t}\nabla\bar{\ell}(w_{t-1},z_{t},\bar{z})$

7: Suffer loss $\hat{\ell}(w^{\prime}_{t},z_{t},\hat{z})$

8: Update model
$w_{t}=w^{\prime}_{t}-\gamma_{t}\nabla\hat{\ell}(w^{\prime}_{t},z_{t},\hat{z})$

9: Update the moving average and the random example,
$\displaystyle\bar{z}\leftarrow\frac{z*(t-1)+z_{t}}{t}$
$\displaystyle\hat{z}\leftarrow z_{t}\text{ if Bernoulli(p) = 1}$

10: end for

11: return $w_{T}$

## Regret Analysis

The regret incurred by the online algorithm relative to the optimal hypothesis
in the original space $\mathcal{H}$, denoted as
$w^{*}=\arg\min_{w\in\mathcal{H}}\sum_{t=2}^{T}L_{t}(w)$, when applied to a
sequence of $T$ examples, is formally expressed as:

$\displaystyle\mathcal{R}_{w^{*},T}=\sum_{t=2}^{T}L_{t}(w_{t-1})-\sum_{t=2}^{T}L_{t}(w^{*}),$
(5)

In the context of this equation, the local all-pairs loss $L_{t}(\cdot)$ is
defined according to Equation 1. To delve deeper, we can decompose the regret
represented in Equation 5 by introducing the best-in-class hypothesis within
the approximated space $\bar{\mathcal{H}}$, denoted as
$\bar{w}^{*}=\arg\min_{w\in\bar{\mathcal{H}}}\sum_{t=2}^{T}L_{t}(w)$. This
decomposition takes the form:

$\displaystyle\mathcal{R}_{w^{*},T}=$
$\displaystyle\underbrace{\sum_{t=2}^{T}L_{t}(w_{t-1})-\sum_{t=2}^{T}L_{t}(\bar{w}^{*})}_{T_{1}}$
$\displaystyle+\underbrace{\sum_{t=2}^{T}L_{t}(\bar{w}^{*})-\sum_{t=2}^{T}L_{t}(w^{*})}_{T_{2}}$
(6)

In our analysis, we provide a bound for $T_{1}$ in Theorem 2, followed by an
error bound for $T_{2}$ in Theorem 3. By combining these bounds, we derive the
regret bound for Algorithm 1 as follows:

######  Theorem 1.

Let $\\{z_{t}\in\mathcal{Z}\\}_{t=1}^{T}$ be sequentially accessed by
Algorithm 1. Let $D$ denote the number of random Fourier features in the
kernel mapping from the original space $\mathcal{X}\subset\mathbb{R}^{d}$. Let
$\eta$ be first step size, $\gamma=O(\Gamma M_{t}\eta)$ the second step size,
and $G$ the Lipschitz constant. Then the regret bound compared to
$w^{*}=\arg\min_{w\in\mathcal{H}}\sum_{t=2}^{T}L_{t}(w)$ is bounded with a
probability of at least $1-2^{8}\left(\frac{\sigma
Diam(\mathcal{X})}{\epsilon}\right)\exp(-D\epsilon^{2}/(4d+8))$, as follows:

$\displaystyle\mathcal{R}_{w^{,}T}$ $\displaystyle\leq
GT\|w^{*}\|_{1}\epsilon+\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+2\gamma G^{2}T$ (7)

where, $\epsilon$ denotes the kernel approximation error, $\sigma$ signifies
the kernel width, $\Gamma\geq tr(cov[z_{t}])$, and $M_{t}:M(w_{t})$ stands for
the smoothness parameter of $\ell$ in the neighborhood of $w_{t}$.

######  Remark 1.

Selecting $\eta=\frac{1}{\sqrt{T}}$ and $\epsilon=\frac{1}{\sqrt{T}}$ results
in a sublinear bound, which is optimal. Notably, $D=O(T)$ in general, yet it
can be as low as $\sqrt{T}\log{T}$ for specific kernels. Furthermore, since
$\gamma=O(\Gamma M_{t}\eta)$, the regret depends on the variance of data and
the smoothness of the loss.

### Regret in the Approximated Space $\bar{\mathcal{H}}$

Prior approaches (e.g., (Kar et al. 2013), (Zhao et al. 2011)) aimed to attain
uniform convergence of buffer loss often relied on the assumption of iid
examples. These methods also require maintaining uniform samples from
historical data through a Bernoulli process with diminishing probabilities,
ultimately approaching zero. Moreover, achieving such uniform convergence
necessitated a disentanglement between the model and the buffer to ensure the
validity of expectation. For instance, while (Zhao et al. 2011) employed
Reservoir sampling to enforce an $iid$ buffer, they overlooked the interplay
between the model and the data. Similarly, (Kar et al. 2013) utilized modified
reservoir sampling and Symmetrization of Expectations to address these
challenges. However, their approach demanded a buffer size of $O(\sqrt{T})$,
which proves impractical for larger-scale scenarios. In contrast, (Yang et al.
2021) accomplished sublinear regret with a buffer size of one, focusing on
local loss. Nevertheless, their assumption of a fully $iid$ training set does
not align with real-world online situations. Our forthcoming theorem obviates
the necessity for the $iid$ assumption by harnessing both moving averages and
individual random points.

As previously explained, in each iteration, a single moving average is
employed to ascertain the gradient direction, which establishes a connection
to the all-pairs loss without relying on the $iid$ assumption. This is then
followed by a step involving the gradient of a random example. Notably, this
sample should not be uniformly chosen from historical data. If we denote the
moving average based loss as $\bar{L}_{t}(w_{t-1})$, then using the convexity
of $\ell(\cdot)$ and the linearity of models in $\bar{\mathcal{H}}$, we have
$\bar{L}_{t}(w_{t-1})\leq L(w_{t-1})$ or
$\bar{L}_{t}(w_{t-1})=L(w_{t-1})-g(t)$ where $g(t)\geq 0$ is Jensen’s gap. The
gap is first bounded in lemma 1 using the smoothness constant and the variance
of the data.

######  Lemma 1.

Given that the data $\\{z_{1},\dots,z_{t}\\}$ is linearly separable and have
bounded total variance, i.e. $\Gamma\geq tr(cov[z_{t}])$ and a loss function
with $M_{t}:M(w_{t})$ smoothness parameters, the Jensen’s gap can be bounded
as follow,

$\displaystyle g(t)\leq\frac{1}{2}\Gamma M_{t}$ (8)

Note that since the model is linear in the space $\bar{\mathcal{H}}$, the
smoothness assumption holds for loss gradient w.r.t. $z$. Therefore we can
provide the following bound in theorem 2.

######  Theorem 2.

Assume that assumptions 1,2, and 3 hold. Let $[w_{t}]_{i=1}^{T}$ be the
sequence of models returned by running the algorithm 1 for $T$ times using the
online sequence of data. Then if
$\bar{w}^{*}=\arg\min_{w\in\bar{\mathcal{H}}}\sum_{t=2}^{T}L_{t}(w)$, we have,

$\displaystyle\sum_{t=2}^{T}L_{t}(w_{t-1})$
$\displaystyle\leq\sum_{t=2}^{T}L_{t}(\bar{w}^{*})+\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+2\gamma
G^{2}T$ (9)

Note that $\Gamma\geq tr(cov[z_{t}])$ where $cov[z_{t}]$ is the covariance of
the received examples, and $M_{t}$ is the smoothness parameter of the loss
function.

### Approximation Error for Random Features

The kernel associated with the space $\mathcal{H}$ is a function defined as
$k:\mathcal{X}^{2}\times\mathcal{X}^{2}\mapsto\mathbb{R}^{+}$ with a shorthand
$k_{(x,x^{\prime})}(\cdot):=k((x,x^{\prime}),(\cdot,\cdot))$ and can be
constructed given any uni-variate kernel $\mathcal{G}$ for any
$x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime}\in\mathcal{X}$ as follow,

$\displaystyle
k(x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime})=\mathcal{G}(x_{1},x_{1}^{\prime})+\mathcal{G}(x_{2},x_{2}^{\prime})$
(10)
$\displaystyle-\mathcal{G}(x_{1},x_{2}^{\prime})-\mathcal{G}(x_{2},x_{1}^{\prime})=\langle\mathcal{G}_{x_{1}}-\mathcal{G}_{x_{2}},\mathcal{G}_{x_{1}^{\prime}}-\mathcal{G}_{x_{2}^{\prime}}\rangle_{\mathcal{G}}$

The pairwise kernel $k$ defined earlier demonstrates positive semi-
definiteness over $\mathcal{X}^{2}$. This property aligns it with the
attributes of a Mercer kernel, as indicated in references like (Ying and Zhou
2015), when $\mathcal{G}$ exhibits such characteristics over $\mathcal{X}$. We
further assume there exist a lower dimensional mapping $r:\mathcal{X}\mapsto
R^{D}$, such that, $\mathcal{G}_{x}(\cdot)\approx r(x)^{T}r(\cdot)$

The literature extensively explores the accuracy of approximating the
pointwise kernel $\mathcal{G}$ using random Fourier features, as discussed in
references such as (Rahimi and Recht 2007), (Bach 2017), and (Li 2022). For
example if the kernel is the Gaussian kernel, i.e.,
$\mathcal{G}(x-x^{\prime})=\frac{1}{\sqrt{2\pi}d}\exp{-(x-x^{\prime})^{T}[diag(\sigma)](x-x^{\prime})}$,
we have $p(u)$ from Fourier transform to be
$p(u)\propto\mathcal{N}(\mathbf{0},diag(\sigma))$ where
$\sigma\in\mathbb{R}^{d}$ is the kernel width. Thus kernel $\mathcal{G}$ can
be given by a Monte Carlo estimate as follow,

$\displaystyle\hat{\mathcal{G}}(x,x^{\prime})\approx\frac{2}{D}\sum_{i=1}^{D/2}\cos(u_{i}^{T}(x-x^{\prime}))=\langle
r(x),r(x^{\prime})\rangle$ (11)

Where $u_{i}\sim\mathcal{N}({0},diag(\sigma))$, and
$r(x):=\frac{1}{\sqrt{D/2}}[cos(u_{i}^{T}x),sin(u_{i}^{T},x)]_{i=1}^{D/2}$.
The quality of the estimation is given by (Rahimi and Recht 2007) to be
dependent on the number of random features $D$ used to approximate the Fourier
transform, which can be extended to pairwise kernels as illustrated in the
following corroallary.

######  Corollary 1.

Given $x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime}\in\mathcal{X}$, and pairwise
kernel $k$ defined on $\mathcal{X}^{2}\times\mathcal{X}^{2}$, the random
Fourier estimation of the kernel has mistake bounded with probability at least
$1-8\exp{\frac{-D^{2}\epsilon}{2}}$ as follow,

$\displaystyle\
|\hat{k}_{(x_{1},x_{2})}(x_{1}^{\prime},x_{2}^{\prime})-{k}_{(x_{1},x_{2})}(x_{1}^{\prime},x_{2}^{\prime})|\leq\epsilon$

The following theorem bounds the random Fourier error in equation (Regret
Analysis) using the results of corollary 1.

######  Theorem 3.

Given a pairwise Mercer kernel
$k_{(x,x^{\prime})}:=k((x,x^{\prime}),(\cdot,\cdot))$ defined on
$\mathcal{X}^{2}\times\mathcal{X}^{2}$. Let $\ell(w,z,z^{\prime})$ be a convex
loss that is Lipschits continuous with constant ${G}$. Then for any
$w^{*}=\sum_{i,j\neq i}^{T}a_{i,j}^{*}k_{(x_{i},x_{j})}$, and random Fourier
features number $D$ we have the following with probability at least
$1-2^{8}\left(\frac{\sigma
Diam(\mathcal{X})}{\epsilon}\right)\exp(-D\epsilon^{2}/(4d+8))$,

$\displaystyle\sum_{t=2}^{T}L_{t}(\bar{w}^{*})-\sum_{t=2}^{T}L_{t}(w^{*})\leq
GT\|w^{*}\|_{1}\epsilon$ (12)

where $\|w^{*}\|_{1}=\sum_{i,j\neq i}^{T}|a^{*}_{i,j}|$.

######  Remark 2.

A trade-off between $\|w^{*}\|_{1}$ and $\epsilon$ exists when independent
support vectors are in the kernel space ($|a^{*}_{i,j}|>0$). In such cases, a
larger error requires a smaller $\epsilon$, implying a larger $D$.
Regularization can set insignificant support vectors’ $a^{*}_{i,j}$ to zero, a
subject for potential future research.

## Related Work

While online gradient descent exhibits a time complexity of $O(T^{2})$
(Boissier et al. 2016; Kar et al. 2013; Gao et al. 2013), this remains
impractical for large-scale problems. For example, in AUC maximization,
pairing a data point $x_{t},y_{t}=1$ with all preceding
${x_{t^{\prime}},y_{t^{\prime}}=-1|1<t^{\prime}<t-1}$ to ascertain the true
loss poses a key challenge. This approach necessitates computing gradients
$\nabla{L}t(f{t-1})$ for all received training examples, growing linearly with
$t$.

Buffer-based methods have been explored in the literature to manage this
challenge. (Zhao et al. 2011) introduced buffers $B_{+}$ and $B_{-}$, each of
size $N_{+}$ and $N_{-}$, considering only points within the buffers at each
step $t$. Reservoir Sampling updates the buffers while ensuring uniform
sampling. Theorem 1 in (Zhao et al. 2011) establishes sublinear regret linked
to buffer sizes. However, this work is limited to AUC with linear models and
does not explicitly address the buffer size’s impact on generalization error.

In contrast, (Kar et al. 2013) proposed modified reservoir sampling (RSx) by
performing $s$ Bernoulli processes with probability $1/t$ to replace each
buffer point with the received data point. This approach ensures $s$ data
points are iid samples from the stream and independent from $w$, with the
generalization bound depending on Rademacher complexity and regret.
Nevertheless, optimal generalization requires $s=O(\sqrt{T})$ for buffer loss.

Another parallel research addressing the high variance in the buffer has been
proposed recently (AlQuabeh, Mukhoty, and Gu 2023), coinciding with the
publication of our paper. Recently, (Yang et al. 2021) suggested optimal
generalization with a buffer size of $s=1$, assuming $iid$ data and
independent buffer-data relations. Their work aligns with ours, but we remove
the $iid$ assumption for convergence proof and consider non-linear model
spaces.

## Experiments

We perform experiments on several real-world datasets and compare our proposed
algorithm to both offline and online pairwise algorithms. Specifically, the
proposed method is compared with different pairwise algorithms on the AUC
maximization task with the squared function as the surrogate loss.

### Experimental Setup

Dataset | AOGD | SPAM-NET | OGD | Sparse Kernel | Projection++ | Kar
---|---|---|---|---|---|---
diabetes | 81.91$\pm$0.48 | 82.03$\pm$0.32 | 82.53$\pm$0.31 | 82.64$\pm$0.37 | 77.92$\pm$1.44 | 79.85$\pm$0.28
ijcnn1 |  92.32$\pm$0.77 | 87.01$\pm$0.10 | 89.46$\pm$1.22 | 71.13$\pm$0.59 | 92.20$\pm$0.27 | 83.44$\pm$1.21
a9a |  90.03$\pm$0.41 | 89.95$\pm$0.42 | 90.01$\pm$0.42 | 84.20$\pm$0.17 | 84.42$\pm$0.33 | 77.93$\pm$1.55
mnist |  96.98$\pm$0.38 | 96.57$\pm$0.54 | 89.56$\pm$0.34 | 95.21$\pm$0.15 | 93.82$\pm$0.15 | 89.16$\pm$0.14
rcv1 |  99.38$\pm$0.48 | 98.13$\pm$0.15 | 99.05$\pm$0.57 | 96.26$\pm$0.35 | 94.54 $\pm$0.36 | 97.78 $\pm$0.64
usps |  92.76$\pm$0.48 | 89.12$\pm$0.88 | 89.88$\pm$0.47 | 91.25$\pm$0.84 | 90.14 $\pm$0.22 | 89.58$\pm$0.25
german |  80.75$\pm$0.48 | 78.55$\pm$2.75 | 79.25$\pm$0.42 | 80.11$\pm$0.44 | 78.44 $\pm$0.66 | 77.24$\pm$0.52
Reg. | $l_{2}$ | $l_{1}+l_{2}$ | $l_{2}$ | $l_{2}$ | $l_{2}$ | $l_{2}$
Table 2: AUC maximization results (average $\pm$ standard error)$\times
10^{2}$ using different batch and online algorithms on different datasets

Figure 1: The AUC vs. time comparison of the algorithms in different datasets
showing superior performance of the proposed method.

#### Compared Algorithms.

The compared algorithms including offline and online settings, as illustrated
below,

  * •

SPAM-NET (Reddi et al. 2016) is an online algorithm for AUC with square loss
that is converted to a saddle point problem with non-smooth regularization.

  * •

OGD (Yang et al. 2021), the most similar to our algorithm but with a linear
model, that uses the last point every iteration.

  * •

Sparse Kernel (Kakkar et al. 2017) is an offline algorithm for AUC
maximization that uses the kernel trick.

  * •

Projection ++ (Hu et al. 2015) is an online algorithm with an adaptive support
vector set.

  * •

Kar (Kar et al. 2013) is online algorithm with buffer of size $s$.

#### Implementation.

The datasets of investigation are available on LIBSVM website (Chang and Lin
2011). The experiments are validated for all algorithms by performing a grid
search on the hyperparameters with three-fold cross-validation. For instance,
in all algorithms the step size $\eta\in 2^{[-8:-1]}$, and the regularization
parameters $\lambda\in 10^{[-8:-1]}$. All the algorithms have been run five
times on different folds on a CPU of 4 GHz speed and 16 GB memory.

### Experimental Results and Analysis

Our results on maximizing the area under the curve (AUC) using a squared loss
function verify the efficacy of our random Fourier pairwise online gradient
descent procedure. Table 2 shows that when compared to online and offline,
linear and nonlinear pairwise learning algorithms, our algorithm achieves
improved AUC performance on large-scale datasets.

Figure 1 displays the AUC performance relative to CPU utilization. Results of
ijcnn1 and mnist show that compared to other state-of-the-art online pairwise
learning methods, our algorithm not only has better AUC performance but also
has a faster convergence rate. As compared to the two most popular kernel
algorithms, Sparse kernel, and Projection++, our technique provides superior
AUC performance across the board. Projection++ is impractical for large-scale
datasets because of its hyper-parameters and the requirement to regularly
update a support vector set.

## Conclusion

In summary, our work introduces a lightweight online kernelized pairwise
learning algorithm applicable to both linear and non-linear models. The
algorithm efficiently computes the online gradient using a moving average and
random examples, with the kernel function approximated through random Fourier
mapping. Achieving a gradient complexity of $O(T)$ for linear models and
$O(\frac{D}{d}T)$ for nonlinear models, our approach showcases computational
efficiency. The incorporation of multiple moving averages, a potential avenue
for future research, could further enhance adaptability and robustness.
Notably, we address the challenges of non-iid data by dual evaluation of
gradients, leveraging random Fourier features for efficient kernel estimation
with a sublinear error bound. This integration of kernelization,
approximation, and dynamic averaging transcends linear constraints,
effectively tackling non-iid data challenges and mitigating kernel complexity.

## References

  * AlQuabeh and Abdurahimov (2022) AlQuabeh, H.; and Abdurahimov, A. 2022.  Pairwise Learning via Stagewise Training in Proximal Setting.  _arXiv preprint arXiv:2208.04075_.
  * AlQuabeh, Mukhoty, and Gu (2023) AlQuabeh, H.; Mukhoty, B.; and Gu, B. 2023.  Variance Reduced Online Gradient Descent for Kernelized Pairwise Learning with Limited Memory.  _arXiv preprint arXiv:2310.06483_.
  * Bach (2017) Bach, F. 2017.  On the equivalence between kernel quadrature rules and random feature expansions.  _The Journal of Machine Learning Research_ , 18(1): 714–751.
  * Boissier et al. (2016) Boissier, M.; Lyu, S.; Ying, Y.; and Zhou, D.-X. 2016.  Fast convergence of online pairwise learning algorithms.  In _Artificial Intelligence and Statistics_ , 204–212. PMLR.
  * Chang and Lin (2011) Chang, C.-C.; and Lin, C.-J. 2011.  LIBSVM: A library for support vector machines.  _ACM Transactions on Intelligent Systems and Technology_ , 2: 27:1–27:27.  Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.
  * Du et al. (2016) Du, C.; Du, C.; Long, G.; He, Q.; and Li, Y. 2016.  Online Bayesian Multiple Kernel Bipartite Ranking.  In _UAI_.
  * Gao et al. (2013) Gao, W.; Jin, R.; Zhu, S.; and Zhou, Z.-H. 2013.  One-pass AUC optimization.  In _International conference on machine learning_ , 906–914. PMLR.
  * Gu, Huo, and Huang (2019) Gu, B.; Huo, Z.; and Huang, H. 2019.  Scalable and efficient pairwise learning to achieve statistical accuracy.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, 3697–3704.
  * Hu et al. (2015) Hu, J.; Yang, H.; King, I.; Lyu, M. R.; and So, A. M.-C. 2015.  Kernelized online imbalanced learning with fixed budgets.  In _Twenty-Ninth AAAI Conference on Artificial Intelligence_.
  * Kakkar et al. (2017) Kakkar, V.; Shevade, S.; Sundararajan, S.; and Garg, D. 2017.  A sparse nonlinear classifier design using AUC optimization.  In _Proceedings of the 2017 SIAM International Conference on Data Mining_ , 291–299. SIAM.
  * Kar et al. (2013) Kar, P.; Sriperumbudur, B.; Jain, P.; and Karnick, H. 2013.  On the Generalization Ability of Online Learning Algorithms for Pairwise Loss Functions.  In Dasgupta, S.; and McAllester, D., eds., _Proceedings of the 30th International Conference on Machine Learning_ , volume 28 of _Proceedings of Machine Learning Research_ , 441–449. Atlanta, Georgia, USA: PMLR.
  * Kulis et al. (2012) Kulis, B.; et al. 2012.  Metric learning: A survey.  _Foundations and trends in machine learning_ , 5(4): 287–364.
  * Li (2022) Li, Z. 2022.  Sharp Analysis of Random Fourier Features in Classification.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, 7444–7452.
  * Lin et al. (2017) Lin, J.; Lei, Y.; Zhang, B.; and Zhou, D.-X. 2017.  Online pairwise learning algorithms with convex loss functions.  _Information Sciences_ , 406: 57–70.
  * Natole, Ying, and Lyu (2018) Natole, M., Jr.; Ying, Y.; and Lyu, S. 2018.  Stochastic Proximal Algorithms for AUC Maximization.  In Dy, J.; and Krause, A., eds., _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , 3710–3719. PMLR.
  * Rahimi and Recht (2007) Rahimi, A.; and Recht, B. 2007.  Random features for large-scale kernel machines.  _Advances in neural information processing systems_ , 20.
  * Reddi et al. (2016) Reddi, S. J.; Hefny, A.; Sra, S.; Poczos, B.; and Smola, A. 2016.  Stochastic variance reduction for nonconvex optimization.  In _International conference on machine learning_ , 314–323. PMLR.
  * Wang et al. (2012) Wang, Y.; Khardon, R.; Pechyony, D.; and Jones, R. 2012.  Generalization bounds for online learning algorithms with pairwise loss functions.  In _Conference on Learning Theory_ , 13–1. JMLR Workshop and Conference Proceedings.
  * Yang et al. (2021) Yang, Z.; Lei, Y.; Wang, P.; Yang, T.; and Ying, Y. 2021.  Simple Stochastic and Online Gradient Descent Algorithms for Pairwise Learning.  _Advances in Neural Information Processing Systems_ , 34.
  * Ying, Wen, and Lyu (2016) Ying, Y.; Wen, L.; and Lyu, S. 2016.  Stochastic online AUC maximization.  _Advances in neural information processing systems_ , 29: 451–459.
  * Ying and Zhou (2015) Ying, Y.; and Zhou, D.-X. 2015.  Online pairwise learning algorithms with kernels.  _arXiv preprint arXiv:1502.07229_.
  * Zhang et al. (2016) Zhang, Z.; Shi, Q.; McAuley, J.; Wei, W.; Zhang, Y.; and Van Den Hengel, A. 2016.  Pairwise matching through max-weight bipartite belief propagation.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 1202–1210.
  * Zhao et al. (2011) Zhao, P.; Hoi, S. C.; Jin, R.; and YANG, T. 2011.  Online AUC maximization.  _Proceedings of the 28th International Conference on Machine Learning ICML 2011:_.

## Appendix

## Proof of Theorem 2

Given the linearity of the model, expressed as $f_{w}(x)=\langle w,x\rangle$,
the convexity of the loss function in the weight variable $w$ implies
convexity in the input variable. In other words, the linearity of the model
suggests that if the loss function is convex with respect to the weight
variable $w$, it will also exhibit convexity concerning the input variable.
This connection further reinforces the relationship between the convexity of
the loss function, linearity of the model, and convexity properties in both
the weight and input variables. Therefore, given the convexity property of the
loss function $\ell$, upon the arrival of a new example $z_{t}$, the loss
based on the average, denoted as $\bar{z}:=(\mathbb{E}[x|y\neq
y_{t}],y)=\left(\sum_{i=1}^{t-1}\frac{x_{i}}{t-1},y_{t-1}|{y_{t-1}\neq
y_{t}}\right)$, can be linked to the local all-pairs loss using Jensen’s
inequality as follow,

$\ell(w,z_{t},\bar{z})\leq\frac{1}{t-1}\sum_{t=1}^{t-1}\ell(w,z_{t},z_{i})$
(13)

Alternatively, representing the average-based loss as $\bar{L}(w)$ and the
cumulative loss as $L(w)=\frac{1}{t-1}\sum_{t=1}^{t-1}\ell(w,z_{t},z_{i})$, we
introduce a slack variable at $w=w_{t-1}$ as follows,

$\bar{L}_{t}(w{t-1})=L(w_{t-1})-g(t)$ (14)

Here, $g(t)\geq 0$ denotes Jensen’s gap. This relation establishes a
mathematical connection between the average-based loss and the all-pairs loss
through the introduction of this slack variable. Additionally, it’s noteworthy
that this gap can be bounded by leveraging the smoothness of the loss function
and considering the variance inherent in the received variables, as
demonstrated in Lemma 2.

######  Lemma 2.

Given that the data $\\{z_{1},\dots,z_{t}\\}$ is linearly separable and have
bounded total variance, i.e. $\Gamma\geq tr(cov[z_{t}])$ and a loss function
with $M_{t}:M(w_{t})$ smoothness parameters at $w=w_{t}$, the Jensen’s gap can
be bounded as follow,

$\displaystyle g(t)\leq\frac{1}{2}\Gamma M_{t}$ (15)

###### Proof.

Given that the loss function $\ell(w_{t-1})$ is $M_{t}$-smooth, and
considering the linearity of the model in the approximated space
$\bar{\mathcal{H}}$, expressed as $f_{w}(x)=\langle w,x\rangle$, it follows
that if the function is smooth in $w$, it is also smooth in $x$ (or $z$ in our
context). Let $Q:={z_{1},\dots,z_{t-1}}$ represent the set of historical
examples. Define $\bar{z}:=(\mathbb{E}[x|y\neq y_{t}],y)$ where the
expectation is taken with respect to the uniform probability over the set $Q$.
Utilizing the smoothness property of $\ell$ from assumption 2 and considering
expectations,

$\displaystyle\ell(w_{t-1},z_{t},z)$
$\displaystyle\leq\ell(w_{t-1},z_{t},\bar{z})+\nabla_{z}\ell(w_{t-1},z_{t},\bar{z})^{T}(z-\bar{z})+\frac{M_{t}}{2}(z-\bar{z})^{T}(z-\bar{z})$
$\displaystyle\mathbb{E}_{z}\ell(w_{t-1},z_{t},z)-\ell(w_{t-1},z_{t},\bar{z})$
$\displaystyle\leq\mathbb{E}_{z}\nabla_{z}\ell(w_{t-1},z_{t},\bar{z})^{T}(z-\bar{z})+\mathbb{E}_{z}\frac{M_{t}}{2}(z-\bar{z})^{T}(z-\bar{z})$
$\displaystyle L(w_{t-1})-\bar{L}(w_{t-1})$
$\displaystyle\leq\mathbb{E}_{z}\nabla_{z}\ell(w_{t-1},z_{t},\bar{z})^{T}(z-\bar{z})+\mathbb{E}_{z}\frac{M_{t}}{2}(z-\bar{z})^{T}(z-\bar{z})$
$\displaystyle=\frac{M_{t}}{2}\mathbb{E}_{z}\|z-\bar{z}\|^{2}=\frac{M_{t}}{2}\mathbb{E}_{z}\sum_{i}^{d}(z_{i}-\bar{z}_{i})^{2}$
$\displaystyle=\frac{M_{t}}{2}\sum_{i}^{d}\mathbb{E}_{z}(z_{i}-\bar{z}_{i})^{2}=\frac{M_{t}}{2}\sum_{i}^{d}Var[z_{i}]=\frac{M_{t}}{2}tr(cov[z])$
(16)

Replace the expression for $g(t)$ in order to derive the outcomes. ∎

It’s important to highlight that $tr(\text{cov}[z])$ denotes the trace of the
covariance matrix related to the variable $z$ (or the feature $x$ of $z$). We
will colloquially refer to this as the total variance. Next we present the
following lemma, which provides an upper bound on the buffer loss (based on
the average), before we establish the regret bounds for the all-pairs loss in
Theorem 2.

######  Lemma 3.

Assume that assumptions 1 and 3 hold. Let $[w_{t}]_{i=1}^{T}$ be the sequence
of functions returned by running the algorithm 1 for $T$ times using online
sequence of data. Then if
$\bar{w}^{*}=\arg\min_{w\in\bar{\mathcal{H}}}\sum_{t=2}^{T}{L}_{t}(w)$, we
have the following:

$\displaystyle\sum_{t=2}^{T}\bar{L}_{t}(w_{t-1})\leq\sum_{t=2}^{T}\bar{L}_{t}(\bar{w}^{*})+\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+\gamma
G^{2}T-\sum_{t=2}^{T}\frac{\gamma}{\eta}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})$
(17)

###### Proof.

Let $\bar{\ell}_{t}(\cdot)=\ell(\cdot,z_{t},\bar{z})$ be the loss based on the
average of history examples, where $\bar{z}:=(\mathbb{E}x|{y\neq y_{t},y})$
and the expectation is w.r.t. uniform probability over previous examples , and
$\bar{v}_{t}\in\partial\bar{\ell}_{t}(w_{t-1})$ then the first update is
$w^{\prime}_{t}:=w_{t-1}-\eta_{t}\bar{v}$. Let
$\hat{\ell}_{t}(\cdot)=\ell(\cdot,z_{t},z)$ where $z$ is a random sample from
$\\{z_{1},\dots,z_{t-1}\\}$ . Let
$\hat{v}_{t}\in\partial\hat{\ell}_{t}(w^{\prime}_{t})$. If we take the
distance of two subsequent models to the optimal model we have,

$\displaystyle\|w_{t}-\bar{w}^{*}\|^{2}-\|w_{t-1}-\bar{w}^{*}\|^{2}=\|w_{t-1}-\eta_{t}(\bar{v}_{t}+\frac{\gamma_{t}}{\eta_{t}}\hat{v}_{t})-\bar{w}^{*}\|^{2}-\|w_{t-1}-\bar{w}^{*}\|^{2}$
$\displaystyle=\|w_{t-1}-\bar{w}^{*}\|^{2}-2\eta_{t}(\bar{v}_{t}+\hat{v}_{t})^{T}(w_{t-1}-\bar{w}^{*})+\eta_{t}^{2}\|\bar{v}_{t}+\hat{v}_{t}\|^{2}-\|w_{t-1}-\bar{w}^{*}\|^{2}$
$\displaystyle=-2\eta_{t}(\bar{v}_{t}+\frac{\gamma_{t}}{\eta_{t}}\hat{v}_{t})^{T}(w_{t-1}-\bar{w}^{*})+\eta_{t}^{2}\|\bar{v}_{t}+\frac{\gamma_{t}}{\eta_{t}}\hat{v}_{t}\|^{2}$
$\displaystyle\overset{(a)}{\leq}-2\gamma_{t}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})-2\eta_{t}(\bar{\ell}_{t}(w_{t-1})-\bar{\ell}_{t}(\bar{w}^{*}))+\eta_{t}^{2}\|\bar{v}_{t}+\frac{\gamma_{t}}{\eta_{t}}\hat{v}_{t}\|^{2}$
$\displaystyle\overset{(b)}{\leq}-2\gamma_{t}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})-2\eta_{t}(\bar{\ell}_{t}(w_{t-1})-\bar{\ell}_{t}(\bar{w}^{*}))+2\gamma_{t}\eta_{t}G^{2}$

Where the inequality (a) implements assumption 3, i.e.
$-v_{t}^{T}(w_{t-1}-\bar{w}^{*})\leq-(g_{t}(w_{t-1})-g_{t}(\bar{w}^{*}))$. And
inequality (b) implements assumption 1, i.e. $\|\nabla\ell(\cdot)\|\leq G$.
Setting the step size $\eta_{t}=\eta$ and $\gamma_{t}=\gamma$ for all $t$, we
can get

$\displaystyle\bar{\ell}(w_{t-1},z_{t},\bar{z})-\bar{\ell}(\bar{w}^{*},z_{t},\bar{z})$
$\displaystyle\quad\leq\frac{\|w_{t-1}-\bar{w}^{*}\|^{2}-\|w_{t}-\bar{w}^{*}\|^{2}}{2\eta}+\gamma
G^{2}-\frac{\gamma}{\eta}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})$ (18)

Finally if we sum from $t=2$ to $t=T$ and set $w_{1}=0$.

$\displaystyle\sum_{t=2}^{T}\ell(w_{t-1},z_{t},\bar{z})-\sum_{t=2}^{T}\ell(\bar{w}^{*},z_{t},\bar{z})$
$\displaystyle\leq\frac{\|\bar{w}^{*}\|^{2}-\|w_{T}-\bar{w}^{*}\|^{2}}{2\eta}+\sum_{t=2}^{T}\gamma
G^{2}-\sum_{t=2}^{T}\frac{\gamma}{\eta}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})$
$\displaystyle\leq\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+\gamma
G^{2}T-\sum_{t=2}^{T}\frac{\gamma}{\eta}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})$
(19)

Substitute the definition of $\bar{L}_{t}(\cdot)=\ell(\cdot,z_{t},\bar{z})$ to
have the results. This completes the proof. ∎

Now we are ready to prove theorem 2 using the previous two lemmas.

######  Theorem 4 (Theorem 2 restated).

Assume that assumptions 1,2, and 3 hold. Let $[w_{t}]_{i=1}^{T}$ be the
sequence of models returned by running the algorithm for $T$ times using
online sequence of data. Then if
$\bar{w}^{*}=\arg\min_{w\in\bar{\mathcal{H}}}\sum_{t=2}^{T}{L}_{t}(w)$, we
have the following,

$\displaystyle\sum_{t=2}^{T}L_{t}(w_{t-1})$
$\displaystyle\leq\sum_{t=2}^{T}L_{t}(\bar{w}^{*})+\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+2\gamma
G^{2}T$ (20)

###### Proof.

Starting from lemma 2 final’s inequality,

$\displaystyle\sum_{t=2}^{T}\bar{L}_{t}(w_{t-1})$
$\displaystyle\underset{}{\leq}\sum_{t=2}^{T}\bar{L}_{t}(\bar{w}^{*})+\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+\gamma
G^{2}T-\sum_{t=2}^{T}\frac{\gamma}{\eta}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})$
$\displaystyle\overset{(a)}{\leq}\sum_{t=2}^{T}{L}_{t}(\bar{w}^{*})+\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+\gamma
G^{2}T-\sum_{t=2}^{T}\frac{\gamma}{\eta}\hat{v}_{t}^{T}(w_{t-1}-\bar{w}^{*})$
$\displaystyle\sum_{t=2}^{T}{L}_{t}(w_{t-1})-\sum_{t=2}^{T}{L}_{t}(\bar{w}^{*})$
$\displaystyle\overset{(b)}{\leq}\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+\gamma
G^{2}T+\sum_{t=2}^{T}\frac{M_{t}}{2}\Gamma-\frac{\gamma}{\eta}\nabla\ell(w_{t}^{\prime},z_{t},\hat{z})_{t}^{T}(w_{t-1}-\bar{w}^{*})$
$\displaystyle\overset{(c)}{=}\frac{\|\bar{w}^{*}\|^{2}}{2\eta}$
$\displaystyle+\gamma
G^{2}T+\sum_{t=2}^{T}\frac{M_{t}}{2}\Gamma-\frac{\gamma}{\eta}\nabla\ell(w_{t}^{\prime},z_{t},\hat{z})^{T}(w^{\prime}_{t}-\bar{w}^{*}+\eta\nabla\ell(w_{t-1},z_{t},\bar{z}))$
$\displaystyle\overset{(d)}{\leq}\frac{\|\bar{w}^{*}\|^{2}}{2\eta}$
$\displaystyle+\gamma
G^{2}T+\sum_{t=2}^{T}\frac{M_{t}}{2}\Gamma-\frac{\gamma}{2\eta
M}\|\nabla\ell(w_{t}^{\prime},z_{t},\hat{z})\|^{2}$
$\displaystyle-\gamma\nabla\ell(w_{t}^{\prime},z_{t},\hat{z})^{T}\nabla\ell(w_{t-1},z_{t},\bar{z})$
$\displaystyle\overset{(e)}{\leq}\frac{\|\bar{w}^{*}\|^{2}}{2\eta}+2\gamma
G^{2}T+\sum_{t=2}^{T}\frac{M_{t}}{2}\Gamma-\frac{\gamma}{2\eta
M}\|\nabla\ell(w_{t}^{\prime},z_{t},\hat{z})\|^{2}$ (21)

where inequality (a) implements the convexity of the loss, inequality (b)
implements lemma 2, equality (c) implements update
$w^{\prime}_{t}+\eta\nabla\ell(w_{t-1},z_{t},\bar{z})=w_{t-1}$, inequality (d)
implements the smoothness of the loss, and (e) is the cauchy inequality.
Finally choosing the second step size to be $\gamma=O(\Gamma M_{t}\eta)$ makes
the last two terms term cancel out, for any random example $\hat{z}$. Note
that $\Gamma\geq tr(cov[z_{t}])$ where $cov[z_{t}]$ is the covariance of the
received examples, and $M_{t}$ is the smoothness parameter of the loss
function. ∎

## Proof of Theorem 3

######  Theorem 5 (Theorem 3 restated).

Given a pairwise Mercer kernel
$k_{(x,x^{\prime})}:=k((x,x^{\prime}),(\cdot,\cdot))$ defined on
$\mathcal{X}^{2}\times\mathcal{X}^{2}$. Let $\ell(w,z,z^{\prime})$ be convex
loss that is Lipschits continuous with constant $G$. Then for any
$w^{*}=\sum_{i,j\neq i}^{T}a_{i,j}^{*}k_{(x_{i},x_{j})}$, and random Fourier
features number $D$ we have the following with probability at least
$1-2^{8}\left(\frac{\sigma
Diam(\mathcal{X})}{\epsilon}\right)\exp(-D\epsilon^{2}/(4d+8))$,

$\displaystyle\sum_{t=2}^{T}L_{t}(\bar{w}^{*})-\sum_{t=2}^{T}L_{t}(w^{*})\leq
GT\|w^{*}\|_{1}\epsilon$ (22)

where $\|w^{*}\|_{1}=\sum_{i,j\neq i}^{T}|a^{*}_{i,j}|$, and
$D=\Omega\left(\frac{d}{\epsilon^{2}}\log\frac{\sigma
diam(\mathcal{X})}{\epsilon}\right)$.

The proof depends on the main claim of (Rahimi and Recht 2007), we rewrite the
results here.

######  Lemma 4 (Claim 1 (Rahimi and Recht 2007)).

Given $x,x^{\prime}\in\mathcal{X}$, and kernel $\mathcal{G}$ defined on
$\mathcal{X\times X}$, we have the following bound on the error of the
approximation $\hat{\mathcal{G}}$ using $D$ random Fourier features.

$\displaystyle\mathbb{P}\left(|\hat{\mathcal{G}}(x,x^{\prime})-\mathcal{G}(x,x^{\prime})|\geq\epsilon\right)\leq
2\exp\left(\frac{-D\epsilon^{2}}{2}\right)$ (23)

The pairwise kernel defined as,

$\displaystyle
k(x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime})=\mathcal{G}(x_{1},x_{1}^{\prime})+\mathcal{G}(x_{2},x_{2}^{\prime})$
(24)
$\displaystyle-\mathcal{G}(x_{1},x_{2}^{\prime})-\mathcal{G}(x_{2},x_{1}^{\prime})=\langle\mathcal{G}_{x_{1}}-\mathcal{G}_{x_{2}},\mathcal{G}_{x_{1}^{\prime}}-\mathcal{G}_{x_{2}^{\prime}}\rangle_{\mathcal{G}}$

The following Corollary bounds the error of the pairwise kernel defined in
equation (24) using main theorem in (Rahimi and Recht 2007).

######  Corollary 2.

Given $x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime}\in\mathcal{X}$, and pairwise
kernel $k$ defined on $\mathcal{X}^{2}\times\mathcal{X}^{2}$, the random
Fourier estimation of the kernel has an error bounded with probability at
least $1-8\exp{\frac{-D^{2}\epsilon}{2}}$ as follow,

$\displaystyle\sup_{x_{i},x_{t}\in\mathcal{X}}|\hat{k}_{(x_{1},x_{2})}(x_{1}^{\prime},x_{2}^{\prime})-{k}_{(x_{1},x_{2})}(x_{1}^{\prime},x_{2}^{\prime})|\leq\epsilon$

Now we are ready to prove theorem 3, where we start by difference of the loss
between the model $\bar{w}^{*}$ and $w^{*}$,

###### Proof.

$\displaystyle\sum_{t=2}^{T}L_{t}(\bar{w}^{*})-\sum_{t=2}^{T}L_{t}(w^{*})$
$\displaystyle=\sum_{t=2}^{T}\frac{1}{t-1}\sum_{i=1}^{t-1}\ell(\bar{w}^{*},z_{t},z_{i})-\sum_{t=2}^{T}\frac{1}{t-1}\sum_{i=1}^{t-1}\ell({w}^{*},z_{t},z_{i})$
$\displaystyle\overset{a}{\leq}\sum_{t=2}^{T}\frac{1}{t-1}\sum_{i=1}^{t-1}|\ell(\bar{w}^{*},z_{t},z_{i})-\ell({w}^{*},z_{t},z_{i})|$
$\displaystyle\overset{b}{\leq}\sum_{t=2}^{T}\frac{1}{t-1}\sum_{i=1}^{t-1}G\|\bar{w}^{*}-{w}^{*}\|_{L^{2}(\rho)}$
(25)

where inequality (a) implements the triangle inequality, and inequality (b)
applies assumption 1 and $\|\cdot\|_{L^{2}(\rho)}\geq\|\cdot\|_{2}$.
Using the Representer theorem and assumption 4, we assume that the space
$\mathcal{H}$ and $\bar{\mathcal{H}}$ are dense in $L^{2}(\rho)$ (the space of
square integrable function under the probability distribution of the received
examples), then we can approximate any function in $\mathcal{H}$ by a function
in $L^{2}(\rho)$, i.e. without loss of generality we assume that
$\bar{w}^{*}=\sum_{j=1,k\neq j}^{t-1}\alpha_{j,k}^{*}\hat{k}_{z_{t},z_{i}}$,
then we have,

$\|\bar{w}^{*}-{w}^{*}\|_{L^{2}(\rho)}=\|\sum_{j=1,k\neq
j}^{t-1}\alpha_{j,k}^{*}(\hat{k}_{z_{t},z_{i}}-{k}_{z_{t},z_{i}})\|_{L^{2}(\rho)}$
(26)

and using the triangle inequality we have,

$\displaystyle\sum_{t=2}^{T}L_{t}(\bar{w}^{*})-\sum_{t=2}^{T}L_{t}(w^{*})$
$\displaystyle{\leq}\sum_{t=2}^{T}\sup_{x_{i},x_{t}\in\mathcal{X}}G\sum_{j=1,k\neq
j}^{t-1}\|\alpha_{j,k}^{*}(\hat{k}_{z_{t},z_{i}}-{k}_{z_{t},z_{i}})\|_{L^{2}(\rho)}$
$\displaystyle\overset{c}{\leq}\sum_{t=2}^{T}G\sum_{j=1,k\neq
j}^{t-1}|\alpha_{j,k}^{*}|\epsilon\leq GT\|w^{*}\|_{1}\epsilon$ (27)

where inequality (c) implements corollary 2, and $\|w^{*}\|_{1}=\sum_{i,j\neq
i}^{T}|a^{*}_{i,j}|$.
This completes the proof. ∎


# Leptoquark-induced CLFV decays with a light SM-singlet scalar

Bibhabasu De<EMAIL_ADDRESS>Department of Physics, The ICFAI University
Tripura, Kamalghat-799210, India

###### Abstract

The Standard Model (SM), if augmented with a light SM-singlet scalar $\phi$
and a TeV-scale scalar leptoquark (LQ) $S_{1}$, 2-body charged lepton flavor
violating (CLFV) decay channels can be accessed with $\phi$ as one of the
final states, where the leading order effective interactions between $\phi$
and the SM fields arise at one-loop level. Further, in the presence of
$S_{1}$, $\phi$ can mediate 3-body CLFV processes with either two photons or
two gluons in the final states. Thus, the model predicts an exotic 3-body CLFV
channel: $\ell_{A}\to\ell_{B}gg$, which can be tested/constrained only through
some future high-energy experiments looking for di-gluon signals from a
leptonic decay.

##  I Introduction

At the perturbative regime, Standard Model possesses an accidental symmetry
represented by the global gauge group $\mathcal{G}_{LB}\supset
U(1)_{B+L}\times U(1)_{B-L}\times U(1)_{L_{\mu}-L_{\tau}}\times
U(1)_{L_{\mu}+L_{\tau}-2L_{e}}$ Heeck:2016xwg which ensures the conservation
of baryon number $B$ and the individual lepton numbers $L_{A}$
($A=e,\,\mu,\,\tau$). $L=\sum_{A=e,\,\mu,\,\tau}L_{A}$ denotes the total
lepton number. However, at the non-perturbative scale, the $B+L$ part of the
symmetry group is broken by 6 units tHooft:1976rip . At the same time, the
neutrino oscillations strongly establish that the lepton flavor can’t be a
protected symmetry, resulting in a broken $U(1)_{L_{\mu}-L_{\tau}}\times
U(1)_{L_{\mu}+L_{\tau}-2L_{e}}$ in the neutrino sector. The observation
definitely motivates a search for similar flavor violations in the charged
sector. However, for the charged leptons, $\mathcal{G}_{LB}$ is still a
conserved symmetry as ensured by the Glashow-Iliopoulos-Maiani (GIM) mechanism
Glashow:1970gm — any CLFV process induced by the non-zero neutrino masses and
a non-trivial lepton mixing matrix will be too suppressed to detect. For
example, the CLFV decay $\ell_{A}\to\ell_{B}\gamma$ induced through the Dirac
neutrinos at one-loop level results in
BR$(\ell_{A}\to\ell_{B}\gamma)<10^{-53}$ for all the possible channels
Petcov:1976ff . Therefore, any observation of charged lepton flavor violation
can be marked as a significant signature of Beyond Standard Model (BSM)
contribution, as well as physics beyond neutrino oscillations.

Experimental searches for CLFV processes have been going on for a long time
Hincks:1948vr ; PhysRev.74.1364 , with the current sensitivity reaching
$\sim\mathcal{O}(10^{-13})$ for a few flavor violating channels MEG:2016leq ;
SINDRUM:1987nra ; SINDRUMII:2006dvw and being significantly improved for the
others. However, in the absence of any positive signal, interest has grown in
searching for CLFV decays involving low-energy BSM particles. For example,
Belle-II presents the best current upper limit on $\tau\to e\phi$ and
$\tau\to\mu\phi$ at 95% C.L. Belle-II:2022heu , where $\phi$ is an invisible
spin-0 boson in the mass range $0-1.6$ GeV. Similar bounds for $\mu\to e\phi$
are given by Refs. TWIST:2014ymv ; Jodidio:1986mz . However, if $\phi$ decays
to visible particles inside the detector, more stringent constraints could be
placed from the non-observation of 3-body CLFV processes.

Various models have been proposed to address the current bounds on
$\ell_{A}\to\ell_{B}\phi$, with $\phi$ being either a light
scalar/pseudoscalar PhysRevLett.48.11 ; Wilczek:1982rv ; Grinstein:1985rt ;
Berezhiani:1989fp ; Feng:1997tn ; Hirsch:2009ee ; Jaeckel:2013uva ;
Celis:2014iua ; Celis:2014jua ; Galon:2016bka ; Calibbi:2016hwq ; Ema:2016ops
; Bjorkeroth:2018dzu ; Bauer:2019gfk ; Heeck:2019guh ; Cornella:2019uxs ;
Calibbi:2020jvd ; Escribano:2020wua or a gauge boson Farzan:2015hkd ;
Heeck:2016xkh ; Farzan:2016wym ; Ibarra:2021xyk , associated with some
spontaneously broken symmetry. This paper considers $\phi$ to be a generic
real SM-singlet scalar where the flavor-violating interactions arise at one-
loop level in the presence of a scalar LQ $S_{1}$. LQs (for a review, see Ref.
Dorsner:2016wpm ) appear naturally in the grand unified theories (GUT)
Pati:1974yy ; Georgi:1974sy ; Georgi:1974my and work as an excellent BSM
candidate to induce flavor violation. Though in the simplest GUT models LQs
acquire mass at a scale $\Lambda_{\rm GUT}\gg\mathcal{O}(1)$ TeV, there exist
GUT formulations that can ensure the stability of proton with a TeV-scale
scalar LQ Dorsner:2004jj ; Aydemir:2019ynb ; Dorsner:2024seb . Thus, the
present paper, while extending the SM with a TeV-scale scalar LQ to describe
the interactions between $\phi$ and the SM fields, effectively portrays a UV-
complete theory. Ref. Mandal:2019gff has already elaborated on various CLFV
observables arising in the presence of a TeV-scale scalar LQ and used the
experimental upper limits to constrain the New Physics (NP) couplings.
However, interesting consequences may appear if one considers a light SM-
singlet scalar within the $S_{1}$-extension of SM. For example, with
$M_{\phi}<(m_{\mu}-m_{e})$, $\phi$ can be produced through the on-shell decays
of $\tau$ and $\mu$ leading to the CLFV processes $\tau\to\ell_{B}\phi$
[$\ell_{B}=e,\,\mu$] and $\mu\to e\phi$, respectively. In practice, the term
light will be used in this paper to refer to $M_{\phi}<2m_{e}$, i.e., a sub-
MeV gauge-singlet scalar. In this mass regime, it can be easily ensured that
the di-lepton and di-quark decay modes of $\phi$ De:2024tbo are kinematically
forbidden. However, $\phi$ can always decay to two photons at one-loop level,
leading to 3-body CLFV channels of the form $\ell_{A}\to\ell_{B}\gamma\gamma$.
Thus, the experimental bounds for $\ell_{A}\to\ell_{B}\gamma\gamma$ can be
used to check the consistency of the parameter space allowed through the
$\ell_{A}\to\ell_{B}\phi$ searches. Further, LQs being colored particles,
within this proposed BSM formulation $\phi$ can also decay to two gluons. This
results in a rarely discussed unconstrained CLFV decay process:
$\ell_{A}\to\ell_{B}gg$. Indeed, this is a significant phenomenological
outcome of this model as it predicts a set of unexplored CLFV observables.
Therefore, any future search aiming for $\ell_{A}\to\ell_{B}gg$ will be
crucial to test/falsify the model.

The rest of the paper is arranged as follows. Sec. II describes the NP
interactions at the TeV scale, which connect $\phi$ with the SM fields. The
analytical results for $\ell_{A}\to\ell_{B}\phi$ have been established in Sec.
III. Sec. IV presents a parameter space consistent with all the existing CLFV
constraints that can appear in this chosen framework and, hence, numerically
analyze the detection prospects of $\ell_{A}\to\ell_{B}\phi$ within the
considered mass regime. The possible decay modes of $\phi$ and the
corresponding $\phi$-mediated 3-body CLFV processes have been discussed in
Sec. V. Finally, the work has been concluded in Sec. VI.

##  II The Model: A Simple Extension of the SM

The model considers a GUT-motivated minimal extension of the SM where NP
interactions arise at a scale $\Lambda_{\rm NP}\sim\mathcal{O}(1)$ TeV. The
augmented particle spectrum includes a scalar LQ $S_{1}$ which transforms as
($\bar{\mathbf{3}}$, 1, 1/3) under the SM gauge group $\mathcal{G}_{\rm
SM}\supset SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$. Further, a real SM-
singlet scalar $\phi$ can be proposed to exist in the sub-MeV mass regime such
that its decays to the SM fields are either kinematically forbidden or highly
suppressed. Defining the electromagnetic (EM) charge as $Q_{\rm EM}=T_{3}+Y$,
the complete particle spectrum of the present framework has been listed in
Table 1.

Fields | Generations | $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$
---|---|---
$L_{L}=(\nu_{L}\quad\ell_{L})^{T}$ | 3 |  (1, 2, -1/2)
$\ell_{R}=(e_{R},\,\mu_{R},\,\tau_{R})$ | 3 |  (1, 1, -1)
$Q_{L}=(u_{L}\quad d_{L})^{T}$ | 3 |  (3, 2, 1/6)
$U_{R}=(u_{R},\,c_{R},\,t_{R})$ | 3 |  (3, 1, 2/3)
$D_{R}=(d_{R},\,s_{R},\,b_{R})$ | 3 |  (3, 1, -1/3)
$H=(H^{+}\quad H^{0})^{T}$ | 1 |  (1, 2, 1/2)
$S_{1}$ | 1 |  ($\bar{\mathbf{3}}$, 1, 1/3)
$\phi$ | 1 |  (1, 1, 0)
Table 1: Fields and their transformations under $\mathcal{G}_{\rm SM}$.

The NP interactions can be described through

$\displaystyle\mathcal{L}_{\rm NP}=~{}$
$\displaystyle(\mathcal{D}^{\mu}S_{1})^{\dagger}(\mathcal{D}_{\mu}S_{1})+\frac{1}{2}(\partial^{\mu}\phi)(\partial_{\mu}\phi)$
$\displaystyle-\Big{[}Y_{L}^{ij}(\bar{Q}_{L}^{Cia}\epsilon^{ab}L_{L}^{jb})S_{1}+Y_{R}^{ij}(\bar{U}_{R}^{Ci}\ell_{R}^{j})S_{1}+{\rm
h.c.}\Big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\mathbb{V}(H,\,S_{1},\,\phi),$
(1)

where $Y_{L}$ and $Y_{R}$ are completely arbitrary $3\times 3$ Yukawa matrices
in the flavor basis. The superscript $C$ defines the charge-conjugated states;
$\\{a,b\\}$ and $\\{i,j\\}$ stand for the $SU(2)_{L}$ and flavor indices,
respectively. $\epsilon^{ab}=(i\sigma_{2})^{ab}$, $\sigma_{k}$ ($k=1,\,2,\,3$)
being the Pauli matrices. The color indices have been suppressed for
simplicity. Note that, the gauge interactions of $S_{1}$ can easily be
obtained with an explicit definition of the covariant derivative
$\mathcal{D}^{\mu}$. These interactions can be useful to produce $\phi$ at the
colliders Bhaskar:2020kdr ; De:2024tbo .

The freedom to rotate equal-isospin fermion fields in the flavor basis,
permits one to assume the charged lepton and down-type quark Yukawas to be
diagonal so that the transformation from flavor to mass basis is given by
$u_{L}\to\left(\mathbf{V}_{\rm CKM}^{\dagger}\right)u_{L}$ and
$\nu_{L}\to\left(\mathbf{U}_{\rm PMNS}\right)\nu_{L}$. Here $\mathbf{V}_{\rm
CKM}$ and $\mathbf{U}_{\rm PMNS}$ represent the CKM and PMNS matrices,
respectively. Therefore, in the physical basis,

$\displaystyle\mathcal{L}_{\rm NP}$
$\displaystyle=(\mathcal{D}^{\mu}S_{1})^{\dagger}(\mathcal{D}_{\mu}S_{1})+\frac{1}{2}(\partial^{\mu}\phi)(\partial_{\mu}\phi)$
$\displaystyle-\Big{[}\left\\{\bar{u}_{L}^{Ci}\left(\mathbf{V}^{*}_{\rm
CKM}Y_{L}\right)^{ij}\ell_{L}^{j}\right\\}S_{1}-\left\\{\bar{d}_{L}^{Ci}\left(Y_{L}\mathbf{U}_{\rm
PMNS}\right)^{ij}\nu_{L}^{j}\right\\}S_{1}$
$\displaystyle\qquad\qquad+Y_{R}^{ij}(\bar{U}_{R}^{Ci}\ell_{R}^{j})S_{1}+{\rm
h.c.}\Bigg{]}-\mathbb{V}(H,\,S_{1},\,\phi)$
$\displaystyle=(\mathcal{D}^{\mu}S_{1})^{\dagger}(\mathcal{D}_{\mu}S_{1})+\frac{1}{2}(\partial^{\mu}\phi)(\partial_{\mu}\phi)-\Big{[}\left(\bar{u}_{L}^{Ci}\xi_{L}^{ij}\ell_{L}^{j}\right)S_{1}$
$\displaystyle-\left(\bar{d}_{L}^{Ci}\chi^{ij}\nu_{L}^{j}\right)S_{1}+\xi_{R}^{ij}(\bar{U}_{R}^{Ci}\ell_{R}^{j})S_{1}+{\rm
h.c.}\Bigg{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\mathbb{V}(H,\,S_{1},\,\phi),$
(2)

where, $\xi_{L}^{ij}=\left(\mathbf{V}^{*}_{\rm CKM}Y_{L}\right)^{ij}$,
$\chi^{ij}=\left(Y_{L}\mathbf{U}_{\rm PMNS}\right)^{ij}$, and
$\xi_{R}^{ij}=Y_{R}^{ij}$. However, the neutrino sector doesn’t play any role
in the present analysis. The scalar potential can be cast as,

$\displaystyle\mathbb{V}$
$\displaystyle(H,\,S_{1},\,\phi)=\lambda_{HS_{1}}\left(H^{\dagger}H\right)\left(S^{\dagger}_{1}S_{1}\right)+\beta\phi\left(S^{\dagger}_{1}S_{1}\right)+\delta\phi(H^{\dagger}H)$
$\displaystyle+\lambda_{H\phi}(H^{\dagger}H)\phi^{2}+\lambda_{\phi
S_{1}}\phi^{2}\left(S^{\dagger}_{1}S_{1}\right)+\frac{1}{2}M^{2}_{\phi}\phi^{2}+M^{2}_{S_{1}}\left(S^{\dagger}_{1}S_{1}\right).$
(3)

However, one can assume $\lambda_{H\phi}$ and $\delta$ to be negligible so
that the presence of $\phi$ doesn’t affect the SM scalar sector. Further, by
setting $\lambda_{HS_{1}}\to 0$, one can ensure that the trilinear Higgs
coupling assumes its SM-predicted value and doesn’t introduce any additional
constraint to the parameter space. Note that, $\beta$ is a mass dimensional
coupling which should, in principle, represent the highest scale of the theory
and hence must be $\sim\mathcal{O}(1)$ TeV. After electroweak symmetry
breaking (EWSB), the physical masses can be obtained via,

$\displaystyle H=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\
v+h\end{pmatrix},\quad\quad\phi=\phi,$ (4)

where $h$ denotes the physical Higgs field, and $v=246$ GeV defines the
electroweak vacuum expectation value (VEV). However, in the limit
$\lambda_{HS_{1}}\to 0$ and $\lambda_{H\phi}\to 0$, $M_{S_{1}}$ and $M_{\phi}$
represent the physical masses of $S_{1}$ and $\phi$, respectively.

##  III Lepton Flavor Violation: $\boldsymbol{\Large\ell_{A}\to\ell_{B}+\phi}$

The proposed model, which extends the SM with a scalar LQ $S_{1}$, possesses a
rich BSM phenomenology in the lepton sector. For example, the discrepancies in
muon and electron $(g-2)$, semi-leptonic decays, non-observation of CLFV
processes including $\ell_{A}\to\ell_{B}\gamma$,
$\ell_{A}\to\bar{\ell}_{B}\ell_{B}\ell_{B}$, $h\to\bar{\ell}_{A}\ell_{B}$, and
$\mu-e$ conversion can be explained with $M_{S_{1}}\sim\mathcal{O}(1)$ TeV.
However, in the presence of an additional light SM-singlet scalar, new CLFV
channels can be introduced within the considered framework where the SM
leptons ($\tau$ and $\mu$) can decay to $\phi$ and a lighter lepton through
the effective vertex shown in Fig. 1. To be specific, for $M_{\phi}\leq 1$
MeV, $\tau\to e\phi$, $\tau\to\mu\phi$, and $\mu\to e\phi$ CLFV modes are
kinematically accessible.

Figure 1: Leading order contribution to $\ell_{A}\to\ell_{B}\phi$. Here, $k$
denotes the loop momentum, and $p_{A}=-(p_{B}+p_{\phi})$ defines the momentum
conservation.

The effective coupling corresponding to Fig. 1 is given by,

$\displaystyle-i$ $\displaystyle
Y_{\ell_{A}\ell_{B}\phi}=\sum_{i=\,u,c,t}(\xi_{L}^{Ai})^{*}\xi_{R}^{iB}\beta
m_{i}\int\frac{d^{4}k}{(2\pi)^{4}}\Bigg{[}\frac{(k^{2}+m_{i}^{2})}{(k^{2}-m_{i}^{2})^{2}}\times$
$\displaystyle\qquad\qquad\frac{1}{\\{(k+p_{B})^{2}-M_{S_{1}}^{2}\\}\\{(k+p_{B}+p_{\phi})^{2}-M_{S_{1}}^{2}\\}}\Bigg{]}$
$\displaystyle=$
$\displaystyle\sum_{i=\,u,c,t}(\xi_{L}^{Ai})^{*}\xi_{R}^{iB}\beta
m_{i}\int\frac{d^{4}k}{(2\pi)^{4}}\Bigg{[}$
$\displaystyle\frac{1}{(k^{2}-m_{i}^{2})\\{(k+p_{B})^{2}-M_{S_{1}}^{2}\\}\\{(k+p_{B}+p_{\phi})^{2}-M_{S_{1}}^{2}\\}}$
$\displaystyle+\frac{2m_{i}^{2}}{(k^{2}-m_{i}^{2})^{2}\\{(k+p_{B})^{2}-M_{S_{1}}^{2}\\}\\{(k+p_{B}+p_{\phi})^{2}-M_{S_{1}}^{2}\\}}\Bigg{]}$
$\displaystyle=$
$\displaystyle\left(\frac{i\,\beta}{16\pi^{2}}\right)\sum_{i=\,u,c,t}\Omega^{AB}_{i}\,m_{i}\,\Bigg{[}\mathcal{J}_{1}+2m_{i}^{2}\times\mathcal{J}_{2}\Bigg{]},$
(5)

where $\Omega^{AB}_{i}=(\xi_{L}^{Ai})^{*}\xi_{R}^{iB}$, and

$\displaystyle\mathcal{J}_{1}=$ $\displaystyle\
-C_{0}\left(m_{B}^{2},\,M_{\phi}^{2},\,m_{A}^{2},\,m_{i}^{2},\,M_{S_{1}}^{2},\,M_{S_{1}}^{2}\right),$
$\displaystyle\mathcal{J}_{2}=$ $\displaystyle\
D_{0}\left(0,\,m_{B}^{2},\,M_{\phi}^{2},\,m_{A}^{2},\,m_{B}^{2},\,m_{A}^{2},\,m_{i}^{2},\,m_{i}^{2},\,M_{S_{1}}^{2},\,M_{S_{1}}^{2}\right).$
(6)

Eq. (5) has been obtained assuming an on-shell decay of $\ell_{A}$ so that
$p_{A}^{2}=(p_{B}+p_{\phi})^{2}=m_{A}^{2}$. Here, $C_{0}$ and $D_{0}$ denote
the standard Passarino-Veltman functions Romao:2020 ; Hahn:1998yk for the
scalar 3-point and 4-point one-loop integrals, respectively.

In the rest-frame of $\ell_{A}$, the decay width for $\ell_{A}\to\ell_{B}\phi$
can be calculated as,

$\displaystyle\Gamma(\ell_{A}\to\ell_{B}\phi)=$
$\displaystyle~{}\frac{\mathcal{N}_{C}\times|Y_{\ell_{A}\ell_{B}\phi}|^{2}}{16\pi}\times
m_{A}\times\Big{[}(1+\rho_{B})^{2}-\rho_{\phi}^{2}\Big{]}$
$\displaystyle\times\Bigg{[}\Big{\\{}1-(\rho_{\phi}+\rho_{B})^{2}\Big{\\}}\Big{\\{}1-(\rho_{\phi}-\rho_{B})^{2}\Big{\\}}\Bigg{]}^{1/2}\,,$
(7)

where $\rho_{B}=m_{B}/m_{A}$, and $\rho_{\phi}=M_{\phi}/m_{A}$.
$\mathcal{N}_{C}=3$ stands for the color degeneracy factor. However, to detect
any signal from $\ell_{A}\to\ell_{B}\phi$ decay,
$\ell_{A}\to\ell_{B}\bar{\nu}_{B}\nu_{A}$ acts as the strongest background.
Therefore, to parametrize the detectability of $\ell_{A}\to\ell_{B}\phi$ over
$\ell_{A}\to\ell_{B}\bar{\nu}_{B}\nu_{A}$, a variable $R_{AB}$ can be defined
as,

$\displaystyle R_{AB}=\frac{{\rm BR}(\ell_{A}\to\ell_{B}\phi)}{{\rm
BR}(\ell_{A}\to\ell_{B}\bar{\nu}_{B}\nu_{A})},$ (8)

where BR symbolizes the branching ratio for a particular decay mode. In the
subsequent analysis, $R_{AB}$ will be considered as a good observable to
constrain the parameter space.

##  IV Numerical Analysis

As stated in Sec. III, the $S_{1}$-extension of SM can accommodate several BSM
observations and possibilities resulting in strong bounds on the NP couplings
$\xi^{ij}_{L,R}$. Thus, to analyze the $\ell_{A}\to\ell_{B}\phi$ processes,
one must choose a parameter space consistent with all the existing leptonic
constraints. Following Ref. Mandal:2019gff , Table 2 enlists the most
stringent upper limits on $|\Omega^{AB}_{i}|^{2}$ relevant for the present
study.

CLFV | $|\Omega^{AB}_{i}|^{2}$ | Experimental | Considered
---|---|---|---
Processes |  |  Bounds $\times\left(M_{S_{1}}/{\rm TeV}\right)^{4}$ | Values
| $|\Omega^{\tau e}_{u}|^{2}$ |  $<\,3.0\times 10^{-6}$ ParticleDataGroup:2018ovx | $10^{-6}$
$\tau\to e\phi$ | $|\Omega^{\tau e}_{c}|^{2}$ |  $<\,2.1\times 10^{-4}$ BaBar:2009hkt | $10^{-6}$
| $|\Omega^{\tau e}_{t}|^{2}$ |  $<\,3.8\times 10^{-7}$ BaBar:2009hkt | $[10^{-12},\,10^{-6}]$
| $|\Omega^{\tau\mu}_{u}|^{2}$ |  $<\,5.7\times 10^{-6}$ ParticleDataGroup:2018ovx | $10^{-6}$
$\tau\to\mu\phi$ | $|\Omega^{\tau\mu}_{c}|^{2}$ |  $<\,2.8\times 10^{-4}$ BaBar:2009hkt | $10^{-6}$
| $|\Omega^{\tau\mu}_{t}|^{2}$ |  $<\,5.0\times 10^{-7}$ BaBar:2009hkt | $[10^{-12},\,10^{-6}]$
| $|\Omega^{\mu e}_{u}|^{2}$ |  $<\,6.7\times 10^{-12}$ SINDRUMII:2006dvw | $10^{-12}$
$\mu\to e\phi$ | $|\Omega^{\mu e}_{c}|^{2}$ |  $<\,1.7\times 10^{-12}$ MEG:2016leq | $[10^{-16},\,10^{-12}]$
| $|\Omega^{\mu e}_{t}|^{2}$ |  $<\,3.0\times 10^{-15}$ MEG:2016leq | $10^{-17}$
Table 2: Experimental constraints on the NP couplings and their considered
values used for the computations.

The first column of Table 2 specifies a particular $\ell_{A}\to\ell_{B}\phi$
process, while the third column represents the experimental bounds on the
associated NP couplings arising through the other CLFV processes. However, the
bounds drastically deteriorate with increasing LQ mass. In general, the
couplings considered for the present analysis assume $M_{S_{1}}\geq 2$ TeV
ParticleDataGroup:2022pth and have been listed in the last column of Table 2.
As an illustrative choice, $\beta=5$ TeV will be used for the analysis.

###  IV.1 $\boldsymbol{\tau}$-Sector

The $\tau e\phi$ and $\tau\mu\phi$ effective couplings can be obtained from
Eq. (5) with $m_{A}=m_{\tau}$ and $m_{B}=m_{e},\,m_{\mu}$, respectively. Fig.
2 shows their variations with increasing LQ mass. All the NP couplings
$\Omega^{\tau B}_{i}$ have been fixed at their respective upper limits, i.e.,
$10^{-3}$, while $M_{\phi}$ has been randomly varied between 1 keV and 1 MeV
to generate Fig. 2.

Figure 2: Variation of $|Y_{\tau e\phi}|$ (red) and $|Y_{\tau\mu\phi}|$
(green) as a function of $M_{S_{1}}$ with $M_{\phi}\in[1~{}{\rm keV},1~{}{\rm
MeV}]$.

However, the results suggest that the variation of $M_{\phi}$ in the
considered range negligibly affects $|Y_{\tau\ell_{B}\phi}|$
($\ell_{B}=e,\,\mu$). Further, the lighter lepton masses (i.e., $m_{e}$ and
$m_{\mu}$) being ignorable with respect to $\Lambda_{\rm NP}$, $|Y_{\tau
e\phi}|$ and $|Y_{\tau\mu\phi}|$ overlap.

((a))  (a)

((b))  (b)

Figure 3: Variation of (a) $R_{\tau e}$, and (b) $R_{\tau\mu}$ as a function
of $M_{\phi}$ for $M_{S_{1}}=2$ TeV (red), 3 TeV (sky), and 4 TeV (yellow).
The solid black lines denote the observed upper limits (at 95% C.L.) from
Belle II collaboration Belle-II:2022heu .

Fig. 3 (a) and 3 (b) display the variation of $R_{\tau e}$ and $R_{\tau\mu}$
as a function of $M_{\phi}$, respectively. For Fig. 3 (a), $\Omega^{\tau
e}_{t}$ has been randomly varied within the range $[10^{-6},\,10^{-3}]$,
whereas one has to vary $\Omega^{\tau\mu}_{t}$ to generate the scattered plot
in Fig. 3 (b). Currently, the most stringent bounds on $R_{\tau\ell_{B}}$ are
set by Belle II collaboration Belle-II:2022heu as shown by the solid black
lines in Fig. 3. The results indicate a significant parameter space where the
non-observation of $\tau\to e\phi$ and $\tau\to\mu\phi$ can be explained along
with the other CLFV processes. Moreover, the allowed parameter space improves
with a heavier $S_{1}$.

###  IV.2 $\boldsymbol{\mu}$-Sector

$\phi$ being in the sub-MeV mass regime, $\mu$ can also decay to $\phi$ and
$e$ through the effective $Y_{\mu e\phi}$ coupling. Fig. 4 depicts its
variation as a function of $M_{S_{1}}$. As before, all the NP couplings have
been fixed at their allowed upper bounds with $M_{\phi}$ varying randomly
within the range 1 keV to 1 MeV.

Figure 4: Variation of $|Y_{\mu e\phi}|$ as a function of $M_{S_{1}}$ with
$M_{\phi}\in[1~{}{\rm keV},1~{}{\rm MeV}]$.

Note that, $|Y_{\mu e\phi}/Y_{\tau\ell_{B}\phi}|\sim 10^{-5}$. The large
suppression for $Y_{\mu e\phi}$ stems from the NP couplings $\Omega^{\mu
e}_{i}$ ($i=u,\,c,\,t$), which are tightly constrained through the
experimental bounds on $\mu\to e\gamma$ MEG:2016leq and $\mu-e$ conversion on
gold nuclei SINDRUMII:2006dvw .

Figure 5: Variation of $R_{\mu e}$ as a function of $M_{\phi}$ for
$M_{S_{1}}=2$ TeV (red), 3 TeV (sky), and 4 TeV (yellow). The black solid line
marks the observed upper limit (at 90% C.L.) from Jodidio et al.
Jodidio:1986mz .

Fig. 5 presents the $M_{\phi}$ dependence of $R_{\mu e}$ for three different
LQ masses. The plot has been generated through a random variation of
$\Omega^{\mu e}_{c}$ within the range $[10^{-8},\,10^{-6}]$, whereas
$\Omega^{\mu e}_{u}$ and $\Omega^{\mu e}_{t}$ have been fixed at $10^{-6}$ and
$3.16\times 10^{-9}$, respectively.

TWIST collaboration sets the most stringent bound on $R_{\mu e}$, when $\phi$
is a massive invisible boson in the range $13-80$ MeV TWIST:2014ymv . However,
in the sub-MeV mass regime, $\phi$ can be treated as effectively massless, and
hence, the best experimental upper limit can be read as $R_{\mu e}<2.6\times
10^{-6}$ (Jodidio et al. Jodidio:1986mz ). Thus, in Fig. 5, the portion below
the black solid line defines the allowed parameter space for $\mu\to e\phi$.

##  V Decay of $\boldsymbol{\phi}$ and Other CLFV Constraints

As indicated in the Introduction, Refs. Jodidio:1986mz ; TWIST:2014ymv ;
Belle-II:2022heu set the best constraints for the considered 2-body CLFV
decays if $\phi$ is an invisible light boson, i.e., either $\phi$ decays to
invisible particles (e.g., dark matter, neutrinos) or features a significantly
large decay-length. In the former case, only a sharp energy peak of electron
or muon can be searched over the continuous Michel spectrum. However, if
$\phi$ decays to visible particles, there is always a chance to have better
constraints from the 3-body CLFV processes.

In the proposed framework, the gauge-singlet scalar $\phi$ produced through
the on-shell decays of $\tau$ and/or $\mu$ can further decay via di-photon and
di-gluon channels with the leading order contributions arising at one-loop
level [see Fig. 6]. Note that, for $M_{\phi}<2m_{e}$, other decay modes of
$\phi$ are kinematically forbidden.

Figure 6: Leading order contributions to $\phi\to\gamma\gamma\,(gg)$
processes.

Thus, the total decay width of $\phi$ can be formulated as,

$\displaystyle\Gamma_{\phi}=\Gamma(\phi\to gg)+\Gamma(\phi\to\gamma\gamma),$
(9)

where the partial decay widths are given by Dorsner:2016wpm ; Bhaskar:2020kdr
,

$\displaystyle\Gamma(\phi\rightarrow gg)$ $\displaystyle\ =\ \frac{G_{\rm
F}\alpha_{S}^{2}M_{\phi}^{3}}{64\sqrt{2}\pi^{3}}\left|\frac{\beta
v}{2M_{S_{1}}^{2}}\mathcal{P}\left(\frac{M_{\phi}^{2}}{4M_{S_{1}}^{2}}\right)\right|^{2},$
$\displaystyle\Gamma(\phi\rightarrow\gamma\gamma)$ $\displaystyle\ =\
\frac{G_{\rm F}\alpha_{\rm
EM}^{2}M_{\phi}^{3}}{128\sqrt{2}\pi^{3}}\left|\frac{\beta
v}{6M_{S_{1}}^{2}}\mathcal{P}\left(\frac{M_{\phi}^{2}}{4M_{S_{1}}^{2}}\right)\right|^{2}.$
(10)

Here, $G_{\rm F}$ is the Fermi constant, $\alpha_{S}$ and $\alpha_{\rm EM}$
define the strong and electromagnetic coupling constants, respectively. The
function $\mathcal{P}$ can be defined as,

$\displaystyle\mathcal{P}(x)=-\frac{[x-\psi(x)]}{x^{2}}\,,$ (11)

where,

$\displaystyle\psi(x)$ $\displaystyle=\Bigg{\\{}\begin{array}[]{cc}{\rm
Arcsin}^{2}(\sqrt{x}),&x\leq 1\\\ -\frac{1}{4}\left[{\rm
ln}\left(\frac{1+\sqrt{1-x^{-1}}}{1-\sqrt{1-x^{-1}}}\right)-i\pi\right]^{2},&x>1\,.\end{array}$
(14)

Thus, the proposed framework opens up 3-body CLFV decay routes of the form
$\ell_{A}\to\ell_{B}\gamma\gamma$ and $\ell_{A}\to\ell_{B}gg$, which can
possibly be used to check the validity of the available parameter space.
However, the former decay channel, where the heavier lepton decays to a
lighter one and two photons, is well-known in the literature and has also been
explored through experiments, whereas the latter is mostly unconstrained and
rarely discussed in theory.

$\Gamma_{\phi}$ being extremely suppressed
$\left(\leq\mathcal{O}(10^{-22})~{}{\rm GeV}\right)$, one can use the narrow
width approximation to study the $\phi$-mediated 3-body CLFV decays Cordero-
Cid:2005vca . Therefore, the corresponding branching ratios can be cast as,

$\displaystyle{\rm BR}(\ell_{A}\to\ell_{B}XX)\simeq{\rm
BR}(\ell_{A}\to\ell_{B}\phi)\times{\rm BR}(\phi\to XX),$ (15)

where $X$ represents either photon or gluon.

###  V.1 $\boldsymbol{\ell_{A}\to\ell_{B}\gamma\gamma}$

((a))  (a)

((b))  (b)

((c))  (c)

Figure 7: Variation of (a) BR$(\tau\to e\gamma\gamma)$, (b)
BR$(\tau\to\mu\gamma\gamma)$, and (c) BR$(\mu\to e\gamma\gamma)$ as a function
of $M_{\phi}$ for $M_{S_{1}}=2$ TeV (red), 3 TeV (sky), and 4 TeV (yellow).
The black solid lines represent the corresponding experimental sensitivities
(at 90% C.L.).

This type of 3-body CLFV processes has been searched for a long time
(particularly for $\mu\to e\gamma\gamma$) through various experiments, and
till now, no flavor violation has been observed. Thus, the null results have
set certain upper limits on BR($\ell_{A}\to\ell_{B}\gamma\gamma$), which may
further propagate into the $(\ell_{A}\to\ell_{B}\phi)$-allowed parameter space
as more stringent constraints. Table 3 presents the current experimental
bounds on BR($\ell_{A}\to\ell_{B}\gamma\gamma$).

CLFV Observables | Upper Limits | Experiments
---|---|---
BR($\mu\to e\gamma\gamma$) | $7.2\times 10^{-11}$ |  Crystal Box Bolton:1988af
BR($\tau\to e\gamma\gamma$) | $2.5\times 10^{-4}$ |  Bryman et al. Bryman:2021ilc
BR($\tau\to\mu\gamma\gamma$) | $1.5\times 10^{-4}$ |  ATLAS Angelozzi:2017oeg
Table 3: Experimental bounds on $\ell_{A}\to\ell_{B}\gamma\gamma$ processes.

Clearly, the experimental bounds for BR$(\tau\to\ell_{B}\gamma\gamma)$
[$\ell_{B}=e,\,\mu$] are less stringent compared to those on
$R_{\tau\ell_{B}}$ from Belle-II. It can be easily understood from Fig. 7 (a)
and 7 (b) that the 3-body decay doesn’t add any additional constraint in the
$\tau$-sector, and effectively the entire parameter space is allowed with the
existing upper limits. Therefore, Fig. 3 alone displays the actual permitted
region where the production prospects of $\phi$ via $\tau$-decays can be
studied. Fig. 7 (c) depicts the $M_{\phi}$ dependence of BR$(\mu\to
e\gamma\gamma)$ for different values of $M_{S_{1}}$. The currently available
experimental sensitivity restricts BR$(\mu\to e\gamma\gamma)<7.2\times
10^{-11}$ Bolton:1988af and, thus, introduces the strongest bound on the
considered parameter space. With all the $(\mu\to e\phi)$-specific NP
couplings (i.e., $|\Omega^{\mu e}_{i}|$) being fixed at their allowed
values/ranges as listed in Table 2, the results show that a sub-MeV SM-singlet
scalar can only be produced through the flavor violating muon decay if
$M_{S_{1}}>2$ TeV.

###  V.2 $\boldsymbol{\ell_{A}\to\ell_{B}gg}$

((a))  (a)

((b))  (b)

((c))  (c)

Figure 8: Variation of (a) BR$(\tau\to e\,gg)$, (b) BR$(\tau\to\mu\,gg)$, and
(c) BR$(\mu\to e\,gg)$ as a function of $M_{\phi}$ for $M_{S_{1}}=2$ TeV
(red), 3 TeV (sky), and 4 TeV (yellow).

The proposed extension of SM can accommodate an exotic 3-body CLFV channel
where the heavier charged leptons can decay to the lighter ones and two
gluons. To the best of the author’s knowledge, very little has been explored
about this particular type of CLFV decays, as they can never be probed through
low-energy experiments. Colliders using deep inelastic scattering can only be
an option to search for $\ell_{A}\to\ell_{B}gg$ AbdulKhalek:2022hcn . For
example, Refs. Takeuchi:2017btl ; Cirigliano:2021img have discussed the
gluonic contribution to the $\tau$-specific CLFV processes in an Electron-
Ion/Proton Collider.

Therefore, these decays might lead to completely new observables for probing
the lepton flavor violation in the charged sector. Fig. 8 shows the model
predictions for $\tau\to egg$, $\tau\to\mu gg$, and $\mu\to egg$. Note that,
${\rm BR}(\ell_{A}\to\ell_{B}gg)/{\rm
BR}(\ell_{A}\to\ell_{B}\gamma\gamma)\sim\mathcal{O}(10^{3})$ as $\phi$
dominantly decays through the di-gluon channel. Therefore, in principle,
$\ell_{A}\to\ell_{B}gg$ should have a better detection prospect compared to
$\ell_{A}\to\ell_{B}\gamma\gamma$ if the background is substantially reduced
ATL-PHYS-PUB-2017-017 ; CMS-DP-2017-027 ; Andrews:2019faz . Though, currently,
there is no direct experimental bound on these processes, detectors looking
for di-gluon signals from CLFV decays may test/falsify the present predictions
in the future.

##  VI Conclusion

By the next decade, experiments searching for CLFV processes are going to play
a vital role in discovering physics beyond the SM. The chances will definitely
be enhanced with the increasing number of CLFV observables. Moreover, the
observation will also be significant in understanding the possible NP theory
associated with such a positive signal. This paper has extended the SM with a
light SM-singlet scalar $\phi$ with mass $M_{\phi}<2m_{e}$ and a TeV-scale
scalar LQ $S_{1}$. There is a wide range of theories where such a generic
gauge-singlet scalar can be associated with a spontaneously broken global
symmetry. The effective couplings of $\phi$ with the SM fields can be obtained
at one-loop level in the presence of $S_{1}$. Phenomenologically, scalar LQs
are well-motivated hypothetical particles that generate a large variety of
CLFV processes. Thus, the non-observation of any CLFV process with the
existing experimental sensitivities results in stringent constraints on the NP
couplings at the quark-LQ-lepton vertices. Considering such a pre-constrained
parameter space, the present paper has studied the discovery prospects of
$\tau\to e\phi$, $\tau\to\mu\phi$, and $\mu\to e\phi$. Though for
illustrations, 1 keV $\leq M_{\phi}\leq$ 1 MeV has been chosen, the results
are equally valid for any generic scalar lighter than $2m_{e}$. $S_{1}$ being
an EM-charged color-triplet scalar, $\phi$ can decay through di-photon and di-
gluon channels. This induces the possibility of 3-body CLFV decays
$\ell_{A}\to\ell_{B}\gamma\gamma$ and $\ell_{A}\to\ell_{B}gg$, where $\phi$
acts as the mediator. Though for the $\tau$-sector, the experimental bounds on
$\tau\to\ell_{B}\gamma\gamma$ [$\ell_{B}=e,\,\mu$] are less stringent compared
to those on the $\tau\to\ell_{B}\phi$, for $\mu$-sector BR$(\mu\to
e\gamma\gamma)$ becomes crucial to constrain the parameter space. The most
remarkable feature of this model is the possibility to predict CLFV
observables corresponding to the 3-body decay $\ell_{A}\to\ell_{B}gg$, for all
the possible channels. No direct experimental limit exists for such leptonic
decays with two gluons in the final states. Using narrow-width approximation,
the model has predicted for BR$(\tau\to egg)$, BR$(\tau\to\mu gg)$, and
BR$(\mu\to egg)$ within the allowed parameter space. Therefore, a future
experiment searching for $\ell_{A}\to\ell_{B}gg$ will be significant to
test/constrain the model or may lead to a discovery.

## References

  * (1) J. Heeck, _Interpretation of Lepton Flavor Violation_ , _Phys. Rev. D_ 95 (2017) 015022, [1610.07623].
  * (2) G. ’t Hooft, _Symmetry Breaking Through Bell-Jackiw Anomalies_ , _Phys. Rev. Lett._ 37 (1976) 8–11.
  * (3) S. L. Glashow, J. Iliopoulos and L. Maiani, _Weak Interactions with Lepton-Hadron Symmetry_ , _Phys. Rev. D_ 2 (1970) 1285–1292.
  * (4) S. T. Petcov, _The Processes $\mu\rightarrow e+\gamma,\mu\rightarrow e+\overline{e},\nu^{\prime}\rightarrow\nu+\gamma$ in the Weinberg-Salam Model with Neutrino Mixing_, _Sov. J. Nucl. Phys._ 25 (1977) 340. [Erratum: Sov.J.Nucl.Phys. 25, 698 (1977), Erratum: Yad.Fiz. 25, 1336 (1977)].
  * (5) E. P. Hincks and B. Pontecorvo, _Search for gamma-radiation in the 2.2-microsecond meson decay process_ , _Phys. Rev._ 73 (1948) 257–258.
  * (6) R. D. Sard and E. J. Althaus, _A search for delayed photons from stopped sea level cosmic-ray mesons_ , _Phys. Rev._ 74 (Nov, 1948) 1364–1371.
  * (7) MEG collaboration, A. M. Baldini et al., _Search for the lepton flavour violating decay $\mu^{+}\rightarrow\mathrm{e}^{+}\gamma$ with the full dataset of the MEG experiment_, _Eur. Phys. J. C_ 76 (2016) 434, [1605.05081].
  * (8) SINDRUM collaboration, U. Bellgardt et al., _Search for the Decay $\mu^{+}\to e^{+}e^{+}e^{-}$_, _Nucl. Phys. B_ 299 (1988) 1–6.
  * (9) SINDRUM II collaboration, W. H. Bertl et al., _A Search for muon to electron conversion in muonic gold_ , _Eur. Phys. J. C_ 47 (2006) 337–346.
  * (10) Belle-II collaboration, I. Adachi et al., _Search for Lepton-Flavor-Violating $\tau$ Decays to a Lepton and an Invisible Boson at Belle II_, _Phys. Rev. Lett._ 130 (2023) 181803, [2212.03634].
  * (11) TWIST collaboration, R. Bayes et al., _Search for two body muon decay signals_ , _Phys. Rev. D_ 91 (2015) 052020, [1409.0638].
  * (12) A. Jodidio et al., _Search for Right-Handed Currents in Muon Decay_ , _Phys. Rev. D_ 34 (1986) 1967. [Erratum: Phys.Rev.D 37, 237 (1988)].
  * (13) A. Davidson and K. C. Wali, _Minimal flavor unification via multigenerational peccei-quinn symmetry_ , _Phys. Rev. Lett._ 48 (Jan, 1982) 11–14.
  * (14) F. Wilczek, _Axions and Family Symmetry Breaking_ , _Phys. Rev. Lett._ 49 (1982) 1549–1552.
  * (15) B. Grinstein, J. Preskill and M. B. Wise, _Neutrino Masses and Family Symmetry_ , _Phys. Lett. B_ 159 (1985) 57–61.
  * (16) Z. G. Berezhiani and M. Y. Khlopov, _Cosmology of Spontaneously Broken Gauge Family Symmetry_ , _Z. Phys. C_ 49 (1991) 73–78.
  * (17) J. L. Feng, T. Moroi, H. Murayama and E. Schnapka, _Third generation familons, b factories, and neutrino cosmology_ , _Phys. Rev. D_ 57 (1998) 5875–5892, [hep-ph/9709411].
  * (18) M. Hirsch, A. Vicente, J. Meyer and W. Porod, _Majoron emission in muon and tau decays revisited_ , _Phys. Rev. D_ 79 (2009) 055023, [0902.0525]. [Erratum: Phys.Rev.D 79, 079901 (2009)].
  * (19) J. Jaeckel, _A Family of WISPy Dark Matter Candidates_ , _Phys. Lett. B_ 732 (2014) 1–7, [1311.0880].
  * (20) A. Celis, J. Fuentes-Martin and H. Serodio, _An invisible axion model with controlled FCNCs at tree level_ , _Phys. Lett. B_ 741 (2015) 117–123, [1410.6217].
  * (21) A. Celis, J. Fuentes-Martín and H. Serôdio, _A class of invisible axion models with FCNCs at tree level_ , _JHEP_ 12 (2014) 167, [1410.6218].
  * (22) I. Galon, A. Kwa and P. Tanedo, _Lepton-Flavor Violating Mediators_ , _JHEP_ 03 (2017) 064, [1610.08060].
  * (23) L. Calibbi, F. Goertz, D. Redigolo, R. Ziegler and J. Zupan, _Minimal axion model from flavor_ , _Phys. Rev. D_ 95 (2017) 095009, [1612.08040].
  * (24) Y. Ema, K. Hamaguchi, T. Moroi and K. Nakayama, _Flaxion: a minimal extension to solve puzzles in the standard model_ , _JHEP_ 01 (2017) 096, [1612.05492].
  * (25) F. Björkeroth, E. J. Chun and S. F. King, _Flavourful Axion Phenomenology_ , _JHEP_ 08 (2018) 117, [1806.00660].
  * (26) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _Axionlike Particles, Lepton-Flavor Violation, and a New Explanation of $a_{\mu}$ and $a_{e}$_, _Phys. Rev. Lett._ 124 (2020) 211803, [1908.00008].
  * (27) J. Heeck and H. H. Patel, _Majoron at two loops_ , _Phys. Rev. D_ 100 (2019) 095015, [1909.02029].
  * (28) C. Cornella, P. Paradisi and O. Sumensari, _Hunting for ALPs with Lepton Flavor Violation_ , _JHEP_ 01 (2020) 158, [1911.06279].
  * (29) L. Calibbi, D. Redigolo, R. Ziegler and J. Zupan, _Looking forward to lepton-flavor-violating ALPs_ , _JHEP_ 09 (2021) 173, [2006.04795].
  * (30) P. Escribano and A. Vicente, _Ultralight scalars in leptonic observables_ , _JHEP_ 03 (2021) 240, [2008.01099].
  * (31) Y. Farzan and I. M. Shoemaker, _Lepton Flavor Violating Non-Standard Interactions via Light Mediators_ , _JHEP_ 07 (2016) 033, [1512.09147].
  * (32) J. Heeck, _Lepton flavor violation with light vector bosons_ , _Phys. Lett. B_ 758 (2016) 101–105, [1602.03810].
  * (33) Y. Farzan and J. Heeck, _Neutrinophilic nonstandard interactions_ , _Phys. Rev. D_ 94 (2016) 053010, [1607.07616].
  * (34) A. Ibarra, M. Marín and P. Roig, _Flavor violating muon decay into an electron and a light gauge boson_ , _Phys. Lett. B_ 827 (2022) 136933, [2110.03737].
  * (35) I. Doršner, S. Fajfer, A. Greljo, J. F. Kamenik and N. Košnik, _Physics of leptoquarks in precision experiments and at particle colliders_ , _Phys. Rept._ 641 (2016) 1–68, [1603.04993].
  * (36) J. C. Pati and A. Salam, _Lepton Number as the Fourth Color_ , _Phys. Rev. D_ 10 (1974) 275–289. [Erratum: Phys.Rev.D 11, 703–703 (1975)].
  * (37) H. Georgi and S. L. Glashow, _Unity of All Elementary Particle Forces_ , _Phys. Rev. Lett._ 32 (1974) 438–441.
  * (38) H. Georgi, _The State of the Art—Gauge Theories_ , _AIP Conf. Proc._ 23 (1975) 575–582.
  * (39) I. Dorsner and P. Fileviez Perez, _Could we rotate proton decay away?_ , _Phys. Lett. B_ 606 (2005) 367–370, [hep-ph/0409190].
  * (40) U. Aydemir, T. Mandal and S. Mitra, _Addressing the ${\mathbf{R}_{D^{(*)}}}$ anomalies with an ${\mathbf{S}_{1}}$ leptoquark from $\mathbf{SO(10)}$ grand unification_, _Phys. Rev. D_ 101 (2020) 015011, [1902.08108].
  * (41) I. Doršner and S. Saad, _Is Doublet-Triplet Splitting Necessary?_ , 2404.09021.
  * (42) R. Mandal and A. Pich, _Constraints on scalar leptoquarks from lepton and kaon physics_ , _JHEP_ 12 (2019) 089, [1908.11155].
  * (43) B. De, _Direct production of SM-singlet scalars at the muon collider_ , _Phys. Lett. B_ 852 (2024) 138634, [2401.08101].
  * (44) A. Bhaskar, D. Das, B. De and S. Mitra, _Enhancing scalar productions with leptoquarks at the LHC_ , _Phys. Rev. D_ 102 (2020) 035002, [2002.12571].
  * (45) J. C. Rom$\tilde{a}$o, “Advanced Quantum Field Theory.” https://porthos.tecnico.ulisboa.pt/Public/textos/tca.pdf, 2020.
  * (46) T. Hahn and M. Perez-Victoria, _Automatized one loop calculations in four-dimensions and D-dimensions_ , _Comput. Phys. Commun._ 118 (1999) 153–165, [hep-ph/9807565].
  * (47) Particle Data Group collaboration, M. Tanabashi et al., _Review of Particle Physics_ , _Phys. Rev. D_ 98 (2018) 030001.
  * (48) BaBar collaboration, B. Aubert et al., _Searches for Lepton Flavor Violation in the Decays $\tau^{\pm}\to e^{\pm}\gamma$ and $\tau^{\pm}\to\mu^{\pm}\gamma$_, _Phys. Rev. Lett._ 104 (2010) 021802, [0908.2381].
  * (49) Particle Data Group collaboration, R. L. Workman et al., _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (50) A. Cordero-Cid, G. Tavares-Velasco and J. J. Toscano, _Implications of a very light pseudoscalar boson on lepton flavor violation_ , _Phys. Rev. D_ 72 (2005) 117701, [hep-ph/0511331].
  * (51) R. D. Bolton et al., _Search for Rare Muon Decays with the Crystal Box Detector_ , _Phys. Rev. D_ 38 (1988) 2077.
  * (52) D. A. Bryman, S. Ito and R. Shrock, _Upper limits on branching ratios of the lepton-flavor-violating decays $\tau\to\ell\gamma\gamma$ and $\tau\to\ell X$_, _Phys. Rev. D_ 104 (2021) 075032, [2106.02451].
  * (53) I. Angelozzi, _In pursuit of lepton flavour violation : A search for the $\tau\to\mu\gamma\gamma$ decay with ATLAS at $\sqrt{s}$ =8 TeV_.  PhD thesis, U. Amsterdam, IHEF, 2017.
  * (54) R. Abdul Khalek et al., _Snowmass 2021 White Paper: Electron Ion Collider for High Energy Physics_ , 2203.13199.
  * (55) M. Takeuchi, Y. Uesaka and M. Yamanaka, _Higgs mediated CLFV processes $\mu$N ( eN )→ $\tau$X via gluon operators_, _Phys. Lett. B_ 772 (2017) 279–282, [1705.01059].
  * (56) V. Cirigliano, K. Fuyuto, C. Lee, E. Mereghetti and B. Yan, _Charged Lepton Flavor Violation at the EIC_ , _JHEP_ 03 (2021) 256, [2102.06176].
  * (57) ATLAS collaboration, “Quark versus Gluon Jet Tagging Using Jet Images with the ATLAS Detector.” https://cds.cern.ch/record/2275641, 2017.
  * (58) CMS collaboration, “New Developments for Jet Substructure Reconstruction in CMS.” https://cds.cern.ch/record/2275226, 2017.
  * (59) M. Andrews, J. Alison, S. An, P. Bryant, B. Burkle, S. Gleyzer et al., _End-to-end jet classification of quarks and gluons with the CMS Open Data_ , _Nucl. Instrum. Meth. A_ 977 (2020) 164304, [1902.08276].