#  $(g-2)_{\mu}$ in the 2HDM and slightly beyond – an updated view

P. M. Ferreira  Instituto Superior de Engenharia de Lisboa, Instituto
Politécnico de Lisboa 1959-007 Lisboa, Portugal  Centro de Física Teórica e
Computacional, Faculdade de Ciências, Universidade de Lisboa, Campo Grande,
Edifício C8 1749-016 Lisboa, Portugal  B. L. Gonçalves  Departamento de Física
and CFTP, Instituto Superior Técnico, Universidade de Lisboa, Lisboa, Portugal
Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade
de Lisboa, Campo Grande, Edifício C8 1749-016 Lisboa, Portugal  F. R. Joaquim
Departamento de Física and CFTP, Instituto Superior Técnico, Universidade de
Lisboa, Lisboa, Portugal  Marc Sher  High Energy Theory Group, William & Mary,
Williamsburg, VA 23187, USA

###### Abstract

The recent measurement of the muon $g-2$ anomaly continues to defy a Standard
Model explanation, but can be accommodated within the framework of two Higgs
doublet models. If one further includes extra fermion content in the form of a
generation of vector-like leptons, the allowed parameter range that explains
the anomaly is even further extended, and thorny issues of perturbativity and
clashes with $B$-decay constraints avoided. We will show how the muon magnetic
moment anomaly can be fit within these models, under the assumption that the
vector-like leptons do not mix with the muon. We update previous analyses and
include all theoretical and experimental constraints, including searches for
extra scalars. It is shown that the inclusion of vector-like fermions allows
the lepton-specific and muon-specific models to perform much better in fitting
the muon’s $g-2$. These fits do require the Yukawa coupling between the Higgs
and the vector-like leptons to be large, causing potential problems with
perturbativity and unitarity, and thus models in which the vector-like leptons
mix with the muon may be preferred.

##  I Introduction

Recently, the Muon g-2 Collaboration at Fermilab reported new results results
from Run 1 of their experiment measuring the anomalous magnetic moment of the
muon $a_{\mu}$. Prior to this announcement, the discrepancy between the
experimental measurement $a_{\mu}^{\rm exp}$ Bennett:2006fi  and the Standard
Model (SM) theoretical prediction $a_{\mu}^{\rm SM}$ Blum:2018mom ;
Keshavarzi:2018mgv ; Davier:2019can ; Aoyama:2020ynm  was

$\Delta a_{\mu}^{\rm exp}=a_{\mu}^{\rm exp}-a_{\mu}^{\rm SM}=(279\pm 76)\times
10^{-11}\quad(3.7\sigma),$ (1)

while the new combined result is results

$\Delta a_{\mu}^{\rm exp}=(251\pm 59)\times 10^{-11}\qquad(4.2\sigma).$ (2)

There are hundreds of papers with new-physics explanations for the
$(g-2)_{\mu}$ anomaly such as supersymmetric models, left-right symmetric
models, scotogenic models, 331 models, $L_{\mu}-L_{\tau}$ models, seesaw
models, the Zee-Babu model as well as two-Higgs doublet models (2HDMs); an
extensive review can be found in Ref. Lindner:2016bgg . In this paper, we will
focus on 2HDMs (for a review, see Ref. Branco:2011iw ), and discuss the
implications of the new result from the Muon g-2 Collaboration.

In the 2HDM it is possible to ensure that tree-level flavor changing neutral
currents mediated by scalars do not exist, by imposing a discrete $Z_{2}$
symmetry on the model. There are four such versions of the 2HDM, referred to
as the type-I, type-II, type-X (sometimes called lepton-specific) and type-Y
(sometimes called flipped) models. In the type-II and type-X models, the
coupling of the muon to the heavy Higgs bosons is enhanced by a factor of
$\tan\beta$, which is the ratio of the two vacuum expectation values, but in
type-II models, the $Q=-1/3$ quark couplings also get the same enhancement,
leading to possible problems with radiative $B$-meson decays. In fact,
explaining the $(g-2)_{\mu}$ anomaly without affecting $B$-decays is one of
the motivations for studies of the type-X model. On the other hand, in the
type-I and type-Y models, the couplings of the muon to heavy Higgs bosons are
suppressed by $\tan\beta$ and, thus, these models are not favored for
explaining the $(g-2)_{\mu}$ discrepancy. For simplicity, we will only
consider models without tree-level flavor-changing neutral currents and with
CP conservation in the Higgs sector; for a discussion of the effects of CP
violation in the Higgs sector on $(g-2)_{\mu}$, see Ref. Keus:2017ioh .

One of the earlier $(g-2)_{\mu}$ studies in 2HDMs following the discovery of
the light Higgs boson was the work of Broggio et al. Broggio:2014mna . They
restricted their analysis to the “alignment limit” ($\cos(\beta-\alpha)=0$) in
which the tree-level couplings of the lighest 2HDM scalar are identical to
those of the SM Higgs. As noted above, they found that only the type-II and
type-X models could account for the $(g-2)_{\mu}$ discrepancy. In both models,
a very light pseudoscalar Higgs mass ($m_{A}$) is required, typically below
$100$ GeV, as is a relatively large $\tan\beta$, typically greater than $60$.
Another analysis of the type-X 2HDM with low-mass pseudoscalars is that of
Wang:2014sda , where it was shown that $m_{A}$ could be as low as 10 GeV. With
values for $m_{A}$ as in these references, unitarity and electroweak precision
constraints then force the charged Higgs mass ($m_{H^{+}}$) to be less than
about $200$ GeV. This is a problem within the type-II model, since radiative
$B$ decays, $\Delta m_{B_{s}}$ and the hadronic $Z\rightarrow\bar{b}b$
branching ratio force $m_{H^{+}}\gtrsim 600$ GeV Arbey:2017gmh . As a result,
the type-X model is favored. A subsequent analysis in Ref. Abe:2015oca
focusing explicitly on the type-X model considered all experimental
constraints. It was noted that bounds from leptonic $\tau$ decay restricted
the parameter-space further and that the discrepancy in $(g-2)_{\mu}$ could be
explained at the $2\sigma$ level with $m_{A}$ between $10$ and $30$ GeV,
$m_{H^{+}}$ between $200$ and $350$ GeV and $\tan\beta$ between $30$ and $50$.
Shortly thereafter, Ref. Chun:2016hzs  included more recent data from lepton
universality tests and found bounds that were somewhat weaker, but in general
agreement with Abe:2015oca . In the same token, studies have been performed in
Refs. Chun:2015hsa ; Ilisie:2015tra , and a more recent work including limits
from Higgs decays to $AA$ Wang:2018hnw  also found similar restrictions on the
parameter-space.

Another attempt to explain the $(g-2)_{\mu}$ discrepancy involved adding
vector-like leptons (VLL). It was shown Dermisek:2013gta ; Falkowski:2013jya
that this alternative works if the VLLs mix with the muon. However, in the SM
this not only alters the Higgs dimuon branching ratio but also affects
excessively the diphoton branching ratio and is thus phenomenologically
unacceptable. Recently, the addition of VLLs to the type-II and type-X models
was considered Frank:2020smf  and it was shown that the parameter-space can be
significantly expanded even without mixing the VLLs with muons, since the VLLs
contribute to $\Delta a_{\mu}$ in two-loop Barr-Zee diagrams. In this case,
the parameter-space of the type-X model is substantially widened and the type-
II model is not completely excluded. Even more recently, Dermisek, Hermanek
and McGinnis Dermisek:2020cod  explored the type-II model with VLLs, but now
allowed mixing with the muon. Here, the effects on the muon coupling to the
weak vector and scalar bosons must be considered. They showed that the extra
Higgs bosons and VLLs could be extremely heavy and still explain the
$(g-2)_{\mu}$ discrepancy without major effects on the dimuon SM Higgs decay.
A more detailed and comprehensive paper by Dermisek et al. has just appeared
Dermisek:2021ajd . None of the above papers deal with the smaller discrepancy
for the $(g-2)$ of the electron. A model explaining both discrepancies was
proposed by Chun and Mondal Chun:2020uzw  in which VLLs mix with the muon and
the electron. It turns out that also in such case a light $(<100\,{\rm GeV})$
pseudoscalar and large $\tan\beta$ are needed.

The belief that only the standard four 2HDMs can avoid flavor-changing neutral
currents at tree level was shown to be incorrect by Abe, Sato and Yagyu
Abe:2017jqo . They proposed a muon-specific ($\mu{\rm Spec}$) model in which
one doublet couples to the muon and the other couples to all other fermions.
This was implemented with one of the usual $Z_{2}$ symmetries, complemented
with a muon number conservation symmetry (which is a global $U(1)$, although
they referred to it as an overall $Z_{4}$ symmetry). This symmetry cannot be
implemented in the quark sector since it eliminates CKM mixing, but neutrino
mixing can easily be produced in a heavy Majorana sector. They studied
$(g-2)_{\mu}$ in the model and found that the discrepancy could be solved with
very large $\tan\beta$ values of approximately $1000$ (but they also showed
that perturbation theory was still valid). A study of the model by PF and MS
Ferreira:2020ukv  showed that this was fine-tuned and, leaving aside any
explanation of the $(g-2)_{\mu}$ discrepancy, studied other properties of the
model, including substantial effects on the dimuon decay of the Higgs.

Finally, one must take into account that there are many alternative attempts
to explain $(g-2)_{\mu}$, even in the context of 2HDMs. For instance, the
model of Han:2018znu  has a lepton-specific inert doublet with $\mu$ Yukawa
couplings having an opposite sign to those of $e$ and $\tau$; models with
tree-level FCNC were considered in this context in Crivellin:2015hha ;
Crivellin:2019dun ; in Arnan:2019uhr  the impact of a fourth generation of
vector-like fermions and an extra scalar was studied; contributions from
additional gauge scalar singlets were investigated in vonBuddenbrock:2016rmr ;
Sabatta:2019nfg ; an extremely light (sub-GeV) scalar was considered in
Jana:2020pxx ; a 2HDM study of $(g-2)_{\mu}$ considering the effects of warped
space was undertaken in Megias:2017dzd ; and a 2HDM complemented with a fourth
generation of fermions and a $Z^{\prime}$ gauge boson from an extra $U(1)_{X}$
symmetry was considered in Chen:2019wbk . As for non-2HDM explanations, one
has, for instance, an addition of a leptophilic scalar to the SM in
Chen:2021rnl ; UV complete models with a $L_{\mu}-L_{\tau}$ symmetry were
considered in Crivellin:2018qmi ; a simultaneous study of dark matter and
$(g-2)_{\mu}$ with extra scalars and vector-like fermions was undertaken in
Jana:2020joi ; a generic analysis of radiative leptonic-mass generation and
its impact on $(g-2)_{\mu}$ can be found in Yin:2021yqy ; and the interplay
between new physics contributions and the SM effective field theory for
$(g-2)_{\mu}$ is discussed in Fajfer:2021cxa . On a different strand, the
experimental implications of the muon anomalous magnetic moment at a future
muon collider were studied in Yin:2020afe .

All of the above analyses relied on the old experimental result Bennett:2006fi
and one would expect the new result from the Muon g-2 Collaboration results
to have a small effect on those studies. In this paper, we update some of the
2HDM results to include the new combined experimental value for $\Delta
a_{\mu}^{\rm exp}$ given in (2). After discussing the impact of the new
measurement in the context of the type-X and type-II models, we will focus on
the 2HDM extended with VLLs which do not mix with the muon. Adding such mixing
increases the number of parameters, and many of the relevant formulae are in
the recent work of Dermisek et al. Dermisek:2021ajd . In addition, we also
look at the $\mu{\rm Spec}$ model, which can suffer substantial changes in the
Higgs dimuon decay, and analyse the $(g-2)_{\mu}$ discrepancy when we add
VLLs.

##  II Type-II, type-X and $\bm{\mu{\rm Spec}}$ 2HDMs

