# Sequential Sample Average Majorization–Minimization

Gersende Fort, Florence Forbes, Hien Duy Nguyen

► To cite this version:

Gersende Fort, Florence Forbes, Hien Duy Nguyen. Sequential Sample Average Majorization–Minimization. 2024. hal-04607609

HAL Id: hal-04607609

https://hal.science/hal-04607609

Preprint submitted on 10 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sequential Sample Average Majorization–Minimization

Gersende Fort, Florence Forbes, Hien Duy Nguyen

► To cite this version:

Gersende Fort, Florence Forbes, Hien Duy Nguyen. Sequential Sample Average Majorization–Minimization. Statistics and Computing, 2025, 36, pp.31. 10.1007/s11222-025-10780-x. hal-04607609v2

HAL Id: hal-04607609

https://hal.science/hal-04607609v2

Submitted on 6 Sep 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Sequential Sample Average Majorization–Minimization

Gersende Fort<fnref n="1" />, Florence Forbes<fnref n="2" /> and Hien Duy Nguyen<fnref n="3" />

June 10, 2024

**Abstract**

Many statistical inference and machine learning methods rely on the ability to optimize an expectation functional, whose explicit form is intractable. The typical method for conducting such optimization is to approximate the expected value problem by a size-*N* sample average, often referred to as sample average approximation (SAA) or M-estimation. When the solution to the SAA problem cannot be obtained in closed form, the majorization-minimization (MM) algorithm framework constitutes a broad class of incremental optimization solutions, relying on the iterative construction of surrogates, known as majorizers, of the original problem. The ability to solve an SAA problem depends on the availability of all *N* observations, contemporaneously, which is difficult when *N* is large or data are observed as a stream. We propose a stochastic MM algorithm that solves the expected value problem via iterative SAA majorizer constructions using sequential subsets of data, which we call *Sequential Sample Average Majorization-Minimization (SAM2)*. Compared to previous stochastic MM algorithm variants, our method permit an extended definition of majorizers, and does not rely on convexity and smoothness assumptions or make functional restrictions on the class of problems and majorizers. We develop a theory of stochastic convergence for *SAM2*, made possible via the presentation of a novel double array uniform strong law of large numbers. Examples of *SAM2* algorithms are given along with a numerical demonstration of *SAM2* to the quantile regression problem.

<footnote n="1">¹CNRS and Institut de Mathématiques de Toulouse, UMR5219; Université de Toulouse; UPS, F-31062 Toulouse Cedex 9, France.</footnote>

<footnote n="2">²Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK 38000 Grenoble, France.</footnote>

<footnote n="3">³School of Computing, Engineering and Mathematical Sciences, La Trobe Univ., Bundoora 3086, Victoria Australia; Institute of Mathematics for Industry, Kyushu Univ., Nishi Ward, Fukuoka 819-0395, Japan.</footnote>

<footnote>All authors thank the ARC grant DP230100905 and the Wombat Inria Associated Team Project. Part of this work (GF) is funded by the Fondation Simone and Cino Del Duca under the project OpSiMorE.</footnote>
Sequential Sample Average Majorization–Minimization

Gersende Fort<fnref n="1" />,<fnref n="2" />, Florence Forbes<fnref n="3" /> and Hien Duy Nguyen<fnref n="4" />,<fnref n="5" />

September 6, 2025

**Abstract**

Many statistical inference and machine learning methods rely on the ability to optimize an expectation functional, whose explicit form is intractable. The typical method for conducting such optimization is to approximate the expected value problem by a size-$N$ sample average, often referred to as sample average approximation (SAA) or M-estimation. When the solution to the SAA problem cannot be obtained in closed form, the Majorization–Minimization (MM) algorithm framework constitutes a broad class of incremental optimization solutions, relying on the iterative construction of surrogates, known as majorizers, of the original problem. The ability to solve an SAA problem depends on the availability of all $N$ observations, contemporaneously, which is difficult when $N$ is large or data are observed as a stream. We propose a stochastic MM algorithm that solves the expected value problem via iterative SAA majorizer constructions using sequential subsets of data, which we call *Sequential Sample Average Majorization–Minimization (SAM2)*. Compared to previous stochastic MM algorithm variants, our method permit an extended definition of majorizers, and does not rely on convexity assumptions, smoothness assumptions, or restrictions on functional classes for objectives and majorizers. We develop a theory of stochastic convergence for *SAM2*, made possible via the presentation of a novel double array uniform strong law of large numbers. Examples of *SAM2* algorithms are given along with a numerical demonstration of *SAM2* to quantile regression problems, in the regular and sparse parameter settings, including both convex and non-convex objective functions.

<footnote n="1">¹IMT, Université de Toulouse, CNRS, Toulouse, France.</footnote>

<footnote n="2">²LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France.</footnote>

<footnote n="3">³Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK 38000 Grenoble, France.</footnote>

<footnote n="4">⁴School of Computing, Engineering and Mathematical Sciences, La Trobe Univ., Bundoora 3086, Victoria Australia.</footnote>

<footnote n="5">⁵Institute of Mathematics for Industry, Kyushu Univ., Nishi Ward, Fukuoka 819-0395, Japan.</footnote>

<footnote>All authors thank the ARC grant DP230100905 and the Wombat Inria Associated Team Project. Part of this work (GF) is funded by the Fondation Simone and Cino Del Duca under the project OpSiMorE.</footnote>
# 1 Introduction

A common task that often arises when conducting statistical inference or machine learning is to solve the problem of obtaining

$$ \mathcal{L}_{F^*} := \operatorname{argmin}_{\theta \in \mathbb{T}} F(\theta), \qquad (1) $$

where $\mathbb{T}$ is the so-called parameter space, with typical element $\theta$, and $F: \mathbb{T} \to \mathbb{R}$ is an objective function, defined at each $\theta \in \mathbb{T}$ by

$$ F(\theta) := \mathbb{E}[f(\theta, X)], $$

the expectation of $f(\theta, X)$, for some random $X$ in some space $\mathcal{X}$. Typically, the underlying probability measure of $X$ is unknown and thus $F$ cannot be accessed, directly. However, given access to a sample of independent and identically distributed (i.i.d.) replicates of $X: \{X^i, 1 \le i \le N\}$, one can approximate problem (1) by the realizable problem of obtaining

$$ \mathcal{L}_N := \operatorname{argmin}_{\theta \in \mathbb{T}} \frac{1}{N} \sum_{i=1}^{N} f(\theta, X_i). \qquad (2) $$

The approximation of (1) by (2) has appeared throughout the literature under many names, including extremum estimation (Amemiya, 1985; Gourieroux & Monfort, 1995), M-estimation (Serfling, 1980; van der Vaart & Wellner, 2023), minimum contrast estimation (Dacunha-Castelle & Duflo, 1986; Bickel & Doksum, 2015), in statistics, empirical risk minimization (Vapnik, 1998; Vidyasagar, 2003; Shalev-Shwartz & Ben-David, 2014), in machine learning, and sample average approximation (SAA; Bonnans 2019; Shapiro et al. 2021; Cui & Pang 2022) in optimization theory.

In general, (2) is not assumed to be solvable in closed form, and thus iterative optimization routines are required in practice. Given various assumptions regarding the smoothness or convexity of $f(\cdot, x)$ ($x \in \mathcal{X}$), numerous methods are available for solving (2), including Newton's algorithm and (sub-)gradient descent (see, e.g., Polyak 2021), and derivative-free methods, such as the Nelder-Mead algorithm and pattern search (see, e.g., Audet & Hare 2017), among other techniques. However, such techniques often require the parameter of the problem be vector-valued, which is unnatural in many learning settings, where $\mathbb{T}$ consists of elements that take value in the probability simplex or in the positive definite matrices, as in the problem of estimating finite mixtures of Gaussian distributions, for example (cf. McLachlan & Peel 2000).

**EM and MM algorithms** The Expectation-Maximization (EM) approach of Dempster et al. (1977) (see also McLachlan & Krishnan 2008) provides a framework for conducting optimization of probabilistic models in potentially non-vector spaces, via the construction of numerical surrogates for $f$, whose optimization generates a stable sequence of parameter estimates that monotonically improve the objective values. The EM approach can be viewed as a special case of the so-called Majorization-Minimization (MM) algorithms (Hunter & Lange, 2004; Lange, 2013, 2016), which
# 1 Introduction

A common task that often arises when conducting statistical inference or machine learning is obtaining

$$ \mathcal{L}_{F^*} := \operatorname{argmin}_{\theta \in \mathbb{T}} F(\theta), \qquad (1) $$

where $\mathbb{T}$ is the so-called parameter space, with typical element $\theta$, and $F: \mathbb{T} \to \mathbb{R}$ is an objective function, defined at each $\theta \in \mathbb{T}$ by

$$ F(\theta) := \mathbb{E}[f(\theta, X)], $$

the expectation of $f(\theta, X)$, for some random variable $X$ in some space $\mathcal{X}$. Typically, the underlying probability measure of $X$ is unknown and thus $F$ cannot be accessed, directly. However, given access to a sample of independent and identically distributed (i.i.d.) replicates of $X$: $\{X^i, 1 \le i \le N\}$, problem of obtaining

$$ \mathcal{L}_N := \operatorname{argmin}_{\theta \in \mathbb{T}} \frac{1}{N} \sum_{i=1}^{N} f(\theta, X^i). \qquad (2) $$

The approximation of (1) by (2) has appeared throughout the literature under many names, including extremum estimation (Amemiya, 1985; Gourieroux & Monfort, 1995), M-estimation (Serfling, 1980; van der Vaart & Wellner, 2023), minimum contrast estimation (Dacunha-Castelle & Duflo, 1986; Bickel & Doksum, 2015), in statistics; empirical risk minimization (Vapnik, 1998; Vidyasagar, 2003; Shalev-Shwartz & Ben-David, 2014), in machine learning; and sample average approximation (SAA; Bonnans 2019; Shapiro et al. 2021; Cui & Pang 2022) in optimization theory.

In general, (2) is not assumed to be solvable in closed form, and thus iterative optimization routines are required in practice. Given various assumptions regarding the smoothness or convexity of $f(\cdot, x)$ ($x \in \mathcal{X}$), numerous methods are available for solving (2), including Newton's algorithm and (sub-)gradient descent (see, e.g., Polyak 2021), and derivative-free methods, such as the Nelder-Mead algorithm and pattern search (see, e.g., Audet & Hare 2017), among other techniques. However, such techniques often require the parameter of the problem be vector-valued, which is unnatural in many learning settings, where $\mathbb{T}$ consists of elements that take value in the probability simplex or in the positive definite matrices, as in the problem of estimating finite mixtures of Gaussian distributions, for example (cf. McLachlan & Peel 2000).

**EM and MM algorithms** The Expectation-Maximization (EM) approach of Dempster et al. (1977) (see also McLachlan & Krishnan 2008) provides a framework for conducting optimization of probabilistic models in potentially non-vector spaces, via the construction of numerical surrogates for $f$, whose optimization generates a stable sequence of parameter estimates that monotonically improve the objective values. The EM approach can be viewed as a special case of the so-called Majorization-Minimization (MM) algorithms (Hunter & Lange, 2004; Lange, 2013, 2016), which extend the applicability of the EM algorithm beyond probabilistic models, to more general surro-
extend the applicability of the EM algorithm beyond probabilistic models, to more general surrogate constructions for $f$, to include any function of the form $g: \mathbb{T} \times \mathcal{X} \times \mathbb{T} \rightarrow \mathbb{R}$, fulfilling criterion $f(\theta, x) \le g(\theta, x; v)$, for any $\theta, v \in \mathbb{T}$ and $x \in \mathcal{X}$, with equality whenever $\theta = v$. This permits the inclusion of gradient descent, proximal algorithms and proximal gradient descent, quadratic approximation algorithm (Bohning & Lindsay, 1988), and the convex concave procedure (Yuille & Rangarajan, 2003), among other methods, within the MM algorithm framework. Furthermore, the framework has been extended to permit coordinate and block-wise iterations (Meng & Rubin, 1993; Razaviyayn et al., 2013; Chalvidal et al., 2023), stochastic surrogates (Celeux & Diebolt, 1992; Delyon et al., 1999; Fort & Moulines, 2003), and decentralized computation (Dieuleveut et al., 2021; Cadoni et al., 2016).

Beyond the applicability, inclusivity, and extendability, the MM framework has a strong unified theoretical foundation, with numerous available global convergence and convergence rates results available. Generally applicable results in this direction include the global convergence theory of Vaida (2005), Lange (2013), Byrne (2014), Lange (2016), Lange et al. (2021), and Cui & Pang (2022), as well as the block-wise results of Razaviyayn et al. (2013), and the convergence rates of Mairal (2015), Chouzenoux & Pesquet (2016), and Hong et al. (2017).

**Stochastic algorithms** As described, the MM algorithms above solve the sample problem (2), for fixed sample size $N$ to produce an estimate for the solution of the expectation problem (1). Via the general SAA theory of Shapiro et al. (2021), it can be shown that the set of solutions of $\mathcal{L}_N$ converges in set deviation to $\mathcal{L}_{F^*}$, almost surely, as $N$ gets large. However, the process of solving repeated SAAs as $N$ increases assumes the availability of the entirety of the sequence $\{X^i, 1 \le i \le N\}$, for each $N$, which is infeasible in practice when data are often accessed as a stream, and when computer memory is so as to make storage of the entire sequence of data impossible for large $N$.

A more feasible scenario is that one has limited access to the data set, via only the subset $\{X^{t,i}, i = 1, \dots, N_t\}$ at time $t$, where $\{N_t, t \ge 1\}$ is a sequence of positive integers denoting sample sizes. Such situations arise in the so-called online or iterative algorithms setting, and approaches for solving (1) in such situations are most notably exemplified by the popular stochastic gradient descent algorithms and their variants and refinements, as studied in Shalev-Shwartz & Ben-David (2014), Lan (2020), Lin et al. (2020), Shapiro et al. (2021), among an ever expanding body of literature. As per their deterministic variants, these algorithms require parameters be vector-valued, among other restrictions, and thus are not universal in their applicability. These limitations are similarly shared by stochastic Newton algorithms such as those of Schraudolph et al. (2007), Byrd et al. (2016), and Meyn (2022).

To resolve these limitations, one can again turn to stochastic variants of EM and MM algorithms, which have been comprehensively studied over the years. Examples of such works include the pioneering work of Cappé & Moulines (2009), whose online EM algorithm forms the basis of
gate constructions for $f$, to include any function of the form $g: \mathbb{T} \times \mathcal{X} \times \mathbb{T} \rightarrow \mathbb{R}$, fulfilling criterion $f(\theta, x) \le g(\theta, x; v)$, for any $\theta, v \in \mathbb{T}$ and $x \in \mathcal{X}$, with equality whenever $\theta = v$. This permits the inclusion of gradient descent, proximal algorithms and proximal gradient descent, quadratic approximation algorithm (Bohning & Lindsay, 1988), and the convex-concave procedure (Yuille & Rangarajan, 2003), among other methods, within the MM algorithm framework. Furthermore, the framework has been extended to permit coordinate and block-wise iterations (Meng & Rubin, 1993; Razaviyayn et al., 2013; Chalvidal et al., 2023), stochastic surrogates (Celeux & Diebolt, 1992; Delyon et al., 1999; Fort & Moulines, 2003), and decentralized computation (Dieuleveut et al., 2021; Cadoni et al., 2016).

Beyond the applicability, inclusivity, and extendability, the MM framework has a strong unified theoretical foundation, with numerous available global convergence and convergence rates results available. Generally applicable results in this direction include the global convergence theory of Vaida (2005), Lange (2013), Byrne (2014), Lange (2016), Lange et al. (2021), and Cui & Pang (2022), as well as the block-wise results of Razaviyayn et al. (2013), and the convergence rates of Mairal (2015), Chouzenoux & Pesquet (2016), and Hong et al. (2017).

**Stochastic algorithms** As described, the MM algorithms above solve the sample problem (2), for fixed sample size $N$ to produce an estimate for the solution of the expectation problem (1). Via the general SAA theory of Shapiro et al. (2021), it can be shown that the set of solutions of $\mathcal{L}_N$ converges in set deviation to $\mathcal{L}_{F^*}$, almost surely, as $N$ gets large. However, the process of solving repeated SAAs as $N$ increases assumes the availability of the entirety of the sequence $\{X^i, 1 \le i \le N\}$, for each $N$, which is infeasible in practice when data are often accessed as a stream, and when computer memory is so as to make storage of the entire sequence of data impossible for large $N$.

A more feasible scenario is that one has limited access to the data set, via only the subset $\{X^{t,i}, 1 \le i \le N_t\}$ at time $t$, where $\{N_t, t \ge 1\}$ is a sequence of positive integers denoting sample sizes. Such situations arise in the so-called online or iterative algorithms setting, and approaches for solving (1) in such situations are most notably exemplified by the popular stochastic gradient descent algorithms and their variants and refinements, as studied in Shalev-Shwartz & Ben-David (2014), Lan (2020), Lin et al. (2020), Shapiro et al. (2021), among an ever expanding body of literature. As per their deterministic variants, these algorithms require parameters to be vector-valued, among other restrictions, and thus are not universal in their applicability. These limitations are similarly shared by stochastic Newton algorithms such as those of Schraudolph et al. (2007), Byrd et al. (2016), and Meyn (2022).

To resolve these limitations, one can again turn to stochastic variants of EM and MM algorithms, which have been comprehensively studied over the years. Examples of such works include the pioneering work of Cappé & Moulines (2009), whose online EM algorithm forms the basis of the lineage of works including those of Karimi et al. (2019a), Karimi et al. (2019b), Nguyen et al.
the lineage of works including those of Karimi et al. (2019b), Karimi et al. (2019a), Dieuleveut et al. (2021), Karimi & Li (2021), Fort et al. (2021a), and Fort et al. (2021b) regarding the online, stochastic, and mini-batch estimation of probabilistic models. The online EM framework of Cappé & Moulines (2009) was also extended to optimization of general models via the online MM extension of Nguyen et al. (2022). In each of these works, the surrogates of $f$ are required to take a restrictive linearized form, which is required for obtaining theoretical guarantees.

Next, we note the original works of Mairal (2013) and Razaviyayn et al. (2016), who laid the foundation for a convexity-based approaches to theory, where the convergence of the algorithm is obtained via strong convexity and Lipschitz smoothness assumptions on the surrogate of $f$, with $f$ required to be convex or strongly convex to obtain the convergence rates. Other works following this direction include the works of Liu et al. (2018), Liu et al. (2019), Zhang et al. (2019), Mokhtari & Koppel (2020), Chouzenoux & Fest (2022), Karimi et al. (2022), and Lupu & Necoara (2023).

Finally, we note the approach of Cui & Pang (2022) who exchange the convexity and Lipschitz smoothness of the surrogates with the inclusion of a proximal term. This approach has been further developed in Liu et al. (2022) and Liu & Pang (2023).

**Current work and contributions** Taking the approach of Cui & Pang (2022) as a starting point, we derive a novel approach for solving problem (1) via a stream of data, using an MM algorithm approach. We call our method *Sequential Sample Average Majorization-Minimization*, or **SAM2** for brevity.

Like the method of Cui & Pang (2022), our approach solves an iterative sequence of optimization problems, characterized by surrogate functions of $f$, using subsamples $\{X^{t,i}, i = 1, \dots, N_t\}$ at each iteration $t$. However, unlike Cui & Pang (2022), our method does not require a proximal term be included as part of the surrogate function, and thus generalizes the approach, since an MM surrogate functions are closed under the addition of a proximal term, as per Cui & Pang (2022), or a Bregman proximal term, as per Lange (2016), Rossignol et al. (2022), and Khanh Hien et al. (2022). As such **SAM2** can be viewed as a direct generalization of the Cui & Pang (2022) method. We develop a theory of the almost-sure convergence of **SAM2**, when the number of iterations tends to infinity. Our theory makes no use of convexity and Lipschitz smoothness assumptions nor restrictions on the form of the surrogate functions and thus provides the least restrictive framework for constructing stochastic MM algorithms, to the best of our knowledge.

In addition to less restrictive assumptions on MM surrogate functions, our approach permits an expanded definition of such surrogates, in that one may take $v \in \mathbb{U}$ instead of $v \in \mathbb{T}$, for some arbitrary set $\mathbb{U}$. This permits for a broader interpretation of the MM algorithm framework, following in the spirit of Mairal (2015) who consider the allowance for surrogates to depend on iteration $t$, and Naderi et al. (2019) who relax the requirement that the surrogate satisfies $f(\theta, X) = g(\theta, X; v)$ for some $v$. In particular, the expanded definition allows us to solve problem (1), with $f$ only implicitly defined via the condition that for every $\theta \in \mathbb{T}$, $F(\theta) \le \mathbb{E}[g(\theta, X; v)]$ for all $v \in \mathbb{U}$,
(2020), Dieuleveut et al. (2021), Karimi & Li (2021), Fort et al. (2021a), and Fort et al. (2021b) regarding the online, stochastic, and minibatch estimation of probabilistic models. The online EM framework of Cappé & Moulines (2009) was also extended to optimization of general models via the online MM extension of Nguyen et al. (2022). In each of these works, the surrogates of f are required to take a restrictive linearized form, which is required for obtaining theoretical guarantees.

Next, we note the original works of Mairal (2013) and Razaviyayn et al. (2016), who laid the foundation for a convexity-based approaches to theory, where the convergence of the algorithm is obtained via strong convexity and Lipschitz smoothness assumptions on the surrogate of f, with f required to be convex or strongly convex to obtain the convergence rates. Other works following this direction include the works of Liu et al. (2018), Liu et al. (2019), Zhang et al. (2019), Mokhtari & Koppel (2020), Chouzenoux & Fest (2022), Karimi et al. (2022), and Lupu & Necoara (2023).

Finally, we note the approach of Cui & Pang (2022) who exchange the convexity and Lipschitz smoothness of the surrogates with the inclusion of a proximal term. This approach has been further developed in Liu et al. (2022) and Liu & Pang (2023).

**Current work and contributions** Taking the approach of Cui & Pang (2022) as a starting point, we derive a novel approach for solving problem (1) via a stream of data, using an MM algorithm approach. We call our method *Sequential Sample Average Majorization-Minimization*, or SAM2 for brevity.

Like the method of Cui & Pang (2022), our approach solves an iterative sequence of optimization problems, characterized by surrogate functions of f, using subsamples {$X^{t,i}, 1 \le i \le N_t$} at each iteration t. However, unlike Cui & Pang (2022), our method does not require a proximal term be included as part of the surrogate function, and thus generalizes the approach, since MM surrogate functions are closed under the addition of a proximal term, as per Cui & Pang (2022), or a Bregman proximal term, as per Lange (2016), Rossignol et al. (2022), and Khanh Hien et al. (2022). As such SAM2 can be viewed as a direct generalization of the Cui & Pang (2022) method. We develop a theory of the almost-sure convergence of SAM2, when the number of iterations tends to infinity. Our theory makes no use of convexity and Lipschitz smoothness assumptions nor restrictions on the form of the surrogate functions and thus provides the least restrictive framework for constructing stochastic MM algorithms, to the best of our knowledge.

In addition to less restrictive assumptions on MM surrogate functions, our approach permits an expanded definition of such surrogates, in that one may take $v \in \mathbb{U}$ instead of $v \in \mathbb{T}$, for some arbitrary set $\mathbb{U}$. This permits for a broader interpretation of the MM algorithm framework, following in the spirit of Mairal (2015) who considers the allowance for surrogates to depend on iteration $t$, and Naderi et al. (2019) who relaxes the requirement that the surrogate satisfies $f(\theta, \cdot) = g(\theta, \cdot; v)$ for some $v$ dependent upon $\theta$. In particular, the expanded definition allows us to solve problem (1), with $f$ only implicitly defined via the condition that for every $\theta \in \mathbb{T}$, $F(\theta) \le \mathbb{E}[g(\theta, X; v)]$ for all $v \in \mathbb{U}$, and there exists a $v \in \mathbb{U}$ such that $F(\theta) = \mathbb{E}[g(\theta, X; v)]$.
and there exists a $v \in \mathbb{U}$ such that $F(\theta) = \mathbb{E}[g(\theta, X; v)]$.

As a byproduct of proving the global convergence of **SAM2**, we obtain a uniform strong law of large numbers for double arrays, which is novel, to the best of our knowledge. This strong law provides a useful result with verifiable regularity conditions for functions indexed by compact Euclidean sets, which serves as an alternative to the more generic theory of Ziegler (2001). The result can also be viewed within the context of Andrews (1992), who posit a double array extension of some uniform strong laws, but do not make explicit the required conditions.

Beyond our general convergence results, we also provide detailed investigations of an online mirror descent algorithm and an online proximal-gradient algorithm as examples of special cases of the **SAM2** framework. These examples are complemented by numerical results regarding the application of **SAM2** to quantile regression, and in particular, the least absolute deviation (LAD) problem. We benchmark our **SAM2** algorithm against Stochastic Subgradient descent (SSG), where we demonstrate that the performance of SSG is sensitive to user calibrated step size schedule, which is avoided by **SAM2**.

The remainder of the manuscript is organized as follows. The **SAM2** algorithm framework is described in Section 2. Asymptotic convergence analysis of **SAM2** is provided in Section 3. Technical descriptions of online mirror descent and proximal-gradient algorithms as examples are provided in Section 4. Numerical illustrations of **SAM2** to the LAD problem appears in Section 5. And finally, proofs and technical results are reported in Section 6.

**Notations.** Throughout the paper, vectors are column-vectors. $\langle a, b \rangle$ denotes the dot product in $\mathbb{R}^d$ and $\| \cdot \|$ is the associated norm. $\partial f$ is the subdifferential of a function $f$ and $\nabla f$ denotes the gradient of a differentiable function $f$. The theorems, lemmas, corollaries, propositions and examples share the same counter; while the algorithms have a separate counter.

## 2 The SAM2 algorithm

We consider the following optimization problem on a compact subset $\mathbb{T}$ of $\mathbb{R}^d$

$$ \operatorname{argmin}_{\theta \in \mathbb{T}} F(\theta), \qquad (3) $$

when the Majorization–Minimization (MM) framework applies:

