# On the role of ${J/\psi}$ production in electron-ion collisions

Zexuan Chu  Key Laboratory of Nuclear Physics and Ion-beam Application (MOE)
and Institute of Modern Physics, Fudan University, Shanghai 200433, China
Jinhui Chen<EMAIL_ADDRESS>Key Laboratory of Nuclear Physics and
Ion-beam Application (MOE) and Institute of Modern Physics, Fudan University,
Shanghai 200433, China  Shanghai Research Center for Theoretical Nuclear
Physics, NSFC and Fudan University, Shanghai 200438, China  Xiang-Peng Wang
<EMAIL_ADDRESS>Technical University of Munich, TUM School of Natural
Sciences, Physics Department, James-Franck-Strasse 1, 85748 Garching, Germany
Hongxi Xing<EMAIL_ADDRESS>Key Laboratory of Atomic and Subatomic
Structure and Quantum Control (MOE), Guangdong Basic Research Center of
Excellence for Structure and Fundamental Interactions of Matter, Institute of
Quantum Matter, South China Normal University, Guangzhou 510006, China
Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Guangdong Provincial
Key Laboratory of Nuclear Science, Southern Nuclear Science Computing Center,
South China Normal University, Guangzhou 510006, China  Southern Center for
Nuclear-Science Theory (SCNT), Institute of Modern Physics, Chinese Academy of
Sciences, Huizhou 516000, China

###### Abstract

Within the framework of non-relativistic QCD (NRQCD) effective field theory,
we study the leptoproduction of $J/{\psi}$ at next-to-leading order in
perturbative QCD for both unpolarized and polarized electron-ion collisions.
We demonstrate that the $J/{\psi}$-tagged deep inelastic scattering in the
future Electron-Ion Collider can be served as a golden channel for the reasons
including constraining NRQCD long distance matrix elements, probing the
nuclear gluon distribution functions, as well as investigating the gluon
helicity distribution inside a longitudinal polarized proton.

##  I Introduction

${J/\psi}$ production has been the focus of much theoretical and experimental
interest in the past decades. Understanding its production mechanism and using
it as a probe to QCD matter are among the most active and challenging subjects
in particle and nuclear physics Rothkopf:2019ipj ; Chapon:2020heu . Among
several extensively used models and effective field theories for the
description of heavy quarknoium production, the most sophisticated one in high
energy collisions is provided by non-relativistic QCD (NRQCD) effective field
theory Bodwin:1994jh , in which the ${J/\psi}$ hadronization mechanism is
encoded in the non-perturbative long distance matrix elements (LDMEs). In the
past decades, through explicit perturbative calculation of ${J/\psi}$
production at next-to-leading order (NLO) in electron-positron, proton-proton
and electron-proton collisions within the NRQCD factorization framework,
significant progress has been made to extract the LDMEs from world data, see
e.g. Bodwin:2014gia ; Butenschoen:2011yh ; Chao:2012iv ; Gong:2012ug ;
Bodwin:2015iua ; Brambilla:2022rjd ; Brambilla:2022ayc . However, no consensus
has been reached yet to identify the universalities of these LDMEs and to
understand the polarization of the produced ${J/\psi}$ Lansberg:2019adr ;
Chen:2021tmf .

On the other hand, taking advantages of the fundamental properties of heavy
quarkonium and its production mechanism, it has been proposed that heavy
quarkonium can be served as a sensitive probe of the fundamental properties of
nuclear matter. Remarkable examples including the identification of quark
gluon plasma Matsui:1986dk ; Tang:2020ame ; Wang:2022fwq  and probing the
gluon jet quenching effect in heavy ion collisions Zhang:2022rby ;
Zhang:2024owr , searching for signatures of gluon saturation in
ultraperipheral collisions CMS:2023snh , as well as the study of short range
correlation in nucleus Xu:2019wso . Moreover, with the advent of electron-ion
colliders Accardi:2012qut ; Anderle:2021wcy , new proposals of using
${J/\psi}$ production in lepton-proton collisions to probe the nucleon
structure have been put forward recently, see e.g. DAlesio:2021yws ;
Boer:2021ehu ; Qiu:2020xum ; Chen:2023hvu ; Flore:2020jau ; Bacchetta:2018ivt
; Copeland:2023qed ; liu:2023unl . All these proposals are driven by the
prominent feature that ${J/\psi}$ in deep inelastic scattering (DIS) is
dominated by virtual photon-gluon fusion process, which can be inferred from
the leading order diagram as shown in Fig. 1(a). Such a unique feature drives
the expectation that ${J/\psi}$ production is more sensitive to the gluon
content of nucleon and nucleus, comparing to the generally recognized
observables like dijet production Boer:2016fqd ; Xing:2012ii , which suffers
contamination from virtual photon-quark channel.

Similar feature also exists in ${J/\psi}$ photoproduction that has been
measured by HERA ZEUS:2002src ; H1:2010udv ; H1:2002voc  and complete NLO
calculations are performed in Refs. Butenschoen:2009zy ; Artoisenet:2009xh ;
Chang:2009uj . In this paper, we emphasize the exceptional role of
${J/\psi}$-tagged DIS in future Electron-Ion Collider (EIC) Accardi:2012qut ;
AbdulKhalek:2021gbh . In particular, through explicit calculation at NLO for
${J/\psi}$ production in unpolarized electron-proton ($ep$) collisions, we
explore the novel opportunities to constrain the LDMEs
$\langle{\mathcal{O}}^{{J/\psi}}({{}^{1}{S}_{0}^{[8]}})\rangle$ and
$\langle{\mathcal{O}}^{{J/\psi}}({{}^{3}{P}_{0}^{[8]}})\rangle$ considering
the high luminosity of EIC, which could provide new insights into the
${J/\psi}$ production mechanism. In addition, we present numerical result for
the nuclear modification factor at NLO for electron-nucleus ($eA$) collisions,
we find out that the ${J/\psi}$-tagged DIS could be an excellent probe to
gluon collinear distribution of nucleus. Furthermore, we extend such a
calculation to polarized $ep$ collision at NLO for the first time. Numerical
results on double longitudinal spin asymmetry are presented to test the
sensitivity to nucleon helicity distributions. The results presented in this
paper strongly motivate the measurement of ${J/\psi}$-tagged DIS at the future
EIC.

##  II NLO calculation for ${J/\psi}$ production at EIC

We start with prompt ${J/\psi}$ production in $ep$ collision,

$\displaystyle e(\ell)+P(p)\to e(\ell^{\prime})+{J/\psi}+X,$ (1)

where $\ell,p$ and $\ell^{\prime}$ are, respectively, the momentum of the
incoming electron, proton and the outgoing electron. In particular, we focus
on the ${J/\psi}$-tagged DIS event, where the phase space for the prompt
${J/\psi}$ is integrated out. Within the framework of NRQCD and leading twist
collinear factorization, the differential cross section for the
${J/\psi}$-tagged DIS can be written in the following form Fleming:1997fq

$\displaystyle\frac{d\sigma}{dQ^{2}dy}=$ $\displaystyle\sum_{i}\sum_{n}\int
dxf_{i/p}(x,\mu_{F}^{2})C_{i}(n)\langle{\mathcal{O}}^{{J/\psi}}{(n)}\rangle,$
(2)

where $Q^{2}=-q^{2}$ with $q=\ell-\ell^{\prime}$ stands for the momentum of
the virtual photon, $y=p\cdot q/p\cdot\ell$, $f_{i/p}$ is the proton parton
distribution functions (PDFs), $x$ is the momentum fraction of proton carried
by parton $i$ and $\mu_{F}$ is the factorization scale. $C_{i}(n)$ is the
short-distance coefficient with initial parton $i$ and final Fock state
$c\bar{c}[n]$, $\langle{\mathcal{O}}^{{J/\psi}}{(n)}\rangle$ is the LDME
representing the hadronization of $c\bar{c}[n]$ Fock state to the final
observed ${J/\psi}$, which has definite power counting of $v$, the relative
velocity of the heavy quark and anti-quark in the quarkonium rest frame. Up to
$\mathcal{O}(v^{4})$ in $v$ expansion, without considering the $v^{2}$
suppressed relativistic corrections for the color-singlet channel, the
${{}^{3}{S}_{1}^{[1]}},{{}^{1}{S}_{0}^{[8]}},{{}^{3}{S}_{1}^{[8]}},{{}^{3}{P}_{J}^{[8]}}$
channels will contribute to the factorization formula in Eq. (2). For the
${{}^{3}{P}_{J}^{[8]}}$ channel, we use the relation
$\langle{\mathcal{O}}^{{J/\psi}}({{}^{3}{P}_{J}^{[8]}})\rangle=(2J+1)\langle{\mathcal{O}}^{{J/\psi}}({{}^{3}{P}_{0}^{[8]}})\rangle$
to convert all the P-wave LDMEs to
$\langle{\mathcal{O}}^{{J/\psi}}({{}^{3}{P}_{0}^{[8]}})\rangle$, so that
$C_{i}({{}^{3}{P}_{J}^{[8]}})$ represents the P-wave short-distance
coefficients summed over $J=0,1,2$.

Figure 1: Representative Feynman diagrams for LO (a), real correction from
gluon (b) and quark (c) channel at NLO, and virtual correction (d) at NLO for
${J/\psi}$ production in DIS.

As shown in Fig. 1(a), the ${J/\psi}$-tagged DIS at LO is solely determined by
gluon channel. Therefore, only the color octet intermediate states
$c\bar{c}({{}^{1}{S}_{0}^{[8]}})$ and $c\bar{c}({{}^{3}{P}_{J}^{[8]}})$
contribute and the corresponding short distance coefficients can be expressed
as Fleming:1997fq

$\displaystyle C_{g}^{(0)}({{}^{1}{S}_{0}^{[8]}})=$
$\displaystyle\frac{1+(1-y)^{2}}{y}C_{0},$ (3) $\displaystyle
C_{g}^{(0)}({{}^{3}{P}_{J}^{[8]}})=$
$\displaystyle\Bigg{[}\frac{1+(1-y)^{2}}{y}\frac{3Q^{2}/m^{2}+28}{Q^{2}+4m^{2}}$
$\displaystyle\,\,\,-\frac{32yQ^{2}}{(Q^{2}+4m^{2})^{2}}\Bigg{]}C_{0},$ (4)

where
$C_{0}=\frac{2\pi^{2}\alpha_{s}\alpha^{2}e_{c}^{2}}{Q^{2}(Q^{2}+4m^{2})m}\delta(xyS-Q^{2}-4m^{2})$
with $m$ the charm quark mass and $S=(\ell+p)^{2}$, $e_{c}$ is the charm quark
electrical charge, $\alpha_{s}$ and $\alpha$ are, respectively, the QCD and
QED coupling constant.

According to the LO result, one would expect that the ${J/\psi}$-tagged DIS
can be uniquely used to determine the gluon distribution in proton. However,
such a unique feature could be contaminated by the quark channel, which
appears at higher orders. In order to test such an effect and check the
convergence of higher order corrections, we perform the first complete NLO
calculation for ${J/\psi}$-tagged DIS, that is at
$\mathcal{O}(\alpha^{2}\alpha_{s}^{2})$. At NLO, in addition to the real and
virtual corrections for ${{}^{1}{S}_{0}^{[8]}}$ and ${{}^{3}{P}_{J}^{[8]}}$
channels from quark and gluon channels as shown in the representative diagrams
in Fig. 1 (b-d), one also needs to include contributions from both
${}^{3}S_{1}^{[1]}$ and ${}^{3}S_{1}^{[8]}$ channels, where the former only
exists in gluon channel. In our calculation, the Feynman amplitudes are
generated using FeynArts Hahn:2000kx , the color and Dirac algebra are
performed with the help of FeynCalc Mertig:1990an , which results in 4
diagrams for quark channel and 8 diagrams for gluon channel in real
correction, and 17 diagrams for virtual correction including 6 counter term
diagrams.