In the 2HDM the scalar-fermion Yukawa interactions can be expressed in the
compact form:

$\displaystyle\mathcal{L}=$
$\displaystyle\dfrac{\sqrt{2}}{v}H^{+}\left\\{\bar{u}\left[\xi_{u}M_{u}^{\dagger}VP_{L}-\xi_{d}VM_{d}P_{R}\right]d-\xi_{\ell}\bar{\nu}M_{\ell}P_{R}\ell\right\\}$
$\displaystyle-\dfrac{1}{v}\sum_{k,f}y_{f}^{k}S_{k}^{0}\bar{f}M_{f}P_{R}f+{\rm
H.c.}\,,$ (3)

where $S_{k}^{0}=\\{h,H,A\\}$ are the neutral-scalar mass eigenstates ($h$ is
the SM Higgs with mass $m_{h}=125.1$ GeV Zyla:2020zbs ), and $f=u,d,\ell$
denote any SM charged-fermion type with mass matrix $M_{f}$. Following the
notation of Ref. Pich:2009sp , the couplings $y_{f,\ell}^{k}$ are

$\displaystyle y_{d,\ell}^{k}$
$\displaystyle=\mathcal{R}_{k1}+(\mathcal{R}_{k2}+i\,\mathcal{R}_{k3})\,\xi_{d,\ell}\,,$
(4) $\displaystyle y_{u}^{k}$
$\displaystyle=\mathcal{R}_{k1}+(\mathcal{R}_{k2}-i\,\mathcal{R}_{k3})\,\xi_{u}^{\ast}\,,$

where $\mathcal{R}$ is the orthogonal matrix which relates the scalar weak
states with the mass eigenstates $S_{k}^{0}$. For a CP-conserving potential,
$A$ does not mix with $h$ and $H$ implying
$\mathcal{R}_{3j}=\mathcal{R}_{j3}=0$ for $j\neq 3$ and $\mathcal{R}_{33}=1$.
On the other hand, $h-H$ mixing is determined by
$\mathcal{R}_{11}=-\mathcal{R}_{22}=\sin(\beta-\alpha)$ and
$\mathcal{R}_{12}=\mathcal{R}_{21}=\cos(\beta-\alpha)$, with
$\beta-\alpha=\pi/2$ in the alignment limit. For the type-II, type-X and
$\mu{\rm Spec}$ 2HDM, the $\xi_{f}$ couplings are

$\displaystyle{\rm Type\,II:}$
$\displaystyle\;\xi_{d,\ell}=-\tan\beta\;,\;\xi_{u}=\cot\beta\,,$ (5)
$\displaystyle{\rm Type\,X:}$
$\displaystyle\;\xi_{\ell}=-\tan\beta\;\;\;\,,\;\xi_{u,d}=\cot\beta\,,$
$\displaystyle\mu{\rm Spec:}$
$\displaystyle\;\xi_{\mu}=-\tan\beta\;\,\,\,,\;\xi_{u,d,\tau,e}=\cot\beta\,.$

In the alignment limit, the tree-level couplings of the SM-like scalar state,
$h$, are therefore identical to those of the SM Higgs boson – the agreement of
observed Higgs production and decays to that of the SM makes this assumption
quite reasonable.

---
Figure 1: $\Delta a_{\mu}$ allowed regions in the $(m_{A},\tan\beta)$ plane
for the type-X (top) and type-II (bottom) 2HDMs. The green, yellow and grey
shaded regions were obtained taking the new $\Delta a_{\mu}$ interval given in
Eq. (2) at 1,2 and $3\sigma$, respectively. For comparison, we also show the
lines which delimit the same regions when the old result in Eq. (1) is
considered.

In the early post-Higgs discovery study of the $(g-2)_{\mu}$ anomaly of Ref.
Broggio:2014mna , the type-II and type-X models were studied in the alignment
limit. The relevant free parameters in the model are $\tan\beta$, the
pseudoscalar mass $m_{A}$, the charged Higgs mass $m_{H^{+}}$, the heavy
scalar mass, $m_{H}$, and the $Z_{2}$ soft-breaking parameter $m^{2}_{12}$.
Broggio at al. considered constraints on these parameters due to electroweak
precision tests, vacuum stability and perturbativity. These imposed bounds on
the heavy Higgs masses, and they showed that for a pseudoscalar mass of less
than $100$ GeV, the charged Higgs mass could not exceed $200$ GeV. This was
not substantially affected by the value of $\tan\beta$ or by deviations from
the $\cos(\beta-\alpha)=0$ assumption. They then calculated the constraints
from the $(g-2)_{\mu}$ results (they also included constraints from the $g-2$
of the electron, but these are very weak).

We now revisit these findings in light of the new Muon g-2 Collaboration
result considering all one-loop and two-loop Barr-Zee (BZ) contributions to
$\Delta a_{\mu}$ as given in Ref. Ilisie:2015tra . 111Keeping the fermion BZ
diagrams only, as done by Broggio et. al., we were able to reproduce their
results for the type-II and type-X 2HDM in the alignment limit. In Fig. 1 we
show the allowed regions in the $m_{A}$ and $\tan\beta$ plane corresponding to
the $\Delta a_{\mu}$ intervals in Eq. (2) at the $1\sigma$, $2\sigma$ and
$3\sigma$ levels (green, yellow and grey shaded regions) for the type-X (top
panel) and type-II (bottom panel) cases. The solid, dashed and dash-dotted
contours delimit the analogue regions when the old result in Eq. (2) is
considered. As noted earlier, radiative $B$-decays favor $m_{H^{+}}>600$ GeV
in the type-II model Arbey:2017gmh . Thus, we take $m_{H^{+}}=600$ GeV as
reference value in the bottom panel of Fig. 1, while for type X a lighter
$H^{+}$ is considered (changing $m_{H^{+}}$ does not have a significant impact
on $\Delta a_{\mu}$). As pointed out in the Introduction, such small values of
$m_{A}$ in the type-II model are in conflict with unitarity and electroweak
precision constraints. On the other hand, type X is still allowed. The results
show that, as expected, the impact of the new $(g-2)_{\mu}$ result is
marginal. The $\mu$Spec 2HDM requires extreme fine-tuning and values of
$\tan\beta$ of $O(1000)$ in order to accommodate the $(g-2)_{\mu}$
discrepancy, and we do not show a plot in this case. It is worth mentioning
that for $\mu$Spec the BZ contributions to $\Delta a_{\mu}$ are not
significant as they are for type II and type X since both the $b$-quark and
$\tau$-lepton couplings are not $\tan\beta$ enhanced.

##  III 2HDM with vector-like Leptons

Let us now consider a 2HDM extension with vector-like leptons
$\chi_{L,R}=(N\,L^{-})_{L,R}^{T}\sim(\bm{2},-1/2)$ and
$E_{L,R}\sim(\bm{1},-1)$, with $N_{L,R}$ being neutral and
$L^{-}_{L,R},E_{L,R}$ charged. If taken within the context of the muon-
specific 2HDM, the VLLs will have no quantum numbers under the underlying muon
number conservation symmetry of that model. We choose to have no mixing
between the VLLs and the usual leptons – this can be achieved, for instance,
by imposing an extra $Z_{2}$ symmetry on the model, under which all VLL fields
have charge -1, and all other fields have charge +1. The VLL Yukawa and mass
terms are

$\displaystyle-\mathcal{L}_{\rm VLL}=$ $\displaystyle
m_{L}\overline{\chi_{L}}\chi_{R}+m_{E}\overline{E_{L}}E_{R}+\lambda_{L}\overline{\chi_{R}}\Phi_{1}E_{L}$
$\displaystyle+\lambda_{R}\overline{\chi_{L}}\Phi_{1}E_{R}+{\rm H.c.}\,,$ (6)

which, after electroweak symmetry breaking, lead to the following heavy
charged-lepton mass matrix

$\displaystyle\mathcal{M}=\begin{pmatrix}m_{L}&\lambda_{R}v\cos\beta/\sqrt{2}\\\
\lambda_{L}^{\ast}v\cos\beta/\sqrt{2}&m_{E}\end{pmatrix}\,,$ (7)

in the $(L^{-}\;E)_{L,R}^{T}$ basis, while for the new neutral-lepton mass
$M_{N}=|m_{L}|$. The above matrix can be diagonalized as
$U_{L}^{\dagger}\mathcal{M}U_{R}={\rm diag}(M_{1},M_{2})$, where $M_{1,2}$ are
the masses of the $L_{1,2}^{-}$ heavy charged-lepton physical states and
$U_{L,R}$ are $2\times 2$ unitary matrices. As for the new charged-current
interactions, we have

$\displaystyle\\!\\!\\!\\!\mathcal{L}=\dfrac{g}{\sqrt{2}}\bar{N}\gamma^{\mu}\left[(U_{L})_{1a}P_{L}+(U_{R})_{1a}P_{R}\right]L_{a}^{-}W_{\mu}^{+}+{\rm
H.c.},$ (8)

where $g=e/\sin\theta_{W}$ is the SU(2)L gauge coupling. The interactions with
the neutral and charged physical scalars are

$\displaystyle-\mathcal{L}=$
$\displaystyle\dfrac{\sqrt{2}}{v}H^{+}\bar{N}\left[\,\xi_{a}M_{a}(U_{R})_{2a}P_{R}+\xi_{N}M_{N}(U_{L})_{2a}P_{L}\right]L_{a}^{-}$
$\displaystyle+\dfrac{1}{v}\sum_{k,a,b}M_{a}\,y_{ab}^{k}\overline{L^{-}_{a}}P_{R}L^{-}_{b}S_{k}^{0}+{\rm
H.c.}\,,$ (9)

with $S_{k}^{0}=\\{h,H,A\\}$, $a,b=1,2$ and

$\displaystyle\xi_{a}$
$\displaystyle=\xi_{\ell}\frac{\mathcal{M}_{12}}{M_{a}}\;,\;\xi_{N}=\xi_{\ell}\frac{\mathcal{M}_{21}^{\ast}}{M_{N}}\;,\;\xi_{\ell}=-\tan\beta\,,$
(10) $\displaystyle y_{ab}^{k}$
$\displaystyle=(\mathcal{R}_{k1}-\xi_{\ell}\mathcal{R}_{k2})X_{ab}^{+}+i\,\xi_{\ell}\,\mathcal{R}_{k3}X_{ab}^{-}\,.$

For simplicity, we will restrict ourselves to the case of real $\mathcal{M}$,
for which $U_{L,R}$ are orthogonal and parametrized by two mixing angles
$\theta_{L,R}$ with $(U_{L,R})_{11}=(U_{L,R})_{22}=c_{L,R}$ and
$(U_{L,R})_{12}=-(U_{L,R})_{21}=s_{L,R}$ (the notation $\cos\theta_{L,R}\equiv
c_{L,R}$, $\sin\theta_{L,R}\equiv s_{L,R}$ has been used). Under this
assumption, and taking into account that the only relevant VLL couplings are
the flavor-conserving ones ($a=b$), we have:

$\displaystyle X_{aa}^{+}$
$\displaystyle=c_{R}^{2}s_{L}^{2}+c_{L}^{2}s_{R}^{2}-\dfrac{1}{2}s_{2L}s_{2R}\left(\dfrac{M_{1}}{M_{2}}\right)^{(-1)^{a}}\,,$
(11) $\displaystyle X_{aa}^{-}$
$\displaystyle=\dfrac{(-1)^{a+1}}{2}(c_{2L}-c_{2R})\,,$

with $c_{2L,2R}\equiv\cos(2\theta_{L,R})$ and
$s_{2L,2R}\equiv\sin(2\theta_{L,R})$.

In the absence of VLL mixing with the muon, there are no new one-loop
contributions to $\Delta a_{\mu}$ besides those already present in the 2HDM.
Still, BZ diagrams involving the new charged and neutral leptons must be taken
into account (see Fig. 2).

Figure 2: BZ contributions to $\Delta a_{\mu}$ involving the new charged and
neutral leptons $L_{a}^{-}$ and $N$, respectively.