A1. There exists a Borel set $\mathbb{U}$ of $\mathbb{R}^{d'}$ and a measurable function $G: \mathbb{T} \times \mathbb{U} \to \mathbb{R}$ such that

(a) For all $v \in \mathbb{U}$, $F(\cdot) \le G(\cdot; v)$ on $\mathbb{T}$.

(b) For all $\theta \in \mathbb{T}$, there exists $v \in \mathbb{U}$ such that $F(\theta) = G(\theta; v)$.

For all $\theta \in \mathbb{T}$, we define

$$ \mathbb{U}[\theta] := \{v \in \mathbb{U} : F(\theta) = G(\theta; v)\}; $$
As a byproduct of proving the global convergence of SAM2, we obtain a uniform strong law of large numbers for double arrays, also known as triangular arrays, which is novel, to the best of our knowledge. This uniform double array strong law allows for simultaneous control of stochastic variation of sample mean functional evaluations with respect to both the parameter values of the function (cf. Shapiro et al., 2021, Thm. 9.60) and the different sequences of data that the function is computed on (cf. Hu et al., 1989), providing a useful theorem with verifiable regularity conditions for functions indexed by compact Euclidean sets, and serves as an alternative to the more technical theory of Ziegler (2001). The result can also be viewed within the context of Andrews (1992), who posit a double array extension of some uniform strong laws, but did not make explicit the required conditions.

Beyond our general convergence results, we also provide detailed investigations of a version of the online mirror descent algorithm and a version of the online proximal-gradient algorithm as examples of special cases of the SAM2 framework. These examples are complemented by numerical results regarding the application of SAM2 to quantile regression, and in particular, some least absolute deviation (LAD) problems. We benchmark our SAM2 algorithm against Stochastic Subgradient descent (SSG), where we demonstrate that the performance of SSG is sensitive to user calibrated step size schedule, which is avoided by SAM2. SAM2 is further compared with a Stochastic Mirror Descent (SMD) algorithm, which is less sensitive to step sizes than SSG and performs well in the convex case, while being outperformed by SAM2 in the non-convex setting.

The remainder of the manuscript is organized as follows. The SAM2 algorithm framework is described in Section 2. Asymptotic convergence analysis of SAM2 is provided in Section 3. Technical descriptions of online mirror descent and proximal-gradient algorithms as examples are provided in Section 4. Numerical illustrations of SAM2 to the LAD problem appears in Section 5. And finally, proofs and technical results are reported in Section 6.

**Notations.** For a positive integer $N$, $[N]$ denotes $\{1, \dots, N\}$. Throughout the paper, vectors are column-vectors. $\langle a, b \rangle$ denotes the dot product in $\mathbb{R}^d$ and $\|\cdot\|$ is the associated norm. $\partial f$ is the subdifferential of a function $f$ and $\nabla f$ denotes the gradient of a differentiable function $f$. The theorems, lemmas, corollaries, propositions and examples share the same counter; while the algorithms have a separate counter.

## 2 The SAM2 algorithm

We consider the following optimization problem on a compact subset $\mathcal{T}$ of $\mathbb{R}^d$

$$ \operatorname{argmin}_{\theta \in \mathcal{T}} F(\theta), \qquad (3) $$

when the Majorization–Minimization (MM) framework applies:
under A1-b, this set is not empty. We consider the cases when G is an intractable expectation, but stochastic oracles exist.

A2. (a) Let $(\Omega, \mathcal{A}, \mathbb{P})$ be a probability space and $(\mathbb{X}, \mathcal{X})$ be a measurable space. There exists a measurable function $g: \mathbb{T} \times \mathbb{X} \times \mathbb{U} \to \mathbb{R}$ such that $G(\theta, v) := \mathbb{E}[g(\theta, X; v)];$ and for all $\theta \in \mathbb{T}$ and $v \in \mathbb{U}$, $\mathbb{E}[|g(\theta, X; v)|] < \infty$.

(b) A stream of random variables $\{X^{t,i}, t \ge 1, 1 \le i \le N_t\}$, defined on $(\Omega, \mathcal{A}, \mathbb{P})$ and i.i.d. with the same distribution as $X$, is available.

Therefore, a stochastic oracle can be defined as follows, computed from $N$ examples:

$$G(\theta; v) \approx \frac{1}{N} \sum_{i=1}^{N} g(\theta, X^i; v).$$

**Example 1.** Consider the minimization of a function F on $\mathbb{T}$, a compact subset of $\mathbb{R}^d$. Assume that $F(\theta) := \mathbb{E}[f(\langle\theta, X\rangle)]$ where $f: \mathbb{R} \to \mathbb{R}$ is a convex function such that $\mathbb{E}[|f(\langle\theta, X\rangle)|] < \infty$ for all $\theta \in \tilde{\mathbb{T}} \supset \mathbb{T}$.

By writing $\langle\theta, x\rangle = \langle v, x\rangle + \sum_{j=1}^d (\theta_j - v_j)x_j$, it holds since $f$ is convex (see (Lange, 2013, Section 8.3))

$$f(\langle\theta, x\rangle) \le g(\theta, x; v) := \frac{1}{d} \sum_{j=1}^{d} f(\langle v, x\rangle + d(\theta_j - v_j)x_j), \quad v \in \mathbb{T}.$$

We assume that $\tilde{\mathbb{T}}$ is large enough so that $\mathbb{E}[|f(\langle v, X\rangle + d(\theta_j - v_j)X_j)|] < \infty$ for all $\theta, v \in \mathbb{T}$.

Set $\mathbb{U} := \mathbb{T}$ and $G(\theta; v) := \mathbb{E}[g(\theta, X; v)]$. The condition on $\tilde{\mathbb{T}}$ implies that A2 is verified. For any $v \in \mathbb{T}$, we have $g(\theta, x; v) - f(\langle\theta, x\rangle) \ge 0$, thus implying that A1-a holds. Finally $g(\theta, x; \theta) = f(\langle\theta, x\rangle)$ thus showing that A1-b holds and $\mathbb{U}[\theta] \supset \{\theta\}$. If $f$ is strictly convex, then $\mathbb{U}[\theta] = \{\theta\}$. This majorizer is particularly useful for converting a $d$-dimensional optimization problem to iterations that solve $d$ one-dimensional problems instead. Often, the one-dimensional problems have either computationally efficient or closed form solutions, when the $d$-dimensional problem lacks either.

**Example 2.** Consider the maximization of a positive function F on a compact subset $\mathbb{T}$ of $\mathbb{R}^d$, when $F(\theta) := \int f(\theta, x)\nu(dx)$ where $\nu$ is a $\sigma$-finite positive measure on $(\mathbb{X}, \mathcal{X})$ and $f \ge 0$. Equivalently, we minimize $-\log F$. We assume that $F < \infty$ on $\mathbb{T}$.

The MM framework holds for the minimization of $-\log F$ since the Jensen's inequality implies that for any $\tau \in \mathbb{T}$

$$-\log F(\theta) \le G(\theta; \tau) := -\log F(\tau) - \int \log \left( \frac{f(\theta, X)}{f(\tau, X)} \right) \frac{f(\tau, X)}{F(\tau)} \nu(dx).$$

Here, $\mathbb{U} := \mathbb{T}$ and A1-a holds. We have $G(\theta; \theta) = -\log F(\theta)$ for all $\theta \in \mathbb{T}$, so that A1-b holds and $\mathbb{U}[\theta] \supset \{\theta\}$. Set $g(\theta, x; \tau) := -\log F(\tau) - \log (f(\theta, x)/f(\tau, x)) f(\tau, x)/F(\tau)$. Then A2 holds under the integrability condition of $x \mapsto \log (f(\theta, x)/f(\tau, x)) f(\tau, x)$ with respect to $\nu$, for all $\theta, \tau \in \mathbb{T}$.
A1. There exists a Borel set $\mathbb{U}$ of $\mathbb{R}^d$ and a measurable function $G: \mathbb{T} \times \mathbb{U} \to \mathbb{R}$ such that

(a) For all $v \in \mathbb{U}$, $F(v) \le G(v)$ on $\mathbb{T}$.

(b) For all $\theta \in \mathbb{T}$, there exists $v \in \mathbb{U}$ such that $F(\theta) = G(\theta; v)$.

For all $\theta \in \mathbb{T}$, we define

$$ \mathbb{U}[\theta] := \{v \in \mathbb{U} : F(\theta) = G(\theta; v)\}; $$

under A1-b, this set is not empty. We consider the cases when $G$ is an intractable expectation, but stochastic oracles exist.

A2. (a) Let $(\Omega, \mathcal{A}, \mathbb{P})$ be a probability space and $(\mathbb{X}, \mathcal{X})$ be a measurable space. There exists a measurable function $g: \mathbb{T} \times \mathbb{X} \times \mathbb{U} \to \mathbb{R}$ such that $G(\theta, v) := \mathbb{E}[g(\theta, X; v)];$ and for all $\theta \in \mathbb{T}$ and $v \in \mathbb{U}$, $\mathbb{E}[|g(\theta, X; v)|] < \infty$.

(b) A stream of random variables $\{X^{t,i}, t \ge 1, 1 \le i \le N_t\}$, defined on $(\Omega, \mathcal{A}, \mathbb{P})$ and i.i.d. with the same distribution as $X$, is available.

Therefore, a stochastic oracle can be defined as follows, computed from $N$ examples:

$$ G(\theta; v) \approx \frac{1}{N} \sum_{i=1}^{N} g(\theta, X^i; v). $$

**Example 1.** Consider the minimization of a function $F$ on $\mathbb{T}$, a compact subset of $\mathbb{R}^d$. Assume that $F(\theta) := \mathbb{E}[f(\langle\theta, X\rangle)]$ where $f: \mathbb{R} \to \mathbb{R}$ is a convex function such that $\mathbb{E}[|f(\langle\theta, X\rangle)|] < \infty$ for all $\theta \in \tilde{\mathbb{T}} \supset \mathbb{T}$.

By writing $\langle\theta, x\rangle = \langle v, x\rangle + \sum_{j=1}^d (\theta_j - v_j)x_j$, it holds since $f$ is convex (see (Lange, 2013, Section 8.3))

$$ f(\langle\theta, x\rangle) \le g(\theta, x; v) := \frac{1}{d} \sum_{j=1}^{d} f(\langle v, x\rangle + d(\theta_j - v_j)x_j), \quad v \in \mathbb{T}. $$

We assume that $\tilde{\mathbb{T}}$ is large enough so that $\mathbb{E}[|f(\langle v, X\rangle + d(\theta_j - v_j)X_j)|] < \infty$ for all $\theta, v \in \mathbb{T}$.

Set $\mathbb{U} := \mathbb{T}$ and $G(\theta; v) := \mathbb{E}[g(\theta, X; v)]$. The condition on $\tilde{\mathbb{T}}$ implies that A2 is verified. For any $v \in \mathbb{T}$, we have $g(\theta, x; v) - f(\langle\theta, x\rangle) \ge 0$, thus implying that A1-a holds. Finally $g(\theta, x; \theta) = f(\langle\theta, x\rangle)$ thus showing that A1-b holds and $\mathbb{U}[\theta] \supset \{\theta\}$. If $f$ is strictly convex, then $\mathbb{U}[\theta] = \{\theta\}$. This majorizer is particularly useful for converting a $d$-dimensional optimization problem to iterations that solve $d$ one-dimensional problems instead. Often, the one-dimensional problems have either computationally efficient or closed form solutions, when the $d$-dimensional problem lacks either.

**Example 2.** Consider the maximization of a positive function $F$ on a compact subset $\mathbb{T}$ of $\mathbb{R}^d$, when $F(\theta) := \int f(\theta, x)\nu(dx)$ where $\nu$ is a $\sigma$-finite positive measure on $(\mathbb{X}, \mathcal{X})$ and $f \ge 0$. Equivalently, we minimize $-\log F$. We assume that $F < \infty$ on $\mathbb{T}$.
We note that this majorization scheme is exactly that which is used to derive the famous expectation-maximization algorithms of Dempster et al. (1977) (see, also McLachlan & Krishnan, 2008), as per (Lange, 2013, Sec. 9.8) and (Razaviyayn et al., 2013, Sec. 8.5). In that case, f is the complete data likelihood.

**Example 3.** Consider the minimization of a function $F: \mathbb{R}^d \to (-\infty, +\infty]$, on a compact convex subset $\mathbb{T}$ of $\mathbb{R}^d$, when $F(\theta) := \mathbb{E}[f(\theta, X)];$ we assume that $\mathbb{E}[|f(\theta, X)|] < \infty$ for all $\theta \in \mathbb{T}$ so that $\mathbb{T} \subset \text{dom}(F)$. Let $b: \mathbb{R}^d \to (-\infty, +\infty]$ be a lower semicontinuous convex function, such that $\mathbb{T} \subset \text{dom}(b)$, continuously differentiable and strictly convex on an open neighborhood $\mathbb{S}$ of $\mathbb{T}$. Define the Bregman divergence associated to $b$

$$D_b(\theta; \tau) := b(\theta) - b(\tau) - \langle \nabla b(\tau), \theta - \tau \rangle, \quad \theta \in \mathbb{R}^d, \tau \in \mathbb{T}. \qquad (4)$$

It holds for any $\theta \in \mathbb{R}^d$, $x \in \mathbb{X}$, $\tau \in \mathbb{T}$ and $\gamma > 0$

$$f(\theta, x) \le g(\theta, x; (\tau, \gamma)) := f(\theta, x) + \frac{1}{\gamma} D_b(\theta; \tau) + \iota_{\mathbb{T}}(\theta),$$

where $\iota_{\mathbb{T}}: \mathbb{R}^d \to (-\infty, +\infty]$ is the characteristic function of $\mathbb{T}$, defined as $\iota_{\mathbb{T}}(\theta) = 0$ if $\theta \in \mathbb{T}$ and $\iota_{\mathbb{T}}(\theta) = +\infty$, otherwise. Set $\mathbb{U} := \mathbb{T} \times \mathbb{R}_{>0}$, and $G(\theta; \upsilon) := \mathbb{E}[g(\theta, X; \upsilon)]$. Then the condition A2 is verified. For any $\upsilon := (\tau, \gamma) \in \mathbb{U}$, we have $g(\cdot, x; \upsilon) - f(\cdot, x) \ge 0$ on $\mathbb{T}$ since $b$ is convex, so that A 1-a holds. In addition, for all $\theta \in \mathbb{T}$ and $\gamma > 0$, $g(\theta, x; (\theta, \gamma)) = f(\theta, x)$ and $g(\theta, x; (\tau, \gamma)) \ne f(\theta, x)$ when $\tau \ne \theta$ since $b$ is strictly convex; hence, A 1-b holds and $\mathbb{U}[\theta] := \{(\theta, \gamma), \gamma > 0\}$.

**Non uniqueness of the majorizing functions.** If $(\theta, \upsilon) \mapsto G(\theta; \upsilon)$ and $(\theta, x, \upsilon) \mapsto g(\theta, x; \upsilon)$ satisfy A1 and A2, then the functions $(\theta, \upsilon, \tau) \mapsto G(\theta; \upsilon) + \varphi(D_b(\varsigma(\theta); \varsigma(\tau)))$ and $(\theta, x, \upsilon, \tau) \mapsto \tilde{g}(\theta, x; \upsilon, \tau) := g(\theta, x; \upsilon) + \varphi(D_b(\varsigma(\theta); \varsigma(\tau)))$ also satisfy A1 and A2 as soon as

(i) $\varphi: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$ and $\varsigma: \mathbb{T} \to \mathbb{T}$ are measurable functions, and $\varphi(0) = 0$.

(ii) $D_b(\cdot; \cdot)$ is the Bregman divergence (see (4)) associated to a strictly convex and continuously differentiable function $b: \mathbb{S} \to \mathbb{R}$ defined on a neighborhood $\mathbb{S}$ of $\mathbb{T}$.

Note indeed that for all $\theta \in \mathbb{T}$ and $x \in \mathbb{X}$, it holds

$$\begin{cases} \forall (\upsilon, \tau) \in \mathbb{U} \times \mathbb{T}: g(\theta, x; \upsilon) \le \tilde{g}(\theta, x; \upsilon, \tau) \\ \forall \upsilon \in \mathbb{U}[\theta]: g(\theta, x; \upsilon) = \tilde{g}(\theta, x; \upsilon, \theta). \end{cases}$$

The introduction of such a Bregman term in the majorization mechanism is the essence of many proximal-type algorithms such as the **Proximal point algorithm** (see e.g. (Polyak, 2021, Section 3.5)) and the **Mirror Descent algorithm** (see Section 4.1) among examples. The motivation is essentially to make the minimization of the majorizing function easier.
The MM framework holds for the minimization of $-\log F$ since the Jensen's inequality implies that for any $\tau \in \mathbb{T}$

$$-\log F(\theta) \leq G(\theta; \tau) := -\log F(\tau) - \int \log \left( \frac{f(\theta, X)}{f(\tau, X)} \right) \frac{f(\tau, X)}{F(\tau)} \nu(dx).$$

Here, $\mathbb{U} := \mathbb{T}$ and $\mathbb{A}$ 1-a holds. We have $G(\theta; \theta) = -\log F(\theta)$ for all $\theta \in \mathbb{T}$, so that $\mathbb{A}$ 1-b holds and $\mathbb{U}[\theta] \supseteq \{\theta\}$. Set $g(\theta, x; \tau) := -\log F(\tau) - \log (f(\theta, x)/f(\tau, x)) f(\tau, x)/F(\tau)$. Then $\mathbb{A}$ 2 holds under the integrability condition of $x \mapsto \log (f(\theta, x)/f(\tau, x)) f(\tau, x)$ with respect to $\nu$, for all $\theta, \tau \in \mathbb{T}$.

We note that this majorization scheme is exactly that which is used to derive the famous Expectation-Maximization algorithm of Dempster et al. (1977) (see, also McLachlan & Krishnan, 2008), as per Lange (2013, Sec. 9.8) and Razaviyayn et al. (2013, Sec. 8.5). In that case, $f$ is the so-called complete data likelihood.

**Example 3.** Consider the minimization of a function $F: \mathbb{R}^d \to (-\infty, +\infty]$ on a compact convex subset $\mathbb{T}$ of $\mathbb{R}^d$, when $F(\theta) := \mathbb{E}[f(\theta, X)];$ we assume that $\mathbb{E}[|f(\theta, X)|] < \infty$ for all $\theta \in \mathbb{T}$ so that $\mathbb{T} \subseteq \text{dom}(F)$. Let $\mathfrak{b}: \mathbb{R}^d \to (-\infty, +\infty]$ be a lower semicontinuous convex function, such that $\mathbb{T} \subset \text{dom}(\mathfrak{b})$, continuously differentiable and strictly convex on an open neighborhood $\mathbb{S}$ of $\mathbb{T}$. Define the Bregman divergence associated to $\mathfrak{b}$

$$D_{\mathfrak{b}}(\theta; \tau) := \mathfrak{b}(\theta) - \mathfrak{b}(\tau) - \langle \nabla \mathfrak{b}(\tau), \theta - \tau \rangle, \quad \theta \in \mathbb{R}^d, \tau \in \mathbb{T}. \qquad (4)$$

It holds for any $\theta \in \mathbb{R}^d$, $x \in \mathbb{X}$, $\tau \in \mathbb{T}$ and $\gamma > 0$

$$f(\theta, x) \leq g(\theta, x; (\tau, \gamma)) := f(\theta, x) + \frac{1}{\gamma} D_{\mathfrak{b}}(\theta; \tau) + \iota_{\mathbb{T}}(\theta),$$

where $\iota_{\mathbb{T}}: \mathbb{R}^d \to (-\infty, +\infty]$ is the characteristic function of $\mathbb{T}$, defined as $\iota_{\mathbb{T}}(\theta) = 0$ if $\theta \in \mathbb{T}$ and $\iota_{\mathbb{T}}(\theta) = +\infty$, otherwise. Set $\mathbb{U} := \mathbb{T} \times \mathbb{R}_{>0}$, and $G(\theta; \upsilon) := \mathbb{E}[g(\theta, X; \upsilon)]$. Then the condition A 2 is verified. For any $\upsilon := (\tau, \gamma) \in \mathbb{U}$, we have $g(\cdot, x; \upsilon) - f(\cdot, x) \geq 0$ on $\mathbb{T}$ since $\mathfrak{b}$ is convex, so that $\mathbb{A}$ 1-a holds. In addition, for all $\theta \in \mathbb{T}$ and $\gamma > 0$, $g(\theta, x; (\theta, \gamma)) = f(\theta, x)$ and $g(\theta, x; (\tau, \gamma)) \neq f(\theta, x)$ when $\tau \neq \theta$ since $\mathfrak{b}$ is strictly convex; hence, $\mathbb{A}$ 1-b holds and $\mathbb{U}[\theta] := \{(\theta, \gamma), \gamma > 0\}$.

**Non uniqueness of the majorizing functions.** If $(\theta, \upsilon) \mapsto G(\theta; \upsilon)$ and $(\theta, x, \upsilon) \mapsto g(\theta, x; \upsilon)$ satisfy A 1 and A 2, then the functions $(\theta, \upsilon, \tau) \mapsto G(\theta; \upsilon) + \varphi(D_{\mathfrak{b}}(\varsigma(\theta); \varsigma(\tau)))$ and $(\theta, x, \upsilon, \tau) \mapsto \tilde{g}(\theta, x; \upsilon, \tau) := g(\theta, x; \upsilon) + \varphi(D_{\mathfrak{b}}(\varsigma(\theta); \varsigma(\tau)))$ also satisfy A 1 and A 2 whenever

(i) $\varphi: \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}$ and $\varsigma: \mathbb{T} \to \mathbb{T}$ are measurable functions, and $\varphi(0) = 0$.

(ii) $D_{\mathfrak{b}}(\cdot; \cdot)$ is the Bregman divergence (see (4)) associated to a strictly convex and continuously differentiable function $\mathfrak{b}: \mathbb{S} \to \mathbb{R}$ defined on a neighborhood $\mathbb{S}$ of $\mathbb{T}$.
# Population MM and SAM2 algorithms

When $\operatorname{argmin}_{\mathbb{T}} G(\cdot; v)$ is not empty for all $v$ in the set $\mathbb{U}$ (see e.g. Proposition 11-i for sufficient conditions), then the Population MM algorithm given by Algorithm 1 defines a $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 1\}$ satisfying a descent property: $F(\theta^{t+1}) \le F(\theta^t)$ for all $t \ge 0$. We have indeed by A1 and

**Algorithm 1** Population MM

**Require:** an initial value $\theta^0 \in \mathbb{T}$

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 0\}$.

1: **for** $t = 0, \dots, \mathbf{do}$

2: Choose $v^t \in \mathbb{U}[\theta^t]$.

3: Compute $\theta^{t+1} \in \operatorname{Argmin}_{\theta \in \mathbb{T}} G(\theta; v^t)$.

4: **end for**

by definition of $\theta^{t+1}$: $F(\theta^{t+1}) \le G(\theta^{t+1}; v^t) \le G(\theta^t; v^t) = F(\theta^t)$.

When $G$ does not have a closed form expression but stochastic oracles $g(\theta, X; v)$ for $G(\theta; v)$ are available, a stochastic Population MM can be designed. We consider in this paper the Sequential Sample Average Majorization-Minimization algorithm (SAM2) given by Algorithm 2. Note that under the compactness property of $\mathbb{T}$, the set of minimizers in Equation (5) of Algorithm 2 is not empty as soon as for any $x \in \mathbb{X}$ and $v \in \mathbb{T}$, the function $\theta \mapsto g(\theta, x; v)$ is lower semicontinuous on $\mathbb{T}$.

**Algorithm 2** The SAM2 algorithm

**Require:** a sequence $\{N_t, t \ge 1\}$ of positive integers, an initial value $\theta^0 \in \mathbb{T}$

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 0\}$.

1: **for** $t = 0, \dots, \mathbf{do}$

2: Sample a minibatch $\{X^{t+1,i}, i = 1, \dots, N_{t+1}\}$ of size $N_{t+1}$.

3: Choose $v^t \in \mathbb{U}[\theta^t]$.

4: Compute

$$ \theta^{t+1} \in \operatorname{Argmin}_{\theta \in \mathbb{T}} \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta, X^{t+1,i}; v^t). \qquad (5) $$

5: **end for**

*Remark.* In the preceding discussion, we note that $G(\cdot; v)$ for $v \in \mathbb{U}[\theta]$, is a local approximation of $F$ in a neighborhood of $\theta \in \mathbb{T}$ (see A1), whereas the relationship between $G$ and $g$ from A2 is used to estimate $G$ via Monte Carlo simulation. Putting this together, at each iteration of Algorithm 2,

$$ \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta, X^{t+1,i}; v^t) $$

is used as an estimate of the approximation $G(\cdot; v^t)$ of $F$. The forthcoming analysis of Algorithm 2 focuses around the control of the numerical and stochastic errors regarding this joint approximation
Note indeed that for all $\theta \in \mathbb{T}$ and $x \in \mathbb{X}$, it holds that

$$ \left\{ \begin{array}{l} \forall (v, \tau) \in \mathbb{U} \times \mathbb{T} : g(\theta, x; v) \le \tilde{g}(\theta, x; v, \tau), \\ \forall v \in \mathbb{U}[\theta] : g(\theta, x; v) = \tilde{g}(\theta, x; v, \theta). \end{array} \right. $$

## Population MM and SAM2 algorithms

When $\text{argmin}_{\mathbb{T}}G(\cdot; v)$ is not empty for all $v$ in the set $\mathbb{U}$ (see e.g. Proposition 12-i for sufficient conditions), then the **Population MM** algorithm given by Algorithm 1 defines a $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 1\}$ satisfying a descent property: $F(\theta^{t+1}) \le F(\theta^t)$ for all $t \ge 0$. We have indeed by A1 and

**Algorithm 1** Population MM algorithm

**Require:** an initial value $\theta^0 \in \mathbb{T}$.

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 0\}$.

1: **for** $t \ge 0$ **do**

2: Choose $v^t \in \mathbb{U}[\theta^t]$.

3: Compute $\theta^{t+1} \in \text{Argmin}_{\theta \in \mathbb{T}} G(\theta; v^t)$.

4: **end for**

by definition of $\theta^{t+1}$: $F(\theta^{t+1}) \le G(\theta^{t+1}; v^t) \le G(\theta^t; v^t) = F(\theta^t)$.

When $G$ does not have a closed form expression but stochastic oracles $g(\theta, X; v)$ for $G(\theta; v)$ are available, a stochastic **Population MM** can be designed. We consider in this paper the *Sequential Sample Average Majorization-Minimization* algorithm (SAM2) given by Algorithm 2. Note that under the compactness property of $\mathbb{T}$, the set of minimizers in Equation (5) of Algorithm 2 is not empty when, for any $x \in \mathbb{X}$ and $v \in \mathbb{T}$, the function $\theta \mapsto g(\theta, x; v)$ is lower semi-continuous on $\mathbb{T}$.

**Algorithm 2** The SAM2 algorithm

**Require:** a sequence $\{N_t, t \ge 1\}$ of positive integers, an initial value $\theta^0 \in \mathbb{T}$.

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 0\}$.

1: **for** $t \ge 0$ **do**

2: Sample a minibatch $\{X^{t+1,i}, 1 \le i \le N_{t+1}\}$ of size $N_{t+1}$.

3: Choose $v^t \in \mathbb{U}[\theta^t]$.

4: Compute

$$ \theta^{t+1} \in \text{Argmin}_{\theta \in \mathbb{T}} \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta, X^{t+1,i}; v^t). \qquad (5) $$

5: **end for**

*Remark.* In the preceding discussion, we note that $G(\cdot; v)$ for $v \in \mathbb{U}[\theta]$, is a local approximation of $F$ in a neighborhood of $\theta \in \mathbb{T}$ (see A1), whereas the relationship between $G$ and $g$ from A2 is used to estimate $G$ via Monte Carlo simulation. Putting this together, at each iteration of Algorithm 2,

$$ \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta, X^{t+1,i}; v^t) $$
and estimation process. Such requirement to control both types of errors is not unique to the SAM2 algorithm, as similar analysis is required for gradient-based algorithms when the objective function F is an intractable expectation: gradient-based algorithms operate directly on the sample gradient of F. Unlike such algorithms, however, computation of (5) does not always require information regarding the gradient of F, which is an advantage when such objects are difficult or costly to compute.

# 3 Asymptotic convergence analysis of SAM2

This section is devoted to the almost-sure convergence of SAM2 when the number of iterations tends to infinity, under the additional assumption that each majorizing function $G(\cdot; v)$ possesses a unique minimizer denoted by $T(v)$ (see A3-Item b). Regularity conditions on $g$ and the objective function F are also assumed (see A3-Item a,Item c) which imply continuity properties of $G$ and $T$ (see Proposition 11 and Proposition 12). The proof consists in comparing one iteration of SAM2 to one iteration of Population MM and proving that under conditions on the stochastic perturbations, both algorithms have the same limiting set.

We start with discussing the fixed points of Population MM and how this set of points is related to the minimizers of F (see Theorem 4). Then, a general result on the convergence of perturbed iterative scheme is derived (see Theorem 5), whose main assumptions are the existence of a Lyapunov function for the iterated map, and the perturbation is vanishing in a sense related to this Lyapunov function; this convergence theorem covers general iterative scheme and is not restricted to the case of the MM one. In the MM framework, the objective function F plays the role of the Lyapunov function and Proposition 7 shows that the vanishing perturbation condition is satisfied as soon as the perturbation is vanishing in a sense related to the majorizing functions G. Combining Theorem 5, Proposition 7 and a condition on the batch size (see A4) yields the convergence result for SAM2 (see Theorem 9).

The convergence analysis is derived under the following regularity conditions on F and g.