We use dimensional regularization ($D=4-2\epsilon$) to regularize infrared and
collinear divergences in real correction, and infrared, collinear and
ultraviolet (UV) ones in virtual correction. In particular, HVBM scheme
tHooft:1972tcz  is implemented to deal with $\gamma_{5}$. In real correction,
we integrate over the whole phase space of the outgoing particles except for
the outgoing electron, which results in soft ($\sim 1/\epsilon$), collinear
($\sim 1/\epsilon$) and soft-collinear ($\sim 1/\epsilon^{2}$) divergences.
Regarding virtual corrections, we use the packages Apart Feng:2012iq  and FIRE
Smirnov:2014hma  for the reduction of loop integrals. The derived UV
divergences are removed through renormalization procedure. Notice that we
adopt a mixed renormalization schemeKlasen:2004az , in which the
renormalization constants of the heavy quark mass ($Z_{m}$), heavy quark field
($Z_{2}$) and gluon field ($Z_{3}$) are defined in the on-shell scheme, while
the renormalization constant of the strong coupling ($Z_{g}$) is defined in
$\rm\overline{MS}$ scheme, thus the self-energy diagrams for the external
fields are excluded in the computation. The soft, collinear and soft-collinear
divergences are cancelled in the sum of real and virtual corrections or
absorbed into the redefinition of the nucleon/nucleus PDF
$f_{g/p}(x,\mu_{F}^{2})$. Since we focus on the DIS event with $Q^{2}>0$ in
our study, there will be no divergence resulting from the QED splitting
$e\rightarrow e\gamma$ Qiu:2020xum .

Likewise, one can also calculate the cross section for ${J/\psi}$-tagged DIS
at NLO in longitudinal polarized collisions. As we expected, one encounters
the same singularities as in unpolarized case, except for the redefinition of
longitudinal polarized PDFs $\Delta f_{i/p}(x,\mu_{F}^{2})$. After the
regularization of all singularities, we are left with the finite part, which
results into the following differential cross section at NLO,

$\displaystyle\frac{d\Delta\sigma}{dQ^{2}dy}=$
$\displaystyle\sum_{i}\sum_{n}\int dx\Delta f_{i/p}(x,\mu_{F}^{2})\Delta
C_{i}(n)\langle{\mathcal{O}}^{{J/\psi}}{(n)}\rangle.$ (5)

In this calculation, we obtain for the first time the fully analytical
expressions for the short distance coefficients $C_{i}(n)$ for unpolarized and
$\Delta C_{i}(n)$ for longitudinal polarized cases at NLO 111Refer to an
ancillary file attached to the arXiv submission for the full analytical
expressions of $C_{i}(n)$ and $\Delta C_{i}(n)$.. Such analytical results will
be very efficient in the global analysis of LDMEs and PDFs for nucleon/nucleus
when experimental data are available.

##  III ${J/\psi}$ as a probe to nonperturbative functions

In this section, we show phenomenological results by employing the obtained
NLO analytical expressions under future EIC designs. We choose a center-of-
mass energy of $\sqrt{s}=100$ GeV to illustrate the power of ${J/\psi}$ as a
probe to nonperturbative functions. In our numerical calculation, the heavy
quark mass $m=1.55$ GeV and the number of light quark flavor $n_{f}=3$, the
factorization scale $\mu_{F}$ and the renormalization scale $\mu_{R}$ are
chosen as $\mu_{F}=\mu_{R}=\sqrt{Q^{2}+4m^{2}}$.

###  III.1 Sensitivity to NRQCD LDMEs

We first present the dependence of ${J/\psi}$-tagged DIS cross section on the
NRQCD LDMEs, which are determined through global analysis of world data at NLO
accuracy within the framework of NRQCD. With different fitting strategies,
various sets of ${J/\psi}$ LDMEs are obtained, challenges the universality of
these nonperturbative functions. In this work, we adopt the LDMEs from five
groups: Bodwin et al. Bodwin:2014gia , Butenschoen et al. Butenschoen:2011yh ,
Chao et al. Chao:2012iv , Gong et al. Gong:2012ug  and Brambilla et al
Brambilla:2022rjd ; Brambilla:2022ayc . The optimized values are summarized in
Table 1.

Table 1: ${J/\psi}$ NRQCD LDMEs from five different groups. | $\langle{\mathcal{O}}(^{3}S_{1}^{[1]})\rangle$ | $\langle{\mathcal{O}}(^{1}S_{0}^{[8]})\rangle$ | $\langle{\mathcal{O}}(^{3}S_{1}^{[8]})\rangle$ | $\langle{\mathcal{O}}(^{3}P_{0}^{[8]})\rangle$
---|---|---|---|---
|  GeV3 |  $10^{-2}$ GeV3 |  $10^{-2}$ GeV3 |  $10^{-2}$ GeV5
Bodwin | 0 | 9.9 | 1.1 | 1.1
Butenschoen | 1.32 | 3.04 | 0.16 | $-0.91$
Chao | 1.16 | 8.9 | 0.30 | 1.26
Gong | 1.16 | 9.7 | $-0.46$ | $-2.14$
Brambilla | 1.18 | $-0.63$ | 1.4 | 5.83

The LDMEs from Bodwin et al., Chao et al. and Gong et al. give
${}^{1}S_{0}^{[8]}$ dominance scenario and thus naturally conform the almost
unpolarized $J/\psi$ hadron production data at large $p_{T}$. However, the
${}^{1}S_{0}^{[8]}$ dominance scenario contradicts with the $\eta_{c}$ hadron
production data under the assumptions of heavy quark spin symmetry
(HQSS)Butenschoen:2014dra ; Han:2014jya ; Zhang:2014ybe ; Bodwin:1994jh . The
LDMEs from Butenschoen et.al. lead to the $J/\psi$ polarization puzzle problem
Brambilla:2010cs . Although, the LDMEs from Brambilla et.al do not have the
above issues, $\langle{\mathcal{O}}^{{J/\psi}}({{}^{1}{S}_{0}^{[8]}})\rangle$
from Brambilla et.al is poorly constrained and suffers from large
uncertainties.

Figure 2: NLO differential cross section as a function of $Q^{2}$ for
${J/\psi}$-tagged DIS process.

In Fig.2, we plot the ${J/\psi}$-tagged DIS differential cross section as a
function of $Q^{2}$ using CT18NLO Hou:2019efy  for proton PDFs and the five
sets of LDMEs shown in Table 1. We find that the LDMEs from five groups lead
to tremendous differences. We have checked that the differential cross
sections are dominated by ${{}^{1}{S}_{0}^{[8]}}$ and ${{}^{3}{P}_{J}^{[8]}}$
channels. The values of the LDMEs
$\langle{\mathcal{O}}^{J/\psi}(^{1}S_{0}^{[8]})\rangle,\langle{\mathcal{O}}^{J/\psi}(^{3}P_{0}^{[8]})\rangle$
from the groups of Bodwin and Chao are very close, thus leads to almost
identical predictions, which are an order of magnitude larger compared with
the predictions obtained using the LDMEs from Butenschoen et al. and
significantly larger than those derived using the LDMEs from Gong et al.. The
LDMEs from Brambilla et al. give similar predictions with Bodwin et al. and
Chao et al. at low $Q^{2}$ but gradually converge towards those using the
LDMEs from Gong et al. at large $Q^{2}$, which leads to a distinct shape of
$Q^{2}$ dependence. The drastic difference of ${J/\psi}$-tagged DIS
differential cross section evaluated with the LDMEs from five different
groups, clearly demonstrates that this observable is a very sensitive probe of
NRQCD LDMEs, especially for
$\langle{\mathcal{O}}^{J/\psi}(^{1}S_{0}^{[8]})\rangle$ and
$\langle{\mathcal{O}}^{J/\psi}(^{3}P_{0}^{[8]})\rangle$.

###  III.2 Sensitivity to nuclear PDFs

The high luminosity of electron-nucleus program at the EIC, combined with the
gluon-dominated ${J/\psi}$-tagged DIS process, promises profound insights into
nuclear gluon distribution. To delve deeper into the role of the
${J/\psi}$-tagged DIS process at the EIC in elucidating nuclear partonic
structure, we investigate the nuclear modification factor defined as follows,

$\displaystyle
R_{eA}=\left.\frac{d\sigma^{eA}}{dx_{B}}\right/\frac{d\sigma^{ep}}{dx_{B}},$
(6)

in which we have adopted EPPS21 parametrization Eskola:2021nhw  for nuclear
PDFs and CT18NLO for proton PDFs. In the top panel of Fig. 3, we present the
nuclear modification factor plotted against the Bjorken variable
$x_{B}=Q^{2}/p\cdot q$. We observe that the five sets of LDMEs yield nearly
identical results for $R_{eA}$, indicating the insensitivity of LDME choices
and substantially reducing theoretical uncertainty arising from LDME
variations. In the middle panel of Fig. 3, we use the LDMEs from Bodwin et.al.
as an example to show the normalized (by central values) uncertainty of the
predicted $R_{eA}$. The blue band represents normalized uncertainty resulting
from the nuclear PDFs uncertainty using EPPS21 parametrization, from which we
can see that the relative uncertainty increases as $x_{B}$ getting larger and
becomes significant when $x_{B}$ approaching $0.6$. The normalized uncertainty
coming from the statistical error, which is based on the best fit of EPPS21
and the expected EIC integrated luminosity 10 $fb^{-1}/A$ for $eA$ collisions
AbdulKhalek:2021gbh  and the detector response simulation Li:2022kwn , is
shown in the vertical red bars. It is obvious that the uncertainty coming from
the estimated statistical error is negligible comparing to those resulting
from the EPPS21 uncertainty. It is important to note that the $eA$
differential cross section is dominated by the gluon channel in all region of
$x_{B}$ considered as it was shown in the bottom panel of Fig. 3, where
$R_{g}$ is defined as the ratio of ${J/\psi}$ tagged differential cross
sections in $eA$ collision without and with contributions from quark channels.
The fact that the gluon channel gives $100\%$ to $110\%$ of the differential
cross section calculated at $\mathcal{O}(\alpha^{2}\alpha_{s}^{2})$, which
means the quark channels only give small negative contribution at
$\mathcal{O}(\alpha^{2}\alpha_{s}^{2})$, strongly indicates that the
${J/\psi}$-tagged DIS process can be served as a clean probe to access the
nuclear gluon distribution at the future EIC.

Figure 3: The nuclear modification factor as a function of $x_{B}$ for
leptoproduction of ${J/\psi}$ in $ePb$ collision at $\sqrt{s}=100$ GeV.

###  III.3 Sensitivity to nucleon helicity distribution

We extend our NLO calculation of leptoproduction of ${J/\psi}$ from
unpolarized to polarized $ep$ collisions, aiming at probing the proton gluon
helicity distribution, which is of significant importance in understanding the
configuration of proton spin Ji:2020ena .