We have computed both contributions and, for diagram (a), our result agrees
with that of Ref. Ilisie:2015tra , doing the appropriate replacements
according to the interactions shown in (9). Namely,

$\displaystyle\Delta a_{\mu}^{\rm(a)}=$
$\displaystyle\dfrac{\alpha\,m_{\mu}^{2}}{4\pi^{3}v^{2}}\sum_{k,a}Q_{a}^{2}\Biggl{[}\text{Re}\left(y_{aa}^{k}\right)\text{Re}\left(y_{\mu}^{k}\right)\mathcal{F}^{(1)}\left(\dfrac{M_{a}^{2}}{M_{k}^{2}}\right)\Biggr{.}$
$\displaystyle+\,\Biggl{.}\text{Im}\left(y_{aa}^{k}\right)\text{Im}\left(y_{\mu}^{k}\right)\mathcal{F}^{(2)}\left(\dfrac{M_{a}^{2}}{M_{k}^{2}}\right)\Biggr{]}\,,$
(12)

with

$\mathcal{F}^{(1,2)}(r)=\dfrac{r}{2}\int_{0}^{1}dx\dfrac{\mathcal{N}^{(1,2)}}{r-x(1-x)}\text{ln}\left[\dfrac{r}{x(1-x)}\right]\,,$
(13)

and $\mathcal{N}^{(1)}=2x(1-x)-1$ and $\mathcal{N}^{(2)}=1$. Regarding diagram
(b), caution must be taken since, contrarily to what happens in the 2HDM,
charged-current interactions are now left- and right-handed – see Eq. (8). The
result in this case is

$\displaystyle\Delta a_{\mu}^{\rm(b)}=$
$\displaystyle\dfrac{\alpha\,m_{\mu}^{2}}{32\pi^{3}s_{W}^{2}v^{2}\left(M^{2}_{H^{+}}-M_{W}^{2}\right)}\sum_{a}\int_{0}^{1}dx\,Q_{a}(1-x)$
$\displaystyle\times\left[\mathcal{G}\left(\dfrac{M_{N}^{2}}{M^{2}_{H^{+}}},\dfrac{M_{a}^{2}}{M^{2}_{H^{+}}}\right)-\mathcal{G}\left(\dfrac{M_{N}^{2}}{M^{2}_{W}},\dfrac{M_{a}^{2}}{M^{2}_{W}}\right)\right]$
$\displaystyle\times\bigl{(}\omega_{L}+\omega_{R}\bigr{)}\,,$ (14)

where the function $\mathcal{G}$ is

$\mathcal{G}(r_{1},r_{2})=\dfrac{\text{ln}\left[\dfrac{r_{1}x+r_{2}(1-x)}{x(1-x)}\right]}{x(1-x)-r_{1}x-r_{2}(1-x)}\,.$
(15)

The term proportional to $\omega_{L(R)}$ takes the left (right)-handed
component of the charged-current interactions with the new leptons, such that:

$\displaystyle\omega_{L}$
$\displaystyle=\text{Re}\left[(U_{L})_{1a}^{*}\,\xi_{a}M_{a}(U_{R})_{2a}\,\xi_{\mu}^{*}\right]M_{a}\left[(1-x)-(1-x)^{2}\right]$
(16)
$\displaystyle-\text{Re}\left[(U_{L})_{1a}^{*}\,\xi_{N}M_{N}(U_{L})_{2a}\,\xi_{\mu}^{*}\right]M_{N}x(1+x)\,,$
$\displaystyle\omega_{R}$
$\displaystyle=-\text{Re}\left[(U_{R})_{1a}^{*}\,\xi_{N}M_{N}(U_{L})_{2a}\,\xi_{\mu}^{*}\right]M_{a}\left[(1-x)+(1-x)^{2}\right]$
$\displaystyle+\text{Re}\left[(U_{R})_{1a}^{*}\,\xi_{a}M_{a}(U_{R})_{2a}\,\xi_{\mu}^{*}\right]M_{N}\,x(1-x)\,.$

The first term is analogous to that for quarks presented in Ref.
Ilisie:2015tra , while the second has the same structure as the general terms
obtained in Ref. Yang:2009zzh .

While fitting the experimental values of $(g-2)_{\mu}$ within the context of a
2HDM, we must make sure that all experimental and theoretical constraints of
that model are satisfied. We have already mentioned $B$-physics constraints –
particularly relevant for the type-II model, wherein $b\rightarrow s\gamma$
results force the charged Higgs mass to be very high – and bounded from below
and unitarity bounds, which limit the values of the scalar potential’s quartic
couplings (see Branco:2011iw  for the explicit expressions). Equally important
is to consider electroweak precision constraints, in the form of the Peskin-
Takeuchi parameters $S$, $T$ and $U$ Peskin:1990zt . Their explicit expression
for the 2HDM can be found in numerous references (see for instance eqs. (12)
and (13) of Kanemura:2011sj ), but if one considers vector-like leptons
present as well, their contributions to the $T$ parameter must also be taken
into account (see Eq. (4.5) of Chun:2020uzw ) – the upshot of their inclusion
in a 2HDM fit is that it usually is always possible, for VLL masses $M_{1}$
and $M_{2}$ below $\sim$ 500 GeV, to find a vector-like neutrino mass $m_{N}$
to satisfy the current constraints on $T$. Indeed, these electroweak precision
constraints, as is usual in the vanilla 2HDM, tend to reduce the splittings
between the extra (scalar and VLL) masses in the theory.

With or without VLLs, LHC precision constraints on the properties of the 125
GeV $h$ scalar are, with the exception of the diphoton branching ratio,
trivially satisfied within the alignment limit, since $h$’s tree-level
couplings become identical to the SM Higgs when $\cos(\beta-\alpha)=0$. The
presence of a charged scalars and two charged VLLs, however, contribute to the
diphoton decay width. The decay amplitude $\mathcal{A}$ for
$h\rightarrow\gamma\gamma$ includes therefore a contribution identical to the
SM, and another from the charged and VLL sector, given by

$\mathcal{A}=\mathcal{A}_{\rm
SM}-\frac{\Lambda}{2m^{2}_{H^{+}}}A^{H}_{0}(\tau_{+})+y^{h}_{11}A^{H}_{1/2}(\tau_{1})+y^{h}_{22}A^{H}_{1/2}(\tau_{2})\,,$
(17)

where $\tau_{X}=m_{h}^{2}/(4m^{2}_{X})$,
$\Lambda=3m^{2}_{h}-2m^{2}_{H^{+}}-4m^{2}_{12}/s_{2\beta}$ and $A^{H}_{0}$,
$A^{H}_{1/2}$ are the well-known scalar and fermionic form factors of the
Higgs (see for instance Djouadi:2005gj ). Depending on the charged and VLL
masses and couplings, their contributions to $h\rightarrow\gamma\gamma$ can
actually cancel each other, and will be less important for higher masses.

---
Figure 3: $(m_{A},\tan\beta)$ allowed regions for the type-X, type-II and
$\mu$Spec 2HDM with VLLs (from top to bottom panel). See text for details.

To illustrate the main features of the VLL models we are interested in, we
focus on two benchmark cases (B1 and B2) within the type-II, type-X and
$\mu$Spec 2HDM with heavy charged-lepton masses:

$\displaystyle{\rm B}_{1}:M_{2}=2.5M_{1}=250\,{\rm GeV}\,,$ $\displaystyle{\rm
B}_{2}:M_{2}=1.5M_{1}=150\,{\rm GeV}\,,$

taking $s_{L}(s_{R})=0.5(0.4)$. The neutral-lepton mass $M_{N}$ turns out to
be fixed since $M_{N}=|m_{L}|=|M_{1}c_{L}c_{R}+M_{2}s_{L}s_{R}|$, yielding
$M_{N}\simeq 129\,(109)$ GeV for B1 (B2). The scalar masses $m_{H^{+}}$ and
$m_{H}$ are the same as those considered in Fig. 1 (for $\mu$Spec we adopt the
same setting as for type X) and, again, $\cos(\beta-\alpha)=0$.222The impact
on $\Delta a_{\mu}$ of considering $\cos(\beta-\alpha)\neq 0$ within the
experimentally-allowed limits is marginal. For instance, taking
$\cos(\beta-\alpha)=0.05$ the change in the results is almost imperceptible.
In Fig. 3 we show the $(m_{A},\tan\beta)$ allowed regions in green and yellow
(blue and red) for B1 (B2) at 1 and $2\sigma$, respectively. Above the solid
(dashed) horizontal line, $\lambda_{\rm max}\equiv{\rm
max}\\{|\lambda_{L}|,|\lambda_{R}|\\}>4\pi$ for B1 (B2).

The results show that in all cases the ranges of $m_{A}$ can be substantially
enlarged with respect to Fig. 1, while simultaneously shifting $\tan\beta$ to
lower values. This effect is more pronounced for B1 since $M_{2}/M_{1}$ is
larger than for B2, and so are the $L_{2}L_{2}S^{0}_{k}$ couplings – see Eqs.
(10) and (11). The main drawback of this scenario is that the larger VLL
coupling with the Higgs doublet $\Phi_{1}$ is required to be close or above
the perturbative limit, i.e. $\lambda_{\rm max}=4\pi$. This can be easily
understood noting that the coupling modifiers $\xi_{a,N}$ in Eq. (10) scale as
$\lambda v\cos\beta/(\sqrt{2}M)$ (with respect to the pure 2HDM case) where
$\lambda$ and $M$ stand for a generic VLL coupling and mass, respectively.
Requiring that, at least, $\xi_{a,N}\simeq\tan\beta$ implies $\lambda
v\cos\beta/(\sqrt{2}M)\simeq 1$. Thus, for $M\gtrsim 100$ GeV and large
$\tan\beta$, $\lambda$ must be large to avoid spoiling the $\tan\beta$
enhancement required to explain the $\Delta a_{\mu}$ discrepancy in 2HDMs. For
the $\mu$Spec model, the effect of adding VLLs allows to drastically lower
$\tan\beta\gtrsim\mathcal{O}(1000)$ down to $\tan\beta\sim 40-50$ for
$m_{A}\gtrsim 200$ GeV.

Since the benchmark analysis shows type X is favoured for agreement with the
muon anomaly, we have performed a general parameter space scan of the 2HDM for
that model, considering extra scalars and VLLs (charged and neutral) to have
masses ranging from 100 to 1000 GeV, as well as all possible values for the
angles $\theta_{L}$ and $\theta_{R}$. We only accepted combinations of
parameters for which the model satisfies unitarity; boundedness from below;
electroweak precision constraints; $b\rightarrow s\gamma$ constraints; and
predicts a Higgs diphoton branching ratio within 10% of the SM values.
Finally, since extra scalars are required to explain the $(g-2)_{\mu}$
anomaly, we need to verify whether current LHC searches would not have
discovered those particles already. The $B$-decay constraints we mentioned
take care of the charged Higgs; because we are working in the alignment limit,
searches for the heavier CP-even scalar $H$ in $ZZ$ or $WW$ are automatically
satisfied; indeed, the channel which we need to worry the most would be
searches for $H$ and specially $A$ decaying into tau pairs, having been
produced via gluon fusion. Tau decays can be $\tan\beta$-enhanced and produce
strong signal rates already excluded by LHC collaborations Aaboud:2017sjh . We
verified that, after all other constraints had been taken into account, our
scan included only a small fraction of points above current experimental
limits on ditau searches for extra scalars, and excluded those points from our
$(g-2)_{\mu}$ calculations. The result of that scan show that the type X
behaviour shown in Fig. 3 can be reproduced by a wealth of different
combinations of parameters; the values of allowed $\tan\beta$ are even more
general, with, for instance, $\tan\beta\simeq 25$ permitted for $m_{A}=400$
GeV. However, and also confirming the previous conclusions, though unitarity
and perturbativity can be satisfied within the scalar sector, agreement with
$(g-2)_{\mu}$ requires large VLL yukawas – in particular, though $\lambda_{L}$
can be found small, $|\lambda_{R}|$ is always found to be above $\sim$ 20, and
thus non-perturbative.