A3. a) For any $x \in \mathcal{X}$, the function $(\theta, v) \mapsto g(\theta, x; v)$ is continuous on $\mathbb{T} \times \mathcal{U}$. In addition, there exists $p > 1$ such that $\mathbb{E}[\sup_{(\theta,v) \in \mathbb{T} \times \mathcal{U}} |g(\theta, X; v)|^p] < \infty$.

b) For all $v \in \mathcal{U}$, the set $\text{Argmin}_{\mathbb{T}} G(\cdot; v)$ has exactly one element, denoted by $T(v)$.

c) The function $\theta \mapsto F(\theta)$ is continuous on $\mathbb{T}$.

Under A3-a, a minimizer $\theta^{t+1}$ in (5) exists: SAM2 is well defined. A3-a also implies that $\theta \mapsto G(\theta; v)$ is continuous for all $v \in \mathcal{U}$ (see Proposition 11); since $\mathbb{T}$ is compact, the set $\text{Argmin}_{\mathbb{T}} G(\cdot; v)$ has at least one element. Finally, A3-a and A3-b imply that $T$ is a continuous point-to-point map from $\mathcal{U}$ to $\mathbb{T}$ (see Proposition 12).
is used as an estimate of the approximation $G(\cdot; v^t)$ of F. The forthcoming analysis of Algorithm 2 focuses around the control of the numerical and stochastic errors regarding this joint approximation and estimation process. Such requirement to control both types of errors is not unique to the **SAM2** algorithm, as similar analysis is required for gradient-based algorithms when the objective function F is an intractable expectation: gradient-based algorithms operate directly on the sample gradient of F. Unlike such algorithms, however, computation of (5) does not always require information regarding the gradient of F, which is an advantage when such objects are difficult or costly to compute.

# 3 Asymptotic convergence analysis of SAM2

This section is devoted to the almost-sure convergence of **SAM2** when the number of iterations tends to infinity, under the additional assumption that each majorizing function $G(\cdot; v)$ possesses a unique minimizer denoted by $T(v)$ (see A3-Item b). Regularity conditions on $g$ and the objective function F are also assumed (see A3-Item a, Item c) which imply continuity properties of $G$ and $T$ (see Proposition 12 and Proposition 13). The proof consists in comparing one iteration of **SAM2** to one iteration of **Population MM** and proving that under conditions on the stochastic perturbations, both algorithms have the same limiting set.

We start with discussing the fixed points of **Population MM** and how this set of points is related to the minimizers of F (see Theorem 4). Then, a general result on the convergence of perturbed iterative scheme is derived (see Theorem 5), whose main assumptions are the existence of a *Lyapunov function* for the iterated map, and the perturbation is vanishing in a sense related to this Lyapunov function. This convergence theorem covers general iterative scheme and is not restricted to the case of the **MM** algorithm. In the **MM** framework, the objective function F plays the role of the Lyapunov function and Proposition 7 shows that the vanishing perturbation condition is satisfied when the perturbation is vanishing in a sense related to the majorizing functions $G$. Combining Theorem 5, Proposition 7, and a condition on the batch size (see A 4) yields the convergence result for **SAM2** (see Theorem 9).

The convergence analysis is derived under the following regularity conditions on F and g.

A3. a) For any $x \in \mathcal{X}$, the function $(\theta, v) \mapsto g(\theta, x; v)$ is continuous on $\mathbb{T} \times \mathcal{U}$. In addition, there exists $p > 1$ such that $\mathbb{E}[\sup_{(\theta,v) \in \mathbb{T} \times \mathcal{U}} |g(\theta, X; v)|^p] < \infty$.

b) For all $v \in \mathcal{U}$, the set $\text{Argmin}_{\mathbb{T}} G(\cdot; v)$ has exactly one element, denoted by $T(v)$.

c) The function $\theta \mapsto F(\theta)$ is continuous on $\mathbb{T}$.

Under A 3-a, a minimizer $\theta^{t+1}$ in (5) exists: i.e., **SAM2** is well defined. A 3-a also implies that $\theta \mapsto G(\theta; v)$ is continuous for all $v \in \mathcal{U}$ (see Proposition 12); since $\mathbb{T}$ is compact, the set
**Example 3** (to follow). Assume in addition that (i) $\mathbf{F}$ is a lower semicontinuous convex function on $\mathbb{R}^d$; (ii) for all $x \in \mathbb{X}$, $\theta \mapsto \mathbf{f}(\theta, x)$ is continuous on $\mathbb{T}$ and there exists $p > 1$ such that $\mathbb{E}[\sup_{\theta \in \mathbb{T}} |\mathbf{f}(\theta, X)|^p] < \infty$; (iii) $\gamma$ is restricted to a compact ball in $\mathbb{R}_{>0}$: $\gamma \in [g_-, g+]$.

The dominated convergence theorem implies that $\mathbf{F}$ is continuous on $\mathbb{T}$ so that A 3-c holds. Since $\mathbf{b}$ is continuous and $\gamma$ is lower bounded from zero, then A 3-a also holds. In addition, for all $(\tau, \gamma) \in \mathbb{U}$, $\theta \mapsto \mathbf{G}(\theta; (\tau, \gamma))$ is strictly convex on $\mathbb{T}$ since $\mathbf{b}$ is strictly convex and $\mathbf{F}$ is convex; hence, $\mathrm{Argmin}_{\theta \in \mathbb{T}} \mathbf{G}(\theta; (\tau, \gamma))$ is not empty and has exactly one element.

Theorem 4 provides sufficient conditions for $\mathrm{Argmin}_{\mathbb{T}}\mathbf{F}$ not to be empty and included in the set $\mathcal{L}$, defined by

$$ \mathcal{L} := \{ \theta \in \mathbb{T} : \forall v \in \mathbb{U}[\theta], \forall \tau \in \mathrm{Argmin}_{\mathbb{T}} \mathbf{G}(\cdot; v), \mathbf{F}(\tau) = \mathbf{F}(\theta) \}. \quad (6) $$

$\mathcal{L}$ contains points in $\mathbb{T}$ such that, starting from such a point, one iteration of **Population MM** can not decrease the value of the objective function. Theorem 4 shows that it is included in the subset of $\mathbb{T}$ consisting of points such that, starting from such a point, one iteration of **Population MM** does not decrease the value of the majorizing function whatever it is. Such a property follows from A1, which implies that for all $\theta \in \mathbb{T}$, $v \in \mathbb{U}[\theta]$ and $\tau \in \mathrm{Argmin}_{\mathbb{T}} \mathbf{G}(\cdot; v)$, it holds

$$ \mathbf{F}(\tau) \le \min_{\mathbb{T}} \mathbf{G}(\cdot; v) \le \mathbf{G}(\theta; v) = \mathbf{F}(\theta). \quad (7) $$

When $\mathrm{Argmin}_{\mathbb{T}} \mathbf{G}(\cdot; v)$ has exactly one element denoted by $\mathbb{T}(v)$, then Theorem 4 also shows that

$$ \mathcal{L} = \{ \theta \in \mathbb{T} : \forall v \in \mathbb{U}[\theta], \mathbf{F}(\mathbb{T}(v)) = \mathbf{F}(\theta) = \min_{\mathbb{T}} \mathbf{G}(\cdot; v) \} = \{ \theta \in \mathbb{T} : \forall v \in \mathbb{U}[\theta], \mathbb{T}(v) = \theta \}, \quad (8) $$

so that $\mathcal{L}$ is the set of the fixed points of **Population MM** (see Algorithm 1). The proof of Theorem 4 is given in Section 6.2.

**Theorem 4.** i) Assume A 3-a. For all $v \in \mathbb{U}$, the set $\mathrm{Argmin}_{\mathbb{T}} \mathbf{G}(\cdot; v)$ is not empty.

ii) Assume A 1, A 3-a and A 3-c. The set $\mathrm{Argmin}_{\mathbb{T}} \mathbf{F}$ is not empty and $\mathrm{Argmin}_{\mathbb{T}} \mathbf{F} \subset \mathcal{L} \subset \{ \theta \in \mathbb{T} : \theta \in \mathrm{Argmin}_{\mathbb{T}} \mathbf{G}(\cdot; v) \text{ for all } v \in \mathbb{U}[\theta] \}$.

iii) Assume A 1 and A 3. Then the equalities (8) hold.

**Example 3** (to follow). Under the stated assumptions, A 1 and A 3 are verified. Let us identify the set $\mathcal{L}$.

Let $\theta \in \mathbb{T}$ and $v = (\theta, \gamma) \in \mathbb{U}[\theta]$; for $v := (\theta, \gamma)$, since $\theta \mapsto \mathbf{G}(\theta; v)$ is a lower semicontinuous proper convex function then it holds (see e.g. (Bauschke & Combettes, 2011, Theorem 16.3)): $0 \in \partial[\mathbf{G}(\cdot; v)](\mathbb{T}(v))$. By (Bauschke & Combettes, 2011, Corollary 16.48 and Proposition 16.6), we have

$$ \partial[\mathbf{G}(\cdot; v)](\mathbb{T}(v)) = \gamma \partial[\mathbf{F} + \iota_{\mathbb{T}}](\mathbb{T}(v)) + \nabla \mathbf{b}(\mathbb{T}(v)) - \nabla \mathbf{b}(\theta). $$
Argmin<sub>T</sub>G(·; v) has at least one element. Finally, A 3-a and A 3-b imply that T is a continuous point-to-point map from U to T (see Proposition 13).

**Example 3** (to follow). Assume in addition that (i) F is a lower semi-continuous convex function on ℝ<sup>d</sup>; (ii) for all x ∈ X, θ ↦→ f(θ, x) is continuous on T and there exists p > 1 such that E[sup<sub>θ∈T</sub> |f(θ, X)|<sup>p</sup>] < ∞; (iii) γ is restricted to a compact ball in ℝ<sub>>0</sub>: γ ∈ [g<sub>-</sub>, g<sub>+</sub>].

The dominated convergence theorem implies that F is continuous on T so that A 3-c holds. Since **b** is continuous and γ is bounded away from zero, A 3-a also holds. In addition, for all (τ, γ) ∈ U, θ ↦→ G(θ; (τ, γ)) is strictly convex on T since **b** is strictly convex and F is convex; hence, Argmin<sub>θ∈T</sub>G(θ; (τ, γ)) is not empty and has exactly one element.

Theorem 4 provides sufficient conditions for Argmin<sub>T</sub>F not to be empty and included in the set

$$
\mathcal{L} := \{ \theta \in \mathbb{T} : \forall v \in \mathbb{U}[\theta], \forall \tau \in \operatorname{Argmin}_{\mathbb{T}} \mathbf{G}(\cdot; v), \mathbf{F}(\tau) = \mathbf{F}(\theta) \}. \quad (6)
$$

$\mathcal{L}$ contains points in $\mathbb{T}$ such that, starting from such a point, one iteration of **Population MM** can not decrease the value of the objective function. Theorem 4 shows that it is included in the subset of $\mathbb{T}$ consisting of points such that, starting from such a point, one iteration of **Population MM** does not decrease the value of the majorizing function whatever it is. Such a property follows from A1, which implies that for all $\theta \in \mathbb{T}$, $v \in \mathbb{U}[\theta]$ and $\tau \in \operatorname{Argmin}_{\mathbb{T}} \mathbf{G}(\cdot; v)$, it holds that

$$
\mathbf{F}(\tau) \le \min_{\mathbb{T}} \mathbf{G}(\cdot; v) \le \mathbf{G}(\theta; v) = \mathbf{F}(\theta). \quad (7)
$$

When Argmin<sub>T</sub>G(·; v) has exactly one element denoted by T(v), then Theorem 4 also shows that

$$
\mathcal{L} = \{ \theta \in \mathbb{T} : \forall v \in \mathbb{U}[\theta], \mathbf{F}(\mathbb{T}(v)) = \mathbf{F}(\theta) = \min_{\mathbb{T}} \mathbf{G}(\cdot; v) \} = \{ \theta \in \mathbb{T} : \forall v \in \mathbb{U}[\theta], \mathbb{T}(v) = \theta \}, \quad (8)
$$

so that $\mathcal{L}$ is the set of the fixed points of **Population MM** (see Algorithm 1). The proof of Theorem 4 is given in Section 6.2.

**Theorem 4.** i) Assume A 3-a. For all v ∈ U, the set Argmin<sub>T</sub>G(·; v) is not empty.

ii) Assume A 1, A 3-a and A 3-c. The set Argmin<sub>T</sub>F is not empty and Argmin<sub>T</sub>F ⊂ L ⊂ {θ ∈ T : θ ∈ Argmin<sub>T</sub>G(·; v), for all v ∈ U[θ]}.

iii) Assume A 1 and A 3. Then (8) holds.

**Example 3** (to follow). Under the stated assumptions, A 1 and A 3 are verified. Let us identify the set $\mathcal{L}$.

Let θ ∈ T and v = (θ, γ) ∈ U[θ]: for v := (θ, γ), since θ ↦→ G(θ; v) is a lower semi-continuous proper convex function then it holds (see e.g. Bauschke & Combettes, 2011, Thm. 16.3): 0 ∈
Therefore, the condition $T(v) = \theta$ is equivalent to $0 \in \partial[F + \iota_T](\theta)$: the set $\mathcal{L}$ is the set of the minimizers of F on $T$. Note also that we have $\mathcal{L} = \{\theta \in T : \exists v \in U[\theta], T(v) = \theta\}$.

The proof of the asymptotic convergence consists in observing that **SAM2** is a stochastic perturbation of the **Population MM** iterative scheme $\tau^{t+1} = T(v^t)$ for $v^t \in U[\tau^t]$. While the sequence $\{\tau^t, t \ge 0\}$ possesses a Lyapunov function i.e. $F(\tau^{t+1}) \le F(\tau^t)$ for all $t \ge 0$, this is no more the case for the **SAM2** sequence $\{\theta^t, t \ge 0\}$. Indeed, under A3-b, $\theta^{t+1}$ is the minimizer of an approximation of the majorizing function $G(\cdot; v^t)$ so that in general, $F(\theta^{t+1}) \ne F(T(v^t))$. Nevertheless, the following result shows that as soon as the error $|F(\theta^{t+1}) - F(T(v^t))|$ vanishes when $t \to +\infty$, the sequence $\{\theta^t, t \ge 0\}$ inherits the same limiting behavior as the **Population MM** sequence. The proof of Theorem 5 is given in Section 6.3; it is adapted from (Fort & Moulines, 2003, Proposition 9) which addresses the case $U[\theta] = \{\theta\}$ for all $\theta \in T$.

**Theorem 5.** Let $T$ be a compact subset of $\mathbb{R}^d$, $U$ be a subset of $\mathbb{R}^d'$ and let $\mathcal{L} \subseteq T$ be a set such that $\mathcal{L} \cap T$ is compact. Let $T: U \to T$ and $F: T \to \mathbb{R}$ be a continuous function such that

(H-i) for all $\theta \in T$, there exists $U[\theta] \subseteq U$ such that $F(T(v)) \le F(\theta)$ for all $v \in U[\theta]$.

(H-ii) for any compact subset $K$ in $T \setminus \mathcal{L}$, $\inf_{\theta \in K} \inf_{v \in U[\theta]} (F(\theta) - F(T(v))) > 0$.

Let $\{(\theta^t, v^t), t \ge 0\}$ be a $T \times U$-valued sequence such that $v^t \in U[\theta^t]$ for all $t \ge 0$ and

(H-iii) $\lim_t |F(T(v^t)) - F(\theta^{t+1})| = 0$.

Then the sequence $\{F(\theta^t), t \ge 0\}$ converges to a connected component of $F(\mathcal{L} \cap T)$. If $F(\mathcal{L} \cap T)$ has an empty interior, the sequence $\{F(\theta^t), t \ge 0\}$ converges to $F^*$ and the sequence $\{\theta^t, t \ge 0\}$ converges to the set $T \cap \mathcal{L}_{F^*}$ where $\mathcal{L}_{F^*} := \{\theta \in \mathcal{L} : F(\theta) = F^*\}$.

Let us start with discussing H-i and H-ii on two examples.

**Example 6.** Assume A1, A2 and A3, and $U[\theta] = \{\theta\}$ for all $\theta \in T$. We prove the conditions H-i and H-ii are verified with $\mathcal{L} := \{\theta \in T : F(\theta) - F(T(\theta)) = 0\} = \{\theta \in T : T(\theta) = \theta\}$.

A3-b and A3-c imply that $T$ is a point-to-point map and $F$ is a continuous function. Since $F$ and $T$ are continuous on $T$ (see A3-c and Proposition 12), the set $\mathcal{L}$ is closed which implies that $\mathcal{L} \cap T$ is compact. By A1 (see also (7)), $F(T(\theta)) \le F(\theta)$ for all $\theta \in T$, so that H-i holds. By definition of $\mathcal{L}$, $F(\theta) - F(T(\theta)) > 0$ for all $\theta \in T \setminus \mathcal{L}$; since $\theta \mapsto F(\theta) - F(T(\theta))$ is continuous on $T$, the condition H-ii holds.

**Example 3** (to follow). We prove the conditions H-i and H-ii are verified.

H-i follows from A1, which holds true for this example. For H-ii, observe that by assumption $U[\theta] = \{(\theta, \gamma), \gamma \in [g_-, g_+]\}$ and remember that $F$ and $T$ are continuous functions on $T$ and $U$ respectively (see Proposition 12). If $\theta \in T \setminus \mathcal{L}$ then $\theta$ is not a minimizer of $F$ on $T$, and $\theta \ne T(\theta, \gamma)$
$\partial [G(\cdot; v)] (T(v))$. By Bauschke & Combettes (2011, Cor. 16.48 and Prop. 16.6), we have

$$ \partial [G(\cdot; v)] (T(v)) = \gamma \partial [F + \iota_{\mathbb{T}}] (T(v)) + \nabla b(T(v)) - \nabla b(\theta). $$

Therefore, the condition $T(v) = \theta$ is equivalent to $0 \in \partial [F + \iota_{\mathbb{T}}](\theta)$ and the set $\mathcal{L}$ is the set of the minimizers of $F$ on $\mathbb{T}$. Note also that we have $\mathcal{L} = \{\theta \in \mathbb{T} : \exists v \in \mathbb{U}[\theta], T(v) = \theta\}$.

The proof of asymptotic convergence consists in observing that **SAM2** is a stochastic perturbation of the **Population MM** iterative scheme $\tau^{t+1} = T(v^t)$ for $v^t \in \mathbb{U}[\tau^t]$. While the sequence $\{\tau^t, t \ge 0\}$ possesses a Lyapunov function i.e. $F(\tau^{t+1}) \le F(\tau^t)$ for all $t \ge 0$, this is no longer the case for the **SAM2** sequence $\{\theta^t, t \ge 0\}$. Indeed, under A3-b, $\theta^{t+1}$ is the minimizer of an approximation of the majorizing function $G(\cdot; v^t)$ so that in general, $F(\theta^{t+1}) \ne F(T(v^t))$. Nevertheless, the following result shows that whenever the error $|F(\theta^{t+1}) - F(T(v^t))|$ vanishes as $t \to +\infty$, the sequence $\{\theta^t, t \ge 0\}$ inherits the same limiting behavior as the **Population MM** sequence. The proof of Theorem 5 is given in Section 6.3 and is adapted from (Fort & Moulines, 2003, Proposition 9), which addresses the case $\mathbb{U}[\theta] = \{\theta\}$ for all $\theta \in \mathbb{T}$.

**Theorem 5.** Let $\mathbb{T}$ be a compact subset of $\mathbb{R}^d$, $\mathbb{U}$ be a subset of $\mathbb{R}^d$ and let $\mathcal{L} \subseteq \mathbb{T}$ be a set such that $\mathcal{L} \cap \mathbb{T}$ is compact. Let $T: \mathbb{U} \to \mathbb{T}$ and $F: \mathbb{T} \to \mathbb{R}$ be a continuous function such that

(H-i) for all $\theta \in \mathbb{T}$, there exists $\mathbb{U}[\theta] \subseteq \mathbb{U}$ such that $F(T(v)) \le F(\theta)$ for all $v \in \mathbb{U}[\theta]$,

(H-ii) for any compact subset $K$ in $\mathbb{T} \setminus \mathcal{L}$, $\inf_{\theta \in K} \inf_{v \in \mathbb{U}[\theta]} (F(\theta) - F(T(v))) > 0$.

Further let $\{(\theta^t, v^t), t \ge 0\}$ be a $\mathbb{T} \times \mathbb{U}$-valued sequence such that $v^t \in \mathbb{U}[\theta^t]$ for all $t \ge 0$ and assume that

$$ (H-iii) \lim_t |F(T(v^t)) - F(\theta^{t+1})| = 0. $$

Then the sequence $\{F(\theta^t), t \ge 0\}$ converges to a connected component of $F(\mathcal{L} \cap \mathbb{T})$. If $F(\mathcal{L} \cap \mathbb{T})$ has an empty interior, the sequence $\{F(\theta^t), t \ge 0\}$ converges to $F^*$ and the sequence $\{\theta^t, t \ge 0\}$ converges to the set $\mathbb{T} \cap \mathcal{L}_{F^*}$ where $\mathcal{L}_{F^*} := \{\theta \in \mathcal{L} : F(\theta) = F^*\}$.

Let us start with discussing H-i and H-ii in two examples.

**Example 6.** Assume A1, A2 and A3, and $\mathbb{U}[\theta] = \{\theta\}$ for all $\theta \in \mathbb{T}$. We prove the conditions H-i and H-ii are verified with $\mathcal{L} := \{\theta \in \mathbb{T} : F(\theta) - F(T(\theta)) = 0\} = \{\theta \in \mathbb{T} : T(\theta) = \theta\}$.

A3-b and A3-c imply that $\mathbb{T}$ is a point-to-point map and $F$ is a continuous function. Since $F$ and $\mathbb{T}$ are continuous on $\mathbb{T}$ (see A3-c and Proposition 13), the set $\mathcal{L}$ is closed which implies that $\mathcal{L} \cap \mathbb{T}$ is compact. By A1 (see also (7)), $F(T(\theta)) \le F(\theta)$ for all $\theta \in \mathbb{T}$, so that H-i holds. By definition of $\mathcal{L}$, $F(\theta) - F(T(\theta)) > 0$ for all $\theta \in \mathbb{T} \setminus \mathcal{L}$; since $\theta \mapsto F(\theta) - F(T(\theta))$ is continuous on $\mathbb{T}$, the condition H-ii holds.
for all $(\theta, \gamma) \in \mathbb{U}[\theta]$. Then, $(\theta, \gamma) \mapsto F(\theta) - F(\mathbf{T}(\theta, \gamma))$ is a positive continuous function on $\mathcal{K} \times [g_{-}, g_{+}]$ for any compact set $\mathcal{K} \subset \mathbb{T} \setminus \mathcal{L}$. This establishes H-ii.

When $\gamma$ is not assumed lower bounded away from zero, the condition H-ii may not hold. As a counter-example, consider the minimization of the absolute value function $F(\theta) = |\theta|$ on the compact subset $\mathbb{T} := [-a, a]$ of $\mathbb{R}$; choose $\mathfrak{b}(\tau) := \tau^2$. Then $\mathbf{T}(\theta, \gamma) = \text{sign}(\theta)(|\theta| - \gamma)_+$ which implies that $F(\theta) - F(\mathbf{T}(\theta, \gamma)) = |\theta| \wedge \gamma$. Therefore, $\inf_{\theta \in \mathcal{K}} \inf_{\gamma > 0} F(\theta) - F(\mathbf{T}(\theta, \gamma)) = 0$ for any compact subset $\mathcal{K}$ of $[-a, 0) \cup (0, a]$.