Figure 4: The double longitudinal spin asymmetry as a function of $x_{B}$ for
leptoproduction of ${J/\psi}$ in $ep$ collisions at $\sqrt{s}=100$ GeV.

We focus on the double longitudinal spin asymmetry for leptoproduction of
${J/\psi}$ defined as

$\displaystyle
A_{LL}=\frac{d\sigma(+,+)-d\sigma(+,-)}{d\sigma(+,+)+d\sigma(+,-)}=\frac{d\Delta\sigma}{d\sigma},$
(7)

where the first and the second signs ($+,-$) in the parentheses right after
$d\sigma$ represent the helicity of the incoming lepton and proton,
respectively. We apply NNPDF parametrization Nocera:2014gqa  for the
longitudinal polarized proton PDFs.

As we can see from the top panel of Fig. 4, the predicted central value of
$A_{LL}$ with the LDMEs sets from five groups are considerably different,
except for those from Bodwin et al. and Chao et al. due to their close values
of
$\langle{\mathcal{O}}^{J/\psi}(^{1}S_{0}^{[8]})\rangle,\langle{\mathcal{O}}^{J/\psi}(^{3}P_{0}^{[8]})\rangle$.
Once the values of $J/\psi$ LDMEs are fixed, one can use the theory
predictions of $A_{LL}$ to constrain the gluon helicity distribution since
$A_{LL}$ is dominated by gluon channel, as it was shown in the bottom panel of
Fig. 4, where $R_{g}$ is defined as the ratio of $d\Delta\sigma$ evaluated
without and with the contributions from quark channels, respectively. The
middle panel of Fig. 4 shows the relative large uncertainty coming from
lacking of precise knowledge on the gluon helicity distribution. This can be
reduced significantly when we compare the $A_{LL}$ data in the future EIC with
integrated luminosity 100 $fb^{-1}$ for $ep$ collisions AbdulKhalek:2021gbh
with theory predictions.

##  IV Summary

We calculate the leptoproduction of ${J/\psi}$ with five sets of LDMEs at NLO,
under both unpolarized and polarized situations. The tremendous difference
arises from various LDMEs indicates that the ${J/\psi}$-tagged DIS is a very
sensitive probe to solve the problem of the universality of NRQCD LDMEs by
comparing with the future EIC data. Besides, considering both the high
collision luminosity and highly polarized beams in future EIC, as well as the
dominance of gluon channel in ${J/\psi}$-tagged DIS process, this process can
be also served as an unique opportunity to probe the nuclear gluon
distribution and polarized gluon distribution. All these aspects strongly
motivate the measurement of ${J/\psi}$-tagged DIS at the future EIC.

Acknowledgement. This work was supported in part by the Guangdong Major
Project of Basic and Applied Basic Research No. 2020B0301030008 and
2022A1515010683, by the National Natural Science Foundation of China (NSFC)
(Nos. 12025501, 12035007 and 12022512). The work of X.-P. W. is supported by
the DFG (Deutsche Forschungsgemeinschaft, German Research Foundation) Grant
No. BR4058/2-2. X.-P. W. acknowledges support from the DFG cluster of
excellence “ORIGINS” under Germany’s Excellence Strategy - EXC-2094 -
390783311. X.-P. W. acknowledges support from STRONG-2020- European Union’s
Horizon 2020 research and innovation program under grant agreement No. 824093.

## References

  * (1) A. Rothkopf,  Phys. Rept. 858, 1 (2020), arXiv:1912.02253.
  * (2) E. Chapon et al.,  Prog. Part. Nucl. Phys. 122, 103906 (2022), arXiv:2012.14161.
  * (3) G. T. Bodwin, E. Braaten, and G. P. Lepage,  Phys. Rev. D 51, 1125 (1995), arXiv:hep-ph/9407339,  [Erratum: Phys.Rev.D 55, 5853 (1997)].
  * (4) G. T. Bodwin, H. S. Chung, U.-R. Kim, and J. Lee,  Phys. Rev. Lett. 113, 022001 (2014), arXiv:1403.3612.
  * (5) M. Butenschoen and B. A. Kniehl,  Phys. Rev. D84, 051501 (2011), arXiv:1105.0820.
  * (6) K.-T. Chao, Y.-Q. Ma, H.-S. Shao, K. Wang, and Y.-J. Zhang,  Phys. Rev. Lett. 108, 242004 (2012), arXiv:1201.2675.
  * (7) B. Gong, L.-P. Wan, J.-X. Wang, and H.-F. Zhang,  Phys. Rev. Lett. 110, 042002 (2013), arXiv:1205.6682.
  * (8) G. T. Bodwin et al.,  Phys. Rev. D 93, 034041 (2016), arXiv:1509.07904.
  * (9) N. Brambilla, H. S. Chung, A. Vairo, and X.-P. Wang,  Phys. Rev. D 105, L111503 (2022), arXiv:2203.07778.
  * (10) N. Brambilla, H. S. Chung, A. Vairo, and X.-P. Wang,  JHEP 03, 242 (2023), arXiv:2210.17345.
  * (11) J.-P. Lansberg,  Phys. Rept. 889, 1 (2020), arXiv:1903.09185.
  * (12) A.-P. Chen, Y.-Q. Ma, and H. Zhang,  Adv. High Energy Phys. 2022, 7475923 (2022), arXiv:2109.04028.
  * (13) T. Matsui and H. Satz,  Phys. Lett. B 178, 416 (1986).
  * (14) Z. Tang, W. Zha, and Y. Zhang,  Nucl. Sci. Tech. 31, 81 (2020), arXiv:2105.11656.
  * (15) H. Wang and J. Chen,  Nucl. Sci. Tech. 33, 15 (2022).
  * (16) S.-L. Zhang, J. Liao, G.-Y. Qin, E. Wang, and H. Xing,  Sci. Bull. 68, 2003 (2023), arXiv:2208.08323.
  * (17) S.-L. Zhang and H. Xing,  (2024), arXiv:2403.12704.
  * (18) CMS, A. Tumasyan et al.,  Phys. Rev. Lett. 131, 262301 (2023), arXiv:2303.16984.
  * (19) J. Xu and F. Yuan,  Phys. Lett. B 801, 135187 (2020), arXiv:1908.10413.
  * (20) A. Accardi et al.,  Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701.
  * (21) D. P. Anderle et al.,  Front. Phys. (Beijing) 16, 64701 (2021), arXiv:2102.09222.
  * (22) U. D’Alesio, L. Maxia, F. Murgia, C. Pisano, and S. Rajesh,  JHEP 03, 037 (2022), arXiv:2110.07529.
  * (23) D. Boer, C. Pisano, and P. Taels,  Phys. Rev. D 103, 074012 (2021), arXiv:2102.00003.
  * (24) J.-W. Qiu, X.-P. Wang, and H. Xing,  Chin. Phys. Lett. 38, 041201 (2021), arXiv:2005.10832.
  * (25) L. Chen, H. Xing, and S. Yoshida,  Phys. Rev. D 108, 054021 (2023), arXiv:2306.12647.
  * (26) C. Flore, J.-P. Lansberg, H.-S. Shao, and Y. Yedelkina,  Phys. Lett. B 811, 135926 (2020), arXiv:2009.08264.
  * (27) A. Bacchetta, D. Boer, C. Pisano, and P. Taels,  Eur. Phys. J. C 80, 72 (2020), arXiv:1809.02056.
  * (28) M. Copeland, S. Fleming, R. Gupta, R. Hodges, and T. Mehen,  Phys. Rev. D 109, 074032 (2024), arXiv:2310.13737.
  * (29) H. liu, X. xie, and Z. Lu,  Phys. Lett. B 849, 138439 (2024), arXiv:2310.00609.
  * (30) D. Boer, P. J. Mulders, C. Pisano, and J. Zhou,  JHEP 08, 001 (2016), arXiv:1605.07934.
  * (31) H. Xing, Z.-B. Kang, I. Vitev, and E. Wang,  Phys. Rev. D 86, 094010 (2012), arXiv:1206.1826.
  * (32) ZEUS, S. Chekanov et al.,  Eur. Phys. J. C 27, 173 (2003), arXiv:hep-ex/0211011.
  * (33) H1, F. D. Aaron et al.,  Eur. Phys. J. C 68, 401 (2010), arXiv:1002.0234.
  * (34) H1, C. Adloff et al.,  Eur. Phys. J. C 25, 25 (2002), arXiv:hep-ex/0205064.
  * (35) M. Butenschoen and B. A. Kniehl,  Phys. Rev. Lett. 104, 072001 (2010), arXiv:0909.2798.
  * (36) P. Artoisenet, J. M. Campbell, F. Maltoni, and F. Tramontano,  Phys. Rev. Lett. 102, 142001 (2009), arXiv:0901.4352.
  * (37) C.-H. Chang, R. Li, and J.-X. Wang,  Phys. Rev. D 80, 034020 (2009), arXiv:0901.4749.
  * (38) R. Abdul Khalek et al.,  Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419.
  * (39) S. Fleming and T. Mehen,  Phys. Rev. D57, 1846 (1998), arXiv:hep-ph/9707365.
  * (40) T. Hahn,  Comput. Phys. Commun. 140, 418 (2001), arXiv:hep-ph/0012260.
  * (41) R. Mertig, M. Bohm, and A. Denner,  Comput. Phys. Commun. 64, 345 (1991).
  * (42) G. ’t Hooft and M. J. G. Veltman,  Nucl. Phys. B 44, 189 (1972).
  * (43) F. Feng,  Comput. Phys. Commun. 183, 2158 (2012), arXiv:1204.2314.
  * (44) A. V. Smirnov,  Comput. Phys. Commun. 189, 182 (2015), arXiv:1408.2372.
  * (45) M. Klasen, B. A. Kniehl, L. N. Mihaila, and M. Steinhauser,  Phys. Rev. D 71, 014016 (2005), arXiv:hep-ph/0408280.
  * (46) M. Butenschoen, Z.-G. He, and B. A. Kniehl,  Phys. Rev. Lett. 114, 092004 (2015), arXiv:1411.5287.
  * (47) H. Han, Y.-Q. Ma, C. Meng, H.-S. Shao, and K.-T. Chao,  Phys. Rev. Lett. 114, 092005 (2015), arXiv:1411.7350.
  * (48) H.-F. Zhang, Z. Sun, W.-L. Sang, and R. Li,  Phys. Rev. Lett. 114, 092006 (2015), arXiv:1412.0508.
  * (49) N. Brambilla et al.,  Eur. Phys. J. C 71, 1534 (2011), arXiv:1010.5827.
  * (50) T.-J. Hou et al.,  Phys. Rev. D 103, 014013 (2021), arXiv:1912.10053.
  * (51) K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado,  Eur. Phys. J. C 82, 413 (2022), arXiv:2112.12462.
  * (52) X. Li et al.,  Nucl. Instrum. Meth. A 1048, 167956 (2023), arXiv:2207.10356.
  * (53) X. Ji, F. Yuan, and Y. Zhao,  Nature Rev. Phys. 3, 27 (2021), arXiv:2009.01291.
  * (54) NNPDF, E. R. Nocera, R. D. Ball, S. Forte, G. Ridolfi, and J. Rojo,  Nucl. Phys. B887, 276 (2014), arXiv:1406.5539.


# Mini-batch stochastic three-operator splitting for distributed optimization