##  IV Conclusions

We have studied the theory/experiment discrepancy of the muon’s anomalous
magnetic moment in the 2HDM. Without vector-like leptons, the type-I and
type-Y models cannot explain the discrepancy and the type II requires light
pseudoscalars that are in conflict with perturbativity, unitarity and
electroweak precision constraints. The type-X model can accommodate the
discrepancy but require large values of $\tan\beta$ and also very light
pseudoscalars. After revisiting some 2HDM results in light of the new result
reported by the Muon g-2 collaboration at Fermilab, we considered adding
vector-like leptons to the type-II, type-X and muon-specific 2HDMs. We did not
include any mixing between the vector-like leptons and the muon. We have shown
that the parameter space is substantially widened and much larger values of
the pseudoscalar mass are allowed. $\tan\beta$ remains fairly large in all of
these models. One problem of the VLL scenario analysed in this work is that
the VLL Yukawa couplings to the scalar doublet $\Phi_{1}$ must be close to or
above the perturbation theory limit. This can be alleviated by considering $N$
families of vector-like leptons instead of only one (very roughly, this would
reduce the maximum values of these couplings by a factor of $1/N$).
Alternatively, one can include VLL-muon mixing. This has been done in Refs
Dermisek:2013gta ; Dermisek:2021ajd  in some detail, where it was also shown
that the two-loop Barr-Zee diagrams are substantially less important than one-
loop contributions and that smaller values of $\tan\beta$ can be accommodated.

Acknowledgements: This work is supported by Fundação para a Ciência e a
Tecnologia (FCT, Portugal) through the projects UIDB/00777/2020,
UIDP/00777/2020, UIDB/00618/2020, CERN/FIS-PAR/0004/2019,
CERN/FISPAR/0014/2019, PTDC/FIS-PAR/29436/2017 and by HARMONIA project’s
contract UMO-2015/18/M/ST2/00518. The work of B.L.G. is supported by the FCT
grant SFRH/BD/139165/2018. The work of M.S. was supported by the National
Science Foundation under grant PHY-1819575

## References

  * (1) B. Abi et al. (Muon g-2 Collaboration) “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm” Phys. Rev. Lett. 126, 141801 (2021).
  * (2) G. W. Bennett et al. [Muon g-2], Phys. Rev. D 73, 072003 (2006) doi:10.1103/PhysRevD.73.072003 [arXiv:hep-ex/0602035 [hep-ex]].
  * (3) T. Blum et al. [RBC and UKQCD], Phys. Rev. Lett. 121, no.2, 022003 (2018) doi:10.1103/PhysRevLett.121.022003 [arXiv:1801.07224 [hep-lat]].
  * (4) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 97, no.11, 114025 (2018) doi:10.1103/PhysRevD.97.114025 [arXiv:1802.02995 [hep-ph]].
  * (5) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 80, no.3, 241 (2020) [erratum: Eur. Phys. J. C 80, no.5, 410 (2020)] doi:10.1140/epjc/s10052-020-7792-2 [arXiv:1908.00921 [hep-ph]].
  * (6) T. Aoyama, N. Asmussen, M. Benayoun, J. Bijnens, T. Blum, M. Bruno, I. Caprini, C. M. Carloni Calame, M. Cè and G. Colangelo, et al. Phys. Rept. 887, 1-166 (2020) doi:10.1016/j.physrep.2020.07.006 [arXiv:2006.04822 [hep-ph]].
  * (7) M. Lindner, M. Platscher and F. S. Queiroz, Phys. Rept. 731, 1-82 (2018) doi:10.1016/j.physrep.2017.12.001 [arXiv:1610.06587 [hep-ph]].
  * (8) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher and J. P. Silva, Phys. Rept. 516 (2012) 1 doi:10.1016/j.physrep.2012.02.002 [arXiv:1106.0034 [hep-ph]].
  * (9) V. Keus, N. Koivunen and K. Tuominen, JHEP 09, 059 (2018) doi:10.1007/JHEP09(2018)059 [arXiv:1712.09613 [hep-ph]].
  * (10) A. Broggio, E. J. Chun, M. Passera, K. M. Patel and S. K. Vempati, JHEP 11, 058 (2014) doi:10.1007/JHEP11(2014)058 [arXiv:1409.3199 [hep-ph]].
  * (11) L. Wang and X. F. Han, JHEP 05, 039 (2015) doi:10.1007/JHEP05(2015)039 [arXiv:1412.4874 [hep-ph]].
  * (12) A. Arbey, F. Mahmoudi, O. Stal and T. Stefaniak, Eur. Phys. J. C 78, no.3, 182 (2018) doi:10.1140/epjc/s10052-018-5651-1 [arXiv:1706.07414 [hep-ph]].
  * (13) T. Abe, R. Sato and K. Yagyu, JHEP 07, 064 (2015) doi:10.1007/JHEP07(2015)064 [arXiv:1504.07059 [hep-ph]].
  * (14) E. J. Chun and J. Kim, JHEP 07, 110 (2016) doi:10.1007/JHEP07(2016)110 [arXiv:1605.06298 [hep-ph]].
  * (15) E. J. Chun, Z. Kang, M. Takeuchi and Y. L. S. Tsai, JHEP 11, 099 (2015) doi:10.1007/JHEP11(2015)099 [arXiv:1507.08067 [hep-ph]].
  * (16) V. Ilisie, JHEP 04, 077 (2015) doi:10.1007/JHEP04(2015)077 [arXiv:1502.04199 [hep-ph]].
  * (17) L. Wang, J. M. Yang, M. Zhang and Y. Zhang, Phys. Lett. B 788, 519-529 (2019) doi:10.1016/j.physletb.2018.11.045 [arXiv:1809.05857 [hep-ph]].
  * (18) R. Dermisek and A. Raval, Phys. Rev. D 88, 013017 (2013) doi:10.1103/PhysRevD.88.013017 [arXiv:1305.3522 [hep-ph]].
  * (19) A. Falkowski, D. M. Straub and A. Vicente, JHEP 05, 092 (2014) doi:10.1007/JHEP05(2014)092 [arXiv:1312.5329 [hep-ph]].
  * (20) M. Frank and I. Saha, Phys. Rev. D 102, no.11, 115034 (2020) doi:10.1103/PhysRevD.102.115034 [arXiv:2008.11909 [hep-ph]].
  * (21) R. Dermisek, K. Hermanek and N. McGinnis, [arXiv:2011.11812 [hep-ph]].
  * (22) R. Dermisek, K. Hermanek and N. McGinnis, [arXiv:2103.05645 [hep-ph]].
  * (23) E. J. Chun and T. Mondal, JHEP 11, 077 (2020) doi:10.1007/JHEP11(2020)077 [arXiv:2009.08314 [hep-ph]].
  * (24) T. Abe, R. Sato and K. Yagyu, JHEP 07, 012 (2017) doi:10.1007/JHEP07(2017)012 [arXiv:1705.01469 [hep-ph]].
  * (25) P. M. Ferreira and M. Sher, Phys. Rev. D 101, no.9, 095030 (2020) doi:10.1103/PhysRevD.101.095030 [arXiv:2002.01000 [hep-ph]].
  * (26) X. F. Han, T. Li, L. Wang and Y. Zhang, Phys. Rev. D 99, no.9, 095034 (2019) doi:10.1103/PhysRevD.99.095034 [arXiv:1812.02449 [hep-ph]].
  * (27) A. Crivellin, J. Heeck and P. Stoffer, Phys. Rev. Lett. 116, no.8, 081801 (2016) doi:10.1103/PhysRevLett.116.081801 [arXiv:1507.07567 [hep-ph]].
  * (28) A. Crivellin, D. Müller and C. Wiegand, JHEP 06, 119 (2019) doi:10.1007/JHEP06(2019)119 [arXiv:1903.10440 [hep-ph]].
  * (29) P. Arnan, A. Crivellin, M. Fedele and F. Mescia, JHEP 06, 118 (2019) doi:10.1007/JHEP06(2019)118 [arXiv:1904.05890 [hep-ph]].
  * (30) S. von Buddenbrock, N. Chakrabarty, A. S. Cornell, D. Kar, M. Kumar, T. Mandal, B. Mellado, B. Mukhopadhyaya, R. G. Reed and X. Ruan, Eur. Phys. J. C 76, no.10, 580 (2016) doi:10.1140/epjc/s10052-016-4435-8 [arXiv:1606.01674 [hep-ph]].
  * (31) D. Sabatta, A. S. Cornell, A. Goyal, M. Kumar, B. Mellado and X. Ruan, Chin. Phys. C 44, no.6, 063103 (2020) doi:10.1088/1674-1137/44/6/063103 [arXiv:1909.03969 [hep-ph]].
  * (32) S. Jana, V. P. K. and S. Saad, Phys. Rev. D 101, no.11, 115037 (2020) doi:10.1103/PhysRevD.101.115037 [arXiv:2003.03386 [hep-ph]].
  * (33) N. Chen, B. Wang and C. Y. Yao, [arXiv:2102.05619 [hep-ph]].
  * (34) A. Crivellin, M. Hoferichter and P. Schmidt-Wellenburg, Phys. Rev. D 98, no.11, 113002 (2018) doi:10.1103/PhysRevD.98.113002 [arXiv:1807.11484 [hep-ph]].
  * (35) S. Jana, P. K. Vishnu, W. Rodejohann and S. Saad, Phys. Rev. D 102, no.7, 075003 (2020) doi:10.1103/PhysRevD.102.075003 [arXiv:2008.02377 [hep-ph]].
  * (36) W. Yin, [arXiv:2103.14234 [hep-ph]].
  * (37) S. Fajfer, J. F. Kamenik and M. Tammaro, [arXiv:2103.10859 [hep-ph]].
  * (38) W. Yin and M. Yamaguchi, [arXiv:2012.03928 [hep-ph]].
  * (39) E. Megias, M. Quiros and L. Salas, JHEP 05, 016 (2017) doi:10.1007/JHEP05(2017)016 [arXiv:1701.05072 [hep-ph]].
  * (40) C. R. Chen, C. W. Chiang and K. Y. Lin, Phys. Lett. B 795, 22-28 (2019) doi:10.1016/j.physletb.2019.05.023 [arXiv:1902.01001 [hep-ph]].
  * (41) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020) doi:10.1093/ptep/ptaa104
  * (42) A. Pich and P. Tuzon, Phys. Rev. D 80, 091702 (2009) doi:10.1103/PhysRevD.80.091702 [arXiv:0908.1554 [hep-ph]].
  * (43) X. Y. Yang and T. F. Feng, Phys. Lett. B 675, 43-48 (2009) doi:10.1016/j.physletb.2009.03.071
  * (44) M. E. Peskin and T. Takeuchi, Phys. Rev. Lett. 65, 964-967 (1990) doi:10.1103/PhysRevLett.65.964
  * (45) S. Kanemura, Y. Okada, H. Taniguchi and K. Tsumura, Phys. Lett. B 704, 303-307 (2011) doi:10.1016/j.physletb.2011.09.035 [arXiv:1108.3297 [hep-ph]].
  * (46) A. Djouadi, Phys. Rept. 459, 1-241 (2008) doi:10.1016/j.physrep.2007.10.005 [arXiv:hep-ph/0503173 [hep-ph]].
  * (47) M. Aaboud et al. [ATLAS], JHEP 01, 055 (2018) doi:10.1007/JHEP01(2018)055 [arXiv:1709.07242 [hep-ex]].


# Online Saddle Point Tracking with Decision-Dependent Data

Killian Wood and Emiliano Dall’Anese Department of Applied Mathematics,
University of Colorado, Boulder, CO (killian.wood@colorado.edu).Department of
Electrical, Computer, and Energy Engineering, and Department of Applied
Mathematics, University of Colorado, Boulder, CO
(emiliano.dallanese@colorado.edu).