Checking H-iii is specific to each perturbation of Population MM since it relies on the definition of $\theta^{t+1}$ given $\theta^t$ and $v^t$. The following result, whose proof is given in Section 6.4, provides conditions implying that if $\theta'$ is such that $G(\theta'; v) - \min_{\mathbb{T}} G(\cdot; v)$ is small, then $|F(\theta') - F(\mathbf{T}(v))|$ is small.

**Proposition 7.** Assume A 1 and A 3. Assume also that for any $\delta > 0$, there exists $\tilde{\eta}_{\delta} > 0$ such that

$$ (\theta, \theta') \in \mathbb{T} \times \mathbb{T}, v \in \mathbb{U}[\theta] : \|\theta' - \mathbf{T}(v)\| \ge \delta \implies (G(\theta'; v) - G(\mathbf{T}(v); v)) \ge \tilde{\eta}_{\delta}. \quad (9) $$

For any $\epsilon > 0$, there exists $\alpha_{\epsilon} > 0$ such that for any $\theta, \theta' \in \mathbb{T}$ and for all $v \in \mathbb{U}[\theta]$, it holds

$$ |F(\theta') - F(\mathbf{T}(v))| \le \epsilon + \alpha_{\epsilon} \sup_{\mathbb{T}} |F| \text{Diam}(\mathbb{T}) (G(\theta'; v) - G(\mathbf{T}(v); v)), $$

where $\text{Diam}(\mathbb{T})$ denotes the diameter of the compact set $\mathbb{T}$.

**Example 8 (Example 6 to follow).** Let us check the condition (9). Since $\mathbb{U}[\theta] = \{\theta\}$, we have $\theta' - \mathbf{T}(v) = \theta' - \mathbf{T}(\theta)$ and $G(\theta'; v) - G(\mathbf{T}(v); v) = G(\theta'; \theta) - G(\mathbf{T}(\theta); \theta)$. By Proposition 12, $\mathbf{T}$ is a continuous function so that the set $\mathcal{K}_{\delta} := \{(\theta, \theta') \in \mathbb{T} \times \mathbb{T} : \|\theta' - \mathbf{T}(\theta)\| \ge \delta\}$ is a compact subset of $\mathbb{T} \times \mathbb{T}$. $G(\mathbf{T}(\theta); \theta) = \min_{\mathbb{T}} G(\cdot; \theta)$ and the minimizer is unique by A 3-b; therefore, $G(\theta'; \theta) - G(\mathbf{T}(\theta); \theta) > 0$ for all $(\theta', \theta) \in \mathcal{K}_{\delta}$. Finally, $G$ and $\mathbf{T}$ are continuous (see Proposition 11 and Proposition 12). Hence $\tilde{\eta}_{\delta}$ exists.

**Example 3 (to follow).** The proof of (9) is on the same lines as in the Example 8. Set $\mathcal{K}_{\delta} := \{(\theta, \theta', \gamma) \in \mathbb{T} \times \mathbb{T} \times [g_{-}, g_{+}] : \|\theta' - \mathbf{T}(\theta, \gamma)\| \ge \delta\}$. Under the stated assumptions, $\mathbf{T}$ is continuous on $\mathbb{T} \times [g_{-}, g_{+}]$ so $\mathcal{K}_{\delta}$ is compact. In addition, $G(\theta'; v) - G(\mathbf{T}(v); v) > 0$ for all $(\theta', v) \in \mathcal{K}_{\delta}$. Finally, $G$ is a continuous function. Hence $\tilde{\eta}_{\delta}$ exists.

We are now ready to apply Theorem 5 and provide a result on the asymptotic convergence of SAM2. Theorem 9 establishes the almost-sure convergence of the SAM2 sequence to the set $\mathcal{L} = \{\theta \in \mathbb{T} : \mathbf{T}(v) = \theta$ for all $v \in \mathbb{U}[\theta]\}$. In order to satisfy the condition H-iii of Theorem 5, we apply Proposition 7 and provide sufficient conditions for the property $\lim_{t}\{G(\theta^{t+1}; v^t) - G(\mathbf{T}(v^t); v^t)\} = 0$ to hold almost-surely. Upon noting that

$$ 0 \le G(\theta^{t+1}; v^t) - G(\mathbf{T}(v^t); v^t) \le 2 \sup_{(\theta,v) \in \mathbb{T} \times \mathbb{U}} \left| \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta, X^{t+1,i}; v) - G(\theta; v) \right|, $$
**Example 3** (to follow). We prove the conditions H-i and H-ii are verified.

H-i follows from A 1, which holds true for this example. For H-ii, observe that by assumption $U[\theta] = \{(\theta, \gamma), \gamma \in [g_{-}, g_{+}]\}$ and remember that $F$ and $T$ are continuous functions on $T$ and $U$ respectively (see Proposition 13). If $\theta \in T \setminus L$ then $\theta$ is not a minimizer of $F$ on $T$, and $\theta \neq T(\theta, \gamma)$ for all $(\theta, \gamma) \in U[\theta]$. Then, $(\theta, \gamma) \mapsto F(\theta) - F(T(\theta, \gamma))$ is a positive continuous function on $K \times [g_{-}, g_{+}]$ for any compact set $K \subset T \setminus L$. This establishes H-ii.

When $\gamma$ is not assumed bounded away from zero, the condition H-ii may not hold. As a counter-example, consider the minimization of the absolute value function $F(\theta) = |\theta|$ on the compact subset $T := [-a, a]$ of $\mathbb{R}$; choose $b(\tau) := \tau^2$. Then $T(\theta, \gamma) = \text{sign}(\theta)(|\theta| - \gamma)_+$ which implies that $F(\theta) - F(T(\theta, \gamma)) = |\theta| \wedge \gamma$. Therefore, $\inf_{\theta \in K} \inf_{\gamma > 0} F(\theta) - F(T(\theta, \gamma)) = 0$ for any compact subset $K$ of $[-a, 0) \cup (0, a]$.

Checking H-iii is specific to each perturbation of Population MM since it relies on the definition of $\theta^{t+1}$ given $\theta^t$ and $v^t$. The following result, whose proof is given in Section 6.4, provides conditions implying that if $\theta'$ is such that $G(\theta'; v) - \min_T G(\cdot; v)$ is small, then $|F(\theta') - F(T(v))|$ is small.

**Proposition 7.** Assume A 1 and A 3. Assume also that for any $\delta > 0$, there exists $\tilde{\eta}_{\delta} > 0$ such that

$$ (\theta, \theta') \in T \times T, v \in U[\theta] : \|\theta' - T(v)\| \ge \delta \implies (G(\theta'; v) - G(T(v); v)) \ge \tilde{\eta}_{\delta}. \quad (9) $$

Then, for any $\epsilon > 0$, there exists $\alpha_{\epsilon} > 0$ such that for any $\theta, \theta' \in T$ and for all $v \in U[\theta]$, it holds that

$$ |F(\theta') - F(T(v))| \le \epsilon + \alpha_{\epsilon} \sup_{T} |F| \text{Diam}(T) (G(\theta'; v) - G(T(v); v)), $$

where $\text{Diam}(T)$ denotes the diameter of the compact set $T$.

**Example 8** (Example 6 to follow). Let us check condition (9). Since $U[\theta] = \{\theta\}$, we have $\theta' - T(v) = \theta' - T(\theta)$ and $G(\theta'; v) - G(T(v); v) = G(\theta'; \theta) - G(T(\theta); \theta)$. By Proposition 13, $T$ is a continuous function so that the set $K_{\delta} := \{(\theta, \theta') \in T \times T : \|\theta' - T(\theta)\| \ge \delta\}$ is a compact subset of $T \times T$. $G(T(\theta); \theta) = \min_T G(\cdot; \theta)$ and the minimizer is unique by A 3-b; therefore, $G(\theta'; \theta) - G(T(\theta); \theta) > 0$ for all $(\theta', \theta) \in K_{\delta}$. Finally, $G$ and $T$ are continuous (see Proposition 12 and Proposition 13). Hence $\tilde{\eta}_{\delta}$ exists.

**Example 3** (to follow). The proof of (9) is on the same lines as in the Example 8. Set $K_{\delta} := \{(\theta, \theta', \gamma) \in T \times T \times [g_{-}, g_{+}] : \|\theta' - T(\theta, \gamma)\| \ge \delta\}$. Under the stated assumptions, $T$ is continuous on $T \times [g_{-}, g_{+}]$ so $K_{\delta}$ is compact. In addition, $G(\theta'; v) - G(T(v); v) > 0$ for all $(\theta', v) \in K_{\delta}$. Finally, $G$ is a continuous function. Hence $\tilde{\eta}_{\delta}$ exists.

We are now ready to apply Theorem 5 and provide a result on the asymptotic convergence of SAM2. Theorem 9 establishes the almost-sure convergence of the SAM2 sequence to the set $\mathcal{L} = \{\theta \in T : T(v) = \theta$ for all $v \in U[\theta]\}$. In order to satisfy the condition H-iii of Theorem 5, we apply
(see Section 6.6), we use a uniform strong law of large numbers (see Section 6.5) and assume that the batch size $N_t$ increases rapidly enough.

A4. The size $\{N_t, t \ge 1\}$ of the mini-batch #t satisfies $\sum_{t \ge 1} N_t^{-(p-1)(p/2)} < \infty$, where p is given by A3-a.

When $p \in (1, 2]$, the batch size $N_t$ can be set to $N_t \propto t^{1/(p-1)}(\ln t)^{(1+\beta)/(p-1)}$ for $\beta > 0$. When $p \ge 2$, the batch size $N_t$ can be set to $N_t \propto t^{2/p}(\ln t)^{(1+\beta)2/p}$ for $\beta > 0$. A large value of p yields a slow increase of the batch size $N_t$.

**Theorem 9.** Assume A2, A3-a, A3-b, A4 and $\mathbb{U}$ is compact. Let $\{\theta^t, t \ge 0\}$ be the output of the SAM2 algorithm. Then with probability one

$$i) \lim_t \{G(\theta^{t+1}; v^t) - G(T(v^t); v^t)\} = 0.$$

Set

$$\mathcal{L}^+ := \{\theta \in \mathbb{T} : \exists v \in \mathbb{U}[\theta], T(v) = \theta\}.$$

If in addition A1 and A3-c hold, then with probability one,

$$ii) \{F(\theta^t), t \ge 0\} \text{ converges to a connected component of } F(\mathcal{L}^+) .$$

$$iii) \text{ If } F(\mathcal{L}^+) \text{ has an empty interior, then } \{F(\theta^t), t \ge 0\} \text{ converges to } F^* \text{ and } \{\theta^t, t \ge 0\} \text{ converges to } \{\theta \in \mathcal{L}^+ : F(\theta) = F^*\}.$$

The proof of Theorem 9 is given in Section 6.6. When $\mathbb{U}[\theta] = \{\theta\}$ for all $\theta \in \mathbb{T}$, then $\mathcal{L}^+ = \mathcal{L}$ where $\mathcal{L}$ is given by (7). Note also that for Example 3, $\mathcal{L}^+ = \mathcal{L}$ even though $\mathbb{U}[\theta]$ is not reduced to a singleton.

When the majorizing function $G(\cdot; v^t)$ can be exactly evaluated at each iteration, it is known that the convergence of the objective function along the path $\{\tau^t, t \ge 0\}$ of the Population MM algorithm follows from the monotonicity of $t \mapsto F(\tau^t)$ (see e.g. (Lange et al., 2021, Proposition 2.1)). When only oracles of the majorizing function are available, the monotonicity property along the SAM2 sequence does not hold anymore; nevertheless, the sequence $\{F(\theta^t), t \ge 0\}$ may still converge, as shown by Theorem 9.

Since there is no guarantee that for the SAM2 sequence, $\theta^{t+1} \in \{\theta \in \mathbb{T} : F(\theta) \le F(\theta^t)\}$, the stability of the SAM2 sequence can not follow from the compactness of the level sets of the descent function $F$. In this paper, the stability is ensured by forcing a compact-valued sequence (see (5) and the compactness assumption on $\mathbb{T}$). Nevertheless, a self-stabilization by projections on growing compact sets could be explored, by adapting (Fort & Moulines, 2003, Propositions 10 and 11). Such an extension is out of the scope of this paper.

Iterative algorithms having a descent function converge to the set of the so-called no-progress points (see e.g. (Lange et al., 2021, Proposition 2.2)). Under A1 and the assumption $\mathbb{U}[\theta] = \{\theta\}$, we have $F(\tau) \le F(\theta)$ for all $\tau \in \text{Argmin}_{\mathbb{T}} G(\cdot; \theta)$ and the inequality is strict when $\theta \notin \text{Argmin}_{\mathbb{T}} G(\cdot; \theta)$.
Proposition 7 and provide sufficient conditions for the property $\lim_t \{G(\theta^{t+1}; v^t) - G(T(v^t); v^t)\} = 0$ to hold almost surely. Upon noting that

$$0 \le G(\theta^{t+1}; v^t) - G(T(v^t); v^t) \le 2 \sup_{(\theta, v) \in \mathbb{T} \times \mathbb{U}} \left| \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta, X^{t+1, i}; v) - G(\theta; v) \right|,$$

(see Section 6.6), we use a double array uniform strong law of large numbers (see Section 6.5) and assume that the batch size $N_t$ increases rapidly enough.

A4. The sizes $\{N_t, t \ge 1\}$ of the batch samples satisfy $\sum_{t \ge 1} N_t^{-(p-1)\wedge(p/2)} < \infty$, where $p$ is given by A3-a.

When $p \in (1, 2]$, the batch size $N_t$ can be set to $N_t \propto t^{1/(p-1)}(\ln t)^{(1+\beta)/(p-1)}$ for $\beta > 0$. When $p \ge 2$, the batch size $N_t$ can be set to $N_t \propto t^{2/p}(\ln t)^{(1+\beta)2/p}$ for $\beta > 0$. A large value of $p$ yields a slower increase in the batch sizes $N_t$.

**Theorem 9.** Assume A2, A3-a, A3-b, A4 and $\mathbb{U}$ is compact. Let $\{\theta^t, t \ge 0\}$ be the output of the SAM2 algorithm. Then with probability one,

$$i) \lim_t \{G(\theta^{t+1}; v^t) - G(T(v^t); v^t)\} = 0.$$

Set

$$\mathcal{L}^+ := \{\theta \in \mathbb{T} : \exists v \in \mathbb{U}[\theta], T(v) = \theta\}.$$

If in addition A1 and A3-c hold, then with probability one,

ii) $\{F(\theta^t), t \ge 0\}$ converges to a connected component of $F(\mathcal{L}^+)$.

iii) If $F(\mathcal{L}^+)$ has an empty interior, then $\{F(\theta^t), t \ge 0\}$ converges to $F^*$ and $\{\theta^t, t \ge 0\}$ converges to $\{\theta \in \mathcal{L}^+ : F(\theta) = F^*\}$.

The proof of Theorem 9 is given in Section 6.6. When $\mathbb{U}[\theta] = \{\theta\}$ for all $\theta \in \mathbb{T}$, then $\mathcal{L}^+ = \mathcal{L}$, where $\mathcal{L}$ is given by (6). Note also that for Example 3, $\mathcal{L}^+ = \mathcal{L}$ even though $\mathbb{U}[\theta]$ is not reduced to a singleton.

When the majorizing function $G(\cdot; v^t)$ can be exactly evaluated at each iteration, it is known that the convergence of the objective function along the path $\{\tau^t, t \ge 0\}$ of the Population MM algorithm follows from the monotonicity of $t \mapsto F(\tau^t)$ (see e.g. Lange et al., 2021, Proposition 2.1). When only oracles of the majorizing function are available, the monotonicity property along the SAM2 sequence does not hold anymore; nevertheless, the sequence $\{F(\theta^t), t \ge 0\}$ may still converge, as shown by Theorem 9.

Since there is no guarantee that for the SAM2 sequence, $\theta^{t+1} \in \{\theta \in \mathbb{T} : F(\theta) \le F(\theta^t)\}$, the stability of the SAM2 sequence can not follow from the compactness of the level sets of the function $F$. In this paper, the stability is ensured by forcing a compact-valued sequence (see (5) and
Hence, the set of the no-progress points is the set $\mathcal{L}$ given by (6). Theorem 9 provides sufficient conditions for the SAM2 sequence to converge to $\mathcal{L}$ (observe that $\mathcal{L} = \mathcal{L}^+$ when $\mathbb{U}[\theta] = \{\theta\}$ and $\text{Argmin}_{\mathbb{T}}\mathcal{G}(\cdot; \theta)$ has exactly one element).

# 4 SAM2 examples

## 4.1 An Online Mirror Descent Algorithm

The goal is to minimize $F: \mathbb{R}^d \to \mathbb{R}$ on a compact convex set $\mathbb{T}$, via a Mirror Descent approach (see Nemirovskii & Yudin (1983), see also the overview Bubeck (2015)). Consider the case when

MD1. (a) $F$ is continuously differentiable in an open neighborhood $\mathcal{S}$ of $\mathbb{T}$. There exists a measurable function $h: \mathbb{T} \times \mathbb{X} \to \mathbb{R}^d$ such that $\nabla F(\theta) = \mathbb{E}[h(\theta, X)]$ where $\mathbb{E}[\|h(\theta, X)\|] < \infty$ for all $\theta \in \mathbb{T}$. No exact computation of the expectation is available.

(b) A stream of independent random variables with the same distribution as $X$ is available, and $h(\theta, x)$ can be computed for all $(\theta, x) \in \mathbb{T} \times \mathbb{X}$.

Choose a function $\psi$ satisfying

MD2. $\psi: \mathbb{R}^d \to (-\infty, +\infty]$ is a lower semicontinuous convex function, $\mathbb{T} \subset \text{dom}(\psi)$, and it is a continuously differentiable strictly convex function on $\mathcal{S}$. The function $b := \psi - F$ is strictly convex.

Define the Bregman divergences associated to $\psi$ and $b$

$$D_{\psi}(\theta; \tau) := \psi(\theta) - \psi(\tau) - \langle \nabla \psi(\tau), \theta - \tau \rangle, \quad D_{b}(\theta; \tau) := b(\theta) - b(\tau) - \langle \nabla b(\tau), \theta - \tau \rangle.$$

The Online Mirror Descent algorithm is given by Algorithm 3. The update mechanism (10) is

**Algorithm 3** The Online Mirror Descent algorithm

**Require:** a sequence $\{N_t, t \ge 1\}$ of positive integers, an initial value $\theta^0 \in \mathbb{T}$

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 0\}$.

1: **for** $t = 0, \dots, \mathbf{do}$

2: &nbsp;&nbsp;&nbsp;&nbsp;Sample a minibatch $\{X^{t+1,i}, i = 1, \dots, N_{t+1}\}$ of size $N_{t+1}$.

3: &nbsp;&nbsp;&nbsp;&nbsp;Set $\mathcal{H}^{t+1} := N_{t+1}^{-1} \sum_{i=1}^{N_{t+1}} h(\theta^t, X^{t+1,i})$

4: &nbsp;&nbsp;&nbsp;&nbsp;Compute

$$\theta^{t+1} \in \underset{\theta \in \mathbb{T}}{\operatorname{Argmin}} \langle \mathcal{H}^{t+1}, \theta \rangle + D_{\psi}(\theta; \theta^t). \qquad (10)$$

5: **end for**

equivalent to

$$\theta^{t+1} = \underset{\theta \in \mathbb{T}}{\operatorname{Argmin}}_{\theta \in \mathbb{T}} (\psi(\theta) - \psi(\theta^t) - \langle \nabla \psi(\theta^t) - \mathcal{H}^{t+1}, \theta - \theta^t \rangle)$$
the compactness assumption on $\mathbb{T}$). Nevertheless, a self-stabilization by projections on growing compact sets could be explored, by adapting (Fort & Moulines, 2003, Propositions 10 and 11). Such an extension is out of the scope of this paper.

Iterative algorithms having a descent function converge to the set of so-called no-progress points (see e.g. Lange et al., 2021, Proposition 2.2). Under A1 and the assumption that $U[\theta] = \{\theta\}$, we have $F(\tau) \le F(\theta)$ for all $\tau \in \text{Argmin}_{\mathbb{T}}G(\cdot; \theta)$ and the inequality is strict when $\theta \notin \text{Argmin}_{\mathbb{T}}G(\cdot; \theta)$. Hence, the set of the no-progress points is the set $\mathcal{L}$, given by (6). Theorem 9 provides sufficient conditions for the SAM2 sequence to converge to $\mathcal{L}$ (observe that $\mathcal{L} = \mathcal{L}^+$ when $U[\theta] = \{\theta\}$ and $\text{Argmin}_{\mathbb{T}}G(\cdot; \theta)$ has exactly one element).

# 4 SAM2 examples

We show that some online versions of Mirror Descent and Proximal-Gradient are examples of SAM2 under conditions on the objective function F. For each example, we then provide sufficient conditions for verifying the assumptions of Theorem 9 and identifying the set $\mathcal{L}^+$.

## 4.1 An Online Mirror Descent Algorithm

The goal is to minimize $F: \mathbb{R}^d \to \mathbb{R}$ on a compact convex set $\mathbb{T}$, via a Mirror Descent approach (see Nemirovskii & Yudin, 1983, and the overview by Bubeck, 2015). Suppose that the following assumptions hold.

MD1. (a) F is continuously differentiable on an open set $S \supset \mathbb{T}$. There exists a measurable function $h: \mathbb{T} \times \mathbb{X} \to \mathbb{R}^d$ such that $\nabla F(\theta) = \mathbb{E}[h(\theta, X)]$, where $\mathbb{E}[||h(\theta, X)||] < \infty$ for all $\theta \in \mathbb{T}$, and exact computation of the expectation is not available.

(b) A stream of independent random variables with the same distribution as $X$ is available, and $h(\theta, x)$ can be computed for all $(\theta, x) \in \mathbb{T} \times \mathbb{X}$.

Then, choose a function $\psi$ satisfying the following conditions.

MD2. $\psi: \mathbb{R}^d \to (-\infty, +\infty]$ is a lower semi-continuous convex function, $\mathbb{T} \subset \text{dom}(\psi)$, and it is a continuously differentiable strictly convex function on $S$.

MD3. There exists $\gamma_\star > 0$ such that the function $\gamma_\star^{-1}\psi - F$ is convex on $\mathbb{T}$.

Observe that under MD2 and MD3, $b_\gamma(\cdot) := \gamma^{-1}\psi(\cdot) - F(\cdot)$ is convex for all $\gamma \in (0, \gamma_\star]$ and strictly convex for $\gamma \in (0, \gamma_\star)$. When $\psi$ is $\alpha$-strongly convex and $\nabla F$ is $L$-lipschitz, then MD3 is verified with $\gamma_\star := \alpha/L$: for any $\theta, \tau \in \mathbb{T}$,

$$F(\theta) - F(\tau) - \langle \nabla F(\tau), \theta - \tau \rangle \le \frac{L}{2} \|\theta - \tau\|^2 \le \frac{L}{\alpha} (\psi(\theta) - \psi(\tau) - \langle \nabla \psi(\tau), \theta - \tau \rangle),$$
thus showing that it is Mirror Descent applied with a constant step size $\gamma$ (set to one here without loss of generality, since $\psi$ can be replaced with $\psi/\gamma$). Note also that the minimizer exists and is unique since $\psi$ is strictly convex on $\mathbb{T}$.

**Online Mirror Descent is a SAM2 algorithm.** Following the same lines as in (Beck & Teboulle, 2003, section 3) (see also (Bubeck, 2015, Section 4.2.), (Lange et al., 2021, Section 4.2.2)), we prove that under MD1 and MD2, Algorithm 3 is a SAM2 algorithm. Set $\mathcal{U} := \mathbb{T}$ and define

$$G(\theta; v) := F(v) + \psi(\theta) - \psi(v) - \langle \nabla\psi(v) - \nabla F(v), \theta - v \rangle + \iota_{\mathbb{T}}(\theta), \quad \theta \in \mathbb{R}^d, v \in \mathbb{T}; \quad (11)$$

$\iota_{\mathbb{T}}$ is the characteristic function of $\mathbb{T}$. For $\theta, v \in \mathbb{T}$, we have $G(\theta; v) = F(\theta) + D_b(\theta; v) + \iota_{\mathbb{T}}(\theta)$. Therefore, $G(\theta, v) \ge F(\theta)$ for all $\theta, v \in \mathbb{T}$ since $D_b(\theta; v) \ge 0$ under MD2. In addition, $G(\theta; \theta) = F(\theta)$ for all $\theta \in \mathbb{T}$ since $D_b(\theta; \theta) = 0$, and since $b$ is strictly convex, we have $D_b(\theta; v) > 0$ for $\theta \ne v$. Hence, $\mathcal{U}[\theta] = \{\theta\}$. From (11), we have $G(\theta, v) = \mathbb{E}[g(\theta, X; v)]$ where $g$ is given by

$$g(\theta, x; v) := \langle h(v, x), \theta \rangle + D_{\psi}(\theta; v) + C(v) + \iota_{\mathbb{T}}(\theta);$$

$C(v)$ does not depend on $\theta$. Hence, it is readily seen that solving (5) is equivalent to solving (10): Online Mirror Descent is a SAM2 algorithm.

**The limiting set.** Under MD1, MD2 and

MD3. for all $x \in \mathbb{X}$, $\theta \mapsto h(\theta, x)$ is continuous on $\mathbb{T}$ and $\mathbb{E}[\sup_{\theta \in \mathbb{T}} \|h(\theta, X)\|^p] < \infty$ for some $p > 1$ the conditions A1, A2 and A3 are satisfied. Since $\mathcal{U}[\theta] = \{\theta\}$, then $\mathcal{L}^+ = \mathcal{L}$; let us identify the set $\mathcal{L}$ given by (6). Since $\theta' \mapsto G(\theta'; \theta)$ is a lower semicontinuous convex function with domain $\mathbb{T}$, $\mathbb{T}(\theta)$ solves (see e.g. (Bauschke & Combettes, 2011, Theorem 16.3))

$$0 \in \partial[\psi(\cdot) - \langle \nabla\psi(\theta) - \nabla F(\theta), \cdot \rangle + \iota_{\mathbb{T}}(\cdot)](\mathbb{T}(\theta)).$$

Hence, $\mathbb{T}(\theta) = \theta$ is equivalent to $0 \in \nabla F(\theta) + \partial\iota_{\mathbb{T}}(\mathbb{T}(\theta))$ by (Bauschke & Combettes, 2011, Corollary 16.48). Hence $\mathcal{L}$ is the set of the minimizers of $F$ on $\mathbb{T}$.

For the sequence $\{\theta^t, t \ge 0\}$ given by Algorithm 3, Theorem 9 shows that, as soon as the size of the mini batches is chosen so that A4 holds: the sequence $\{F(\theta^t), t \ge 0\}$ converges to $\min_{\mathbb{T}} F$ and the sequence $\{\theta^t, t \ge 0\}$ converges to the set of the minimizers.

When the majorizing function $G$ relies on a Bregman divergence, as it is the case in this example, the convergence analyses usually adopt a different approach which uses the three-point inequality by Chen & Teboulle (1993): the function that plays the role of a descent function is not $F$ but it is a Bregman distance between the current iterate and a minimizer. In that approach, $\epsilon$-approximate stationarity is generally discussed, most often along an averaged path (see e.g. Nemirovski et al. (2009); Lan et al. (2012); Nedić & Lee (2014); Zhang & He (2018);
which implies the convexity inequality for $\mathfrak{b}_\gamma$. Define the Bregman divergences

$$D_{\psi}(\theta; \tau) := \psi(\theta) - \psi(\tau) - \langle \nabla\psi(\tau), \theta - \tau \rangle, \quad D_{\mathfrak{b}_{\gamma}}(\theta; \tau) := \mathfrak{b}_{\gamma}(\theta) - \mathfrak{b}_{\gamma}(\tau) - \langle \nabla\mathfrak{b}_{\gamma}(\tau), \theta - \tau \rangle.$$

An Online Mirror Descent algorithm is given by Algorithm 3. Observe that update mechanism

**Algorithm 3 Online Mirror Descent**

**Require:** a sequence $\{N_t, t \ge 1\}$ of positive integers, an initial value $\theta^0 \in \mathbb{T}$.

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 0\}$.

1: **for** $t \ge 0$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;Sample a batch $\{X^{t+1,i}, 1 \le i \le N_{t+1}\}$ of size $N_{t+1}$.

3: &nbsp;&nbsp;&nbsp;&nbsp;Set $\mathcal{H}^{t+1} := N_{t+1}^{-1} \sum_{i=1}^{N_{t+1}} \mathbf{h}(\theta^t, X^{t+1,i})$.

4: &nbsp;&nbsp;&nbsp;&nbsp;Choose $\gamma_{t+1} \in (0, \gamma_*)$.

5: &nbsp;&nbsp;&nbsp;&nbsp;Compute

$$\theta^{t+1} \in \underset{\theta \in \mathbb{T}}{\operatorname{Argmin}} \langle \mathcal{H}^{t+1}, \theta \rangle + \frac{1}{\gamma_{t+1}} D_{\psi}(\theta; \theta^t). \quad (10)$$

6: **end for**

(10) is equivalent to

$$\theta^{t+1} \in \underset{\theta \in \mathbb{T}}{\operatorname{Argmin}} \left( \frac{1}{\gamma_{t+1}} \psi(\theta) - \frac{1}{\gamma_{t+1}} \psi(\theta^t) - \left\langle \frac{1}{\gamma_{t+1}} \nabla\psi(\theta^t) - \mathcal{H}^{t+1}, \theta - \theta^t \right\rangle \right). \quad (11)$$

The minimizer exists and is unique since $\gamma_{t+1}^{-1}\psi$ is strictly convex on $\mathbb{T}$ by MD2.

**Algorithm 3 is a SAM2 algorithm.** Following the same arguments as in Beck & Teboulle (2003, Sec. 3) (see also Bubeck, 2015, Sec. 4.2 and Lange et al., 2021, Sec. 4.2.2), we prove that under MD1, MD2 and MD3, Algorithm 3 is a SAM2 algorithm. Set $\mathbb{U} := \mathbb{T} \times (0, \gamma_*)$ and for $\theta \in \mathbb{R}^d$, $\nu := (\tau, \gamma) \in \mathbb{U}$, define

$$\mathcal{G}(\theta; \nu) := \mathbf{F}(\tau) + \frac{1}{\gamma}\psi(\theta) - \frac{1}{\gamma}\psi(\tau) - \left\langle \frac{1}{\gamma}\nabla\psi(\tau) - \nabla\mathcal{F}(\tau), \theta - \tau \right\rangle + \iota_{\mathbb{T}}(\theta); \quad (12)$$

$\iota_{\mathbb{T}}$ is the characteristic function of $\mathbb{T}$. For $\theta \in \mathbb{T}$, $\nu = (\tau, \gamma) \in \mathbb{U}$, we have $\mathcal{G}(\theta; \nu) = \mathbf{F}(\theta) + D_{\mathfrak{b}_{\gamma}}(\theta; \tau) + \iota_{\mathbb{T}}(\theta)$. Therefore, $\mathcal{G}(\theta, \nu) \ge \mathbf{F}(\theta)$, for all $\theta \in \mathbb{T}$, $\nu \in \mathbb{U}$ since $D_{\mathfrak{b}_{\gamma}}(\theta; \tau) \ge 0$ under MD3. In addition, $\mathcal{G}(\theta; (\theta, \gamma)) = \mathbf{F}(\theta)$, for all $\theta \in \mathbb{T}$ and $\gamma \in (0, \gamma_*)$ since $D_{\mathfrak{b}_{\gamma}}(\theta; \tau) = 0$. Hence, $\mathbb{U}[\theta] = \{(\tau, \gamma) \in \mathbb{T} \times (0, \gamma_*], D_{\mathfrak{b}_{\gamma}}(\theta; \tau) = 0\} \supseteq \{(\theta, \gamma), \gamma \in (0, \gamma_*)\}$. From (12), we have $\mathcal{G}(\theta, \nu) = \mathbb{E}[\mathbf{g}(\theta, X; \nu)]$, where $\mathbf{g}$ is given by

$$\mathbf{g}(\theta, x; \nu) := \mathbf{F}(\tau) + \langle \mathbf{h}(\tau, x), \theta - \tau \rangle + \frac{1}{\gamma} D_{\psi}(\theta; \tau) + \iota_{\mathbb{T}}(\theta), \quad \nu = (\tau, \gamma). \quad (13)$$