Barbara Franci and Mathias Staudigl1 *We thank Olivier Bilenne for his
contribution to the early stage of this research. M. Staudigl acknowledges
financial support from the FMJH Program PGMO and from the support of EDF. 1The
authors are with the Department of Data Science and Knowledge Engineering,
Maastricht University, P.O. Box 616, NL–6200 MD Maastricht, The Netherlands
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider a network of agents, each with its own private cost consisting of
a sum of two possibly nonsmooth convex functions, one of which is composed
with a linear operator. At every iteration each agent performs local
calculations and can only communicate with its neighbors. The challenging
aspect of our study is that the smooth part of the private cost function is
given as an expected value and agents only have access to this part of the
problem formulation via a heavy-tailed stochastic oracle. To tackle such
sampling-based optimization problems, we propose a stochastic extension of the
triangular pre-conditioned primal-dual algorithm. We demonstrate almost sure
convergence of the scheme and validate the performance of the method via
numerical experiments.

##  I Introduction

We consider a large class of convex optimization problems given by

$\min_{x\in\mathbb{R}^{n}}f(x)+g(x)+h(Lx),$ (1)

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is a convex and continuously
differentiable function, $g:\mathbb{R}^{n}\to(-\infty,\infty]$ and
$h:\mathbb{R}^{m}\to(-\infty,\infty]$ are closed convex and lower semi-
continuous functions, and $L:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is a given
linear map. Such a structure is very general, and and describes many
applications that range from signal processing to machine learning to control
[1, 2, 3]. In many instances of problem (1), the cost functions are
contaminated by stochastic noise. In such settings, we are given a probability
space $(\Omega,\mathcal{A},\mathbb{P})$ carrying a random variable
$\xi:\Omega\to\Xi\subset\mathbb{R}^{d}$, and a measurable function
$F:\mathbb{R}^{n}\times\Xi\to\mathbb{R}$ so that

$f(x)=\mathbb{E}[F(x,\xi)]\qquad\forall x\in\mathbb{R}^{n}.$ (2)

The presence of stochastic uncertainty challenges any direct solution method
for problem (1), since the smooth function $f$ is not directly accessible in
practice, unless the distribution of the random variable is known. Indeed, the
expectation (2) cannot be computed exactly, and instead we need to restore to
simulation-based techniques. A reasonable assumption is that we can draw
samples from the distribution of the random variable $\xi$, and stochastically
approximate the necessary information about the function $f$. Specifically, we
adopt an online gradient-based stochastic approximation method where a
deterministic version of a numerical algorithm, known to solve problem (1) in
its expected-value formulation, is supplied with stochastic estimators of the
gradient of $f$. In this setting, the mechanism to access $f$ via samples of
the law of $\xi$ is usually named a _stochastic oracle_ (SO). The SO outputs
unbiased approximations of the gradient obtained via an average over a batch
of realizations. The larger the number of samples, the smaller the variance of
point estimators, and thus the more precise information we obtain at the cost
of generating a larger amount of i.i.d random variables. This trade-off
between accuracy of estimators and the available simulation budget is what
makes such _mini-batch_ stochastic approximation approaches efficient methods
of choice [4, 5, 6].

###  I-A Main Contributions and relation to the literature

A standard assumption in stochastic optimization is that the approximation
error is uniformly bounded [5]. Instead, the stochastic approximation approach
developed in this paper allows to handle stochastic oracles with potentially
unbounded moments. This is of relevance in primal-dual dynamics, which usually
act on unbounded domains, for which any a-priori variance bound is rather
unnatural.

Besides weaker hypothesis on the noise structure, this paper is the first
stochastic primal-dual algorithm which is even provably convergent in the
multi-agent formulation of problem (1), i.e., when a finite set of agents
$i\in\\{1,\ldots,m\\}$ cooperatively minimize the composite objective function

$\begin{split}\min_{x_{i}\in\mathbb{R}^{p_{i}},1\leq i\leq
m}\sum_{i=1}^{m}f_{i}(x_{i})+g_{i}(x_{i})+h_{i}(L_{i}x_{i}),\\\ \text{s.t.:
}A_{ij}x_{i}+A_{ji}x_{j}=b_{(i,j)}\quad\forall(i,j)\in E\end{split}$ (3)

Problems of this form appear in several application fields. In distributed
model predictive control, $f_{i}$ can represent individual finite-horizon
costs for each agent, $L_{i}$ model linear dynamics of each agent, and $h_{i}$
models state and input constraints [1]. In machine learning, $f_{i}$ would
represent a smooth data fidelity terms, $Lx$ linear restrictions on the
parameters and $h_{i}$ and $g_{i}$ take over the role of statistical penalties
reflecting a-priori structure properties of the parameters to be estimated
[2].

Compared with existing work, this paper makes the following contributions:

  1. (i)

Our algorithm is a stochastic extension of the triangular preconditioned
primal-dual algorithm (TriPD), developed in [7, 8].

  2. (ii)

To the best of our knowledge, our scheme is the first stochastic approximation
method which is able to solve the three-operator splitting problem
characterizing primal-dual pairs for (1).

  3. (iii)

The analysis immediately generalizes to the multi-agent formulation (3), where
distributed iterations, i.e., locally performed without central supervision,
are obtained.

The only related contributions we are aware of are [9, 10]. Both assume a
uniformly bounded SO, which is very restrictive in primal-dual methods,
essentially forcing a-priori compactness on the domain of $g$. For the special
case when $g=0$, stochastic accelerated algorithms for the centralized problem
(1) have been considered in [11], imposing a uniformly bounded noise condition
on the SO. We include the non-smooth term $g$ and allow for heavy-tailed noise
in the SO.

Our analysis is restricted to synchronous versions of distributed optimization
algorithms. It is possible to extend our analysis to the asynchronous case and
block-coordinate descent strategies, acting on general real separable Hilbert
spaces. We will present these extensions in a future paper.

###  I-B Notation and preliminary results

Throughout $\mathcal{H},\mathcal{G}$ are finite dimensional Euclidean spaces.
Their scalar products and associated norms are denoted by
$\langle\cdot,\cdot\rangle$ and $\|\cdot\|$. We denote
$\mathcal{B}(\mathcal{H},\mathcal{G})$ the space of bounded linear operators
from $\mathcal{H}$ to $\mathcal{G}$. The adjoint of
$L\in\mathcal{B}(\mathcal{H},\mathcal{G})$ is denoted by $L^{\ast}$.
$\operatorname{Id}$ denotes the identity operator. We set
$\mathcal{P}_{\alpha}(\mathcal{H})\triangleq\\{\mathbf{Q}\in\mathcal{B}(\mathcal{H})|\mathbf{Q}^{\ast}=\mathbf{Q}\text{
and }\mathbf{Q}\succ\alpha\operatorname{Id}\\}.$ For $\alpha>0$ and
$\mathbf{Q}\in\mathcal{P}_{\alpha}(\mathcal{H})$, we define a scalar product
and norm on $\mathcal{H}$ by $\langle
x,y\rangle_{\mathbf{Q}}\triangleq\langle\mathbf{Q}x,y\rangle$, and
$\|x\|_{\mathbf{Q}}\triangleq\sqrt{\langle\mathbf{Q}x,x\rangle}$ for all
$x,y\in\mathcal{H}$. We let $\mathtt{M}(\mathbf{Q})=\max_{u:\|u\|=1}\langle
u,\mathbf{Q}u\rangle$ and $\mathtt{m}(\mathbf{Q})=\min_{u:\|u\|=1}\langle
u,\mathbf{Q}u\rangle$. For an extended-valued real function $f$, we use
$\operatorname{dom}f=\\{x\in\mathcal{H}|f(x)<\infty\\}$ for its effective
domain. Let $\mathbf{Q}\in\mathcal{P}_{\alpha}(\mathcal{H})$ for some
$\alpha>0$. $\delta_{C}$ is the indicator function of the set $C$, that is,
$\delta_{C}(x)=0$ if $x\in C$ and $\delta_{C}(x)=\infty$ otherwise. The
weighted prox-operator of $f$ is defined as
$\operatorname{Prox}^{\mathbf{Q}}_{f}(x)\triangleq\operatorname*{argmin}_{u\in\mathcal{H}}\left\\{f(u)+\frac{1}{2}\|u-x\|^{2}_{\mathbf{Q}}\right\\}.$
The conjugate of $f:\mathcal{H}\to(-\infty,\infty]$ is
$f^{\ast}(u)=\sup_{x\in\mathcal{H}}\\{\langle x,u\rangle-f(x)\\}.$ We also
need the celebrated Robbins-Siegmund Lemma for the convergence analysis

######  Lemma 1 ([12])

Let $(\Omega,\mathcal{A},(\mathcal{A}_{n})_{n\in\mathbb{N}},\mathbb{P})$ be a
filtered probability space satisfying the usual conditions. For every
$n\in\mathbb{N}$, let $v_{n},\xi_{n},\zeta_{n}$ and $t_{n}$ be non-negative
$\mathcal{A}_{n}$-measurable random variables such that
$\\{\zeta_{n}\\}_{n\in\mathbb{N}}$ and $\\{t_{n}\\}_{n\in\mathbb{N}}$ are
summable and for all $n\in\mathbb{N}$,

$\mathbb{E}[v_{n}|\mathcal{F}_{n}]\leq(1+t_{n})v_{n}+\zeta_{n}-\xi_{n}\qquad\mathbb{P}-\text{a.s.}$
(4)

Then $\\{v_{n}\\}_{n\in\mathbb{N}}$ converges and
$\\{\xi_{n}\\}_{n\in\mathbb{N}}$ is summable $\mathbb{P}$-a.s.

##  II Stochastic Primal-Dual Algorithm

In this section we propose a stochastic primal-dual algorithm for solving (1).
Let $\mathcal{H}\equiv\mathbb{R}^{n}$ and $\mathcal{G}\equiv\mathbb{R}^{m}$
with inner products
$\langle\cdot,\cdot\rangle_{\mathcal{H}},\langle\cdot,\cdot\rangle_{\mathcal{G}}$,
and corresponding norms $\|\cdot\|_{\mathcal{H}},\|\cdot\|_{\mathcal{G}}$.

###### Assumption 1

Throughout the paper the following assumptions shall be in place:

  * (i)

$g:\mathcal{H}\to(-\infty,\infty],h:\mathcal{G}\to(-\infty,\infty]$ are
proper, closed, convex and lower-semi continuous functions.

  * (ii)

$L:\mathcal{H}\to\mathcal{G}$ is a linear mapping with adjoint $L^{\ast}$.

  * (iii)

$f:\mathcal{H}\to\mathbb{R}$ is convex, continuously differentiable and there
exists $\alpha>0,\mathbf{Q}\in\mathcal{P}_{\alpha}(\mathcal{H})$, and
$\beta_{f}>0$ such that $\|\nabla f(x)-\nabla
f(y)\|_{\mathbf{Q}^{-1}}\leq\beta_{f}\|x-y\|_{\mathbf{Q}}$ for all
$x,y\in\mathcal{H}.$

  * (iv)

Let $\Xi\subset\mathbb{R}^{d}$ be a measurable set and
$(\Omega,\mathcal{A},\mathbb{P})$ a probability space. There exists a
measurable function $F:\mathcal{H}\times\Xi\to\mathbb{R}$ such that (2) holds.

  * (v)