###### Abstract

In this work, we consider a time-varying stochastic saddle point problem in
which the objective is revealed sequentially, and the data distribution
depends on the decision variables. Problems of this type express the
distributional dependence via a distributional map, and are known to have two
distinct types of solutions—saddle points and equilibrium points. We
demonstrate that, under suitable conditions, online primal-dual type
algorithms are capable of tracking equilibrium points. In contrast, since
computing closed-form gradient of the objective requires knowledge of the
distributional map, we offer an online stochastic primal-dual algorithm for
tracking equilibrium trajectories. We provide bounds in expectation and in
high probability, with the latter leveraging a sub-Weibull model for the
gradient error. We illustrate our results on an electric vehicle charging
problem where responsiveness to prices follows a location-scale family based
distributional map.

##  1 Introduction

The general goal of stochastic optimization is to find optimal decisions in
systems with parameters dictated by data [22, 29, 36]. In statistical
learning, optimal decisions represent model parameters that best fit a mapping
between feature and label data (see, e.g., [11, 30]). In the context of
optimization of physical and dynamical systems, they may model externalities
or system parameters that are predicted from data and are accompanied by given
error statistics (see, e.g., [2, 4, 15]). To analyze these problems, works
posit that the data distributions are stationary [5]; in modern machine
learning and cyber-physical systems applications this assumption may be
violated when population data shifts in response to previously deployed
decisions, thus making said decisions sub-optimal. Hence the distribution is
inextricably tied to the decision variables.

This work considers the problem of tracking the solution trajectories for
problems of the form:

$\min_{x\in\mathcal{X}_{t}}\max_{y\in\mathcal{Y}_{t}}\left\\{F_{t}(x,y):=\underset{w\sim
D_{t}(x,y)}{\operatorname{\mathbb{E}}}[f_{t}(x,y,w)]\right\\}$ (1)

where $t$ is a time index, $\mathcal{X}_{t}\subseteq\mathbb{R}^{n}$ and
$\mathcal{Y}_{t}\subseteq\mathbb{R}^{m}$ are convex and compact sets capturing
time-varying constraints,
$f_{t}:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{k}\rightarrow\mathbb{R}$
is a strongly-convex-strongly-concave function revealed at time $t$, and
$D_{t}:\mathbb{R}^{d}\times\mathbb{R}^{n}\rightarrow\mathcal{P}(\mathbb{R}^{k})$
is a distributional map that maps decision variables to the set of finite-
first moment probability distributions supported on $\mathbb{R}^{k}$ denoted
by $\mathcal{P}(\mathbb{R}^{k})$. Without loss of generality, we refer to the
support of $w$ as $\mathbb{R}^{k}$ (even if $w$ is matrix valued, our analysis
holds as $w$ is isomorphic to its vectorization over $\mathbb{R}^{k}$).

Examples of problems of the form (1) emerge in cost maximization in
competitive markets, where the (stochastic) demand shifts in response to
prices (see, e.g., [17, 31]), and in applications in adversarial strategic
classification, finance, energy systems, transportation networks, and ride-
sharing—just to mention a few. Focusing on the first example, consider a
competition between two service providers in an area with $n$ distinct regions
for which each provider seeks to maximize their relative revenue, and when the
demand for each provider’s service changes in response to the price variation
set by both providers. This problem can be written as the saddle point problem

$\min_{x\in\mathcal{X}_{t}}\max_{y\in\mathcal{Y}_{t}}\left\\{F_{t}(x,y)=\operatorname{\mathbb{E}}_{(a,b)\sim
D_{t}(x,y)}\|\Gamma_{t}^{1}x\|^{2}-\|\Gamma_{t}^{2}y\|^{2}-\langle
a+c_{t},x\rangle+\langle b+c_{t},y\rangle\right\\}\,,$ (2)

where $x=(x_{i})_{i=1}^{n}$ and $y=(y_{i})_{i=1}^{n}$ are vectors of price
deviations from a nominal value for providers one and two respectively
(components $x_{i}$ and $y_{i}$ are the prices in region $i\in[n]$);
$\Gamma_{t}^{1},\Gamma_{t}^{2}\in\mathbb{R}^{n\times n}$ are the charging rate
utility matrices; $c_{t}\in\mathbb{R}^{n}$ is the location-based utility
vector (i.e., cost of operation); and $a,b\in\mathbb{R}^{n}$ are changes in
demand in each region (in response to price changes) with distributions
$a\overset{d}{=}a^{t}_{0}+A^{t}_{1}x+B^{t}_{1}y$, and
$b\overset{d}{=}b^{t}_{0}+A^{t}_{2}x+B^{t}_{2}y$. Here $a^{t}_{0}$ and
$b^{t}_{0}$ are random variables drawn from zero-mean stationary
distributions.

Classical solutions to (1) are saddle points, which we denote
$z_{t}^{*}=(x_{t}^{*},y_{t}^{*})\in\mathcal{X}_{t}\times\mathcal{Y}_{t}$.
Under appropriate conditions, namely minimax equality, saddle points satisfy

$x_{t}^{*}\in\mathrm{arg}\min_{x\in\mathcal{X}_{t}}\max_{y\in\mathcal{Y}_{t}}F_{t}(x,y),\
\
y_{t}^{*}\in\mathrm{arg}\max_{y\in\mathcal{Y}_{t}}\min_{x\in\mathcal{X}_{t}}F_{t}(x,y).$
(3)

In this setting, saddle points are optimal decisions that effectively
anticipate the distributional shift, and hence are optimal even after the data
distribution has changed in the system. While these are ideal, finding them is
typically computationally intractable. While sufficient conditions for their
existence and uniqueness have been studied, guarantees for convergence to
saddle points are only approximate or require explicit knowledge of a model
for the distributional map [20, 35]. A common heuristic to overcome
distributional shift in general is to repeatedly retrain the optimal decisions
each time the distribution shifts. This amounts to forming a sequence
$\\{z_{t}^{\ell}\\}_{\ell\geq 0}=\\{(x_{t}^{\ell},y_{t}^{\ell})\\}_{\ell\geq
0}$ at each time $t$ defined by

$\displaystyle
x_{t}^{\ell+1}\in\mathrm{arg}\min_{x\in\mathcal{X}_{t}}\max_{y\in\mathcal{Y}_{t}}\underset{w\sim
D_{t}(x_{t}^{\ell},y_{t}^{\ell})}{\operatorname{\mathbb{E}}}[f_{t}(x,y,w)],$
(4) $\displaystyle
y_{t}^{\ell+1}\in\mathrm{arg}\max_{y\in\mathcal{Y}_{t}}\min_{x\in\mathcal{X}_{t}}\underset{w\sim
D_{t}(x_{t}^{\ell},y_{t}^{\ell})}{\operatorname{\mathbb{E}}}[f_{t}(x,y,w)].$

The fixed points of this repeated retraining procedure have been coined
equilibrium points, and are known to exist under mild conditions. In what
follows we provide algorithms capable of tracking the equilibrium point
trajectory $\\{\bar{z}_{t}\\}=\\{\bar{x}_{t},\bar{y}_{t}\\}$ without requiring
that we take the sequences in LABEL:eqn:repeated_retraining to convergence
($\ell\to\infty)$. This will be crucial for our online setting, as we assume
that each time $t$, a new function and distributional map arrive ([3, 6, 8,
13, 28, 34]).

###  1.1 Related Work

Stochastic Saddle Point Problems. Algorithms for computing saddle points can
be loosely catagorized as primal-dual based or proximal based [14, 18, 21, 22,
37]. Some works seek to find approximate saddle points by analyzing a saddle
point gap. We conduct our analysis in a setting in which solutions are known
to be unique, so we simply track them. Our analysis is primarily conducted
through variational analysis [26]. Hence we define the appropriate gradient
maps and demonstrate that solutions to our problems are the solution to the
variational inequalities induced by said gradient maps.

Decision Dependent Distributions. This work is most closely related to the
literature on stochastic optimization with decision dependent distributions,
or its counterpart in learning, performative prediction. The problem of
finding optimal decisions that are robust to decision-dependent data has been
studied extensively, and in many distinct settings: minimization problems [9,
23], saddle-point problems [35], games [20], online [34], time-varying decay
[25]. Relative to the existing work on saddle point problems in the
literature, this work considers problems for which the objective, constraints
and distribution are time-varying and revealed sequentially in time. The work
on games is related, as specific instances of games such as two-player zero-
sum may be cast into a saddle point problem. Saddle point problems however are
however not a strict subset of games as they exist in their own right; arising
from constrained minimization problems, etc. The most obvious inspiration for
this work is that of [35], as the setting in this work is precisely stochastic
saddle point problems with decision-dependent distributions for time-invariant
problems. Relative to this work, the results we present here are extensions to
the online setting where analysis requires handling of additional noise due to
solution drift.

Online Convex Optimization. Relevant to works on online optimization that are
concerned with tracking trajectories (see the representative works [7, 19, 16,
24, 27]) or given comparator sequences [13]; another line of work is concerned
with finding a sequence of decision that minimize a suitable dynamic regret
metric [12]. Our metric is the distance to the solution of (1) at the current
iteration, where we incorporate the drift of the solution trajectory. We
account for the time-variability of the solution by incorporating the solution
drift into our guarantees.

###  1.2 Contributions

Our contributions are as follows.

_(c1)_ The Online Equilibrium Problem. We propose a notion of equilibrium
points for the time-varying saddle-point problem in 1, provide conditions to
guarantee existence and uniqueness, and provide bounds for the distance
between the unique equilibrium points (1).

_(c2)_ Online Algorithms. We demonstrate that primal-dual algorithms, using
the gradients of $f_{t}$, are effective at finding equilibrium points when the
stochastic objective $f_{t}$ is strongly-convex-strongly-concave for any
realization of $w$. First, we demonstrate effective tracking of a conceptual
algorithm using full gradient information. We then demonstrate that a
stochastic algorithm tracks equilibrium points with additional noise due to
estimation. Furthermore, we provide provide expectation bounds and high
probability bounds that hold for each iteration.

_(c3)_ Experiments. We illustrate our results on the electric vehicle charging
problem in (2) by incorporating synthetic demand data from [10]. Here, the
demand changes in response to prices with a location-scale family based
distributional map.

##  2 Equilibrium Points

In this section we define the equilibrium problem, the fixed points of the
repeated retraining heuristic in LABEL:eqn:repeated_retraining, and provide
sufficient conditions for their existence. We start from the definition of
equilibrium points.

######  Definition 1.

(Equilibrium Points) A pair
$(\bar{x}_{t},\bar{y}_{t})\in\mathcal{X}\times\mathcal{Y}$ is an equilibrium
point if:

$\displaystyle\bar{x}_{t}\in\arg\min_{x\in\mathcal{X}_{t}}\left\\{\max_{y\in\mathcal{Y}_{t}}\underset{w\sim
D_{t}(\bar{x}_{t},\bar{y}_{t})}{\operatorname{\mathbb{E}}}[f_{t}(x,y,w)]\right\\},$
(5)
$\displaystyle\bar{y}_{t}\in\arg\max_{y\in\mathcal{Y}_{t}}\left\\{\min_{x\in\mathcal{X}_{t}}\underset{w\sim
D_{t}(\bar{x}_{t},\bar{y}_{t})}{\operatorname{\mathbb{E}}}[f_{t}(x,y,w)]\right\\}.$

Sequences of equilibrium points are defined as
$(\bar{x}_{t},\bar{y}_{t})_{t\in\mathbb{N}}$. $\Box$

In essence, equilibrium points are the solutions to the stationary saddle
point problem that they induce. In this way, they are optimal decisions when
data distribution is in state $D_{t}(\bar{x}_{t},\bar{y}_{t})$ but need not be
optimal otherwise. Existence of these points is contingent on the
distributional function being continuous on the set of probability
distributions, and $f_{t}$ being at least convex-concave.