Hence, it is readily seen from (13) and (11) that solving (5) is equivalent to solving (10); i.e., Online Mirror Descent is a SAM2 algorithm.
Lan (2020); Dragomir et al. (2021); D'Orazio et al. (2023)), from which asymptotic convergence can be obtained (see e.g. (Beck & Teboulle, 2003, Theorem 4.1) in the deterministic case, or Lei & Zhou (2020) for convergence in expectation of Online Mirror Descent). While our approach addresses the almost-sure convergence under a constant step size, but at the price of a mini-batch size $N_t$ increasing at each iteration, the other approach considers the case $N_t = 1$ but the step size is iteration-varying and its choice is crucial.

## 4.2 An Online Proximal-Gradient algorithm

The objective is to minimize a composite function $F := F_s + F_c$ on a compact convex subset $\mathbb{T}$ of $\mathbb{R}^d$, where

PG1. (a) the function $F_s: \mathbb{R}^d \to \mathbb{R}$ is continuously differentiable with $L$-Lipschitz gradient. There exists a measurable function $h: \mathbb{T} \times \mathbb{X} \to \mathbb{R}^d$ such that $\nabla F_s(\theta) = \mathbb{E}[h(\theta, X)]$ where $\mathbb{E}[\|h(\theta, X)\|] < \infty$ for all $\theta \in \mathbb{T}$. No exact computation of the expectation is available.

(b) A stream of independent random variables with the same distribution as $X$ is available, and $h(\theta, x)$ can be computed for all $(\theta, x) \in \mathbb{T} \times \mathbb{X}$.

(c) the function $F_c: \mathbb{R}^d \to (-\infty, +\infty]$ is a continuous convex function, $\mathbb{T} \subset \text{dom}(F_c)$ and for all $\gamma > 0$, the proximal map of $\gamma F_c + \iota_{\mathbb{T}}$ is a point-to-point map with an explicit expression.

$\iota_{\mathbb{T}}$ denotes the characteristic function of the set $\mathbb{T}$; remember that the proximal map is given by

$$ \text{Prox}_{\gamma F_c + \iota_{\mathbb{T}}}(\theta) := \text{Argmin}_{\tau \in \mathbb{T}} \left( \gamma F_c(\tau) + \frac{1}{2} \|\tau - \theta\|^2 \right). $$

The Online Proximal-Gradient algorithm, given by Algorithm 4, is a stochastic perturbation of a forward-backward splitting process in optimization (see Lions & Mercier (1979), see also (Beck & Teboulle, 2009, Section 1.3.3.) for an introduction as a MM technique) in which the forward operator is not explicit and learnt from a stream of oracles (see e.g. (Combettes & Wajs, 2005, Section 3) for a pioneering contribution on inexact proximal-gradient algorithms).

**Algorithm 4** The Online Proximal-Gradient algorithm

**Require:** a sequence $\{N_t, t \ge 1\}$ of positive integers, an initial value $\theta^0 \in \mathbb{T}$, $\epsilon > 0$

**Ensure:** A $\mathbb{T}$-valued sequence $\{\theta^t, t \ge 0\}$.

1: **for** $t = 0, \dots, \mathbf{do}$

2: &nbsp;&nbsp;&nbsp;&nbsp;Sample a minibatch $\{X^{t+1, i}, i = 1, \dots, N_{t+1}\}$ of size $N_{t+1}$.

3: &nbsp;&nbsp;&nbsp;&nbsp;Set $\mathcal{H}^{t+1} := N_{t+1} \sum_{i=1}^{N_{t+1}} h(\theta^t, X^{t+1, i})$

4: &nbsp;&nbsp;&nbsp;&nbsp;Choose $\gamma_{t+1} \in [\epsilon, 1/L]$

5: &nbsp;&nbsp;&nbsp;&nbsp;Compute

$$ \theta^{t+1} = \text{Prox}_{\gamma_{t+1}F_c + \iota_{\mathbb{T}}} (\theta^t - \gamma_{t+1}\mathcal{H}^{t+1}). \quad (12) $$

6: **end for**
**The limiting set when the step sizes are constant.** The convergence result Theorem 9 requires the set $\mathbb{U}$ to be compact. In this paragraph, we consider the case when the step size sequence is constant and equal to $\bar{\gamma} \in (0, \gamma_*)$ (see Section 4.2 for another strategy in a similar context). Make the following additional assumptions.

**MD4.** For all $x \in \mathbb{X}$, $\theta \mapsto h(\theta, x)$ is continuous on $\mathbb{T}$ and $\mathbb{E}[\sup_{\theta \in \mathbb{T}} \|h(\theta, X)\|^p] < \infty$ for some $p > 1$.

**MD5.** The function $\bar{\gamma}^{-1}\psi - F$ is strictly convex.

Under MD 1, MD 2, MD 3 and MD 4, conditions A1, A2 and A3 are satisfied. Under MD 5, $D_{\bar{\gamma}}(\theta, \tau) = 0$ iff $\theta = \tau$. This implies that $\mathbb{U}[\theta] = \{(\theta, \bar{\gamma})\}$ and yields $\mathcal{L}^+ = \mathcal{L}$, where we identify the set $\mathcal{L}$ by (6) and the set $\mathcal{L}^+$ as in Theorem 9.

Since $\theta' \mapsto G(\theta'; \theta)$ is a lower semi-continuous convex function with domain $\mathbb{T}$, $T(\theta)$ solves (cf. Bauschke & Combettes, 2011, Theorem 16.3)

$$0 \in \partial[\psi(\cdot) - \langle \nabla\psi(\theta) - \nabla F(\theta), \cdot \rangle + \iota_{\mathbb{T}}(\cdot)] (T(\theta)).$$

Hence, $T(\theta) = \theta$ is equivalent to $0 \in \nabla F(\theta) + \partial \iota_{\mathbb{T}}(T(\theta))$ by Bauschke & Combettes (2011, Corollary 16.48). Hence, under the additional assumption that $F$ is convex, $\mathcal{L}$ is the set of the minimizers of $F$ on $\mathbb{T}$. Then, for the sequence $\{\theta^t, t \ge 0\}$ given by Algorithm 3, Theorem 9 shows that, when the sizes of the batches are chosen so that A4 holds: the sequence $\{F(\theta^t), t \ge 0\}$ converges to $\min_{\mathbb{T}} F$ and the sequence $\{\theta^t, t \ge 0\}$ converges to the set of minimizers.

**Remark.** When the majorizing function $G$ relies on a Bregman divergence, as in this example, the convergence analyses usually adopt a different approach which uses the three-point inequality of Chen & Teboulle (1993). In such an approach, $\epsilon$-approximate stationarity is generally discussed, most often along an averaged path (see e.g. Nemirovski et al., 2009; Lan et al., 2012; Nedić & Lee, 2014; Lan, 2020; Dragomir et al., 2021; D'Orazio et al., 2023), from which asymptotic convergence can be obtained (see e.g. Beck & Teboulle, 2003, Thm. 4.1, in the deterministic case, or Lei & Zhou, 2020, for convergence in expectation of Online Mirror Descent). Here, we addressed the almost sure convergence under a constant step size, under the requirement that batch sizes $N_t$ increase at each iteration (see e.g. Ghadimi et al., 2016, Thm. 2, for a finite time convergence analysis in the setting $N_t \ne 1$ and $\gamma_t = \bar{\gamma}$).

**The limiting set in the general case.** We now consider the more general case when $\{\gamma_t, t \ge 1\}$ is a deterministic sequence, possibly time varying, taking values in $(0, \gamma_*)$. Here, $\mathbb{U} := \mathbb{T} \times (0, \gamma_*)$ and is not a compact set; hence, Theorem 9 does not apply but the more general result Theorem 5 remains valid. In this paragraph, we assume MD1, MD2, and MD4, together with:

**MD6.** For all $\gamma \in (0, \gamma_*)$, the function $b_\gamma := \gamma^{-1}\psi - F$ is strictly convex

**MD7.** $\psi$ is $\mu$-strongly convex on $\mathbb{T}$.
**Online Proximal-Gradient is a SAM2 algorithm.** Set $U := \mathbb{T} \times [\epsilon, 1/L]$ and for all $\theta \in \mathbb{T}$, set $U[\theta] := \{(\theta, \gamma), \gamma \in [\epsilon, 1/L]\}$. The function

$$G(\theta'; (\theta, \gamma)) := F_s(\theta) + \langle \nabla F_s(\theta), \theta' - \theta \rangle + \frac{1}{2\gamma} \|\theta' - \theta\|^2 + F_c(\theta') + \iota_{\mathbb{T}}(\theta'), \quad \theta' \in \mathbb{R}^d, \theta \in \mathbb{T}, \gamma \in [\epsilon, 1/L]$$

satisfies A1 under PG1. We also have $G(\theta'; (\theta, \gamma)) = \mathbb{E}[g(\theta', X; (\theta, \gamma))]$ where

$$g(\theta', x; (\theta, \gamma)) := F_s(\theta) + \langle h(\theta, x), \theta' - \theta \rangle + \frac{1}{2\gamma} \|\theta' - \theta\|^2 + F_c(\theta') + \iota_{\mathbb{T}}(\theta') \\ = F_c(\theta') + \frac{1}{2\gamma} \|\theta' - \theta + \gamma h(\theta, x)\|^2 + C(\gamma, \theta) + \iota_{\mathbb{T}}(\theta')$$

where $C(\gamma, \theta)$ does not depend on $\theta'$. Therefore, the minimization (5) is equivalent to the computation of the proximal-gradient operator (12), thus showing that Algorithm 4 is a SAM2 algorithm.

**The limiting set.** Under PG1 and

PG2. (a) $F_s$ is convex,

(b) for all $x \in \mathbb{X}$, $\theta \mapsto h(\theta, x)$ is continuous on $\mathbb{T}$ and $\mathbb{E}[\sup_{\theta \in \mathbb{T}} \|h(\theta, X)\|^p] < \infty$ for some $p > 1$

the conditions A1, A2 and A3 are satisfied. Since $\theta' \mapsto G(\theta'; (\theta, \gamma))$ is a (proper) lower semicontinuous strictly convex function possessing a minimizer $T(\theta, \gamma)$, then this minimizer is the unique solution of (see e.g. (Bauschke & Combettes, 2011, Theorem 16.3))

$$0 \in \partial \left[ \langle \nabla F_s(\theta), \cdot \rangle + \frac{1}{2\gamma} \| \cdot - \theta \|^2 + F_c(\cdot) + \iota_{\mathbb{T}}(\cdot) \right] (T(\theta, \gamma)).$$

By (Bauschke & Combettes, 2011, Corollary 16.48), $T(\theta, \gamma) = \theta$, is equivalent to $0 \in \nabla F_s(\theta) + \partial F_c(\theta) + \partial \iota_{\mathbb{T}}(\theta)$. Observe that this property does not depend on $\gamma$ which implies that $\mathcal{L}^+ = \mathcal{L}$. This property also yields

$$\mathcal{L} = \{\theta \in \mathbb{T} : 0 \in \partial [F + \iota_{\mathbb{T}}](\theta)\}.$$

Since $F$ is a lower semicontinuous proper convex function, then $\mathcal{L}$ is the set of the minimizers of $F$ on $\mathbb{T}$.

For the sequence $\{\theta^t, t \ge 0\}$ given by Algorithm 4, Theorem 9 shows that as soon as the size of the mini batches is chosen so that A4 holds: the sequence $\{F(\theta^t), t \ge 0\}$ converges to $\min_{\mathbb{T}} F$ and the sequence $\{\theta^t, t \ge 0\}$ converges to the set of the minimizers of $F$ on $\mathbb{T}$.

In the literature, the convergence analysis of perturbed Proximal-Gradient algorithms usually relies on the Siegmund-Robbins lemma (see (Robbins & Siegmund, 1971, Theorem 1) for a stochastic version assuming the errors are nonnegative; see also Atchadé et al. (2017) for an extension addressing the case of signed errors and deterministic perturbations; and Lai (1989) for nonnegative errors and possibly non converging cumulated errors). Different settings were considered. Let
Under MD6, we have $\mathbb{U}[\theta] = \{(\theta, \gamma), \gamma \in (0, \gamma_*)\}$; this property follows from: $F(\theta) = G(\theta; (\tau, \gamma))$ iff $D_{b,\gamma}(\theta; \tau) = 0$. Now, let us verify the assumptions of Theorem 5.

**Lemma 10.** Assume MD 1, MD 2, MD 6 and MD 7. Then $\mathcal{T}$ is continuous on $\mathbb{U}$. In addition, for all $\gamma \in (0, \gamma_*)$ and $\theta \in \mathbb{T}$,

$$F(\theta) - F(\mathcal{T}(\theta, \gamma)) \geq \frac{\mu}{2\gamma} \|\theta - \mathcal{T}(\theta, \gamma)\|^2. \quad (14)$$

The proof is postponed to Section 6.7. Lemma 10 establishes Item i of Theorem 5. Set

$$\mathcal{L}^+ := \bigcap_{q \in \mathbb{N}_*} \{\theta \in \mathbb{T} : \exists \gamma \in (0, \gamma_*], \|\mathcal{T}(\theta, \gamma) - \theta\|^2/\gamma < 1/q\},$$

and observe that $\mathcal{L} \subseteq \mathcal{L}^+$, where $\mathcal{L}$ is given by Equation (8). Upon noting that on the closed set $\{\theta \in \mathbb{T} : \forall \gamma \in (0, \gamma_*], \|\mathcal{T}(\theta, \gamma) - \theta\|^2/\gamma \geq 1/q\}$, we have $F(\theta) - F(\mathcal{T}(\theta, \gamma)) \geq \mu/(2q)$ for all $\gamma \in (0, \gamma_*)$, Item ii of Theorem 5 is verified. Finally, let us verify Item iii. By Proposition 7, we have to show $\lim_t G(\theta^{t+1}; (\theta^t, \gamma_{t+1})) - G(\mathcal{T}(\theta^t, \gamma_{t+1})); (\theta^t, \gamma_{t+1})) = 0$, almost surely. Following the same arguments as in the proof of Theorem 9 Item i, such a limit follows from $\lim_t \sup_{\theta \in \mathbb{T}} \|N_{t+1}^{-1} \sum_{i=1}^{N_{t+1}} h(\theta, X^{t+1,i}) - \nabla F(\theta)\| = 0$, with probability one. Integrability, continuity assumptions on $h$ (see MD 4), and increasing batch size $\{N_{t+1}, t \geq 0\}$ at a convenient rate (see A4 with $p$ given by MD 4), yield this limit (see Theorem 14). Full details are omitted.

Theorem 5 implies that the sequence $\{\theta^t, t \geq 0\}$ converges to a subset of $\mathcal{L}^+$, where $\mathcal{L}^+$ contains $\mathcal{L}$ (see Equation (8)) and the set of the minimizers of $F$ (see Theorem 4).

## 4.2 An Online Proximal-Gradient algorithm

The objective is to minimize a composite function $F := F_s + F_c$ on a compact convex subset $\mathbb{T}$ of $\mathbb{R}^d$, where the following conditions hold.

PG1. (a) The function $F_s: \mathbb{R}^d \to \mathbb{R}$ is continuously differentiable with $L$-Lipschitz gradient, there exists a measurable function $h: \mathbb{T} \times \mathbb{X} \to \mathbb{R}^d$ such that $\nabla F_s(\theta) = \mathbb{E}[h(\theta, X)]$, where $\mathbb{E}[\|h(\theta, X)\|] < \infty$ for all $\theta \in \mathbb{T}$, and exact computation of the expectation is not available.

(b) A stream of independent random variables with the same distribution as $X$ is available, and $h(\theta, x)$ can be computed for all $(\theta, x) \in \mathbb{T} \times \mathbb{X}$.

(c) The function $F_c: \mathbb{R}^d \to (-\infty, +\infty]$ is a continuous convex function, $\mathbb{T} \subset \text{dom}(F_c)$ and for all $\gamma > 0$, the proximal map of $\gamma F_c + \iota_{\mathbb{T}}$ is a point-to-point map with an explicit expression.

Recall that $\iota_{\mathbb{T}}$ denotes the characteristic function of the set $\mathbb{T}$ and that the proximal map is given by

$$\text{Prox}_{\gamma F_c + \iota_{\mathbb{T}}}(\theta) := \text{Argmin}_{\tau \in \mathbb{T}} \left( \gamma F_c(\tau) + \frac{1}{2} \|\tau - \theta\|^2 \right).$$
us cite as few examples, the finite-sum setting Nitanda (2014); Fort & Moulines (2023); the convex case under summability assumptions on the errors $\mathcal{H}^{t+1} - \nabla F_s(\theta^t)$ which reveal restrictive for the i.i.d. streaming framework considered here Combettes & Pesquet (2015); the convex case under possibly biased errors Atchadé et al. (2017); the convex case when $F_s$ is uniformly convex at a minimizer and the conditional variance of the error is controled by $\|\nabla F_s(\theta^t)\|^2$ and does not vanish Rosasco et al. (2020). The closest result to the one we obtain is (Atchadé et al., 2017, Theorem 6), which covers Algorithm 4: applied with a lower bounded stepsize sequence $\{\gamma_t, t \ge 1\}$, it provides almost-sure convergence of the SAM2 sequence to a point in $\mathcal{L}$ under the assumption $\sum_t N_t^{-1} < \infty$.

# 5 A numerical illustration: Regression under a quantile loss function

The effectiveness of the SAM2 algorithm is demonstrated in the following simulation analysis. We consider a linear model with heavy-tailed noise, for which most existing theories and implementation are no longer applicable. In contrast, we show that SAM2 applies straightforwardly with theoretical guarantees. More specifically, we consider the following linear model,

$$Y = \langle \theta, \overline{W} \rangle + \epsilon \quad (13)$$

where $\theta \in \mathbb{R}^l$ is the unknown regression parameter, $\overline{W} := (1, W)$ and $W \in \mathbb{R}^{l-1}$ are the covariates and $Y \in \mathbb{R}$ is the response variable. The random variables $W$ and $\epsilon$ are defined in $(\Omega, \mathcal{A}, \mathbb{P})$; the additive noise $\epsilon$ is assumed independent of $W$ and $W$ is integrable $\mathbb{E}[\|W\|] < \infty$.

The goal is to learn the regression parameter $\theta$ from streaming data, independent and with the same distribution as $X := (W, Y)$. We are interested in very heavy-tailed noise $\epsilon$, typically with infinite variance, for which the standard squared-loss approach cannot be applied: the quantile regression (QR) loss is commonly used instead. Let $q \in (0, 1)$ be a specified quantile level, the so-called QR problem consists of finding the solution of the following minimization problem

$$\arg\min_{\theta \in \mathbb{R}^l} \mathbb{E} \left[ \rho_q(Y - \langle \theta, \overline{W} \rangle) - \rho_q(Y) \right], \quad \rho_q(u) := (q - 1)_{u<0} u. \quad (14)$$

This objective function is finite for all $\theta \in \mathbb{R}^l$ under the integrability condition on $W$ (see e.g. Lemma 17). Note that the addition of the second term in (14) does not impact the solution of the minimization problem but guarantees the existence of the loss and allows to consider more general settings. In this section, we consider the possibly penalized QR problem and aim to solve:

$$\arg\min_{\theta} F_{\eta}(\theta), \quad F_{\eta}(\theta) := \mathbb{E} \left[ \rho_q(Y - \langle \theta, \overline{W} \rangle) - \rho_q(Y) \right] + \eta \|\theta\|_1, \quad (15)$$

where $\|\theta\|_1$ is the sparsity inducing $L_1$-norm of $\theta$ and $\eta \ge 0$. When $q = 1/2$, we recover the so-called
The Online Proximal-Gradient algorithm, given by Algorithm 4, is a stochastic perturbation of a forward-backward splitting process in optimization (see Lions & Mercier, 1979, and also Beck & Teboulle, 2009, Sec. 1.3.3, for an introduction as an MM algorithm) in which the forward operator is not explicit and is learnt from a stream of oracles (see e.g. Combettes & Wajs, 2005, Sec. 3, for a pioneering contribution on inexact proximal-gradient algorithms).

**Algorithm 4** An Online Proximal-Gradient algorithm

**Require:** a sequence {$N_t, t \ge 1$} of positive integers, an initial value $\theta^0 \in \mathbb{T}$, $\epsilon > 0$.

**Ensure:** A $\mathbb{T}$-valued sequence {$\theta^t, t \ge 0$}.

1: **for** $t \ge 0$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;Sample a batch {$X^{t+1,i}, 1 \le i \le N_{t+1}$} of size $N_{t+1}$.

3: &nbsp;&nbsp;&nbsp;&nbsp;Set $\mathcal{H}^{t+1} := N_{t+1}^{-1} \sum_{i=1}^{N_{t+1}} h(\theta^t, X^{t+1,i})$.

4: &nbsp;&nbsp;&nbsp;&nbsp;Choose $\gamma_{t+1} \in (0, 1/L]$.

5: &nbsp;&nbsp;&nbsp;&nbsp;Compute

$$ \theta^{t+1} = \text{Prox}_{\gamma_{t+1}F_c + \iota_{\mathbb{T}}} (\theta^t - \gamma_{t+1}\mathcal{H}^{t+1}). \quad (15) $$

6: **end for**

**Algorithm 4 is a SAM2 algorithm.** Set $\mathcal{U} := \mathbb{T} \times (0, 1/L]$ and for all $\theta \in \mathbb{T}$, set $\mathcal{U}[\theta] := \{(\theta, \gamma), \gamma \in (0, 1/L]\}$. The function

$$ \mathcal{G}(\theta'; (\theta, \gamma)) := F_s(\theta) + \langle \nabla F_s(\theta), \theta' - \theta \rangle + \frac{1}{2\gamma} \|\theta' - \theta\|^2 + F_c(\theta') + \iota_{\mathbb{T}}(\theta'), \quad \theta' \in \mathbb{R}^d, \theta \in \mathbb{T}, \gamma \in (0, 1/L], $$

satisfies A1 under PG1. We also have $\mathcal{G}(\theta'; (\theta, \gamma)) = \mathbb{E}[g(\theta', X; (\theta, \gamma))],$ where

$$ \begin{aligned} g(\theta', x; (\theta, \gamma)) &:= F_s(\theta) + \langle h(\theta, x), \theta' - \theta \rangle + \frac{1}{2\gamma} \|\theta' - \theta\|^2 + F_c(\theta') + \iota_{\mathbb{T}}(\theta') \\ &= F_c(\theta') + \frac{1}{2\gamma} \|\theta' - \theta + \gamma h(\theta, x)\|^2 + \mathcal{C}(\gamma, \theta) + \iota_{\mathbb{T}}(\theta') \end{aligned} $$

and $\mathcal{C}(\gamma, \theta)$ does not depend on $\theta'$. Therefore, the minimization (5) is equivalent to the computation of the proximal-gradient operator (15), thus showing that Algorithm 4 is a SAM2 algorithm.

**The limiting set in the case** $\mathcal{U} := \mathbb{T} \times [\epsilon, 1/L]$ **for** $\epsilon > 0$. Here, for simplicity, we consider the case the set $\mathcal{U}$ is a compact interval; under this assumption, $\mathcal{U}[\theta] = \{(\theta, \gamma), \gamma \in [\epsilon, 1/L]\}$. The general case when $\mathcal{U} = \mathbb{T} \times (0, 1/L]$ can be addressed, following the same lines as in Section 4.1. Details are left to the interested reader.

Make the following additional assumption.

PG2. (a) $F_s$ is convex, and

(b) for all $x \in \mathbb{X}$, $\theta \mapsto h(\theta, x)$ is continuous on $\mathbb{T}$ and $\mathbb{E}[\sup_{\theta \in \mathbb{T}} \|h(\theta, X)\|^p] < \infty$ for some $p > 1$.
Least Absolute Deviation (LAD) problem. As the QR loss is non smooth, a number of standard algorithms cannot be applied. More specifically, the two main sources of complication in QR are the heavy-tail of the noise and the non-smoothness of the loss. The sometimes called check function $\rho_q$ is non-differentiable at 0. When the noise is not of finite variance, most least-square-based theories or Huber-based robust approaches are not applicable. To handle non-smoothness, most approaches are two-step approaches with first the design of a smooth approximation of the QR loss and then the optimization of this approximation, with the hope that the distance between the estimates and the true minimizers can be controlled and tends to zero. Various smooth approximations have been proposed, using a perturbation of the check function that can be majorized by a quadratic function (Hunter & Lange, 2000), kernel convolution-based smoothing (Jiang & Yu, 2022; Chen et al., 2019), a huberized pinball loss for robust to outliers but essentially Gaussian noise (Ichinose et al., 2023), transformation of the QR loss into a least-square loss on new response variables (Chen et al., 2020), etc. The type of approximations designed generally guides in turn the subsequent optimization technique. For instance Hunter & Lange (2000) use an MM algorithm, while Chen et al. (2019, 2020); Ichinose et al. (2023); Zheng (2011) exploit the differentiability of their approximation to derive stochastic gradient or stochastic Newton-Raphson algorithms. The position of SAM2 appears quite unique in the QR literature, as SAM2 does not require to change the target loss to a smooth approximation while performing the loss optimization with an MM framework. The following proposition shows that we can indeed formulate the QR problem so as to exhibit a convenient majorizer, which is both easy to optimize and satisfies the SAM2 requirements under mild conditions. Note however, that the result below only holds for linear regression functions while some of the previous mentioned methods, e.g. (Hunter & Lange, 2000; Ichinose et al., 2023) can also deal with non-linear regressions. Regarding the addition of a Lasso penalty for variable selection, most approaches including SAM2 handle this extension straightforwardly.

**Proposition 10.** Assume there exists $p_* > 1$ such that $\mathbb{E}[||W||^{p_*}] < \infty$. Then A 1 and A 2-Item a are satisfied with

$$g_{\eta}(\theta, X; \tau) := \frac{1}{\ell} \sum_{j=1}^{\ell} \left( \rho_q(Y - \langle \tau, \overline{W} \rangle + \ell \overline{W}_j \tau_j - \ell \overline{W}_j \theta_j) - \rho_q(Y) + \eta \ell |\theta_j| \right),$$

and $U[\theta] = \{\theta\}$ for all $\theta$. In addition, A 3 is satisfied with $p = p_*$. The optimization step (5) is explicit for all $q \in (0, 1)$ and $\eta \ge 0$.

The proof is in Section 6.7. Whatever the dimension $\ell$ of $\theta$, $g_{\eta}$ is a sum of functions that involve each $\theta_i$ separately and can then be separately minimized over a one-dimensional space. A $\ell$-dimensional problem is then turned into $\ell$ 1-dimensional problems. In Section 6.7, we show that such 1-dimensional minimization can be solved exactly.
Under PG1 and PG2, the conditions A1, A2 and A3 are satisfied. Since $\theta' \mapsto G(\theta'; (\theta, \gamma))$ is a (proper) lower semi-continuous strictly convex function possessing a minimizer $T(\theta, \gamma)$, the minimizer is the unique solution of (see e.g. Bauschke & Combettes, 2011, Thm. 16.3)

$$0 \in \partial \left[ \langle \nabla F_s(\theta), \cdot \rangle + \frac{1}{2\gamma} \| \cdot - \theta \|^2 + F_c(\cdot) + \iota_{\mathbb{T}}(\cdot) \right] (T(\theta, \gamma)).$$

By Bauschke & Combettes (2011, Cor. 16.48), $T(\theta, \gamma) = \theta$ is equivalent to $0 \in \nabla F_s(\theta) + \partial F_c(\theta) + \partial \iota_{\mathbb{T}}(\theta)$. Observe that this property does not depend on $\gamma$, which implies that $\mathcal{L}^+ = \mathcal{L}$. This property also yields that

$$\mathcal{L} = \{ \theta \in \mathbb{T} : 0 \in \partial [F + \iota_{\mathbb{T}}](\theta) \}.$$

Since $F$ is a lower semi-continuous proper convex function, $\mathcal{L}$ is the set of the minimizers of $F$ on $\mathbb{T}$.

For the sequence $\{\theta^t, t \ge 0\}$ given by Algorithm 4, Theorem 9 shows that if the sizes of the batches are chosen so that A4 holds, the sequence $\{F(\theta^t), t \ge 0\}$ converges to $\min_{\mathbb{T}} F$ and the sequence $\{\theta^t, t \ge 0\}$ converges to the set of the minimizers of $F$ on $\mathbb{T}$.