The set of solutions to (1), denoted by $\mathcal{X}^{\ast}$, is nonempty
Moreover, there exists $x\in\operatorname{relint}(\operatorname{dom}g)$ such
that $Lx\in\operatorname{relint}(\operatorname{dom}h)$.

Let $\mathcal{Z}=\mathcal{G}\times\mathcal{H}$ represent the product space
with the inner product $\langle(y_{1},x_{1}),(y_{2},x_{2})\rangle=\langle
y_{1},y_{2}\rangle_{\mathcal{G}}+\langle x_{1},x_{2}\rangle_{\mathcal{H}},$
and associated norm
$\|(y,x)\|=\sqrt{\|y\|^{2}_{\mathcal{G}}+\|x\|^{2}_{\mathcal{H}}}.$ Whenever
clear from the context, we will suppress the ambient space from the inner
products and norms, respectively.

We remark that by the Baillon-Haddad theorem [13, Corollary 18.17] the
$\beta_{f}$-Lipschitz smoothness of the function $f$ is equivalent to the
$\tfrac{1}{\beta_{f}}$-cocoercivity of $\nabla f$, i.e.,

$\langle\nabla f(x)-\nabla
f(x^{\prime}),x-x^{\prime}\rangle\geq\tfrac{1}{\beta_{f}}\|\nabla f(x)-\nabla
f(x^{\prime})\|^{2}.$

for all $x,x^{\prime}\in\mathcal{H}$.

Problem (1) can be equivalently written as a max-min problem

$\max_{y\in\mathcal{G}}\min_{x\in\mathcal{H}}\\{f(x)+g(x)+\langle
Lx,y\rangle-h^{\ast}(y)\\}$

Then, $(\bar{y},\bar{x})\in\mathcal{G}\times\mathcal{H}$ is called a primal-
dual pair if

$\begin{split}0&\in\partial h^{\ast}(\bar{y})-L\bar{x},\\\ 0&\in\partial
g(\bar{x})+\nabla f(\bar{x})+L^{\ast}\bar{y}.\end{split}$

Introduce the maximally monotone operators

$\displaystyle\mathsf{A}(y,x)\triangleq\partial h^{\ast}(y)\times\partial
g(x),$ $\displaystyle\mathsf{M}(y,x)\triangleq[-Lx,L^{\ast}y],$
$\displaystyle\mathsf{C}(y,x)\triangleq[0,\nabla f(x)].$

Then $\bar{z}\triangleq(\bar{y},\bar{x})$ is a primal-dual pair if and only if
$\bar{z}\in\operatorname{Zer}(\mathsf{G})\triangleq\\{z\in\mathcal{Z}|0\in\mathsf{G}(z)\\},$
where
$\mathsf{G}(y,x)\triangleq\mathsf{A}(y,x)+\mathsf{M}(y,x)+\mathsf{C}(y,x).$

###  II-A Triangular pre-conditioned primal-dual algorithm

If it is possible to access the operator $\mathsf{C}$ directly, then an
application of the Triangular pre-conditioned Primal-Dual (TriPD) algorithm of
[7] would be possible. Let $\tau_{1}>0$ and $\tau_{2}>0$ be two positive real
numbers, and
$\Sigma\in\mathcal{P}_{1/\tau_{1}}(\mathcal{G}),\Gamma\in\mathcal{P}_{1/\tau_{2}}(\mathcal{H})$.
TriPD can be compactly written as the fixed point iteration
$z_{k+1}=T(z_{k})$, where $z_{k}\triangleq(y_{k},x_{k})\in\mathcal{Z}$, and
the mapping $T:\mathcal{Z}\to\mathcal{Z}$ is given by

$\displaystyle\Phi_{1}(z)\triangleq\operatorname{Prox}^{\Sigma^{-1}}_{h^{\ast}}(y+\Sigma
Lx),$ (5)
$\displaystyle\Phi_{2}(z)\triangleq\operatorname{Prox}^{\Gamma^{-1}}_{g}(x-\Gamma\nabla
f(x)-\Gamma\Phi_{1}(z))$ $\displaystyle T(z)\triangleq[\Phi_{1}(z)+\Sigma
L(\Phi_{2}(z)-x),\Phi_{2}(z)].$

First, let us report an important connection between the fixed points of the
mapping $T$ and the set of primal-dual solutions to (1).

###### Lemma 2

[7, Equation (15)]. We have
$\operatorname{Fix}(T)\triangleq\\{z\in\mathcal{Z}|T(z)=z\\}=\operatorname{Zer}(\mathsf{G})$.

Before taking care of our expected valued formulation as in (2), let us
introduce the matrices

$\displaystyle\mathbf{P}\triangleq\left(\begin{array}[]{cc}\Sigma^{-1}&\tfrac{1}{2}L\\\
\tfrac{1}{2}L^{\ast}&\Gamma^{-1}\end{array}\right),\mathbf{K}\triangleq\left(\begin{array}[]{cc}0&-\tfrac{1}{2}L\\\
\tfrac{1}{2}L^{\ast}&0\end{array}\right),$
$\displaystyle\mathbf{S}\triangleq\left(\begin{array}[]{cc}\Sigma^{-1}&0\\\
0&\Gamma^{-1}\end{array}\right),\mathbf{R}\triangleq\mathbf{P}+\mathbf{K}.$

These matrices act like step-sizes and pre-conditioners. Allowing for matrix-
valued step sizes makes the distributed algorithm in Section IV a corollary of
the present analysis.

Initialize: $(x_{0},y_{0})\in\mathcal{Z}$,
$\Sigma\in\mathcal{P}_{1/\tau_{1}}(\mathcal{G})$,
$\Gamma\in\mathcal{P}_{1/\tau_{2}}(\mathcal{H})$

Iteration $k$:  $\displaystyle\hat{Y}_{k}$
$\displaystyle=\operatorname{Prox}^{\Sigma^{-1}}_{h^{\ast}}(Y_{k}+\Sigma
LX_{k})$ $\displaystyle X_{k+1}$
$\displaystyle=\operatorname{Prox}_{g}^{\Gamma^{-1}}(x_{k}-\Gamma\mathcal{F}_{k}(X_{k})-\Gamma
L^{\ast}\hat{Y}_{k})$ $\displaystyle Y_{k+1}$
$\displaystyle=\hat{Y}_{k}+\Sigma L(X_{k+1}-X_{k})$

Algorithm 1 Stochastic Triangularly preconditioned primal-dual algorithm
(STriPD)

###  II-B Stochastic TriPD

Since we cannot evaluate $\mathsf{C}(z)$ directly, we let $\mathcal{C}_{k}(z)$
denote its stochastic estimator. This estimator is constructed by a Monte-
Carlo scheme involving mini-batches. Given an i.i.d. sample
$\xi_{1:N}\triangleq\\{\xi^{i}\\}_{i=1}^{N}$ drawn from the law of $\xi$, let

$\bar{w}(x,\xi_{1:N})\triangleq\tfrac{1}{N}\textstyle{\sum_{i=1}^{N}}(\nabla_{x}F(x,\xi^{i})-\nabla
f(x)).$

The approximation take place at each iteration $k$ and a sequence
$\\{N_{k}\\}_{k\in\mathbb{N}}$ (the _batch size_) defines the number of random
variables we need to sample in each iteration.

###### Assumption 2

We are given a sequence $\\{N_{k}\\}_{k\in\mathbb{N}}\subset\mathbb{N}$ such
that $\lim_{k\to\infty}N_{k}=\infty$.

Since complexity questions are beyond our scope, we do not specify the speed
at which batch sizes grow (see, however, [6]). For a given sequence of batch-
sizes, we set $w_{k}(x,\omega)\triangleq\bar{w}(x,\xi_{1:N_{k}}(\omega))$.

###### Assumption 3

For all $k\geq 0$, $\mathcal{C}_{k}(z,\omega)=[0,\mathcal{F}_{k}(x,\omega)]$,
where $\mathcal{F}_{k}(x,\omega)=\nabla f(x)+w_{k}(x,\omega)$ for all
$(x,\omega)\in\mathcal{H}\times\Omega.$

We let $(\mathcal{A}_{k})_{k\geq 0}$ denote the filtration given by
$\mathcal{A}_{0}\triangleq\\{\emptyset,\Omega\\}$, and for $k\geq 1$,
$\mathcal{A}_{k}\triangleq\sigma(\xi_{1:N_{0}},\ldots,\xi_{1:N_{k-1}})$.
Define $\nu(x)\triangleq\|\nabla_{x}F(x,\xi)-\nabla f(x)\|^{2}$, and
$\sigma(x)\triangleq\sqrt{\mathbb{E}[\nu(x)|x]}.$

###### Assumption 4

For all $k\geq 0$ we have
$\mathbb{E}[\mathcal{F}_{k}(X_{k})|\mathcal{A}_{k}]=\nabla f(X_{k})$ a.s..
Moreover, there exists $x^{\ast}\in\mathcal{X}^{\ast}$ and
$\sigma_{0}(x^{\ast})\geq 0,\sigma_{1}>0$ such that

$\sigma(x)\leq\sigma_{0}(x^{\ast})+\sigma_{1}\|x-x^{\ast}\|^{2}\quad\forall
x\in\mathcal{H}.$

This is a heavy-tailed noise assumption which allows us to consider random
perturbation even with unbounded second moment [4, 14].

In analogy to TriPD and in light of the approximation scheme, the equations in
(5) correspond to those in Algorithm 1 and lead to a sequence of random maps
$\\{\mathbb{T}_{k}\\}_{k\geq 0}:\mathcal{Z}\times\Omega\to\mathcal{Z}$
generating a stochastic process $\\{Z_{k}\\}_{k\geq 0}=\\{(Y_{k},X_{k});k\geq
0\\}$ via recursive updates

$Z_{k+1}=\mathbb{T}_{k}(Z_{k}),\quad Z_{0}\in\mathcal{Z}\text{ given.}$

The exact action of this mapping can be described as follows. For
$z=(y,x)\in\mathcal{Z}$, let

$\displaystyle\mathbb{Y}_{k}(z)\triangleq\operatorname{Prox}^{\Sigma^{-1}}_{h^{\ast}}(y+\Sigma
Lx),$
$\displaystyle\mathbb{X}_{k}(z)\triangleq\operatorname{Prox}^{\Gamma^{-1}}_{g}\left(x-\Gamma\mathcal{F}_{k}(x)-\Gamma\mathbb{Y}_{k}(z)\right)$
$\displaystyle\mathbb{T}_{k}(z)\triangleq\left[\mathbb{Y}_{k}(z)+\Sigma
L\left(\mathbb{X}_{k}(z)-x\right),\mathbb{X}_{k}(z)\right].$

Let $\hat{Z}_{k}\triangleq[\mathbb{Y}_{k}(z),\mathbb{X}_{k}(z)]$, so that
$Y_{k+1}=\mathbb{Y}_{k}(Z_{k})+\Sigma_{k}L\left(\mathbb{X}_{k}(z)-X_{k}\right)$
and $X_{k+1}=\mathbb{X}_{k}(Z_{k})$ for all $k\geq 0$. One can verify that

$\displaystyle\hat{Z}_{k}=(\mathbf{R}+\mathsf{A})^{-1}(\mathbf{R}-\mathsf{M}-\mathcal{C}_{k})(Z_{k}),\text{
and}$
$\displaystyle\mathbb{T}_{k}(z)=z+\mathbf{S}^{-1}(\mathbf{R}+\mathsf{M}^{\ast})(\hat{Z}_{k}-z).$