######  Theorem 2.1.

(Existence of Equilibrium Points) Suppose that the following assumptions hold
at time $t\geq 0$:

i) $x\mapsto f_{t}(x,y,w)$ is convex in $x$ for all $y\in\mathcal{Y}_{t}$ and
for any realization of $w$;

ii) $y\mapsto f_{t}(x,y,w)$ is concave in $y$ for all $x\in\mathcal{X}_{t}$
and for any realization of $w$;

iii) $(x,y)\mapsto f_{t}(x,y,w)$ is continuous on
$\mathcal{X}\times\mathcal{Y}$ for any given $w$;

iv) the sets
$\mathcal{X}_{t}\subset\mathbb{R}^{n},\mathcal{Y}_{t}\subset\mathbb{R}^{n}$
are convex compact subsets;

v) the distributional map $D_{t}:\mathcal{Z}\to(\mathcal{P}(M),W_{1})$ is
continuous.

Then the set of equilibrium points is nonempty and compact. $\Box$

The proof follows from the fact that equilibrium points exist for each problem
at time $t$ due to [35, Theorem 2.10]. The proof strategy amounts to
demonstrating that the repeated retraining map satisfies Kakutani’s Fixed
Point Theorem; see, e.g., [1, Corollary 17.55]; it is not provided due to
space limitations.

###  2.1 Theoretical Framework

In light of our discussion on existence of equilibrium points, we outline the
assumptions and some results that will be necessary later in our analysis. For
notational convenience, we will refer to the stacked variable $z=(x,y)$ and
the Cartesian product set of constraints
$\mathcal{Z}=\mathcal{X}\times\mathcal{Y}$. We will rely on the following
assumptions to hold at each time $t$ throughout this work.

######  Assumption 1.

(Strong-Convexity-Strong-Concavity) The function $(x,y)\mapsto f_{t}(x,y,w)$
is continuously differentiable over $\mathbb{R}^{n}\times\mathbb{R}^{m}$ for
any realization of $w$. The function $(x,y)\mapsto f_{t}(x,y,w)$ is
$\gamma$-strongly-convex-strongly-concave, for any realization of $w$; that
is, $f_{t}$ is $\gamma$-strongly-convex in $x$ for all $y\in\mathbb{R}^{m}$
and $\gamma$-strongly-concave in $y$ for all $x\in\mathbb{R}^{n}$. $\Box$

######  Assumption 2.

(Joint Smoothness) The map
$g_{t}(z,w):=(\nabla_{x}f_{t}(z,w),-\nabla_{y}f_{t}(z,w))$ is $L$-Lipschitz in
$z$ and $w$. Namely, $\|g_{t}(z,w)-g_{t}(z^{\prime},w)\|\leq
L\|z-z^{\prime}\|$, $\|g_{t}(z,w)-g_{t}(z,w^{\prime})\|\leq
L\operatorname{d}(w,w^{\prime})$, for any
$z,z^{\prime}\in\mathbb{R}^{n}\times\mathbb{R}^{m}$ and $w,w^{\prime}$
supported on $\mathbb{R}^{k}$, for some $L\geq 0$, where
$\operatorname{d}:\mathbb{R}^{k}\times\mathbb{R}^{k}\rightarrow\mathbb{R}$ is
some chosen metric on $\mathbb{R}^{k}$. $\Box$

######  Assumption 3.

(Lipschitz-Continuous Distributional Map) The distributional maps
$D_{t}:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathcal{P}(M)$ are
$\varepsilon$-Lipschitz. Namely,
$W_{1}(D_{t}(z),D_{t}(z^{\prime}))\leq\varepsilon\|z-z^{\prime}\|$, for any
$z,z^{\prime}\in\mathbb{R}^{n}\times\mathbb{R}^{m}$, where $W_{1}$ is the
Wasserstein-1 distance. $\Box$

######  Assumption 4.

(Compact Sets) The sets $\mathcal{X}_{t}\subset\mathbb{R}^{n}$ and
$\mathcal{Y}_{t}\subset\mathbb{R}^{m}$ are compact and convex. $\Box$

######  Assumption 5.

(Bounded Drift) There exists a $\Delta>0$ such that the equilibrium drift
sequence defined by $\Delta_{t}:=\|\bar{z}_{t+1}-\bar{z}_{t}\|$ is uniformly
bounded by $\Delta$. Namely, $\Delta_{t}\leq\Delta$ for all $t\geq 0$. $\Box$

These assumptions provided are sufficient to guarantee uniqueness of the
equilibrium point, and convergence of primal-dual algorithms in the batch
setting; see [35].

######  Theorem 2.2.

(Equilibrium Point Uniqueness) If Assumptions 1-4 are satisfied such that
$\varepsilon L<\gamma$, then a unique equilibrium point exists.

Proof of this results amounts to showing that the repeated retraining
heuristic in LABEL:eqn:repeated_retraining is a strict contraction and hence
satisfies the Banach-Picard Fixed Point Theorem. For a detailed proof, see
[35].

Given that the data distribution is shifting, it is necessary to characterize
this shift and its effect on the gradient. The key to computing equilibrium
points will be the gradients of $f_{t}$. We note that this is only one term
required to compute the gradients of $F_{t}$, effectively ignoring the
dependence of $D_{t}$ on the decision variables. For now, we will denote the
decoupled gradient map as the function $G_{t}$ defined by

$G_{t}(z;z^{\prime}):=\underset{w\sim
D_{t}(z^{\prime})}{\operatorname{\mathbb{E}}}g_{t}(z,w)=\left(\underset{w\sim
D_{t}(z^{\prime})}{\operatorname{\mathbb{E}}}\nabla_{x}f_{t}(z,w),\underset{w\sim
D_{t}(z^{\prime})}{\operatorname{\mathbb{E}}}-\nabla_{y}f_{t}(z,w)\right)$ (6)

for all $z,z^{\prime}\in\mathbb{R}^{n}\times\mathbb{R}^{m}$. Note that we
refer to this gradient map as “decoupled” as we separate the decision variable
in the stochastic objective and the distributional map. This will allow us to
characterize these behaviors separately.

######  Lemma 2.3.

(Gradient Map Characterization) If Assumptions 1-4 hold, then:

  1. 1.

(Gradient Deviation) For any fixed $\hat{z}\in\mathcal{Z}$, the map $z\mapsto
G_{t}(\hat{z},z)$ is $\varepsilon L$-Lipschitz-continuous. That is,
$\|G_{t}(\hat{z};z)-G_{t}(\hat{z};z^{\prime})\|\leq\varepsilon
L\|z-z^{\prime}\|$, for all $z,z^{\prime}\in\mathcal{Z}$.

  2. 2.

(Strong-Monotonicity) The map $z\mapsto G_{t}(z,z)$ is $(\gamma-\varepsilon
L)$-strongly-monotonic.

  3. 3.

(Lipschitz-Continuity) The map $z\mapsto G_{t}(z,z)$ is $(L+\varepsilon
L)$-Lipschitz Continuous. $\Box$