**Remark.** In the literature, the convergence analysis of perturbed Proximal-Gradient algorithms usually relies on the Siegmund-Robbins lemma (see Robbins & Siegmund, 1971, Thm. 1, for a stochastic version assuming the errors are non-negative, Atchadé et al., 2017, for an extension addressing the case of signed errors and deterministic perturbations; and Lai, 1989, for non-negative errors and possibly non-converging cumulated errors). Some examples of analyses include the finite-sum setting (Nitanda, 2014; Fort & Moulines, 2023); the convex setting under summability assumptions on the errors $\mathcal{H}^{t+1} - \nabla F_s(\theta^t)$ (Combettes & Pesquet, 2015), which are too restrictive for the i.i.d. streaming framework considered in this work; the convex setting under possibly biased errors (Atchadé et al., 2017); and the convex setting when $F_s$ is uniformly convex at a minimizer and the conditional variance of the error is non-vanishing and controlled by $\|\nabla F_s(\theta^t)\|^2$ (Rosasco et al., 2020). The closest result to the one we obtain under the assumption that $\gamma_t \in [\varepsilon, 1/L]$, is Atchadé et al. (2017, Thm. 6), which covers Algorithm 4. Applied with a lower bounded step size sequence $\{\gamma_t, t \ge 1\}$, it provides almost sure convergence of the SAM2 sequence to a point in $\mathcal{L}$ under the assumption $\sum_t N_t^{-1} < \infty$.

# 5 Illustration: Regression under a quantile loss function

The effectiveness of the SAM2 algorithm is demonstrated in the following simulation analysis. We consider a linear model with heavy-tailed noise, for which most existing theories and implementation are no longer applicable. In contrast, we show how SAM2 applies straightforwardly. More
## 5.1 Simulation Setup

Let $X^i := (W^i, Y^i) \in \mathbb{R}^{l-1} \times \mathbb{R}$ for each $i \in [N]$ be the realization of an i.i.d. sample of $(W, Y)$. We consider the following linear model,

$$Y = \langle \theta^{\text{true}}, \bar{W} \rangle + \epsilon \quad (16)$$

where $\bar{W} := (1, W_1, \dots, W_{l-1})$ is a $l$-dimensional covariate vector. The examples $\{W^i, i \in [N]\}$ are sampled independently, from a multivariate normal distribution $N(0, \Sigma)$ with a covariance matrix $\Sigma$ whose entries are $\Sigma_{r,s} := r^{|r-s|}$ for $r, s \in [l-1]$. We set $l = 11$ and $r = 0.9$. The true coefficient is set to

$$\theta^{\text{true}} := 10 \left( \frac{1}{l}, \frac{2}{l}, \dots, \frac{l-1}{l}, 1 \right).$$

The noise distribution is chosen to be Cauchy $\epsilon \sim \text{Cauchy}(0, 1)$, and the random variables $\epsilon^i$ are independent.

The unknown parameter $\theta$ is learnt from the examples $(X^i, Y^i)$, as a minimizer of the criterion (15) in the case $\eta = 0$ and $q = 0.5$. For SAM2 the sequence $\{N_t, t \ge 1\}$ is chosen so as to satisfy assumption A4. For a Gaussian vector $W$, the assumption A3 is satisfied for any $p > 1$ so that the theory claims that a sufficient condition for the convergence of SAM2 is the $N_t$'s increase as slowly as desired. However, we illustrate below that too small constant values are not recommended. Thus, we first illustrate SAM2 with a linearly growing modified so as to start with batch sizes of 100: $N_t = \max(100, t)$. As will be illustrated in Section 5.3, this burnin has no real impact on the algorithm convergence but starting with batch sizes of 100 provides better starts. The number of iterations is set to $T = 1000$. It follows that the total number of data points processed is $N = \sum_{t=1}^T N_t = 505,450$.

As the QR loss is non-smooth, for comparison, we implement a stochastic subgradient (SSG) algorithm using a subgradient $g_q(\theta; (w, y))$ of $\rho_q(y - \langle \theta, \bar{w} \rangle)$ with respect to $\theta$ which is given by

$$g_q(\theta; (w, y)) := -(q - 1_{y-\langle\theta, \bar{w}\rangle<0}) \bar{w}; \quad \bar{w} := (1, w).$$

The SSG algorithm is essentially the subgradient algorithm, but using noisy subgradients and a sequence of positive step sizes $\{\gamma_t, t \ge 1\}$. It uses the standard update specified in Algorithm 5, written with the possibility to use mini batches also denoted by $\{N_t, t \ge 1\}$.

SSG requires the choice of step sizes and mini-batch sizes, while SAM2 requires only to choose a sequence of mini-batch sizes $\{N_t, t \ge 1\}$.

For all algorithms, the iterate is arbitrarily initialized to $\theta^0 := (1, \dots, 1)$ and a Polyak averaging is also performed starting at $T_0 = 500$, for $t > T_0$: $\bar{\theta}^t = \frac{t-1-T_0}{t-T_0}\bar{\theta}^{t-1} + \frac{1}{t-T_0}\theta^t$. For a quantitative performance assessment, we compute the root mean squared error (RMSE) or L2-error at each
specifically, we consider the following linear model,

$$Y = \langle \theta, \overline{W} \rangle + \epsilon, \qquad (16)$$

where $\theta \in \mathbb{R}^l$ is the unknown regression parameter, $\overline{W} := (1, W)$ and $W \in \mathbb{R}^{l-1}$ are the covariates and $Y \in \mathbb{R}$ is the response variable. The random variables $W$ and $\epsilon$ are defined in $(\Omega, \mathcal{A}, \mathbb{P})$; the additive noise $\epsilon$ is assumed independent of $W$ and $W$ is integrable $\mathbb{E}[\|W\|] < \infty$. Recall that $\| \cdot \|$ denotes the $l_2$-norm. Throughout this section, $\| \cdot \|_r$ is the $l_r$-norm, $r \ge 1$.

The goal is to learn the regression parameter $\theta$ from streaming data, that are independent and with the same distribution as $X := (W, Y)$. We are interested in very heavy-tailed noise $\epsilon$, typically with infinite variance, for which the standard squared-loss approach cannot be applied: the quantile regression (QR) loss is commonly used instead. Let $q \in (0, 1)$ be a specified quantile level, the so-called QR problem consists of finding the solution of the following minimization problem

$$\arg \min_{\theta \in \mathbb{R}^l} \mathbb{E} [\rho_q(Y - \langle \theta, \overline{W} \rangle) - \rho_q(Y)], \quad \rho_q(u) := (q - 1_{u<0})u. \qquad (17)$$

This objective function is finite for all $\theta \in \mathbb{R}^l$ under the integrability condition on $W$ (see e.g., Lemma 18). Note that the addition of the second term in (17) does not impact the solution of the minimization problem but guarantees the existence of the loss and allows us to consider more general settings. In this section, we consider the possibly penalized QR problem and aim to solve:

$$\arg \min_{\theta} \mathbf{F}_{\eta}(\theta), \quad \mathbf{F}_{\eta}(\theta) := \mathbb{E} [\rho_q(Y - \langle \theta, \overline{W} \rangle) - \rho_q(Y)] + \eta \mathbf{pen}(\theta), \qquad (18)$$

where $\eta \ge 0$ and $\mathbf{pen}(\theta) := \sum_{j=1}^l \varphi(|\theta_j|)$ is a penalty that regularizes $\theta = (\theta_1, \dots, \theta_l)$, based on a positive map $\varphi: \mathbb{R}_{\ge 0} \to \mathbb{R}_{\ge 0}$, such that $\varphi(0) = 0$. In the case of $\varphi(a) = a$ for each $a \ge 0$, we have the usual $l_1$ penalty: $\mathbf{pen}(\theta) = \|\theta\|_1$. When $q = 1/2$, we recover the so-called Least Absolute Deviation (LAD) problem. As the QR loss is non-smooth, many standard algorithms cannot be applied. Specifically, the two main sources of complication in QR are the heavy tail of the noise and the non-smoothness of the loss. Sometimes called the check function, $\rho_q$ is non-differentiable at 0. When the noise variance is not finite, most least-square-based theories or Huber-based robust approaches are not applicable. To handle non-smoothness, most approaches are two-step approaches that first design a smooth approximation of the QR loss and then conduct the optimization of this approximation, with the hope that the distance between the estimates and the true minimizers can be controlled and tends to zero. Various smooth approximations have been proposed, using a perturbation of the check function, that can be majorized by a quadratic function (Hunter & Lange, 2000), kernel convolution-based smoothing (Jiang & Yu, 2022; Chen et al., 2019), a Huberized pinball loss for robust-to-outliers but essentially Gaussian noise (Ichinose et al., 2023), transformation of the QR loss into a least-square loss on new response variables (Chen et al., 2020),
**Algorithm 5 The SSG algorithm**

**Require:** a sequence {$N_t, t \ge 1$} of positive integers and of positive step sizes {$\gamma_t, t \ge 1$}, an initial value $\theta^0 \in \mathbb{T}$

**Ensure:** A $\mathbb{T}$-valued sequence {$\theta^t, t \ge 0$}.

1: **for** $t = 0, \dots, \mathbf{do}$

2: Sample a minibatch {$X^{t+1,i}, i = 1, \dots, N_{t+1}$} of size $N_{t+1}$.

3: Compute

$$ \theta^{t+1} = \theta^t - \gamma_{t+1} \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g_q(\theta; X^{t+1,i}) \quad (17) $$

4: **end for**

iteration:

$$ \text{RMSE}_t := ||\theta^t - \theta^{\text{true}}||. $$

SAM2 is then compared to SSG using various settings in order to assess their respective robustness to the step size and batch size choices.

## 5.2 Sensitivity to step sizes

Regarding the SSG algorithm, we first consider a standard setting with mini-batches of equal size. To compare both algorithms using the same amount of observations and perform the same number of parameter updates, the batch size is set to $\lfloor N/T \rfloor = 505$ for all iterations except the last one which is augmented with the remaining data points and then includes 955 samples. This later adjustment has very little impact. SSG requires the user to decide on a step size sequence, which has been reported to be impactful and the main practical limitation of gradient algorithms. Three such sequences are thus tested: $\gamma_t = \frac{1}{(t+1)^{0.51}}$, $\gamma_t = \frac{1}{(t+1)^{0.6}}$ and $\gamma_t = \frac{1}{(t+1)^{0.7}}$.

Figure 1 shows, for one simulated data set, the sequences (left) of RMSE obtained with SAM2 and the three SSG tested settings, with a zoom on the last 500 iterations (right). Two representations of the RMSE evolution are reported in Figure 1, one with respect to the number of parameter updates (iterations, in the first line), and the other with respect to the number of processed samples (second line). The lowest RMSE are obtained for SAM2 and SSG with step sizes $\gamma_t = \frac{1}{(t+1)^{0.51}}$, with SAM2 outperforming the later. For $\gamma_t = \frac{1}{(t+1)^{0.7}}$, the SSG blue curve is not visible on the zoomed plot as its RMSE are too large. Polyak averaging performed with a burnin, starting at $T_0 = 500$, is beneficial to the SAM2 RMSE but not for the SSG RMSE. In terms of observed RMSE, the Polyak averaging performance decreases as the burnin phase $T_0$ becomes smaller. Note that the RSME curve of the Polyak averaging sequence is not the Polyak averaging of the RMSE sequence. Overall, although noisier, the SAM2 estimations appear less biased than the SSG ones. This is illustrated on some parameter sequences in Figure 2. The plots illustrate that SAM2 estimates are in general closer to the true values. In addition, the simulations confirm the often reported fact that SSG performance can be quite sensitive to the step size choice, with 2 out of 3 settings showing significantly larger
etc. The type of approximations designed generally guides in turn the subsequent optimization technique. For instance Hunter & Lange (2000) use an MM algorithm, while Chen et al. (2019), Chen et al. (2020), Ichinose et al. (2023), and Zheng (2011) each exploit the differentiability of their approximation to derive stochastic gradient or stochastic Newton-Raphson algorithms. The position of SAM2 appears quite unique in the QR literature, as SAM2 does not require to change the target loss to a smooth approximation while performing the loss optimization with an MM framework. The following proposition shows that we can indeed formulate the QR problem so as to exhibit a convenient majorizer, which is both easy to optimize and satisfies the SAM2 requirements under mild conditions. Note however, that the result below only holds for linear regression functions, while some of the previous mentioned methods (e.g., Hunter & Lange, 2000; Ichinose et al., 2023) can also deal with non-linear regressions. Regarding the addition of a Lasso ($l_1$) penalty for variable selection or sparse regression problems, most approaches, including SAM2, handle this extension straightforwardly but SAM2 can additionally handle non-convex sparsity-inducing norms that can be more efficient than the $l_1$-norm in terms of support recovery, see Section 5.3 below.

**Proposition 11.** Assume that there exists $p_* > 1$ such that $\mathbb{E}[\|W\|^{p_*}] < \infty$. Then A 1 and A 2-Item a are satisfied with

$$g_{\eta}(\theta, X; \tau) := \frac{1}{\ell} \sum_{j=1}^{\ell} \left( \rho_q(Y - \langle \tau, \overline{W} \rangle + \ell \overline{W}_j \tau_j - \ell \overline{W}_j \theta_j) - \rho_q(Y) + \eta \ell \varphi(|\theta_j|) \right),$$

and $\theta \in \mathbb{U}[\theta]$ for all $\theta$. In addition, A 3 is satisfied with $p = p_*$. Then, if $\varphi$ is continuous, concave, and strictly increasing, the optimization step (5) is explicit and can be solved exactly for all $q \in (0, 1)$ and $\eta \ge 0$.

The proof of the result appears in Section 6.8. Whatever the dimension $\ell$ of $\theta$, $g_{\eta}$ is a sum of functions that involve each $\theta_j$, separately, and can then be separately minimized over a one-dimensional space. A $\ell$-dimensional optimization problem is thus turned into $\ell$ 1-dimensional optimization problems. In Section 6.8, we show that such 1-dimensional minimization can be solved exactly.

## 5.1 Simulation Setup

Let $X^i := (W^i, Y^i) \in \mathbb{R}^{\ell-1} \times \mathbb{R}$ for each $i \in [N]$ be the realization of an i.i.d. sample of $(W, Y)$. We consider the following linear model,

$$Y = \langle \theta^{\text{true}}, \overline{W} \rangle + \epsilon \quad (19)$$

where $\overline{W} := (1, W_1, \dots, W_{\ell-1})$ is a $\ell$-dimensional covariate vector. The examples $\{W^i, i \in [N]\}$ are sampled independently, from a multivariate normal distribution $N(0, \Sigma)$ with a covariance matrix
<figure><img src="image_2.png" /><figcaption>Figure 1: LAD example. First line: RMSE over iterations for SAM2 (black), SSG with different step size schedules, $\gamma_t = \frac{1}{(t+1)^{0.51}}$ (red and purple), $\gamma_t = \frac{1}{(t+1)^{0.6}}$ (green) and $\gamma_t = \frac{1}{(t+1)^{0.7}}$ (blue). All SSG sequences are with constant batch sizes except the purple curve which uses the same batches as SAM2. The left plot shows a zoom on the last half iterations, with additional RMSE curves for the Polyak averaging sequences (dashed lines). Second line: same comparison of RMSE but over number of processed samples.</figcaption></figure>

RMSE. In contrast, SAM2 does not require such tuning. It is quite robust to the choice of the $N_t$'s as soon as they are not kept too small for too many iterations. Typically, when $N_t = t$, which lowers the first batch sizes, SAM2 performs similarly as in Figure 1 but with a slower start. We further investigate the sensitivity to batch sizes in the next section.

## 5.3 Sensitivity to batch sizes

First, for SSG a second strategy is considered for the mini-batch size with a varying size $N_t$ so that at each iteration the SSG algorithm uses the exact same data points as SAM2. This setting is tested only for SSG with $\gamma_t = \frac{1}{(t+1)^{0.51}}$, which provided before the best results. For this sequence of step sizes, the purple curve in Figure 1 shows RMSE obtained when using the same mini-batch sizes and data points for both SAM2 and SSG. The stronger impact is seen on the second line of Figure 1. SAM2 and SSG with linearly increasing batch sizes, show better start reaching lower RMSE faster. This is consistent with the intuition that more frequent parameter updating is preferable than finer optimization steps at the early stages of the algorithms. Overall, final SSG RMSE do not seem to be much impacted by this change of batch sizes. For this specific simulation, the constant batch
Σ whose entries are $\Sigma_{k,k'} := c^{|k-k'|}$ for $k, k' \in [\ell - 1]$. The true parameter is set to

$$ \theta^{\text{true}} := 10 \left( \frac{1}{s}, \frac{2}{s}, \dots, \frac{s-1}{s}, 1, 0, 0, \dots, 0 \right), \qquad (20) $$

where $s \in [\ell]$ denotes the sparsity level, here the number of non-zero coefficients. We consider two cases with $\ell = 11$ and $c = 0.9$:

* a non-sparse case: where $s = \ell = 11$,

* a sparse case: where $s = 5$, which corresponds to 5 non-zero coefficients and 6 zero coefficients.

The distribution of the noise $\epsilon$ is chosen to be Cauchy(0, 1), and the random variables $\{\epsilon^i, i \in [N]\}$ are independent of each other and of the random variables $\{W^i, i \in [N]\}$.

The unknown parameter $\theta^{\text{true}}$ is learned from the examples $(X^i, Y^i)$, as a minimizer of criterion (18), in the case $q = 0.5$. We will consider both the case when $\eta = 0$ (Section 5.2) and the case when $\eta > 0$ (Section 5.3).

**The SAM2 algorithm.** Throughout Section 5, we refer to the procedure in Algorithm 2, using the function $g(\cdot) := g_\eta(\cdot)$ given by Proposition 11, as **SAM2**. The sequence $\{N_t, t \ge 1\}$ is chosen so as to satisfy assumption A4: for a Gaussian vector $W$, the assumption A3 is satisfied for any $p > 1$ and thus our theory states that a sufficient condition for the convergence of **SAM2** is for the $N_t$'s to increase as slowly as desired. Thus, except when specified, **SAM2** is applied with a linearly growing sequence $\{N_t, t \ge 1\}$ modified so as to start with a minimal size of 100: $N_t = \max(100, t)$. Throughout this section, this choice is referred to as **linearly increasing batch size**. Since $p$ in A3 can be chosen arbitrarily large, we will also numerically explore how **SAM2** behaves when $\{N_t, t \ge 1\}$ is constant.

**The Stochastic Subgradient algorithm.** As the QR loss is non-smooth, for comparison, we implement a stochastic subgradient (SSG) algorithm using a subgradient $g_q(\theta; (w, y))$ of $\rho_q(y - \langle \theta, \bar{w} \rangle)$ with respect to $\theta$, which is given by

$$ g_q(\theta; (w, y)) := -(q - 1_{y-\langle \theta, \bar{w} \rangle < 0}) \bar{w}; \quad \bar{w} := (1, w). \qquad (21) $$

The SSG algorithm is essentially the subgradient algorithm, but uses noisy subgradients and a sequence of positive step sizes $\{\gamma_t, t \ge 1\}$. It uses the standard update specified in Algorithm 5, presented with the possibility to use batches with size $N_t$ at iteration $\#t$. We will run SSG for different decreasing step sizes $\{\gamma_t, t \ge 1\}$ and with different batch sizes $\{N_t, t \ge 1\}$.
<figure><img src="image_3.png" /><figcaption>Figure 2: LAD example. Sequences of estimates for components $\theta_1, \theta_2, \theta_3$ and $\theta_{11}$, using SAM2 with linear $N_t = \max(100, t)$ (black) and SSG with $\gamma_t = \frac{1}{(t+1)^{0.51}}$ and constant batch sizes $N_t = 505$ (red). Polyak averaging sequences are shown with dashed lines. True parameter values are indicated by the blue horizontal dashed lines.</figcaption></figure>
**Algorithm 5 The SSG algorithm**

**Require:** a sequence {$N_t, t \ge 1$} of positive integers and of positive step sizes {$\gamma_t, t \ge 1$}, an initial value $\theta^0 \in \mathbb{T}$

**Ensure:** A $\mathbb{T}$-valued sequence {$\theta^t, t \ge 0$}.

1: **for** $t \ge 0$ **do**

2: Sample a batch {$X^{t+1,i}, 1 \le i \le N_{t+1}$} of size $N_{t+1}$.

3: Compute

$$ \theta^{t+1} = \theta^t - \gamma_{t+1} \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g_q(\theta^t; X^{t+1,i}) \quad (22) $$

4: **end for**

**The Stochastic Mirror Descent algorithm.** Additionally, we also run an online version of the stochastic mirror descent (SMD) algorithm, as given by Algorithm 3. Mirror map potentials $\psi$, also referred to as distance generating functions, are restricted to strictly convex functions: when considering sparse problems, the $\ell_1$-norm cannot then be used. As per applications to similar sparse estimation problems (see, e.g., Azizan & Hassibi, 2019 and Ilandarideva et al., 2023), we set the mirror map potential to $\psi(\cdot) := \|\cdot\|_r^r/r$ with $1 < r < 2$. The motivation for such a potential is also due to its strict convexity, and the fact that it approximates the $\ell_1$-norm as $r \to 1$. Furthermore, when $\ell \ge 3$ (in our simulations, $\ell = 11$), if $r = 1+1/\log \ell$ (i.e. $r \approx 1.42$ when $\ell = 11$), the mirror map potential can be shown to be strongly convex on any compact set (Ilandarideva et al., 2023, Sec. 3.2). In Azizan & Hassibi (2019), it was demonstrated that SMD with this mirror map potential was more computationally efficient at sparse recovery than standard $\ell_1$ optimization.

The gradient of the map $\psi$ and the inverse gradient can be computed almost everywhere and are given by

$$ \nabla\psi(\theta) = \text{sign}(\theta) |\theta|^{r-1} \quad \text{and} \quad [\nabla\psi]^{-1}(u) = \text{sign}(u) |u|^{\frac{1}{r-1}}. $$

It follows that, at iteration $\#t$, SMD simplifies to the following steps (see (21) for the definition of $g_q$):

$$ \tilde{\theta}^{t+1} = \nabla\psi(\theta^t) - \gamma_{t+1} \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g_q(\theta^t; X^{t+1,i}) \quad \text{and} \quad \theta^{t+1} = [\nabla\psi]^{-1}(\tilde{\theta}^{t+1}). $$

**Description of the analyses.** SAM2 is compared to SSG and SMD in various settings, to assess the respective robustness of the algorithms to the step size and batch size choices.

For all algorithms, the iterate is arbitrarily initialized to $\theta^0 := (1, \dots, 1)$. We note that SSG and SMD both require the choice of step sizes in addition to the batch sizes, while SAM2 requires only to choose a sequence of batch sizes {$N_t, t \ge 1$}. Different strategies for step sizes and batch sizes will be compared: decreasing step sizes of the form $\gamma_t = 1/(t+1)^a$; and constant batch size or
<figure><img src="image_4.png" /><figcaption>Figure 3: LAD example. First line: RMSE over iterations for SAM2 with different batch size sequences, original SAM2 (black) and constant (505) batch size SAM2 (blue). The left plot shows a zoom on the last half iterations, with additional RMSE curves for the Polyak averaging sequences (dashed lines). Second line: same comparison but with respect to the number of samples.</figcaption></figure>

size setting slightly outperforms eventually the varying one but boxplots in Figure 5 show than on average both setting lead to very similar performance.

Symmetrically, Figure 3 shows then the result of two SAM2 runs, where the previous setting is compared with one where the batch sizes are set to a constant $\lfloor N/T \rfloor = 505$, as for the previous SSG. Although this does not satisfy condition A4 on the $N_t$'s, this setting provides better RMSE in the first iterations but is eventually equivalent as the batch sizes increase. The Polyak averaging sequences suggest that the constant batch size version of SAM2 will eventually be outperformed by the increasing size one, as expected. Figure 4 further shows SAM2 sequences with constant batch sizes successively set to constant between 1 and 10. $N_t$ should not be set too low. For $N_t = 1$, the algorithm diverges very quickly. For $N_t = 5$, RMSE values remain very high and noisy, while $N_t = 10$ still provides much larger RMSE than $N_t = 505$.

This first set of simulations suggests that SAM2 has the advantage, over its stochastic gradient SSG counterpart, not to require the choice of a step size sequence, while choosing batch sizes is not problematic. One could argue that the need to use increasing $N_t$'s could be limiting, typically if one does not control the data stream. But note that this is only a sufficient condition and that SAM2 seems to perform well for constant $N_t$'s too especially when combined with Polyak averaging. For SSG the sensitivity to the step size appears not to be as easily compensated by appropriate
linearly increasing batch size $N_t = \max(100, t)$. In each analysis, sample paths of the algorithms are simulated for $T = 1000$ iterations and it follows that the total number of data points processed in each simulation is $N := \sum_{t=1}^{T} N_t = 505,450$. To compare algorithms using the same amount of data, the constant batch size is set to $\lfloor N/T \rfloor = 505$ for all iterations except the last one, which is augmented with the remaining data points and then includes 955 samples. This later adjustment has very little impact.

Polyak averaging (Polyak & Juditsky, 1992) is also performed, starting at $T_0$: for $t > T_0$, set

$$ \bar{\theta}^t := \frac{t-1-T_0}{t-T_0} \bar{\theta}^{t-1} + \frac{1}{t-T_0} \theta^t = \frac{1}{t-T_0} \sum_{k=T_0+1}^{t} \theta^k. $$

We will set $T_0 = 500$. For a quantitative performance assessment, we compute the root mean squared error (RMSE) or $\ell_2$-error; at iteration $t$, defined by

$$ \text{RMSE}_t := ||\theta^t - \theta^{\text{true}}||. $$

Note that the RMSE curve of a Polyak averaging sequence is not the Polyak averaging of the RMSE sequence.

The two simulation setups will be considered in turn: for the non-sparse case, we will use the objective function given by (18) with $\eta = 0$; see Figure 1 to Figure 5. Then, for the sparse case, we will use $\eta > 0$ and a penalty that makes F no longer convex; see Figure 6 to Figure 8.

Figure 1 compares up to 16 different settings of SAM2, SMD and SSG. For each instance, the algorithm is run over 100 independent simulation replications using the same QR model; the boxplot of the RMSE at iteration $\#T$ is displayed. SAM2 is run with linearly increasing batch size, and both the regular sequence $\{\theta^t, t \ge 1\}$ and the Polyak averaging sequence $\{\bar{\theta}^t, t \ge 1 + T_0\}$ are studied. SMD and SSG are run with three different decreasing step sizes $\gamma_t = 1/(t+1)^a$, for $a \in \{0.51, 0.6, 0.7\}$, and with linearly increasing batch size. Again, the RMSE at the end of the regular sequence and the Polyak averaging one are reported. For the case $a = 0.51$, we also analyze the constant batch size strategy. Combinations of settings that did not produce significantly different information have been omitted in the plots.

In Figure 2, we report the RMSE along simulation runs of several algorithms. SAM2 is run with the linearly increasing batch size. SMD is run with the same batch size strategy and the decreasing step size $\gamma_t = 1/(t+1)^{0.6}$. Four different SSG settings are run: with a linearly increasing batch size and the decreasing step size $\gamma_t = 1/(t+1)^{0.51}$; and with a constant batch size but three different decreasing step sizes $\gamma_t = 1/(t+1)^a$ for $a \in \{0.51, 0.6, 0.7\}$. When run with linearly increasing batch size, all the algorithms use the same data points at each iteration.