##  III Convergence analysis

In this section, we present a number of results that lead to the convergence
proof of Algorithm 1 (Theorem 1). We start with a property of the operator
$\mathsf{C}$.

###### Lemma 3

[7, Lemma II.4]. For all $z=(y,x)$, $z^{\prime}=(y^{\prime},x^{\prime})$,
$z^{\prime\prime}=(y^{\prime\prime},x^{\prime\prime})\in\mathcal{Z}$ we have

$\langle\mathsf{C}(z)-\mathsf{C}(z^{\prime}),z^{\prime\prime}-z\rangle\leq\tfrac{\beta_{f}}{4}\|x^{\prime\prime}-x^{\prime}\|^{2}_{\mathbf{Q}}.$

From Assumption 3, we can write

$\mathcal{C}_{k}(z)=\mathsf{C}(z)+\phi_{k+1}(z),$

where $\phi_{k+1}(z,\omega)\triangleq[0,w_{k}(x,\omega)]$. Let
$z^{\ast}\in\operatorname{Zer}(\mathsf{G})$, or
$-(\mathsf{M}+\mathsf{C})(z^{\ast})\in\mathsf{A}(z^{\ast})$. From the
definition of the update $\hat{Z}_{k}$, we deduce
$\mathbf{R}(Z_{k}-\hat{Z}_{k})-(\mathsf{M}+\mathcal{C}_{k})(Z_{k})\in\mathsf{A}(\hat{Z}_{k})$.
By monotonicity of the operator $\mathsf{A}$, this implies

$\displaystyle\langle\mathbf{R}(Z_{k}-\hat{Z}_{k})-\mathsf{M}(Z_{k}-z^{\ast})+\mathsf{C}(z^{\ast})-\mathsf{C}(Z_{k}),\hat{Z}_{k}-z^{\ast}\rangle$
(6) $\displaystyle-\langle\phi_{k+1}(Z_{k}),\hat{Z}_{k}-z^{\ast}\rangle\geq
0.$

Applying Lemma 3 to the points $z=z^{\ast},z^{\prime}=Z_{k}$ and
$z^{\prime\prime}=\hat{Z}_{k}$, we get
$\langle\mathsf{C}(z^{\ast})-\mathsf{C}(Z_{k}),\hat{Z}_{k}-z^{\ast}\rangle\leq\tfrac{\beta_{f}}{4}\|\mathbb{X}_{k}(Z_{k})-X_{k}\|^{2}_{\mathbf{Q}}.$
To reduce notational clutter, we write in the following
$\hat{X}_{k}\equiv\mathbb{X}_{k}(Z_{k})$. Then, (6) leads to the estimate

$\displaystyle 0$
$\displaystyle\leq\langle\mathbf{R}(Z_{k}-\hat{Z}_{k})+\mathsf{M}(z^{\ast}-Z_{k}),\hat{Z}_{k}-z^{\ast}\rangle$
$\displaystyle+\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle+\tfrac{\beta_{f}}{4}\|\hat{X}_{k}-X_{k}\|^{2}_{\mathbf{Q}}$
$\displaystyle=\langle\mathbf{P}(\hat{Z}_{k}-Z_{k}),z^{\ast}-Z_{k}\rangle$
$\displaystyle+\langle(\mathsf{M}-\mathbf{K})(Z_{k}-z^{\ast}),Z_{k}-\hat{Z}_{k}\rangle-\|\hat{Z}_{k}-Z_{k}\|_{\mathbf{P}}^{2}$
$\displaystyle+\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle+\tfrac{\beta_{f}}{4}\|\hat{X}_{k}-X_{k}\|^{2}_{\mathbf{Q}}$

where the last equality uses the skew-symmetry so that
$\langle(\mathsf{M}-\mathbf{K})z,z\rangle=0$ for all $z\in\mathcal{Z}$. Thanks
to the skew-symmetry, we also observe that

$\displaystyle\langle(\mathsf{M}-\mathbf{K})(Z_{k}-z^{\ast}),Z_{k}-\hat{Z}_{k}\rangle=\langle
z^{\ast}-Z_{k},\mathsf{M}^{\ast}(\hat{Z}_{k}-Z_{k})\rangle$
$\displaystyle+\langle(\mathbf{R}-\mathbf{P})(\hat{Z}_{k}-Z_{k}),z^{\ast}-Z_{k}\rangle.$

Therefore,

$\displaystyle 0$
$\displaystyle\leq\langle(\mathbf{R}+\mathsf{M}^{\ast})(\hat{Z}_{k}-Z_{k}),z^{\ast}-Z_{k}\rangle-\|\hat{Z}_{k}-Z_{k}\|^{2}_{\mathbf{P}}$
$\displaystyle+\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle+\tfrac{\beta_{f}}{4}\|\hat{X}_{k}-X_{k}\|^{2}_{\mathbf{Q}}.$

By definition
$\mathbf{S}(\mathbb{T}_{k}(Z_{k})-Z_{k})=(\mathbf{R}+\mathsf{M}^{\ast})(\hat{Z}_{k}-Z_{k})$,
and

$\displaystyle 0$
$\displaystyle\leq\langle\mathbf{S}(\mathbb{T}_{k}(Z_{k})-Z_{k}),z^{\ast}-Z_{k}\rangle-\|\hat{Z}_{k}-Z_{k}\|^{2}_{\mathbf{P}}$
(7)
$\displaystyle+\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle+\tfrac{\beta_{f}}{4}\|\hat{X}_{k}-X_{k}\|^{2}_{\mathbf{Q}}.$

A straightforward, but slightly tedious computation, yields the next result.

###### Lemma 4

Let
$\mathbf{U}\triangleq\left(\begin{array}[]{cc}\Sigma^{-1}&-\tfrac{1}{2}L\\\
-\tfrac{1}{2}L^{\ast}&\Gamma^{-1}-\tfrac{\beta_{f}}{4}\mathbf{Q}\end{array}\right)$.
Then,

$\|\hat{Z}_{k}-Z_{k}\|^{2}_{\mathbf{P}_{k}}-\tfrac{\beta_{f}}{4}\|\hat{X}_{k}-X_{k}\|^{2}_{\mathbf{Q}}=\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{\mathbf{U}}.$

Eq. (7) and Lemma 4 deliver the relation

$\displaystyle 0$
$\displaystyle\leq\langle\mathbf{S}(\mathbb{T}_{k}(Z_{k})-Z_{k}),z^{\ast}-Z_{k}\rangle$
(8)
$\displaystyle+\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle-\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{\mathbf{U}}.$

We now analyze the noise term
$\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle$. Denote by
$\bar{x}_{k}\triangleq\Phi_{2}(Z_{k})$ the evaluation of the deterministic
generator at $Z_{k}$ as in (5). Then,

$\displaystyle\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle$
$\displaystyle=\langle\mathcal{F}(X_{k})-\nabla
f(X_{k}),x^{\ast}-\hat{X}_{k}\rangle$
$\displaystyle=\langle\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k}),x^{\ast}-\bar{x}_{k}\rangle$
$\displaystyle-\langle\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k}),\hat{X}_{k}-\bar{x}_{k}\rangle.$

Using Cauchy-Schwarz, and the non-expansiveness of the proximal operator
$\operatorname{Prox}^{\Gamma^{-1}}_{g}(\cdot)$ with respect to the norm
$\|\cdot\|_{\Gamma}$ [13], we obtain

$\displaystyle-\langle\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k}),\hat{X}_{k}-\bar{x}_{k}\rangle$
$\displaystyle\leq\|\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k})\|_{\Gamma}\cdot\|\Gamma(\mathcal{F}(X_{k})-\nabla
f(X_{k}))\|_{\Gamma^{-1}}$ $\displaystyle=\|\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k})\|^{2}_{\Gamma}.$

Hence

$\displaystyle\langle\phi_{k+1}(Z_{k}),z^{\ast}-\hat{Z}_{k}\rangle$
$\displaystyle\leq\langle\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k}),x^{\ast}-\bar{x}_{k}\rangle$
$\displaystyle+\|\mathcal{F}_{k}(X_{k})-\nabla f(X_{k})\|^{2}_{\Gamma}.$

Let $\mathbb{E}_{k}[\cdot]\triangleq\mathbb{E}[\cdot|\mathcal{A}_{k}]$ and
recall that Assumption 4 implies that
$\mathbb{E}_{k}[\mathcal{F}_{k}(X_{k})]=\nabla f(X_{k})$ holds a.s.. Then, we
obtain

$\mathbb{E}_{k}[\langle\phi_{k+1},z^{\ast}-\hat{Z}_{k}\rangle]\leq\mathbb{E}_{k}[\|\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k})\|^{2}_{\Gamma}].$

Setting $\bar{\nu}_{k}\triangleq\|\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k})\|^{2}_{\Gamma}$, we obtain from (8)

$\begin{split}\mathbb{E}_{k}\left[\langle\mathbf{S}(\mathbb{T}_{k}(Z_{k})-Z_{k}),Z_{k}-z^{\ast}\rangle\right]\leq\mathbb{E}_{k}[\bar{\nu}_{k}]\\\
-\mathbb{E}_{k}\left[\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{\mathbf{U}}\right].\end{split}$
(9)

Since, by definition,

$\displaystyle\|Z_{k+1}-z^{\ast}\|^{2}_{\mathbf{S}}=\|(\mathbb{T}_{k}(Z_{k})-Z_{k})+(Z_{k}-z^{\ast})\|^{2}_{\mathbf{S}}$
$\displaystyle=\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{\mathbf{S}}+2\langle\mathbf{S}(\mathbb{T}_{k}(Z_{k})-Z_{k}),Z_{k}-z^{\ast}\rangle$
$\displaystyle+\|Z_{k}-z^{\ast}\|^{2}_{\mathbf{S}},$

the estimate (9) delivers

$\displaystyle\mathbb{E}_{k}$
$\displaystyle\left[\|Z_{k+1}-z^{\ast}\|^{2}_{\mathbf{S}}\right]\leq\|Z_{k}-z^{\ast}\|^{2}_{\mathbf{S}}$
$\displaystyle-\mathbb{E}_{k}\left[\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{2\mathbf{U}_{k}-\mathbf{S}}\right]+2\mathbb{E}_{k}[\bar{\nu}_{k}]$

If $2\mathbf{U}-\mathbf{S}\succ 0$ and $\\{\mathbb{E}_{k}[\nu_{k}]\\}_{k\geq
0}$ is summable, then it follows that $\\{Z_{k}\\}_{k\geq 0}$ is quasi-Fejér
monotone with respect to $\operatorname{Zer}(\mathsf{G})$ relative to
$\|\cdot\|_{\mathbf{S}}$, i.e.,

$\mathbb{E}_{k}\left[\|Z_{k+1}-z^{\ast}\|^{2}_{\mathbf{S}}\right]\leq\|Z_{k}-z^{\ast}\|^{2}_{\mathbf{S}}+2\mathbb{E}_{k}[\bar{\nu}_{k}].$

The next two results guarantee this.

###### Lemma 5

Suppose that Assumption 4 holds true. For all $k\in\mathbb{N}$, we have

$\mathbb{E}_{k}[\|\mathcal{F}_{k}(X_{k})-\nabla f(X_{k})\|^{2}_{\Gamma}]\leq
a_{k}+b_{k}\|Z_{k}-z^{\ast}\|^{2}_{\mathbf{S}}\qquad\text{a.s.}$