Proof of the Gradient Deviation property follows by combining the properties
allowed from joint smoothness and lipschitz continuity of the distributional
map (Assumptions 2 and 3 respectively. For a detailed proof, we refer the
reader to [35]. Strong monotonicity and Lipschitz continuity of $z\mapsto
G_{t}(z;z)$ then follow immediately. With this lemma, we can effectively deal
with the decoupled gradient map by passing variables into both the $D_{t}$ and
$g_{t}$ simultaneously. Going forward, we will simply write $G_{t}$ to mean
the gradient map given by $z\mapsto G_{t}(z;z)$.

##  3 Online Algorithms

###  3.1 A Conceptual Primal-Dual Algorithm

In this section, we show that if the decoupled gradient map $G_{t}$ is
available, then tracking the equilibrium points is possible using a primal-
dual algorithm. This provides a basis of comparison for our analysis in the
next section where we use a stochastic gradient estimator in place of $G_{t}$.
Indeed, we denote the equilibrium primal-dual algorithmic map by

$\mathcal{G}_{t}(z)=\Pi_{\mathcal{Z}}\left(z-\eta G_{t}(z)\right)$ (7)

so that the algorithm generates the sequence $\\{z_{t}\\}_{t\geq 0}$ defined
by $z_{t+1}=\mathcal{G}_{t}(z_{t})$, $t\in\mathbb{N}$. To proceed, we observe
that equilibrium points are the fixed points of the primal-dual algorithmic
map.

######  Proposition 3.1.

(Fixed Point Characterization) Let Assumptions 1-4 hold and suppose that
$\frac{\varepsilon L}{\gamma}<1$. A point $\bar{z}_{t}\in\mathcal{Z}$ is an
equilibrium point if and only if $\bar{z}_{t}=\mathcal{G}_{t}(\bar{z}_{t})$.
$\Box$

This proposition will allow us to cast our analysis into a fixed point
framework, using the equilibrium points as the fixed points of the
distributional map.

######  Theorem 3.2.

(Primal-Dual Tracking) Suppose that Assumptions 1-4 hold and that
$\frac{\varepsilon L}{\gamma}<1$. Then the sequence
$z_{t+1}=\mathcal{G}_{t}(z_{t})$ satisfies the bound

$\|z_{t}-\bar{z}_{t}\|\leq\alpha^{t}\|z_{0}-\bar{z}_{0}\|+(1-\alpha)^{-1}\Delta$
(8)

for any initial point $z_{0}\in\mathcal{Z}$, and
$\alpha:=\sqrt{1-\eta(\gamma-\varepsilon L)}$ provided that

$\eta<\min\left\\{\frac{1}{\gamma-\varepsilon L},\frac{\gamma-\varepsilon
L}{(1+\varepsilon)^{2}L^{2}}\right\\}$ (9)

Furthermore, $\\{z_{t}\\}_{t\geq 0}$ ultimately tracks the sequence of unique
equilibrium points $\\{\bar{z}_{t}\\}_{t\geq 0}$ in the sense that
$\limsup_{t\to\infty}\|z_{t}-\bar{z}_{t}\|\leq(1-\alpha)^{-1}\Delta$. $\Box$

###### Proof.

It follows from the triangle inequality that
$\|z_{t+1}-\bar{z}_{t+1}\|\leq\|z_{t+1}-\bar{z}_{t}\|+\|\bar{z}_{t}-\bar{z}_{t+1}\|=\|z_{t+1}-\bar{z}_{t}\|+\Delta_{t}$,
and hence we simply need to bound $\|z_{t+1}-\bar{z}_{t}\|$. We observe that

$\displaystyle\|z_{t+1}-\bar{z}_{t}\|^{2}$
$\displaystyle=\|\Pi_{\mathcal{Z}}\left(z_{t}-\eta
G_{t}(z_{t})\right)-\Pi_{\mathcal{Z}}\left(\bar{z}_{t}-\eta
G_{t}(\bar{z}_{t})\right)\|^{2}$
$\displaystyle\leq\|(z_{t}-\bar{z}_{t})-\eta\left(G_{t}(z_{t})-G_{t}(\bar{z}_{t})\right)\|^{2}$
$\displaystyle\leq\|z_{t}-\bar{z}_{t}\|^{2}-2\eta\langle
z_{t}-\bar{z}_{t},G_{t}(z_{t})-G_{t}(\bar{z}_{t})\rangle+\eta^{2}\|G_{t}(z_{t})-G_{t}(\bar{z}_{t})\|^{2}.$

If we denote $\hat{\gamma}=\gamma-\varepsilon L$ and $\hat{L}=L+\varepsilon
L$, then from Lemma 2.3 we have that $G_{t}$ is $\hat{\gamma}$-strongly
monotone and $\hat{L}$-Lipschitz continuous. Combining these facts yields

$\langle
z_{t}-\bar{z}_{t},G_{t}(z_{t})-G_{t}(\bar{z}_{t})\rangle\geq\frac{\hat{\gamma}}{2}\|z_{t}-\bar{z}_{t}\|^{2}+\frac{\hat{\gamma}}{2\hat{L}^{2}}\|G_{t}(z_{t})-G_{t}(\bar{z}_{t})\|^{2}.$

Substituting into the above yields

$\displaystyle\|z_{t+1}-\bar{z}_{t}\|^{2}$
$\displaystyle\leq(1-\eta\hat{\gamma})\|z_{t}-\bar{z}_{t}\|^{2}+\eta\left(\eta-\frac{\hat{\gamma}}{\hat{L}^{2}}\right)\|G_{t}(z_{t})-G_{t}(\bar{z}_{t})\|^{2}\leq(1-\eta\hat{\gamma})\|z_{t}-\bar{z}_{t}\|^{2}$

where the last inequality follows provided that
$\eta\leq\hat{\gamma}/\hat{L}^{2}$. It follows that if $\eta<1/\hat{\gamma}$
as well, then $1-\eta\hat{\gamma}<1$ and the bound in Theorem (8) follows.
Considering the limit supremum of the bound in (8) yields the result. ∎

We note that the noise due to the drift in (8) increases as we decrease the
step size $\eta$. Hence it is impossible to completely remove this disturbance
from the algorithm. This reflects intuition however as very small step sizes
would make it difficult to ever reach the solution trajectory. Meanwhile,
larger step sizes decrease this noise while simultaneously decreasing the rate
at which we overcome the error $z_{t+1}-\bar{z}_{t}$ between successive
iterates. We build on this intuition in our stochastic algorithm. This
concludes our discussion of the conceptual primal-dual algorithm. In the next
section, we demonstrate tracking of a stochastic primal-dual algorithm.

###  3.2 A Stochastic Primal-Dual Algorithm

In previous section, we demonstrated that a conceptual first-order algorithm
is capable of tracking the trajectory of equilibrium point. We say conceptual
because having access to full information in $G_{t}$ requires the ability to
compute the expectation with respect to the distributional map $D_{t}$ at each
algorithmic step—which is of course impractical. Hence, we are concerned with
a more pragmatic setting in which we merely have access to a stochastic
gradient oracle, which we will denote $H_{t}$. We make the implicit assumption
throughout that $H_{t}$ is a function of the stochastic gradient function
$g_{t}$ defined in (6). Such functions are typically of the form

$H_{t}(z)=\begin{cases}g_{t}(z,w_{1}),\quad w_{1}\sim D_{t}(z),\\\
\frac{1}{N}\sum_{i=1}^{N}g_{t}(z,w_{i}),\quad
w_{1},\ldots,w_{N}\overset{i.i.d.}{\sim}D_{t}(z).\\\ \end{cases}$ (10)

Then, given a starting point $z_{0}$, the stochastic primal-dual algorithm
performs the update

$z_{t+1}=\hat{\mathcal{G}}_{t}(z_{t}),\ \text{where}\
\hat{\mathcal{G}}_{t}(z_{t})=\Pi_{\mathcal{Z}}\left(z_{t}-\eta
H_{t}(z_{t})\right)$ (11)

Crucial to our analysis will be providing reasonable assumptions regarding the
quality of the gradient estimator $H_{t}$. The case where
$H_{t}(z)=g_{t}(z,w_{1})$ is particularly appealing in applications such as
competitive markets, strategic classification, etc., where $g_{t}(z,w_{1})$
can be computed using an observation of $w$ (in our example in competitive
markets, we would observe the demands $a$ and $b$).

We are interested in providing results in expectation as well as high-
probability. A common assumption throughout the literature is to use a sub-
Gaussian error model on this gradient error quantity—an observation supported
by the central limit theorem when using a sufficiently large batch size $N$ in
(10). While this may hold in some cases, it has been observed that this
requires a prohibitively large set of data while also assuming the data is of
sufficiently good quality [32]; this has also been observed in works on
stochastic gradient methods such as in [11, 30]. A more mild assumption then
is to assume a larger class of heavy-tailed distributions known as sub-Weibull
distributions, which we formalize in the following.

######  Definition 2.

(Sub-Weibull Random Variable [32]) The distribution of a random variable $\xi$
is sub-Weibull, denoted $\xi\sim\operatorname{subW}(\theta,\nu)$, if there
exists $\theta>0,\nu>0$ such that $\|z\|_{p}\leq\nu p^{\theta}$, for all
$p\geq 1$. $\Box$

######  Assumption 6.

(Stochastic Framework) Denote the gradient error incurred throughout the
stochastic algorithm as $\xi_{t}=H_{t}(z_{t})-G_{t}(z_{t})$. Then there exists
constant $\theta,\nu>0$ and a sequence $\\{\nu_{t}\\}_{t\geq
0}\subseteq\mathbb{R}_{+}$ such that the following hold:

  1. 1.

Sub-Weibull Gradient Error. For each $t\geq 0$, $\|\xi_{t}\|$ is a sub-Weibull
random variable such that
$\|\xi_{t}\|\sim\operatorname{subW}(\theta,\nu_{t})$.

  2. 2.

Bounded Variance Proxies. The sequence of variance proxies
$\\{\nu_{t}\\}_{t\geq 0}$ is bounded by $\nu$. $\Box$

With this assumption, the main convergence result is stated next.

######  Theorem 3.3.

Suppose that Assumptions 1-6 hold and $\frac{\varepsilon L}{\gamma}<1$. If
$\eta$ satisfies the bound in 9 then the following hold:

  1. 1.

Expectation. The sequence $\\{z_{t}\\}_{t\geq 0}$ satisfies the bound in
expectation

$\operatorname{\mathbb{E}}[\|z_{t}-\bar{z}_{t}\|]\leq\alpha^{t}\|z_{0}-\bar{z}_{0}\|+(1-\alpha)^{-1}\Delta+(1-\alpha)^{-1}\eta\nu.$
(12)

for all $t\geq 0$, for any initial point $z_{0}\in\mathcal{Z}$, and
$\alpha:=\sqrt{1-\eta(\gamma-\varepsilon L)}$.

  2. 2.

High Probability. For any $\delta\in(0,1)$, and $t\geq 0$,

$\mathbb{P}\left(\|z_{t}-\bar{z}_{t}\|\leq\alpha^{t}\|z_{0}-\bar{z}_{0}\|+\frac{\Delta}{1-\alpha}+c(\theta)\log^{\theta}\left(\frac{2}{\delta}\right)\frac{\eta\nu}{1-\alpha}\right)\geq
1-\delta\,.$ (13)

with $c(\theta):=\left(\frac{2e}{\theta}\right)^{\theta}$, for any initial
point $z_{0}\in\mathcal{Z}$. $\Box$

We note that the noise terms above are diametrically opposed functions of the
step-size. While the drift term grows larger for small step-size, the gradient
noise decreases for smaller step-size values. This relationship makes removing
the contribution of any one source of perturbation impossible. We also note
that the high-probability bound scales as $\log(\delta^{-1})$, as opposed to
classical bounds derived using Markov’s bound that scale as $\delta^{-1}$.

Before proving the theorem, we provide supporting lemmas that will be used in
the proof.

######  Lemma 3.4.

(Equivalent Characterizations) If $\xi$ is a sub-Weibull random variable with
tail parameter $\theta>0$, then the following characterizations are equivalent
(we recall that $\|z\|_{k}=\operatorname{\mathbb{E}}[|z|^{k}]^{1/k}$):

  1. (c1)

Tail Probability: $\exists\ \nu_{1}>0$ such that
$\mathbb{P}(|z|\geq\epsilon)\leq
2\exp(-\left(\epsilon/\nu_{1})\right)^{1/\theta}$ for all $\epsilon>0$.

  2. (c2)

Moment: $\exists\ \nu_{2}>0$ such that $\|z\|_{k}\leq\nu_{2}k^{\theta}$ for
all $k\geq 1$.

Moreover, if (c2) holds for a given $\nu_{2}>0$, then (c1) holds with
$\nu_{1}=\left(\frac{2e}{\theta}\right)^{\theta}\nu_{2}$. $\Box$

######  Lemma 3.5.

(Sub-Weibull Inclusion) If $\xi\sim\operatorname{subW}(\theta,\nu)$ based on
(c2) and $\theta^{\prime},\nu^{\prime}>0$ such that
$\theta\leq\theta^{\prime}$ and $\nu\leq\nu^{\prime}$ then
$\xi\sim\operatorname{subW}(\theta^{\prime},\nu^{\prime})$. $\Box$

######  Lemma 3.6.

(Sub-Weibull Closure) If $\xi_{1}\sim\operatorname{subW}(\theta_{1},\nu_{1})$,
$\xi_{2}\sim\operatorname{subW}(\theta_{2},\nu_{2})$ are (possibly coupled)
sub-Weibull random variables based on (c2) and $c\in\mathbb{R}$, then the
following hold:

  1. 1.

$\xi_{1}+\xi_{2}\sim\operatorname{subW}(\max\\{\theta_{1},\theta_{2}\\},\nu_{1}+\nu_{2})$;

  2. 2.

$\xi_{1}\xi_{2}\sim\operatorname{subW}(\theta_{1}+\theta_{2},\psi(\theta_{1},\theta_{2})\nu_{1}\nu_{2})$,
$\psi(\theta_{1},\theta_{2}):=(\theta_{1}+\theta_{2})^{\theta_{1}+\theta_{2}}/(\theta_{1}^{\theta_{1}}\theta_{2}^{\theta_{2}})$;

  3. 3.

$c\xi_{1}\sim\operatorname{subW}(\theta_{1},|c|\nu_{1})$. $\Box$

The proofs of these lemmas can be found in [32, 33].

###### Proof.

_of 3.3_. As before, we have that
$\|z_{t+1}-\bar{z}_{t+1}\|\leq\|z_{t+1}-\bar{z}_{t}\|+\Delta_{t}$ where

$\displaystyle\|z_{t+1}-\bar{z}_{t}\|$
$\displaystyle\leq\|(z_{t}-\bar{z}_{t})-\eta(H_{t}(z_{t})-G_{t}(\bar{z}_{t})\|$
$\displaystyle=\|(z_{t}-\bar{z}_{t})-\eta(G_{t}(z_{t})-G_{t}(\bar{z}_{t})-\eta\xi_{t}\|\leq\alpha\|z_{t}-\bar{z}_{t}\|+\eta\|\xi_{t}\|.$

This yields that stochastic recursion
$\|z_{t}-\bar{z}_{t}\|\leq\alpha^{t}\|z_{0}-\bar{z}_{0}\|+\Delta\sum_{i=0}^{t}\alpha^{i}+\eta\sum_{i=0}^{t}\alpha^{i}\|\xi_{t-i}\|$.
Recall that when $\eta$ satisfies the condition in (9), $\alpha<1$. Hence
assuming this fact and taking the expectation of both sides yields

$\operatorname{\mathbb{E}}\|z_{t}-\bar{z}_{t}\|\leq\alpha^{t}\|z_{0}-\bar{z}_{0}\|+\frac{\Delta}{1-\alpha}+\eta\sum_{i=0}^{t}\alpha^{i}\operatorname{\mathbb{E}}\|\xi_{t-i}\|$

so that the result in (12) follows. To prove the result in (13), we denote
$e_{t}=\|z_{t}-\bar{z}_{t}\|$,
$\omega_{t}=\alpha_{t}\|z_{0}-\bar{z}_{0}\|+\Delta(1-\alpha)^{-1}$, and
$\sigma_{t}=\eta\sum_{i=0}^{t}\alpha^{i}\xi_{t-i}$. Observe that, due to our
closure properties,

$\|\sigma_{t}\|_{p}\leq\sum_{i=0}^{t}\alpha^{i}\operatorname{\mathbb{E}}[\|\xi_{t}\|]^{p}]^{1/p}\leq\frac{\eta\nu}{1-\alpha}p^{\theta}$

for any $p\geq 1$ and hence
$\sigma_{t}\sim\operatorname{subW}(\theta,\eta\nu(1-\alpha)^{-1})$. It follows
from Lemma (c1) (c1) that

$\mathbb{P}\left(\sigma_{t}\geq\epsilon\right)\leq
2\exp\left(-\frac{\theta}{2e}\left(\frac{(1-\alpha)\epsilon}{\eta\nu}\right)^{\frac{1}{\theta}}\right),$
(14)

and setting the right hand side above equal to $\delta>0$ yields
$\epsilon=c(\theta)\log^{\theta}\left(\frac{2}{\delta}\right)\eta\nu(1-\alpha)^{-1}$.
Now, observe that our stochastic recursion implies that for any $a>0$,
$\mathbb{P}(\omega_{t}+\sigma_{t}\geq a)\geq\mathbb{P}(e_{t}\geq a)$. It
follows that setting $a=\omega_{t}+\epsilon$ yields
$\mathbb{P}(e_{t}\leq\omega_{t}+\epsilon)\geq\mathbb{P}(\omega_{t}+\sigma_{t}\leq\omega_{t}+\epsilon)=\mathbb{P}(\sigma_{t}\leq\epsilon)\geq
1-\delta$, thus the result follows. ∎

##  4 Numerical Simulations on Electric Vehicle Charging

(a)

(b)

Figure 1: Demand time series visualization: horizontal axis is time of day,
vertical axis is the day of the year between 1 and 365. Brightness indicates
intensity of the demand.

(a)

(b)

Figure 2: Results: in (a) we depict the evolution of the equilibrium points
over the time horizon plotted in absolute value. In (b), we depict the
tracking error for both algorithms.

In this section, we provide a demonstration of an online electric vehicle
charging market, as described in (2), with time series demand data from [10].
The data describes a years of worth of electricity demand with entries for
each minute of the year. Each file represents a different type of charging
station positioned near commercial uses with varying number of ports (2 or 6),
frequency of use (2, 8, or 16 event), and port power output (50, 150, or 350
kW). We randomly allocate each provider with three $8$-event stations and draw
samples from each day of the year. The demand data is normalized by first
subtracting the mean across each minute and dividing by the variance. A
representative example of the raw data is provided in Figure 1, with time in
minutes along the horizontal axis, day of the year along the vertical, and
color intensity representing demand value. The price elasticity is dictated by
the function $h_{t}(p)=(-c(p)/m|t-m|+c(p))$ where $p$ is the station’s port
power and $c(p)$ is given by $c(p)=0.3$ for $p\in\\{50,150\\}$ and $c(p)=0.5$
for $p=350$. The elasticity matrices are then given by
$(A_{1}^{t})_{ij}=-h_{t}(p_{i})\delta_{i,j}$,
$(B_{1}^{t})_{ij}=-h_{t}(p_{i})\delta_{i,j}$ for $i\in[3]$ where $p_{i}$ is
the power of each port at the $i$th station belonging to the provider and
$B_{2}^{t}=-B_{1}^{t}$ and $A_{2}^{t}=-A_{1}^{t}$. For the sake of simplicity,
we consider service providers with unit charging speed utility rates and zero
location-based utility. From this we conclude that for all $t$, $G_{t}$ is
$1$-strongly monotone and $1$-Lipschitz. Hence our results apply provided that
$\eta<1/3$.

We compute the equilibrium points by executing a batch primal-dual algorithm
for 2000 iterations with a step size of $\eta=0.001$. We then run the online
primal-dual and stochastic primal-dual algorithms over each minute of the time
series data and plot the distance to the solutions in Figure 2. We observe
that the primal-dual algorithm is capable of reasonably tracking the
trajectory. The noise incurred by the stochastic algorithm clearly prevents it
from having identical performance, however the trajectory does decrease to an
acceptable level after overcoming transient behavior for approximately 200
time steps.

## References

  * [1] Charalambos D Aliprantis and Kim Border  “Infinite dimensional analysis: A Hitchhiker’s Guide”  Springer, 2006
  * [2] Julian Berberich, Johannes Köhler, Matthias A Müller and Frank Allgöwer  “Data-driven model predictive control with stability and robustness guarantees”  In _IEEE Transactions on Automatic Control_ 66.4 IEEE, 2020, pp. 1702–1717
  * [3] Omar Besbes, Yonatan Gur and Assaf Zeevi  “Non-stationary stochastic optimization”  In _Operations research_ 63.5 INFORMS, 2015, pp. 1227–1244
  * [4] Gianluca Bianchin, Miguel Vaquero, Jorge Cortes and Emiliano Dall’Anese  “Online stochastic optimization for unknown linear systems: Data-driven synthesis and controller analysis”  In _arXiv preprint arXiv:2108.13040_ , 2021
  * [5] John R Birge and Francois Louveaux  “Introduction to stochastic programming”  Springer Science & Business Media, 2011
  * [6] Xuanyu Cao, Junshan Zhang and H Vincent Poor  “Online stochastic optimization with time-varying distributions”  In _IEEE Transactions on Automatic Control_ 66.4 IEEE, 2020, pp. 1840–1847
  * [7] Joshua Cutler, Dmitriy Drusvyatskiy and Zaid Harchaoui  “Stochastic optimization under time drift: iterate averaging, step-decay schedules, and high probability guarantees”  In _Advances in Neural Information Processing Systems_ 34, 2021
  * [8] Emiliano Dall’Anese, Andrea Simonetto, Stephen Becker and Liam Madden  “Optimization and learning with information streams: Time-varying algorithms and applications”  In _IEEE Signal Processing Magazine_ 37.3 IEEE, 2020, pp. 71–83
  * [9] Dmitriy Drusvyatskiy and Lin Xiao  “Stochastic optimization with decision-dependent distributions”  In _Mathematics of Operations Research_ INFORMS, 2022
  * [10] Madeline Gilleran et al.  “Impact of electric vehicle charging on the power demand of retail buildings”  In _Advances in Applied Energy_ 4 Elsevier, 2021, pp. 100062
  * [11] Mert Gürbüzbalaban, Umut Şimşekli and Lingjiong Zhu  “The heavy-tail phenomenon in SGD”  In _International Conference on Machine Learning (ICML)_ 139, 2021, pp. 3964–3975
  * [12] Elad Hazan  “Introduction to online convex optimization”  In _arXiv preprint arXiv:1909.05207_ , 2019
  * [13] Ali Jadbabaie, Alexander Rakhlin, Shahin Shahrampour and Karthik Sridharan  “Online optimization: Competing with dynamic comparators”  In _Artificial Intelligence and Statistics_ , 2015, pp. 398–406  PMLR
  * [14] Jayash Koshal, Angelia Nedić and Uday V Shanbhag  “Multiuser optimization: Distributed algorithms and error analysis”  In _SIAM Journal on Optimization_ 21.3 SIAM, 2011, pp. 1046–1081
  * [15] Dan Li, Dariush Fooladivanda and Sonia Martinez  “Online optimization and learning in uncertain dynamical environments with performance guarantees”  In _arXiv preprint arXiv:2102.09111_ , 2021
  * [16] Liam Madden, Stephen Becker and Emiliano Dall’Anese  “Bounds for the tracking error of first-order online optimization methods”  In _Journal of Optimization Theory and Applications_ 189.2 Springer, 2021, pp. 437–457
  * [17] Chinmay Maheshwari et al.  “Zeroth-Order Methods for Convex-Concave Minmax Problems: Applications to Decision-Dependent Risk Minimization”  In _arXiv preprint arXiv:2106.09082_ , 2021
  * [18] Aryan Mokhtari, Asuman Ozdaglar and Sarath Pattathil  “A unified analysis of extra-gradient and optimistic gradient methods for saddle point problems: Proximal point approach”  In _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 1497–1507  PMLR
  * [19] Aryan Mokhtari, Shahin Shahrampour, Ali Jadbabaie and Alejandro Ribeiro  “Online optimization in dynamic environments: Improved regret rates for strongly convex problems”  In _IEEE Conference on Decision and Control_ , 2016, pp. 7195–7201
  * [20] Adhyyan Narang et al.  “Learning in Stochastic Monotone Games with Decision-Dependent Data”  In _International Conference on Artificial Intelligence and Statistics_ , 2022, pp. 5891–5912  PMLR
  * [21] Arkadi Nemirovski  “Prox-method with rate of convergence $O(1/t)$ for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems”  In _SIAM Journal on Optimization_ 15.1 SIAM, 2004, pp. 229–251
  * [22] Arkadi Nemirovski, Anatoli Juditsky, Guanghui Lan and Alexander Shapiro  “Robust stochastic approximation approach to stochastic programming”  In _SIAM Journal on optimization_ 19.4 SIAM, 2009, pp. 1574–1609
  * [23] Juan Perdomo, Tijana Zrnic, Celestine Mendler-Dünner and Moritz Hardt  “Performative prediction”  In _International Conference on Machine Learning_ , 2020, pp. 7599–7609  PMLR
  * [24] A Yu Popkov  “Gradient methods for nonstationary unconstrained optimization problems”  In _Automation and Remote Control_ 66.6 Springer, 2005, pp. 883–891
  * [25] Mitas Ray, Lillian J Ratliff, Dmitriy Drusvyatskiy and Maryam Fazel  “Decision-dependent risk minimization in geometrically decaying dynamic environments”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , 2022, pp. 8081–8088
  * [26] R Tyrrell Rockafellar and Roger J-B Wets  “Variational analysis”  Springer Science & Business Media, 2009
  * [27] Daniel D. Selvaratnam, Iman Shames, Jonathan H. Manton and Mohammad Zamani  “Numerical Optimisation of Time-Varying Strongly Convex Functions Subject to Time-Varying Constraints”  In _IEEE Conference on Decision and Control_ , 2018, pp. 849–854
  * [28] Iman Shames and Farhad Farokhi  “Online stochastic convex optimization: Wasserstein distance variation”  In _arXiv preprint arXiv:2006.01397_ , 2020
  * [29] Alexander Shapiro and Arkadi Nemirovski  “On complexity of stochastic programming problems”  In _Continuous optimization_ Springer, 2005, pp. 111–146
  * [30] Umut Şimşekli, Levent Sagun and Mert Gürbüzbalaban  “A tail-index analysis of stochastic gradient noise in deep neural networks”  In _International Conference on Machine Learning (ICML)_ 97, 2019, pp. 5827–5837
  * [31] Berkay Turan and Mahnoosh Alizadeh  “Competition in electric autonomous mobility on demand systems”  In _IEEE Transactions on Control of Network Systems_ IEEE, 2021
  * [32] Mariia Vladimirova, Stéphane Girard, Hien Nguyen and Julyan Arbel  “Sub-Weibull distributions: Generalizing sub-Gaussian and sub-Exponential properties to heavier tailed distributions”  In _Stat_ 9.1 Wiley Online Library, 2020, pp. e318
  * [33] Kam Chung Wong, Zifan Li and Ambuj Tewari  “Lasso guarantees for $\beta$-mixing heavy-tailed time series”  In _The Annals of Statistics_ 48.2 Institute of Mathematical Statistics, 2020, pp. 1124–1142
  * [34] Killian Wood, Gianluca Bianchin and Emiliano Dall’Anese  “Online projected gradient descent for stochastic optimization with decision-dependent distributions”  In _IEEE Control Systems Letters_ 6 IEEE, 2021, pp. 1646–1651
  * [35] Killian Wood and Emiliano Dall’Anese  “Stochastic Saddle Point Problems with Decision-Dependent Distributions”  In _arXiv preprint arXiv:2201.02313_ , 2022
  * [36] Junyu Zhang, Mingyi Hong, Mengdi Wang and Shuzhong Zhang  “Generalization bounds for stochastic saddle point problems”  In _International Conference on Artificial Intelligence and Statistics_ , 2021, pp. 568–576  PMLR
  * [37] Xuan Zhang, Necdet Serhat Aybat and Mert Gürbüzbalaban  “Robust accelerated primal-dual methods for computing saddle points”  In _arXiv preprint arXiv:2111.12743_ , 2021