The RMSE evolution $t \mapsto \text{RMSE}_t$ is reported first with respect to the number of parameter updates or equivalently, the number of iterations (left); and second, with respect to the number of processed samples (middle). A focus on the last 500 iterations is also provided (right). For some algorithms,
<figure><img src="image_5.png" /><figcaption>Figure 4: LAD example. RMSE over iterations for SAM2 with constant batch size set from $N_t = 2$ to $N_t = 10$: (a) For $N_t = 1$ only the first 100 iterations are shown due to the explosion of the values, (b) For $N_t = 2$ note the very high RMSE values too, (c) Comparison of the $N_t = 5$ and $N_t = 10$, (d) $N_t = 10$ RMSE sequence and comparison with $N_t = 505$.</figcaption></figure>
the RMSE along the Polyak averaging sequence is compared to the one along the regular sequence.

Figure 3 presents the simulation run of one instance of SAM2, SSG and SMD when estimating the components #1, #2, #3, and #11 of parameter $\theta^{\text{true}}$. More specifically, the evolution of the relative error $t \mapsto (\theta_i^t - \theta_i^{\text{true}})/\theta_i^{\text{true}}$ is displayed. The relative error along the Polyak averaging sequence $\{\bar{\theta}_i^t, t \ge 1 + T_0\}$ is also displayed.

The performance of SAM2 applied with a constant batch size is then investigated. Figure 4 shows the evolution of the RMSE along two paths of SAM2 run with two different strategies for the batch size sequence: the linearly increasing one $N_t = \max(100, t)$ and the constant one $N_t = 505$ (except for the last iteration, where $N_T = 955$), shown with respect to iterations and with respect to number of processed samples. Figure 5 compares different constant batch sizes, $N_t = 1, N_t = 2, N_t = 5$ and $N_t = 10$, by showing the corresponding RMSE evolution along the regular SAM2 sequences.

The last three plots pertain to the non-convex objective function. Figure 6 shows the evolution of the RMSE computed along the regular and Polyak averaged sequences for SAM2, and the Polyak averaged sequences for SMD and SSG. All algorithms are run with the linearly increasing batch size; the step size is equal to $\gamma_t = 1/(t+1)^{0.6}$ for SMD and to $\gamma_t = 1/(t+1)^{0.51}$ for SSG. The RMSE is first computed for the $\ell$ components of the parameter, then for the components #j where $\theta_j^{\text{true}}$ is zero, and finally for the components #j where $\theta_j^{\text{true}}$ is not zero. Results are shown for two data sets to illustrate the algorithms stability. Here $\eta = 0.01$.

Then, the behavior at the end of $T$ iterations is studied: boxplots, over 100 repetitions, of the RMSE at iteration #T are displayed in Figure 7, for the regular algorithm SAM2, SMD and SSG and their Polyak averaging versions; the batch size increases linearly and the step size is $\gamma_t = 1/(t+1)^{0.6}$ for SMD and $\gamma_t = 1/(t+1)^{0.51}$ for SSG.

Similarly, Figure 8 shows iteration #T RMSE boxplots, computed from 100 independent replications, but for a larger $\eta$ set to 0.05. Only the results for SAM2 and SMD are shown, for both the regular sequence and the Polyak averaging one. They are run with linearly increasing batch size, and the step size is $\gamma_t = 1/(t+1)^{0.6}$ for SMD. A simulation run of SAM2 estimates is also displayed and compared to the sparse vector $\theta^{\text{true}}$.

## 5.2 The non-sparse case

We first consider the case $s = \ell = 11$ with strictly positive components for the true parameter $\theta^{\text{true}}$. Thus, parameter support recovery is not a goal and we set $\eta = 0$ in (18).

**Average assessment, after $T$ iterations.** Figure 1 shows that SSG is very sensitive to the step size strategy. For the most favorable step size $\gamma_t = 1/(t+1)^{0.51}$, constant batch size and linearly increasing batch size behave equivalently. Furthermore, Polyak averaging is not beneficial in any
<figure><img src="image_6.png" /><figcaption>Figure 5: LAD example, 100 simulations. Last iteration RMSE boxplots for SAM2 and SSG with different batch and step size schedules. Left plot: From left to right, SAM2, SAM2 with Polyak averaging, SSG and its Polyak averaging, with varying batch size, with constant batch sizes set to 505, both with $\gamma_t = \frac{1}{(t+1)^{0.51}}$, SSG with constant batch sizes set to 505 with $\gamma_t = \frac{1}{(t+1)^{0.6}}$ and $\gamma_t = \frac{1}{(t+1)^{0.7}}$. Right plot: zoom on the 6 best algorithms. Polyak averaging is performed for all algorithms with $T_0 = 500$.</figcaption></figure>

batch sizes.

## 5.4 Average assessment

The comparison is completed by repeating the same experiment for 100 different simulated data sets using the same QR model. Figure 5 shows the boxplots obtained with the SAM2 and SSG algorithms, using the RMSE at the last iteration. More specifically, eight settings are considered: SAM2 with $N_t = \max(100, t)$ and its Polyak averaging, SSG for $\gamma_t = \frac{1}{(t+1)^{0.51}}$ with both varying batch sizes and constant $N_t = 505$, and their respective Polyak averaging. SSG with $N_t = 505$ for $\gamma_t = \frac{1}{(t+1)^{0.6}}$ and $\gamma_t = \frac{1}{(t+1)^{0.7}}$ are also shown.

The same conclusions as in the previous sections hold. The sensitivity of SSG to step sizes is confirmed while SAM2 outperforms the other settings with a simple increasing sequence of $N_t$'s. SAM2 produces noisier sequences which are easily smoothed via Polyak averaging, which further improves the final SAM2 estimates. In contrast, Polyak averaging is not beneficial to SSG in none of the configurations. This is essentially due to the parameter sequences that are much more biased as illustrated in Figure 2. At last, the good performance of SAM2 is at the expense of a higher computation cost due to the resolution of an optimization problem at each iteration, while SSG involves only a straightforward update.
<figure><img src="image_2.png" /><figcaption>Figure 1: LAD example, 100 data sets. **Left plot:** 16 tested configurations with last iteration RMSE boxplots for SAM2 (black), SMD (red) and SSG (purple) with different batch and step size schedules. PA indicates Polyak averaging; postfixes 051, 06 and 07 stand respectively for $\gamma_t = 1/(t+1)^{0.51}$, $\gamma_t = 1/(t+1)^{0.6}$ and $\gamma_t = 1/(t+1)^{0.7}$. Batch sizes are set to $N_t = \max(100, t)$ except when specified by postfix EQ meaning constant batch sizes $N_t = 505$. **Right plot:** focus on the 14 best configurations.</figcaption></figure>

of the configurations. SMD is less sensitive to the decreasing rate of the step size than SSG; a linearly increasing batch size performs better than a constant one. Polyak averaging improves the variability of the SMD final estimates, except in the setting when $\gamma_t = 1/(t+1)^{0.7}$. Polyak averaging also improves the final RMSE for SAM2. The global comparison shows that SAM2 and SMD with slowly vanishing step size, both run with linearly increasing batch size and post-processed with Polyak averaging, are the best strategies. SSG is not able match such performances. The SAM2 and SMD estimations appear less biased than the SSG ones (see also Figure 3).

Based on Figure 1, we now analyze the simulation runs for a selection of the 16 instances. For SMD and SSG, we will consider the case when the batch size increases linearly and $\gamma_t = 1/(t+1)^{0.6}$ and $\gamma_t = 1/(t+1)^{0.51}$, respectively. Since SSG behaved poorly for other choices of step size when $N_t$ increases linearly, we decide to consider these step sizes with constant batch size. The results are displayed on Figure 2 and Figure 3.

**Sensitivity to step sizes.** SSG and SMD require the user to decide on a step size sequence, which has been reported to be impactful and the main practical limitation of gradient algorithms. The simulations confirm that SSG performance can be quite sensitive to the step size choice, with two out of three settings showing significantly larger RMSE.

On Figure 1, we observe that SMD is much less sensitive to the step size choice than SSG. The lowest RMSE are obtained for SAM2, for SMD with step sizes $\gamma_t = 1/(t+1)^{0.6}$, and for SSG with step sizes $\gamma_t = 1/(t+1)^{0.51}$, with SAM2 and SMD outperforming SSG. The optimal choice for SMD differs from that of SSG. This is not surprising and has been reported in the literature; see e.g., Guilmeau et al. (2025, Fig. 7). When compared through the RMSE at iteration #T (see Figure 1), SAM2 and
# 6 Proofs and Technical results

## 6.1 Technical results on F, G and T

**Proposition 11.** i) Assume A3-a. G is continuous on $T \times U$, and $\mathrm{Argmin}_T G(\cdot; v)$ is not empty for any $v \in U$.

ii) Assume A3-c. $\mathrm{Argmin}_T F$ is not empty.

iii) Assume A1. For any $v \in U$, $\min_T F \le \min_T G(\cdot; v)$. In addition, for any $\theta \in T$ and $v \in U[\theta]$, $0 \le G(\theta; v) + \min_T F \le F(\theta) + \min_T G(\cdot; v)$.

*Proof.* Under A3-a, the dominated convergence theorem implies that the function $(\theta, v) \to G(\theta; v)$ is continuous. Since $T$ is compact, $\inf_T G(\cdot; v)$ is not empty for any $v \in U$.

Under A3-c, $\theta \mapsto F(\theta)$ is continuous on the compact set $T$; hence, $\mathrm{Argmin}_T F$ is not empty.

Let $v \in U$. A1-a implies that $\min_T F \le \min_T G(\cdot; v)$. This inequality, combined with the equality in A1-b, yields the second statement. □

**Proposition 12.** Assume A3-a and A3-b. Then

i) $T$ is continuous on $U$.

ii) $(\theta, v) \mapsto G(\theta; v) - G(T(v); v)$ is continuous on $T \times U$.

*Proof.* (i) Let $v^* \in U$ and $\{v^t, t \ge 0\}$ be a $U$-valued sequence such that $\lim_t v^t = v^*$. Then, by A3-b, for all $\theta \in T$

$$G(T(v^t); v^t) \le G(\theta; v^t).$$

By A3-a, $G$ is continuous on $T \times U$ (see Proposition 11); then $\lim_t G(\theta; v^t) = G(\theta; v^*)$. The sequence $\{T(v^t), t \ge 0\}$ is a $T$-valued sequence and there exists a subsequence $\{v^{\rho(t)}, t \ge 0\}$ such that $\lim_t T(v^{\rho(t)})$ exists and is in $T$. Note that, as a subsequence of $\{v^t, t \ge 0\}$, we have $\lim_t v^{\rho(t)} = v^*$. Using again that $G$ is continuous, we have for all $\theta \in T$

$$G(\lim_t T(v^{\rho(t)}); v^*) \le G(\theta; v^*).$$

This inequality implies that $\lim_t T(v^{\rho(t)})$ is a minimizer of $\theta \mapsto G(\theta; v^*)$. Since this minimizer is unique by A3-b, we have $\lim_t T(v^{\rho(t)}) = T(v^*)$. This holds true for any subsequence $\{v^{\rho(t)}, t \ge 0\}$ and concludes the proof.

(ii) It follows from the continuity of $T$ and of $G$ (see Proposition 11). □

## 6.2 Proof of Theorem 4

For all $v \in U$, $\mathrm{Argmin}_T G(\cdot; v)$ is not empty by Proposition 11-i. $\mathrm{Argmin}_T F$ is not empty by Proposition 11-ii.
<figure><img src="image_3.png" /><figcaption>Figure 2: LAD example. RMSE with respect to the number of iterations (left) and with respect to the number of samples (middle). Linearly increasing batch sizes $N_t = \max(100, t)$ are used, except when the postfix "EQ" is indicated, meaning constant batch sizes. **SMD** (red) is run for a step size schedule $\gamma_t = 1/(t+1)^{0.6}$. **SSG** is run with different step size schedules, $\gamma_t = 1/(t+1)^{0.51}$ (blue and **purple**), $\gamma_t = 1/(t+1)^{0.6}$ (green) and $\gamma_t = 1/(t+1)^{0.7}$ (dashed **dark green**). The right plot shows a zoom on the last half iterations, with additional RMSE curves for the Polyak averaging sequences (dashed lines). Note that the y-axis is shown in a log scale transformation to emphasize small differences.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: LAD example. Relative error $t \mapsto (\theta_i^t - \theta_i^{\text{true}})/\theta_i^{\text{true}}$ for estimates of components $i \in \{1, 2, 3, 11\}$, using **SAM2** (black), **SMD** (red) and **SSG** (purple), all with linearly increasing batch sizes $N_t = \max(100, t)$. Step sizes are set to $\gamma_t = 1/(t+1)^{0.6}$ for **SMD** and to $\gamma_t = 1/(t+1)^{0.51}$ for **SSG**. Polyak averaging sequences are shown with dashed lines. The blue horizontal line indicates 0.</figcaption></figure>
*Proof of Item ii.* Let $\theta \in \operatorname{Argmin}_T \mathbf{F}$. From (7), we have for any $v \in \mathbb{U}[\theta]$, $\min_T \mathbf{F} = \mathbf{F}(\theta) \ge \mathbf{F}(\tau)$ for any $\tau \in \operatorname{Argmin}_T \mathbf{G}(\cdot; v)$. This yields $\tau \in \operatorname{Argmin}_T \mathbf{F}$ and $\mathbf{F}(\theta) = \mathbf{F}(\tau)$. Hence $\theta \in \mathcal{L}$.

Let $\theta^* \in \mathcal{L}$. Let $v \in \mathbb{U}[\theta^*]$ and $\tau \in \operatorname{Argmin}_T \mathbf{G}(\cdot; v)$. Since $\theta^* \in \mathcal{L}$, we have $\mathbf{F}(\tau) = \mathbf{F}(\theta^*)$. Combined with (7), this yields $\mathbf{G}(\theta^*; v) = \mathbf{G}(\tau; v)$. Therefore, $\mathbf{G}(\theta^*; v) = \min_T \mathbf{G}(\cdot; v)$ and $\theta^* \in \operatorname{Argmin}_T \mathbf{G}(\cdot; v)$.

*Proof of Item iii.* By definition, we have $\mathcal{L} = \{\theta \in \mathbb{T} : \mathbf{F}(\mathbf{T}(v)) = \mathbf{F}(\theta) \text{ for all } v \in \mathbb{U}[\theta]\}$.

Let $\theta \in \mathbb{T}$ such that $\theta \in \operatorname{Argmin}_T \mathbf{G}(\cdot; v)$ for all $v \in \mathbb{U}[\theta]$. By A3-b, $\operatorname{Argmin}_T \mathbf{G}(\cdot; v) = \{\mathbf{T}(v)\}$ which implies that $\theta = \mathbf{T}(v)$ for all $v \in \mathbb{U}[\theta]$. Hence, $\mathbf{F}(\theta) = \mathbf{F}(\mathbf{T}(v))$ for all $v \in \mathbb{U}[\theta]$ and $\theta \in \mathcal{L}$. The proof is concluded with (7).

## 6.3 Proof of Theorem 5

Since $\mathbf{F}$ is continuous and $\mathcal{L} \cap \mathbb{T}$ is compact, the set $\mathcal{D} := \mathbf{F}(\mathcal{L} \cap \mathbb{T})$ is a compact subset of $\mathbb{R}$. Denote by $\mathcal{D}_\beta$ the $\beta$-neighborhood of the closed set $\mathcal{D}$: $\mathcal{D}_\beta := \{w \in \mathbb{R} : \mathrm{d}(w, \mathcal{D}) < \beta\}$; $\mathrm{d}$ is the Euclidean distance from a point to a closed set. Since $\mathcal{D}$ is compact, we have

$$ \mathcal{D} = \bigcap_{\beta>0} \mathcal{D}_{\beta}. \qquad (18) $$

**Step 1.** Let $\alpha > 0$. Since $\mathcal{D}_\alpha$ is a finite union of disjoint bounded intervals, there exist a positive integer $n_\alpha$ and two increasing real-valued sequences $\{a_\alpha(k), 1 \le k \le n_\alpha\}$ and $\{b_\alpha(k), 1 \le k \le n_\alpha\}$ such that

$$ \mathcal{D}_{\alpha} = \bigcup_{k=1}^{n_{\alpha}} (a_{\alpha}(k), b_{\alpha}(k)). \qquad (19) $$

**Step 2.** Since $\mathbf{F}$ is continuous and $\mathcal{D}_{\alpha/2}$ is open, $\mathbf{F}^{-1}(\mathcal{D}_{\alpha/2})$ is an open covering of $\mathcal{L} \cap \mathbb{T}$. Define

$$ \epsilon_{\alpha} := \inf_{\theta \in \mathbb{T} \setminus \mathbf{F}^{-1}(\mathcal{D}_{\alpha/2}), v \in \mathbb{U}[\theta]} (\mathbf{F}(\theta) - \mathbf{F}(\mathbf{T}(v))); \quad \rho_{\alpha} := \epsilon_{\alpha} \wedge \alpha. \qquad (20) $$

Observe that since $\mathbb{T} \setminus \mathbf{F}^{-1}(\mathcal{D}_{\alpha/2})$ is a compact subset of $\mathbb{T} \setminus \mathcal{L}$, then $\epsilon_\alpha > 0$ (and therefore $\rho_\alpha > 0$) by H-ii.

**Step 3.** We write

$$ \mathbf{F}(\theta^t) - \mathbf{F}(\theta^{t+1}) = \mathbf{F}(\theta^t) - \mathbf{F}(\mathbf{T}(\nu^t)) + \mathbf{F}(\mathbf{T}(\nu^t)) - \mathbf{F}(\theta^{t+1}). \qquad (21) $$

By H-iii, there exists $T_\alpha$ such that

$$ \forall t \ge T_\alpha, \quad |\mathbf{F}(\mathbf{T}(\nu^t)) - \mathbf{F}(\theta^{t+1})| \le \rho_\alpha/2. \qquad (22) $$

Hence, if $t \ge T_\alpha$ and $\theta^t \in \mathbb{T} \setminus \mathbf{F}^{-1}(\mathcal{D}_{\alpha/2})$ then (20)-(22) imply that $\mathbf{F}(\theta^t) - \mathbf{F}(\theta^{t+1}) \ge \epsilon_\alpha - \rho_\alpha/2 \ge \rho_\alpha/2$. Consequently, upon noting that $\inf_T \mathbf{F} > -\infty$, the sequence $\{\mathbf{F}(\theta^t), t \ge 0\}$ is infinitely often in $\mathcal{D}_{\alpha/2}$ and therefore in $\mathcal{D}_\alpha$. This implies that the sequence $\{\mathbf{F}(\theta^t), t \ge 0\}$ is infinitely often in an interval of
SMD perform similarly with a slight advantage for SMD as also confirmed in Figure 2. Figure 3 illustrates that SAM2 and SMD estimates are in general closer to the true values, with SAM2 slightly closer than SMD in this specific example.

**Sensitivity to batch sizes.** Figure 2 (middle) illustrates that SAM2, SMD and SSG with linearly increasing batch sizes, show better start, reaching lower RMSE faster, than SSG run with constant batch size, regardless of the decreasing rate of the step size. Figure 2 also shows that a constant batch size strategy does not help SSG when the step size rapidly vanishes to zero: the sensitivity to the step size appears not to be as easily compensated by appropriate batch sizes. Also from Figure 1, we conclude that the best strategy for the step size is a slow decay ($\gamma_t = 1/(t+1)^{0.51}$), whatever the strategy for the batch size. For such a step size, the conclusions of Figure 2 (left, middle) are that more frequent parameter updating is preferable than finer optimization steps at the early stages of the algorithms. When the number of parameter updates (i.e. iterations) is large, the constant batch size and the linearly increasing strategies are almost equivalent, with a slight advantage to the first prior choice. The constant batch size strategy is not advocated for SMD, as illustrated on Figure 1 excepted if it is combined with a Polyak averaging post-processing.

<figure><img src="image_5.png" /><figcaption>Figure 4: LAD example. RMSE over iterations (left) and over samples (right), for SAM2 with different batch size sequences: linearly increasing batch size (black) and constant batch size ($N_t = 505$) (blue).</figcaption></figure>

Figure 4 shows the result of two SAM2 simulation runs, comparing two strategies for the batch size: the linearly increasing one and the constant one (although this does not satisfy the sufficient condition A4 on the $N_t$'s). This second setting provides better RMSE in the first iterations but is eventually equivalent as the batch sizes increase.

Figure 5 further shows SAM2 sequences with constant batch sizes successively set to constant between 1 and 10. It illustrates that $N_t$ should not be set too small. For $N_t = 1$, the algorithm diverges very quickly. For $N_t = 5$, RMSE values remain very high and noisy, while $N_t = 10$ still provides much larger RMSE than $N_t = 505$. Hence, taking $\{N_t, t \ge 1\}$ to be constant and small is not recommended. However, in the increasing $N_t$ setting, one can hope that starting with small $N_t$ will essentially slow down the estimation at the start without a critical impact if $N_t$ is gradually
(19): set $I_\alpha := (a_\alpha(k_\star), b_\alpha(k_\star))$ where $k_\star$ is defined by $k_\star := \max\{k \in \{1, \dots, n_\alpha\} : \liminf_t F(\theta^t) > a_\alpha(k)\}$.

**Step 4.** Let $t \ge T_\alpha$ such that $F(\theta^t) \in I_\alpha$. We prove by induction that for any $s \ge 0$, $F(\theta^{t+s}) < b_\alpha(k_\star)$. The property holds true for $s = 0$; assume it holds for $s > 0$. Using again (21), we distinguish two cases

* if $F(\theta^{t+s}) \in D_\alpha \setminus D_{\alpha/2}$, then (20)-(22) imply that $F(\theta^{t+s}) - F(\theta^{t+s+1}) \ge \epsilon_\alpha - \rho_\alpha/2 \ge \rho_\alpha/2$. Hence $F(\theta^{t+s+1}) \le F(\theta^{t+s}) - \rho_\alpha/2 < b_\alpha(k_\star) - \rho_\alpha/2 < b_\alpha(k_\star)$.

* if $F(\theta^{t+s}) \in D_{\alpha/2}$, then (21)-(22) and H-i imply that $F(\theta^{t+s}) - F(\theta^{t+s+1}) \ge 0 - \rho_\alpha/2$. Since $F(\theta^{t+s}) < b_\alpha(k_\star) - \alpha/2$ and $\rho_\alpha \le \alpha$, then $F(\theta^{t+s+1}) < b_\alpha(k_\star)$.

This concludes the induction.

**Step 5.** We have $\liminf_t F(\theta^t) > a_\alpha(k_\star)$ and there exists $t$ such that for all $\tau \ge 0$, $F(\theta^{t+\tau}) < b_\alpha(k_\star)$. Hence, $F(\theta^t) \in I_\alpha$ for any $t$ large enough: the set $\mathcal{F}$ of the limiting points of $\{F(\theta^t), t \ge 0\}$ is not empty and is included in the interval $I_\alpha$.

Let $0 < \alpha_1 < \alpha_2$: $I_{\alpha_1} \subseteq I_{\alpha_2}$ so that for any positive sequence $\{\alpha_n, n \ge 0\}$ decreasing to zero, $\mathcal{F} \subset \bigcap_n I_{\alpha_n}$ and $\bigcap_n I_{\alpha_n}$ is an interval. Note also that $\bigcap_n I_{\alpha_n} \subset F(\mathcal{L} \cap \mathbb{T})$ by (18). Therefore, the sequence $\{F(\theta^t), t \ge 0\}$ converges to an interval in $F(\mathcal{L} \cap \mathbb{T})$; this concludes the proof of the first claim.

**Step 6.** Let us prove the last statement. If $F(\mathcal{L} \cap \mathbb{T})$ has an empty interior, there exists $F^*$ such that $\lim_t F(\theta^t) = F^*$. Since $\{\theta^t, t \ge 0\}$ is a compact sequence and $F$ is continuous, then $\{\theta^t, t \ge 0\}$ converges to the set $\{\theta \in \mathbb{T} : F(\theta) = F^*\}$.

Let a converging subsequence $\{\theta^{\rho(t)}, t \ge 0\}$, with limiting value $\theta^* \in \mathbb{T}$; such a sequence exists since $\mathbb{T}$ is compact. The proof is by contradiction: assume that $\theta^* \notin \mathcal{L}$. Since $\mathcal{L} \cap \mathbb{T}$ is closed, there exists $\delta > 0$ such that the compact ball $\mathcal{K}_\delta := \{\theta \in \mathbb{T} : \|\theta - \theta^*\| \le \delta\}$ is in $\mathbb{T} \setminus \mathcal{L}$. The assumption H-ii implies that there exists $\epsilon_\delta > 0$ such that for all $\theta \in \mathcal{K}_\delta$ and for all $v \in \mathbb{U}[\theta]$, $F(\theta) > \epsilon_\delta + F(\mathbb{T}(v))$. Since $\lim_t \theta^{\rho(t)} = \theta^*$, there exists $T_\delta$ such that for all $t \ge T_\delta$, $\theta^{\rho(t)} \in \mathcal{K}_\delta$. Therefore for all $t \ge T_\delta$, $F(\theta^{\rho(t)}) > \epsilon_\delta + F(\mathbb{T}(\nu^{\rho(t)}))$ which yields

$$ \lim_t F(\theta^{\rho(t)}) \ge \epsilon_\delta + \lim_t F(\mathbb{T}(\nu^{\rho(t)})). \quad (23) $$

On the other hand, the assumption H-iii and the result $\lim_t F(\theta^t) = F^*$ imply that

$$ \lim_t F(\mathbb{T}(\nu^{\rho(t)})) = \lim_t F(\theta^{\rho(t)+1}), \quad F^* = \lim_t F(\theta^t) = \lim_t F(\theta^{\rho(t)}) = \lim_t F(\theta^{\rho(t)+1}). \quad (24) $$

The results (23) and (24) imply that $F^* \ge \epsilon_\delta + F^*$ which is a contradiction.
<figure><img src="image_6.png" /><figcaption>Figure 5: LAD example. RMSE over iterations for SAM2 with constant batch size set from $N_t = 2$ to $N_t = 10$. Left: for $N_t = 1$ only the first 100 iterations are shown due to the explosion of the values. Middle: for $N_t = 2$ note the very high RMSE values too. Right: comparison of the $N_t = 5$ (black) and $N_t = 10$ (blue).</figcaption></figure>

increased.

Overall, SAM2 has the advantage over SSG (i) not to require the choice of a step size sequence, (ii) while choosing batch sizes is not problematic. One could argue that the need to use increasing $N_t$'s could be limiting, typically if one does not control the data stream. But note that this is only a sufficient condition and that SAM2 seems to perform well for constant $N_t$'s too.

**Polyak averaging.** The role of Polyak averaging is to reduce the variability of a path. This is exactly what we observe in Figure 3; see also Figure 2 (right). Particularly, we observe that the addition of such post-processing to the regular SAM2 sequence and to the regular SMD sequence (with step size $\gamma_t = 1/(t+1)^{0.6}$), both run with a linearly increasing batch size, drastically reduces the variability of the estimates. The relative error criterion shows a small advantage to SAM2 when compared to SMD (see Figure 3). On the contrary, the RMSE criterion shows a small advantage to SMD. As expected, Polyak averaging does not help when the regular sequence is not convergent; see e.g. the purple curves in Figure 2 (right) and Figure 3.

## 5.3 The sparse case: quantile regression with a non-convex log-penalty

In the previous sections, we showed that SAM2 outperformed SSG in convex settings, while avoiding the need to tune the step sizes (SMD does not suffer from the same step size issue as SSG).

In this section, we now illustrate how SAM2 is able to handle a non-convex objective corresponding to the interesting scenario of sparsity-inducing penalties. While SMD can provide reasonable results in this setting, we show that only SAM2 is able to recover exact zeros of the coefficients vector, and thus outperforms both SSG and SMD in terms of support recovery.

We consider the same QR loss as before but with simulated data corresponding to the sparse parameter scenario. In order to study a non-convex example, we consider the case when $\eta > 0$ in Equation (18) with pen taken to be the log-penalty, which is a well-known non-convex regularization function used to promote sparsity more aggressively than the $\ell_1$-norm, while still being
## 6.4 Proof of Proposition 7

Let $\epsilon > 0$. Since $F$ is continuous (see A3-c) and $T$ is compact: there exists $\eta_\epsilon > 0$ such that

$$ \|\tau - \theta'\| \le \eta_\epsilon \implies |F(\tau) - F(\theta')| \le \epsilon/2. $$