where $a_{k}\triangleq
2\sigma^{2}_{0}(x^{\ast})^{2}\frac{\mathtt{M}(\Gamma)}{N_{k}}$ and
$b_{k}\triangleq
2\sigma_{1}^{2}\frac{\mathtt{M}(\Gamma)}{\mathtt{m}(\Gamma)N_{k}}$.

###### Proof:

We compute

$\displaystyle\mathbb{E}_{k}[\|\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k})\|^{2}_{\Gamma}]\leq\mathtt{M}(\Gamma)\mathbb{E}_{k}[\|\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k})\|^{2}]$
$\displaystyle=\mathtt{M}(\Gamma)\mathbb{E}[\|w_{k}(X_{k})\|^{2}|X_{k}]=\mathtt{M}(\Gamma)\tfrac{\sigma^{2}(X_{k})}{N_{k}}\quad\text{a.s.}$

Using Assumption 4, we get

$\displaystyle\sigma^{2}(X_{k})$ $\displaystyle\leq
2\sigma^{2}_{0}(x^{\ast})+2\sigma_{1}^{2}\|X_{k}-x^{\ast}\|^{2}$
$\displaystyle\leq
2\sigma^{2}_{0}(x^{\ast})+2\sigma_{1}^{2}\tfrac{1}{\mathtt{m}(\Gamma)}\|X_{k}-x^{\ast}\|^{2}_{\Gamma}$
$\displaystyle\leq
2\sigma^{2}_{0}(x^{\ast})+2\sigma_{1}^{2}\tfrac{1}{\mathtt{m}(\Gamma)}\|Z_{k}-z^{\ast}\|^{2}_{\mathbf{S}}$

∎

###### Lemma 6

[8, Lemma 5.1] Assume that

$\mathtt{m}(\Gamma^{-1}-\tfrac{\beta_{f}}{2}\mathbf{Q})>\tfrac{\|L\|^{2}}{\mathtt{m}(\Sigma^{-1})}=\|L\|^{2}\mathtt{M}(\Sigma)$
(10)

Then, $2\mathbf{U}-\mathbf{S}\in\mathcal{P}_{\tau}(\mathcal{Z})$ for some
$\tau>0$.

Initialize: $x^{0}_{i}\in\mathcal{H}_{i},y^{0}_{i}\in\mathcal{G}_{i}$ for
$i=1,\ldots,m$ and $w_{(i,j)}\in\mathcal{W}^{\ast}_{(i,j)}$ for all $(i,j)\in
E$.

Iteration $k$: For all agents $i=1,2,\ldots,m$
$\displaystyle\bar{w}_{(i,j),k}^{i}=\tfrac{1}{2}(w^{i}_{(i,j),k}+w^{j}_{(i,j),k})+\tfrac{\tau_{(i,j)}}{2}(A_{ij}x^{i}_{k}+A_{ji}x^{j}_{k}-b_{(i,j)})\quad\forall
j\in\mathcal{N}_{i},$
$\displaystyle\bar{v}^{i}_{k}=\operatorname{Prox}_{\sigma_{i}h^{\ast}_{i}}(v_{k}^{i}+\sigma_{i}L_{i}x^{i}_{k}),$
$\displaystyle
x^{i}_{k+1}=\operatorname{Prox}_{\gamma_{i}g^{\ast}_{i}}[x^{k}_{i}-\gamma\mathcal{F}^{i}_{k}(x^{k})-\gamma_{i}L^{\ast}_{i}\bar{v}_{i}^{k}-\gamma_{i}\textstyle{\sum_{j\in\mathcal{N}_{i}}}A_{ij}w^{i}_{(i,j),k}]$
$\displaystyle
v^{i}_{k+1}=\bar{v}^{i}_{k}+\tau_{(i,j)}(x^{i}_{k+1}-x^{i}_{k})$
$\displaystyle
w^{i}_{(i,j),k+1}=\bar{w}^{i}_{(i,j),k}+\tau_{i,j}A_{i,j}(x^{i}_{k+1}-x^{i}_{k})\quad\forall
j\in\mathcal{N}_{i}.$

Algorithm 2 Distributed STriPD

We can finally prove the main result of this paper.

###### Theorem 1

Let Assumptions 1-4 hold true. Choose $\Gamma,\Sigma$ such that condition (10)
holds. Then, the stochastic process $\\{Z_{k}\\}_{k\in\mathbb{N}}$ generated
by Algorithm 1 converges a.s. to a random variable
$Z_{\infty}\in\operatorname{Zer}(\mathsf{G})$, i.e., to a solution of (1).

###### Proof:

Using Lemma 5, we get

$\displaystyle\mathbb{E}_{k}\left[\|Z_{k+1}-z^{\ast}\|^{2}_{\mathbf{S}}\right]\leq\|Z_{k}-z^{\ast}\|^{2}_{\mathbf{S}}+2\mathbb{E}_{k}[\bar{\nu}_{k}]$
$\displaystyle-\mathbb{E}_{k}\left[\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{2\mathbf{U}-\mathbf{S}}\right]$
$\displaystyle\leq(1+2b_{k})\|Z_{k}-z^{\ast}\|^{2}_{\mathbf{S}}-\mathbb{E}_{k}[\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{2\mathbf{U}-\mathbf{S}}]+2a_{k}.$

Set $v_{k}=\|Z_{k}-z^{\ast}\|_{\mathbf{S}}^{2}$, $t_{k}=2b_{k}$,
$\zeta_{k}=2a_{k}$ and
$\theta_{k}=\mathbb{E}_{k}[\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{2\mathbf{U}-\mathbf{S}}]$,
and apply Lemma 1 to deduce
$\mathbb{P}\left(\sum_{k\in\mathbb{N}}\theta_{k}<\infty\right)=1$ and that
$v_{k}$ converges to a finite random variable $v_{\infty}$ a.s. Furthermore,
using [15, Proposition 2.3], we know that $\\{Z_{k}\\}_{k\in\mathbb{N}}$ is
almost surely bounded. Hence, there exists a measurable set
$\Omega_{0}\subseteq\Omega$ with $\mathbb{P}(\Omega_{0})=1$ such that
$\lim_{k\to\infty}\|Z_{k}(\omega)-z^{\ast}\|=0$ for all $\omega\in\Omega_{0}$,
Fix such an event $\omega\in\Omega_{0}$ and sequence
$\\{k_{j}\\}\subset\mathbb{N}$ with $k_{j}\uparrow\infty$. Consider a
converging subsequence $\\{Z_{k_{j}}(\omega)\\}_{j\in\mathbb{N}}$ with
$Z_{k_{j}}(\omega)\to Z_{\infty}(\omega)$. To reduce notational clutter, we
omit the relabeling and simply denote $Z_{k}(\omega)$ the converging
subsequence. We have to show that
$Z_{\infty}(\omega)\in\operatorname{Zer}(G)$. Set
$\mu\triangleq\mathtt{m}(2\mathbf{U}-\mathbf{S})\geq\tau>0$. We deduce

$\displaystyle\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{2\mathbf{U}-\mathbf{S}}\geq\mu\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}$
$\displaystyle=\mu\left[\|\mathbb{T}_{k}(Z_{k})-T(Z_{k})\|^{2}+\|T(Z_{k})-Z_{k}\|^{2}\right]$
$\displaystyle-2\mu\langle\mathbb{T}_{k}(Z_{k})-T(Z_{k}),Z_{k}-T(Z_{k})\rangle$
$\displaystyle\geq\tfrac{-\mu(1-\alpha)}{\alpha}\|\mathbb{T}_{k}(Z_{k})-T(Z_{k})\|^{2}+(1-\alpha)\mu\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2},$

where in the last inequality we have used the Fenchel-Young inequality. A
simple computation shows that

$\displaystyle\|\mathbb{T}_{k}(Z_{k})-T_{k}(Z_{k})\|^{2}=\|\Sigma_{k}L(\hat{X}_{k}-\bar{x}_{k})\|^{2}+\|\hat{X}_{k}-\bar{x}_{k}\|^{2}$
$\displaystyle\leq\left(1+\mathtt{M}(L^{\ast}\Sigma^{2}L)\right)\|\hat{X}_{k}-\bar{x}_{k}\|^{2}$
$\displaystyle\leq\tfrac{\left(1+\mathtt{M}(L^{\ast}\Sigma^{2}L)\right)}{\mathtt{m}(\Gamma^{-1})}\|\hat{X}_{k}-\bar{x}_{k}\|^{2}_{\Gamma^{-1}}$
$\displaystyle\leq\tfrac{\left(1+\mathtt{M}(L^{\ast}\Sigma^{2}L)\right)}{\mathtt{m}(\Gamma^{-1})}\|\mathcal{F}_{k}(X_{k})-\nabla
f(X_{k})\|^{2}_{\Gamma}.$

Hence, for $\alpha\in(0,1)$, we get

$\displaystyle\mathbb{E}_{k}[\|\mathbb{T}_{k}(Z_{k})-Z_{k}\|^{2}_{2\mathbf{U}-\mathbf{S}}]\geq(1-\alpha)\mu\|T_{k}(Z_{k}(\omega))-Z_{k}(\omega)\|^{2}$
$\displaystyle-\tfrac{\mu(1-\alpha)}{\alpha}\tfrac{\left(1+\mathtt{M}(L^{\ast}\Sigma^{2}L)\right)}{\mathtt{m}(\Gamma^{-1})}\mathbb{E}_{k}[\bar{\nu}_{k}](\omega).$

By continuity of the mappings $T(\cdot)$, it follows

$0=\lim_{k\to\infty}\|T(Z_{k}(\omega))-Z_{k}(\omega)\|^{2}=\|T(Z_{\infty}(\omega))-Z_{\infty}(\omega)\|.$

Hence, $Z_{\infty}\in\operatorname{Fix}(T)=\operatorname{Zer}(\mathsf{G})$. As
$\omega\in\Omega_{0}$ is arbitrary, the claim follows. ∎

##  IV Distributed Optimization

In this section we consider a network of agents, whose aim is to solve problem
(3) in a cooperative way. Consider an undirected graph $G=(V,E)$ over a vertex
set $V=\\{1,\ldots,m\\}$ with edge set $E\subset V\times V$. Each vertex is
associated with an agent, which is assumed to have a local memory and
computational unit and can only communicate with its neighbors. We define the
neighborhood of agent $i$ as $\mathcal{N}_{i}=\\{j\in V|(i,j)\in E\\}$. Each
agent in the network is characterized by a private cost function
$f_{i}(x_{i})+g_{i}(x_{i})+h_{i}(L_{i}x_{i}),$ defined over the vector space
$\mathcal{H}_{i}\equiv\mathbb{R}^{p_{i}}$. Let
$\mathcal{H}=\prod_{i=1}^{m}\mathcal{H}_{i}$. Furthermore, the decisions of
the agents in the network are subject to affine constraints, coupling the
decisions of neighboring agents.

###### Assumption 5

For each $i=1,\ldots,m$:

  * (i)

For all $j\in\mathcal{N}_{i}$,
$b_{(i,j)}\in\mathcal{E}_{(i,j)}\equiv\mathbb{R}^{l_{ij}}$ and
$A_{ij}\in\mathbb{R}^{p_{i}\times l_{ij}}$;

  * (ii)

$g_{i}:\mathbb{R}^{p_{i}}\to(-\infty,\infty],h_{i}:\mathbb{R}^{q_{i}}\to(-\infty,\infty]$
are proper closed convex and lower semi-continuous functions, and
$L_{i}\in\mathbb{R}^{p_{i}\times q_{i}}$;

  * (iii)

$f_{i}:\mathcal{H}_{i}\to\mathbb{R}$ is a convex, continuously differentiable
and for some $\beta_{i}\geq 0$, $\nabla f_{i}$ is $\beta_{i}$-Lipschitz
continuous in the norm $\|\cdot\|_{\mathbf{Q}_{i}}$;

  * (iv)

The graph $G$ is connected;

  * (v)

The set of solutions of (3) is nonempty and there exists
$x_{i}\in\operatorname{relint}(\operatorname{dom}g_{i})$ such that
$L_{i}x_{i}\in\operatorname{relint}(\operatorname{dom}h_{i})$ and
$A_{ij}x_{i}+A_{ji}x_{j}=b_{(i,j)}$ for $(i,j)\in E$.

Let $\mathcal{W}_{(i,j)}\triangleq\mathbb{R}^{2l_{ij}}$ for all $(i,j)\in E$,
with generic element $w_{ij}=[w_{(i,j)}^{i},w_{(i,j)}^{j}]$. The
interpretation is that $w_{(i,j)}^{i}$ is controlled by agent $i$ and
$w_{(i,j)}^{j}$ is controlled by agent $j$. For each edge $(i,j)\in E$ define
the set

$\mathsf{K}_{(i,j)}\triangleq\\{(w_{1},w_{2})\in\mathcal{W}_{(i,j)}|w_{1}+w_{2}=b_{(i,j)}\\},$

and the linear map $\mathsf{N}_{(i,j)}:\mathcal{H}\to\mathbb{R}^{2l_{ij}}$ by

$\langle\mathsf{N}_{(i,j)}x,w_{ij}\rangle=\langle
A_{ij}x_{i},w_{ij}^{i}\rangle+\langle
A_{ji}x_{j},w_{ij}^{j}\rangle\qquad\forall x\in\mathcal{H}.$

Accordingly, we let $\mathsf{N}:\mathcal{H}\to\prod_{(i,j)\in
E}\mathbb{R}^{2l_{ij}}$ to be the operator defined by
$\langle\mathsf{N}x,w\rangle=\sum_{(i,j)\in E}(\langle
A_{ij}x_{i},w_{(i,j)}^{i}\rangle+\langle A_{ji}x_{j},w_{(i,j)}^{j}\rangle).$
Using these concepts, we can reformulate problem (3) as

$\min_{x_{1},\ldots,x_{m}}\sum_{i=1}^{m}[f_{i}(x_{i})+g_{i}(x_{i})+h_{i}(L_{i}x_{i})]+\sum_{(i,j)\in
E}\delta_{\mathsf{K}_{(i,j)}}(\mathsf{N}_{(i,j)}x).$

Let $\mathsf{K}\triangleq\prod_{(i,j)\in E}\mathsf{K}_{(i,j)}$ and
$L:\mathcal{H}\to\prod_{i}\mathbb{R}^{q_{i}}$ be defined by
$Lx=[L_{1}x_{1},\ldots,L_{m}x_{m}]$. Set $Dx=[Lx,\mathsf{N}x]$. Define the
functions $F(x)\triangleq\sum_{i=1}^{m}f_{i}(x_{i})$,
$G(x)\triangleq\sum_{i=1}^{m}g_{i}(x_{i})$ and $H(Dx)\triangleq
h(Lx)+\delta_{\mathsf{K}}(\mathsf{N}x)$. With this notation, we have converted
problem (3) to problem (1) with the functions $F,G$ and $H$.

As in Section II, the primal-dual optimality conditions can be written in the
compact form as an inclusion problem involving the (maximally monotone)
operators

$\displaystyle\mathsf{A}(v,w,x)=[\partial
h^{\ast}(v),\partial\delta_{\mathsf{C}}^{\ast}(w),\partial g(x)],$
$\displaystyle\mathsf{M}(v,w,x)=[-Lx,-\mathsf{N}x,L^{\ast}v+\mathsf{N}^{\ast}w],$
$\displaystyle\mathsf{C}(v,w,x)=[0,0,\nabla f(x)].$

The dual variable is the pair $y=[v,w]\in\mathcal{G}$ and we can apply
Algorithm 1 directly to solve the distributed optimization problem (3). The
resulting stochastic process $\\{Z_{k}\\}_{k\geq 0}$ decomposes to agent-
specific updates as described in Algorithm 2 (cf. [7]). Given the
identification of the operators characterizing the optimality conditions of a
primal-dual pair, the convergence of the synchronous STriPD (Algorithm 2)
follows immediately from Theorem 1.

(a) Distance from the solution.

(b) Distance from minimum cost.

(c) Feasibility of the iterates.

##  V Numerical Example

We test Algorithm 1 on an economic dispatch problem for power grids, inspired
by [16, 3]. Consider $n$ control areas, each with a generator that supply
power $x_{i}\in\mathbb{R}$ and a local demand $b_{i}\in\mathbb{R}_{+}$ that
should be satisfied. Each generator has security bound of the form
$\underline{x}_{i}\leq x_{i}\leq\bar{x}_{i}$, with
$\underline{x}_{i},\bar{x}_{i}\geq 0$ for all $i=1,\dots,m$. Each area has a
local generation cost $f_{i}(x_{i}):\mathbb{R}\to\mathbb{R}$ so that the
optimization problem is to minimize the overall cost
$f(x)=\sum_{i=1}^{m}f_{i}(x_{i})$, subject to the satisfaction of the demand:

$\begin{cases}\min\limits_{x_{i},i=1,\dots,m}&f(x)=\sum_{i=1}^{m}f_{i}(x_{i})\\\
\quad\;\text{s.t.}&\sum_{i=1}^{m}x_{i}=\sum_{i=1}^{m}b_{i}\\\
&\underline{x}_{i}\leq x_{i}\leq\bar{x}_{i},i=1,\dots,m\end{cases}$ (11)

To write the problem in (11) in the form of problem (3), let us take
$g_{i}(x_{i})=\delta_{X_{i}}(x_{i})$, i.e., the indicator function of the
local constraints $X_{i}\triangleq[\underline{x}_{i},\bar{x}_{i}]$, and
$h(Lx)=\delta_{C}(Lx)$, where $L=\mathbf{1}_{m}$ and
$C=\\{y\in\mathbb{R}^{m}:y=b\\}$ represents the coupling constraints.

Similarly to [3], we consider $m=5$ generators with cost functions
$f_{i}(x_{i})=\mathbb{E}[q_{i}(\xi)x_{i}^{2}+p_{i}x_{i}]$ where the
coefficients are chosen with mean $q=[0.094,0.078,0.105,0.082,0.074]$ and
$p=[1.22,3.41,2.53,4.02,3.17]$. The random variable is drawn according to a
normal distribution with mean $q$. The local bounds on the supply power are
$\underline{x}=[10,8,3.8,5.4,4.2]$ and $\bar{x}=[80,60,40,45,18]$ while the
demands take the values $b=[35,20,25,30,10]$. We run the algorithm 100 times
and plot in a thick blue line the average results; the transparent areas
indicate the minimum and maximum values reached. Figure 1(a) displays the
distance of the iterates from the solution and Fig 1(b) illustrates the
distance of the cost from the optimal value. In Figure 1(c), we show that
asymptotically the constraints are satisfied.

##  VI Conclusions

We propose a stochastic triangular preconditioned primal-dual algorithm
(STriPD) for solving a large family of structured convex optimization
problems, and multi-agent versions thereof. Many extensions of the present
work will be investigated in a more elaborate investigation. Such extensions
will include block-coordinate descent implementations, allowing for
asynchronous updates in the distributed case. Furthermore, we will be
investigating the iteration and oracle complexity of the method in detail.

## References

  * [1] B. Jin, H. Li, W. Yan, and M. Cao, “Distributed model predictive control and optimization for linear systems with global constraints and time-varying communication,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 7, pp. 3393–3400, 2021.
  * [2] T. Hastie, R. Tibshirani, and M. Wainwright, _Statistical learning with sparsity: the lasso and generalizations_. Chapman and Hall/CRC, 2015.
  * [3] H. Li, E. Su, C. Wang, J. Liu, Z. Zheng, Z. Wang, and D. Xia, “A primal-dual forward-backward splitting algorithm for distributed convex optimization,” _IEEE Transactions on Emerging Topics in Computational Intelligence_ , 2021\.
  * [4] A. Iusem, A. Jofré, R. I. Oliveira, and P. Thompson, “Extragradient method with variance reduction for stochastic variational inequalities,” _SIAM Journal on Optimization_ , vol. 27, no. 2, pp. 686–724, 2017.
  * [5] ——, “Variance-based extragradient methods with line search for stochastic variational inequalities,” _SIAM Journal on Optimization_ , vol. 29, no. 1, 2019.
  * [6] R. I. Boţ, P. Mertikopoulos, M. Staudigl, and P. T. Vuong, “Minibatch forward-backward-forward methods for solving stochastic variational inequalities,” _Stochastic Systems_ , vol. 11, no. 2, pp. 112–139, 2021.
  * [7] P. Latafat, N. M. Freris, and P. Patrinos, “A new randomized block-coordinate primal-dual proximal algorithm for distributed optimization,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 10, pp. 4050–4065, 2019.
  * [8] P. Latafat and P. Patrinos, “Asymmetric forward-backward-adjoint splitting for solving monotone inclusions involving three operators,” _Computational Optimization and Applications_ , vol. 68, no. 1, pp. 57–93, 2017.
  * [9] A. Yurtsever, B. C. Vu, and V. Cevher, “Stochastic three-composite convex minimization,” _Advances in Neural Information Processing Systems_ , vol. 29, 2016.
  * [10] R. Zhao and V. Cevher, “Stochastic three-composite convex minimization with a linear operator,” in _Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, A. Storkey and F. Perez-Cruz, Eds., vol. 84. PMLR, 09–11 Apr 2018, pp. 765–774.
  * [11] Y. Chen, G. Lan, and Y. Ouyang, “Optimal primal-dual methods for a class of saddle point problems,” _SIAM Journal on Optimization_ , vol. 24, no. 4, pp. 1779–1814, 2014.
  * [12] H. Robbins and D. Siegmund, _A convergence theorem for non negative almost supermartingales and some applications_. Academic Press, 1971, pp. 233–257.
  * [13] H. H. Bauschke and P. L. Combettes, _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_. Springer - CMS Books in Mathematics, 2016.
  * [14] A. Jofré and P. Thompson, “On variance reduction for stochastic smooth convex optimization with multiplicative noise,” _Mathematical Programming_ , vol. 174, no. 1, pp. 253–292, 2019.
  * [15] P. L. Combettes and J.-C. Pesquet, “Stochastic quasi-fejér block-coordinate fixed point iterations with random sweeping,” _SIAM Journal on Optimization_ , vol. 25, no. 2, pp. 1221–1248, 2015.
  * [16] P. Yi, Y. Hong, and F. Liu, “Initialization-free distributed algorithms for optimal resource allocation with feasibility constraints and application to economic dispatch of power systems,” _Automatica_ , vol. 74, pp. 259–269, 2016.