and $\sup_T |F| < \infty$. It holds

$$ \begin{aligned} |F(T(\nu)) - F(\theta')| &= |F(T(\nu)) - F(\theta')| (1_{\|T(\nu)-\theta'\| \le \eta_\epsilon} + 1_{\|T(\nu)-\theta'\| > \eta_\epsilon}) \\ &\le \frac{\epsilon}{2} + \frac{2}{\eta_\epsilon} \sup_T |F| \|T(\nu) - \theta'\|. \end{aligned} $$

Set $\delta_\epsilon := \epsilon\eta_\epsilon/(4\sup_T |F|)$. Then we write by using (9)

$$ \begin{aligned} \|T(\nu) - \theta'\| &\le \|T(\nu) - \theta'\| 1_{\|T(\nu)-\theta'\| < \delta_\epsilon} + \|T(\nu) - \theta'\| 1_{\|T(\nu)-\theta'\| \ge \delta_\epsilon} \\ &\le \delta_\epsilon + \|T(\nu) - \theta'\| 1_{\|T(\nu)-\theta'\| \ge \delta_\epsilon} \\ &\le \delta_\epsilon + 2 \operatorname{Diam}(T) 1_{G(\theta';\nu)-G(T(\nu);\nu) \ge \tilde{\eta}\delta_\epsilon} \\ &\le \delta_\epsilon + \frac{2}{\tilde{\eta}\delta_\epsilon} \operatorname{Diam}(T) (G(\theta';\nu) - G(T(\nu);\nu)). \end{aligned} $$

As a conclusion, by using the definition of $\delta_\epsilon$, we have

$$ |F(T(\nu)) - F(\theta')| \le \epsilon + \frac{4}{\eta_\epsilon \tilde{\eta}\delta_\epsilon} \sup_T |F| \operatorname{Diam}(T) (G(\theta'; \nu) - G(T(\nu); \nu)). $$

## 6.5 Uniform strong Law of Large Numbers for double arrays

Theorem 13 establishes a uniform strong Law of Large Numbers for double arrays, with the proof following the same lines of argument as that of (Andrews, 1992, Thm. 3), for example.

**Theorem 13.** Assume

B1. i) $V$ is a compact subset of $\mathbb{R}^v$ and $(X, \mathcal{X})$ is a measure space.

ii) $h: V \times X \to \mathbb{R}$ is measurable and there exists $p > 1$ such that $E[\sup_{u \in V} |h(u, X)|^p] < \infty$, where $X$ is a $\mathcal{X}$-valued random variable defined on $(\Omega, \mathcal{A}, \mathbb{P})$.

iii) For all $x \in X$, $u \mapsto h(u, x)$ is continuous on $V$.

iv) The deterministic $\mathbb{Z}_{>0}$-valued sequence $\{N_t, t \ge 1\}$ satisfies $\sum_{t \ge 1} N_t^{-(p-1)\wedge(p/2))} < \infty$.

v) The random variables $\{X^{t,i}, t \ge 1, 1 \le i \le N_t\}$, defined on $(\Omega, \mathcal{A}, \mathbb{P})$, are i.i.d. with the same distribution as $X$.
differentiable almost everywhere. It is a non-convex alternative to the $l_1$-norm, and has been widely used in applications in compressed sensing, sparse regression, and signal processing (see, e.g., Zhao, 2018). For $\theta \in \mathbb{R}^l$, the log-penalty is defined as

$$ \text{pen}(\theta) := \sum_{j=1}^{l} \log \left( 1 + \frac{|\theta_j|}{\epsilon} \right), \qquad (23) $$

where $\epsilon > 0$ is a small constant to control smoothness near zero. As $\epsilon \to 0$, pen($\theta$) approximates the $l_0$-norm more closely. The $l_0$-norm of $\theta$ is the number of non-zero elements in $\theta$ and is therefore a strong sparsity-inducing norm. In contrast, when $\epsilon \to \infty$, pen($\theta$) tends to 0 (no-penalization). In between, as $\epsilon$ increases, the log-penalty becomes more closely resembling of the $l_1$-norm penalty:

$$ \text{pen}(\theta) = \sum_{j=1}^{l} \log \left( 1 + \frac{|\theta_j|}{\epsilon} \right) \approx \sum_{j=1}^{l} \frac{|\theta_j|}{\epsilon} = \frac{1}{\epsilon} \|\theta\|_1. $$

If the weight of the log-penalty is set to $\eta\epsilon$, the $l_1$-norm penalty with weight $\eta$ is recovered.

A property of the log-penalty is that it is locally concave for small and moderate values of $\epsilon$. The key effect is that it penalizes small coefficients $\theta_j$ more, and large ones less, than the standard $l_1$-norm. The log-penalty becomes sharper near 0, promoting stronger sparsity. It is differentiable almost everywhere and due to its separability, its gradient is easy to compute by the fact that for all $j \in [l]$,

$$ \frac{d}{d\theta_j} \log \left( 1 + \frac{|\theta_j|}{\epsilon} \right) = \frac{\text{sign}(\theta_j)}{\epsilon + |\theta_j|}, $$

which is larger near 0 than the $l_1$-norm gradient, and thus shrinks small values more aggressively than the $l_1$-norm.

Figures 6, 7 and 8 illustrate that separable non-convex and locally concave penalties are particularly favorable to the SAM2 algorithm, which is not sensitive to non-convexity and can benefit from the stronger sparsity-inducing effect of these penalties. In contrast, the performance of SSG and SMD deteriorates when dealing with non-convex objectives and the algorithms are much less effective in recovering the parameter support. They also produce estimates that are generally more biased.

We consider an aggressive log-penalty with $\epsilon = 0.001$ in (23), with penalty weight $\eta = 0.01$. Figure 6 shows the effect of the log-penalty on the $\theta$ sequences. As expected, neither SSG nor SMD can handle the non-convex penalty and both result in worse estimations than SAM2. For SMD, the non-zero coefficient estimates are similar to that of SAM2 (Figure 6, right), but the two algorithms differ significantly in their ability to recover the zero coefficients (Figure 6, middle). In particular, SAM2 is able to shrink these coefficients to exactly zero after only a few iterations. The better performance of SAM2 is further illustrated in Figure 7, which shows boxplots for the RMSE at iteration #T. Not only was SAM2 able to recover exactly the zero coefficients (Figure 7, middle),
Then, with probability one,

$$ \limsup_{t \to \infty} \sup_{u \in \mathbb{V}} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} \mathbf{h}(u, X^{t,i}) - \mathbb{E}[\mathbf{h}(u, X)] \right| = 0. $$

Set $\mathbf{H}(u) := \mathbb{E}[\mathbf{h}(u, X)]$. Observe that under B 1-ii and B 1-iii, the dominated convergence theorem implies that $\mathbf{H}$ is continuous on $\mathbb{V}$.

The proof of Theorem 13 is based on Proposition 14: statement i) establishes a strong Law of Large numbers for a fixed $u \in \mathbb{V}$, and statement ii) establishes a strong Law of Large numbers for the modulus of continuity defined for any $\delta > 0$, by

$$ \omega_{\delta}(x) := \sup_{u \in \mathbb{V}, u' \in \mathbb{V}, \|u-u'\| \le \delta} |\mathbf{h}(u, x) - \mathbf{h}(u', x)|. $$

In the terminology of Davidson (2021, Sec. 22.4), statement i) is referred to as almost sure pointwise convergence of the quantity $N_t^{-1} \sum_{i=1}^{N_t} \mathbf{h}(u, X^{t,i}) - \mathbf{H}(u)$ to zero, and statement ii) implies its strong asymptotic uniform stochastic equicontinuity (often shorted to strong stochastic equicontinuity, for brevity), via the Markov's inequality. By Davidson (2021, Thm. 22.8), when $\mathbb{V}$ is compact (and thus totally bounded and separable), taken together, almost sure pointwise convergence and strong stochastic equicontinuity are equivalent to almost sure uniform convergence to zero of $N_t^{-1} \sum_{i=1}^{N_t} \mathbf{h}(\cdot, X^{t,i}) - \mathbf{H}$. Thus, under A6, the conclusions of Theorem 13 and Proposition 14 are almost equivalent.

**Proposition 14.** Assume B 1.

i) For any $u \in \mathbb{V}$, there exists $\Omega_u \in \mathcal{A}$ such that $\mathbb{P}(\Omega_u) = 1$ and on $\Omega_u$

$$ \lim_{t \to \infty} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} \mathbf{h}(u, X^{t,i}) - \mathbf{H}(u) \right| = 0. $$

ii) For any $\delta > 0$, $\mathbb{E}[(\omega_{\delta}(X))^p] < \infty$ and there exists $\bar{\Omega}_{\delta} \in \mathcal{A}$ such that $\mathbb{P}(\bar{\Omega}_{\delta}) = 1$ and on $\bar{\Omega}_{\delta}$

$$ \lim_{t \to \infty} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} \omega_{\delta}(X^{t,i}) - \mathbb{E}[\omega_{\delta}(X)] \right| = 0. $$

Finally, $\lim_{\delta \to 0} \mathbb{E}[\omega_{\delta}(X)] = 0$.

*Proof. Item i.* Let $u \in \mathbb{V}$. Set $Y^{t,i} := \mathbf{h}(u, X^{t,i}) - \mathbf{H}(u)$. We prove that for any $\varepsilon > 0$, $\mathbb{P}(\lim \sup_t A_{\varepsilon}^t) = 0$ where $A_{\varepsilon}^t := \{ N_t^{-1} \left| \sum_{i=1}^{N_t} Y^{t,i} \right| \ge \varepsilon \}$. This will imply that with probability one, $\lim_{t \to \infty} N_t^{-1} \sum_{i=1}^{N_t} Y^{t,i} = 0$ and conclude the proof.
<figure><img src="image_7.png" /><figcaption>Figure 6: LAD example with a non-convex sparsity-inducing log-penalty ($\varepsilon = 0.001$, $\eta = 0.01$). The two lines show two different simulations. RMSE with respect to the number of samples, for all coefficients (left), zero coefficients (middle) and non-zero coefficients (right). All algorithms use linearly increasing batch sizes $N_t = \max(100, t)$. SSG and SMD are run with a step size schedule respectively of $\gamma_t = 1/(t+1)^{0.51}$ and $\gamma_t = 1/(t+1)^{0.6}$. For better viewing, the first iterations are not displayed and the plots start after 130, 200 samples have been processed. PA stands for Polyak averaging. Note that the y-axis is shown in a pseudo-log scale transformation $\log(0.01 + y)$ to emphasize small differences while handling zero values.</figcaption></figure>

but its estimations of the non-zero coefficients were also less biased (Figure 7, right). On average, SMD provides satisfactory RMSE values (Figure 7, left), but with a number of failures essentially due to poor zero coefficients estimations. This is visible in the second row of Figure 6 (middle) and in the numerous outliers in the boxplots of Figure 7 (middle). We emphasize that only SAM2 was able to provide a shrinkage exactly to zero (Figures 6 and 7, middle).

The better performance of SAM2 is particularly clear when the penalty weight is increased from $\eta = 0.01$ to $\eta = 0.05$, making the non-convexity of the objective function F further exaggerated. Here, SMD completely fails, as shown in Figure 8(top), where the RMSE values can be very large. In contrast, SAM2 is still able to handle this case and was still able to accurately estimate all zero and non-zero coefficients, except for the smallest non-zero component $\theta_1^{\text{true}} = 2$, which is wrongly shrunken to 0, as illustrated in Figure 8 (bottom right). As a result of this error regarding $\theta_1^{\text{true}}$, SAM2 provides RMSE values around 2, see Figure 8 (bottom middle).

## 5.4 Computational costs

The good performance of SAM2 can be expected to be at the expense of a higher computation cost due to the resolution of an optimization problem at each iteration, while SSG and SMD generally involve only a straightforward update provided that the (sub)gradient of the objective function
We write, by using the Markov inequality,

$$ \sum_{t \ge 1} \mathbb{P}(A_{\varepsilon}^t) = \sum_{t \ge 1} \mathbb{P} \left( \left| \sum_{i=1}^{N_t} Y^{t,i} \right| \ge \varepsilon \right) \le \varepsilon^{-p} \sum_{t \ge 1} N_t^{-p} \mathbb{E} \left[ \left| \sum_{i=1}^{N_t} Y^{t,i} \right|^p \right] $$

where $p > 1$ is given by B1-ii. The random variables $\{Y^{t,i}, t \ge 1, 1 \le i \le N_t\}$ are i.i.d. and centered under B1-ii and B1-v. Hence

$$ \mathbb{E} \left[ \left| \sum_{i=1}^{N_t} Y^{t,i} \right|^p \right] = \mathbb{E} \left[ \left| \sum_{i=1}^{N_t} Y^{\star,i} \right|^p \right], $$

where $\{Y^{\star,i}, i \ge 1\}$ are i.i.d. random variables with the same distribution as $Y^{1,1}$. $n \mapsto \sum_{i=1}^n Y^{\star,i}$ is a martingale sequence with $L^p$-moment (see B1-ii): by (Hall & Heyde, 1980, Section 2.4 Theorem 2.10), there exists a constant $C_p$ such that for any $t \ge 1$,

$$ \mathbb{E} \left[ \left| \sum_{i=1}^{N_t} Y^{\star,i} \right|^p \right] \le C_p \mathbb{E} \left[ \left| \sum_{i=1}^{N_t} |Y^{\star,i}|^2 \right|^{p/2} \right]. $$

From standard calculations (see e.g. Lemma 15), this yields

$$ \mathbb{E} \left[ \left| \sum_{i=1}^{N_t} Y^{\star,i} \right|^p \right] \le C_p N_t^{(p/2)\vee 1} \mathbb{E} \left[ |Y^{\star,1}|^p \right]. $$

Under B1-iv, the Borel-Cantelli lemma implies that $\mathbb{P}(\limsup_t A_{\varepsilon}^t) = 0$.

*Item ii.* Let $\delta > 0$. We write $\omega_{\delta}(x) \le 2 \sup_{u \in \mathbb{V}} |h(u, x)|$; from B1-ii, we have $\mathbb{E}[(\omega_{\delta}(X))^p] < \infty$.

Set $Y^{t,i} := \omega_{\delta}(X^{t,i}) - \mathbb{E}[\omega_{\delta}(X)]$. Following the same lines as in the proof of Item i, it can be proved that under B1-ii, B1-iv and B1-v, $\lim_t N_t^{-1} \sum_{i=1}^{N_t} Y^{t,i} = 0$ with probability one; details are omitted.

Let us apply the dominated convergence theorem. Let $x \in \mathbb{X}$. By B1-i and B1-iii, for any $\epsilon > 0$, there exists $\eta_{\epsilon,x} > 0$ such that $\|u-u'\| \le \eta_{\epsilon,x}$ implies $|h(u,x)-h(u',x)| \le \epsilon$. Therefore, for any $\epsilon > 0$, there exists $\eta_{\epsilon,x} > 0$ such that for any $\delta \in [0, \eta_{\epsilon,x}]$, $\omega_{\delta}(x) \le \epsilon$. This implies that $\lim_{\delta \to 0} \omega_{\delta}(x) = 0$. In addition, $\omega_{\delta}(X) \le 2 \sup_{u \in \mathbb{V}} |h(u,X)|$ and the RHS is integrable by B1-ii. Therefore, the assumptions of the dominated convergence Theorem are satisfied and $\lim_{\delta \to 0} \mathbb{E}[\omega_{\delta}(X)] = 0$. $\square$

**Proof of Theorem 13.** Define $\mathbb{I}_{\mathbb{Q}} := \{1/q, q \in \mathbb{Z}_{>0}\}$. Since $\mathbb{V}$ is compact (see B1-i), for any $\delta \in \mathbb{I}_{\mathbb{Q}}$, there exists a finite covering of $\mathbb{V}$ by closed balls $\{B_{\delta,1}, \dots, B_{\delta,L_{\delta}}\}$ of radius $\delta$; let us denote by $v_{\delta,1}, \dots, v_{\delta,L_{\delta}}$ the centers of these balls. Set

$$ \Omega^* := \left( \bigcap_{\delta \in \mathbb{I}_{\mathbb{Q}}, \ell \in \{1, \dots, L_{\delta}\}} \Omega_{v_{\delta, \ell}} \right) \cap \left( \bigcap_{\delta \in \mathbb{I}_{\mathbb{Q}}} \bar{\Omega}_{\delta} \right), $$
<figure><img src="image_8.png" /><figcaption>Figure 7: LAD example with a significantly non-convex sparsity-inducing log-penalty ($\varepsilon = 0.001$, $\eta = 0.01$), 100 simulations. Boxplots of the RMSE at iteration $T$ for **SAM2** (black), **SMD** (red) and **SSG** (purple). **SMD** and **SSG** are run with step size schedules respectively set to $\gamma_t = 1/(t+1)^{0.6}$ and $\gamma_t = 1/(t+1)^{0.51}$. Batch sizes are set to $N_t = \max(100, t)$ for all algorithms. **Top row:** RMSE for all (left), zero (middle) and non-zero (right) coefficients. **PA** indicates Polyak averaging. **Bottom row:** focus on the comparison between **SAM2** and **SMD**.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: LAD example with a significantly non-convex sparsity-inducing log-penalty ($\varepsilon = 0.001$, $\eta = 0.05$), 100 simulations. Last iteration RMSE boxplots for **SAM2** (black) and **SMD**. **SMD** is run with a step size schedule set to $\gamma_t = 1/(t+1)^{0.6}$. Batch sizes are set to $N_t = \max(100, t)$ for both algorithms. **First row:** RMSE for all (first column), zero (second column) and non-zero (third column) coefficients. **PA** indicates Polyak averaging. **Second row:** focus on the **SAM2** zero and non-zero coefficients RMSE. The right plot shows the 11 coefficient sequences for one **SAM2** run, with the true $\theta$ coefficients indicated with dotted black lines.</figcaption></figure>
where $\Omega_v$ and $\bar{\Omega}_\delta$ are given respectively by Proposition 14-i and Proposition 14-ii. Since $\Omega^*$ is a countable intersection of sets of probability one, then $\mathbb{P}(\Omega^*) = 1$. We prove that on $\Omega^*$,

$$ \limsup_{t \to \infty} \sup_{u \in \mathbb{V}} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} h(u, X^{t,i}) - H(u) \right| = 0. $$

Let $\epsilon > 0$. By Item ii, there exists $\delta_\epsilon > 0$ such that $|\mathbb{E}[\omega_{\delta_\epsilon}(X)]| \le \epsilon/6$; without loss of generality we can assume that $\delta_\epsilon \in \mathbb{I}_Q$ and we do so. We write

$$ \sup_{u \in \mathbb{V}} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} h(u, X^{t,i}) - H(u) \right| \le \sup_{\ell \in \{1, \dots, L_{\delta_\epsilon}\}} \sup_{u \in B_{\delta_\epsilon, \ell}} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} h(u, X^{t,i}) - H(u) \right|. $$

We now prove that on $\Omega^*$, there exists $T_\epsilon$ such that for all $t \ge T_\epsilon$, each of the $L_{\delta_\epsilon}$ terms in the RHS is upper bounded by $\epsilon$. This will conclude the proof.

Fix $\ell \in \{1, \dots, L_{\delta_\epsilon}\}$. We write $\sup_{u \in B_{\delta_\epsilon, \ell}} |N_t^{-1} \sum_{i=1}^{N_t} h(u, X^{t,i}) - H(u)| \le \sum_{j=1}^2 \mathcal{T}_{j,\epsilon,\ell}(t) + \mathcal{T}_{3,\epsilon,\ell}$ where

$$ \mathcal{T}_{1,\epsilon,\ell}(t) := \sup_{u \in B_{\delta_\epsilon,\ell}} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} h(u, X^{t,i}) - \frac{1}{N_t} \sum_{i=1}^{N_t} h(u_{\delta_\epsilon,\ell}, X^{t,i}) \right|, $$

$$ \mathcal{T}_{2,\epsilon,\ell}(t) := \left| \frac{1}{N_t} \sum_{i=1}^{N_t} h(u_{\delta_\epsilon,\ell}, X^{t,i}) - H(u_{\delta_\epsilon,\ell}) \right|, $$

$$ \mathcal{T}_{3,\epsilon,\ell} := \sup_{u \in B_{\delta_\epsilon,\ell}} |\mathbf{H}(u_{\delta_\epsilon,\ell}) - \mathbf{H}(u)|. $$

It holds by definitions of $\omega_\delta$ and $\delta_\epsilon$,

$$ \mathcal{T}_{1,\epsilon,\ell}(t) \le \frac{1}{N_t} \sum_{i=1}^{N_t} \omega_{\delta_\epsilon}(X^{t,i}) \le \frac{\epsilon}{6} + \left| \frac{1}{N_t} \sum_{i=1}^{N_t} \omega_{\delta_\epsilon}(X^{t,i}) - \mathbb{E}[\omega_{\delta_\epsilon}(X)] \right|; $$

by Item ii, since $\delta_\epsilon \in \mathbb{I}_Q$, on $\Omega^*$ there exists $T_{1,\epsilon}$ such that for any $t \ge T_{1,\epsilon}$, we have $\mathcal{T}_{1,\epsilon,\ell} \le \epsilon/3$.

By Item i, on $\Omega^*$ there exists $T_{2,\epsilon}$ such that for any $t \ge T_{2,\epsilon}$, $\mathcal{T}_{2,\epsilon,\ell}(t) \le \epsilon/3$.

Finally, $\mathcal{T}_{3,\epsilon,\ell} \le \mathbb{E}[\omega_{\delta_\epsilon}(X)] \le \epsilon/6$ by definition of $\delta_\epsilon$.

Set $T_\epsilon := T_{1,\epsilon} \lor T_{2,\epsilon}$. On $\Omega^*$, for any $t \ge T_\epsilon$ and for any $\ell \in \{1, \dots, L_{\delta_\epsilon}\}$,

$$ \sup_{u \in B_{\delta_\epsilon, \ell}} \left| \frac{1}{N_t} \sum_{i=1}^{N_t} h(u, X^{t,i}) - H(u) \right| \le \frac{2\epsilon}{3} + \frac{\epsilon}{6} \le \epsilon. $$
is available. Exact running times depend on many factors, including the objective function, the implementation of the pseudo-code, and the number of iterations. Thus, it is more informative to provide a complexity comparison in our QR example.

For one iteration $t$, with a batch size $N_t$, for a parameter $\theta$ in dimension $\ell$, SAM2 involves solving $\ell$ 1-dimensional optimization problems that simplify to finding a minimum over $N_t$ scalar values. The cost is then of order $O(\ell N_t \log N_t)$ operations, which can be reduced to $O(N_t \log N_t)$ if the $\ell$ optimization problems can be run in parallel. In comparison SMD and SSG only require gradient updates, and thus have iterative costs in $O(\ell N_t)$. As an illustration, without parallelization, for our QR example implemented by a custom R code, in dimension $\ell = 11$, for 1000 iterations and linearly increasing batch sizes, a SAM2 run takes $22.9 \pm 0.6$ seconds on average over 30 repetitions, while SSG and SMD runs take $1.1 \pm 0.1$ seconds, on an Intel MacBook Pro laptop with 2.4 Ghz processor. Thus, SAM2 running times can be further reduced to approximately 2 seconds when handling the $\ell = 11$ optimizations in parallel.

Lastly, we note that unlike constant batch size algorithms, where the per-iteration complexity remains constant as a function of the index $t$, in our increasing batch size setting, the complexity at iteration $t$ depends on $t$, specifically of order $O(N_t \log N_t)$. For our linearly increasing choice of $N_t$, this becomes $O(t \log t)$ per iteration. Consequently, after $T$ iterations, the total computational effort for our increasing batch size scheme is $O(T^2 \log T)$, which is greater than the constant batch size setting, whose total computational effort is $O(T)$.

# 6 Proofs and Technical results

## 6.1 Technical results on F, G and T

**Proposition 12.** i) Under A 3-a, $G$ is continuous on $\mathbb{T} \times \mathbb{U}$, and $\text{Argmin}_{\mathbb{T}}G(\cdot; v)$ is not empty for any $v \in \mathbb{U}$.

ii) Under A 3-c, $\text{Argmin}_{\mathbb{T}}F$ is not empty.

iii) Under A 1, for any $v \in \mathbb{U}$, $\min_{\mathbb{T}} F \le \min_{\mathbb{T}} G(\cdot; v)$. In addition, for any $\theta \in \mathbb{T}$ and $v \in \mathbb{U}[\theta]$, $0 \le G(\theta; v) + \min_{\mathbb{T}} F \le F(\theta) + \min_{\mathbb{T}} G(\cdot; v)$.

*Proof.* (i) Under A 3-a, the dominated convergence theorem implies that the function $(\theta, v) \to G(\theta; v)$ is continuous and since $\mathbb{T}$ is compact, $\text{Argmin}_{\mathbb{T}}G(\cdot; v)$ is not empty for any $v \in \mathbb{U}$. (ii) Under A 3-c, $\theta \mapsto F(\theta)$ is continuous on the compact set $\mathbb{T}$, hence $\text{Argmin}_{\mathbb{T}}F$ is not empty. (iii) Let $v \in \mathbb{U}$. Then, A 1-a implies that $\min_{\mathbb{T}} F \le \min_{\mathbb{T}} G(\cdot; v)$. This inequality, combined with the equality in A 1-b, yields the second statement. $\square$

**Proposition 13.** Assume A 3-a and A 3-b. Then

i) $\mathbb{T}$ is continuous on $\mathbb{U}$.
**Lemma 15.** Let {$Y^i, i \ge 1$} be i.i.d. variables such that $\mathbb{E}[|Y^1|^p] < \infty$ for some $p > 1$. Then,

$$ \mathbb{E}\left[\left|\sum_{i=1}^{N} (Y^i)^2\right|^{p/2}\right] \le N^{(p/2)\nu 1} \mathbb{E}[|Y^1|^p]. $$

*Proof.* First case: $p/2 \le 1$. We write

$$ \mathbb{E}\left[\left|\sum_{i=1}^{N} (Y^i)^2\right|^{p/2}\right] \le \mathbb{E}\left[\left|\sum_{i=1}^{N} |(Y^i)^2|^{p/2}\right|\right] = \mathbb{E}\left[\left|\sum_{i=1}^{N} |Y^i|^p\right|\right] = N \mathbb{E}[|Y^1|^p]. $$

Second case: $p/2 \ge 1$. By the Minkowski inequality, we write

$$ \mathbb{E}\left[\left|\sum_{i=1}^{N} (Y^i)^2\right|^{p/2}\right] = \left(\left\|\sum_{i=1}^{N} (Y^i)^2\right\|_{p/2}\right)^{p/2} \le \left(\sum_{i=1}^{N} \|(Y^i)^2\|_{p/2}\right)^{p/2}. $$

The RHS is equal to

$$ \left(\sum_{i=1}^{N} \mathbb{E}[|Y^i|^p]^{2/p}\right)^{p/2} = \left(\sum_{i=1}^{N} \mathbb{E}[|Y^1|^p]^{2/p}\right)^{p/2} = \left(N \mathbb{E}[|Y^1|^p]^{2/p}\right)^{p/2} = N^{p/2} \mathbb{E}[|Y^1|^p]. $$

This concludes the proof. $\square$

## 6.6 Proof of Theorem 9

### 6.6.1 Proof of Item i

We write

$$ \begin{aligned} & G(\theta^{t+1}; v^t) - G(T(v^t); v^t) \\ & \le G(\theta^{t+1}; v^t) - \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta^{t+1}, X^{t+1,i}; v^t) + \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(T(v^t), X^{t+1,i}; v^t) - G(T(v^t); v^t) \end{aligned} $$

where we used that, by definition of $\theta^{t+1}$,

$$ \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(\theta^{t+1}, X^{t+1,i}; v^t) - \frac{1}{N_{t+1}} \sum_{i=1}^{N_{t+1}} g(T(v^t), X^{t+1,i}; v^t) \le 0. $$