---
abstract: |
    We propose a single time-scale stochastic subgradient method for constrained optimization of a composition of several nonsmooth and nonconvex functions. The functions are assumed to be locally Lipschitz and differentiable in a generalized sense. Only stochastic estimates of the values and generalized derivatives of the functions are used. The method is parameter-free. We prove convergence with probability one of the method, by associating with it a system of differential inclusions and devising a nondifferentiable Lyapunov function for this system. For problems with functions having Lipschitz continuous derivatives, the method finds a point satisfying an optimality measure with error of order $1/\sqrt{N}$, after executing $N$ iterations with constant stepsize.

    *Keywords:*[ Stochastic Composition Optimization, Nonsmooth Optimization, Stochastic Approximation, Risk-Averse Optimization, Stochastic Variational Inequality]{}\
    *AMS*: 90C15, 49J52, 62L20
author:
- 'Andrzej Ruszczyński[^1]'
date: 'January 29, 2020'
title: 'A Stochastic Subgradient Method for Nonsmooth Nonconvex Multi-Level Composition Optimization '
---

Introduction
============

We consider the composition optimization problem $$\label{main_prob}
\min_{x\in X}\; f_1\Big(x,f_2\big(x, \cdots f_{M-1}\big(x,f_M(x)\big)\cdots \big)\Big),$$ where $X\subset {\mathbbm{R}}^n$ is convex and closed, and $f_m:{\mathbbm{R}}^n\times {\mathbbm{R}}^{d_{m+1}}\to {\mathbbm{R}}^{d_m}$, $m=1,\dots,M-1$, and $f_M:{\mathbbm{R}}^n\to {\mathbbm{R}}^{d_M}$ are locally Lipschitz continuous functions, possibly neither convex nor smooth.

The Clarke derivatives of $f_m(\cdot,\cdot)$ are not available; instead, we postulate access to their random estimates. Such situations occur in *stochastic composition optimization*, where $$\label{stoch_prob}
f_m(x,u_{m+1}) = {\mathbbm{E}}\big[ \varphi_m(x,u_{m+1},\xi_m)\big],\quad m =1,\dots,M,$$ in which $\xi_m$ is a random vector, and $\mathbb{E}$ denotes the expected value. Two examples illustrate the relevance of the problem.

\[e:1\] [Given some random loss function $H:{\mathbbm{R}}^n\times \Omega\to {\mathbbm{R}}^n$, we consider the random variable $Z = H(x)$ as an element of a Banach space ${\mathcal{Z}}$ (for example, ${\mathcal{L}}_2(\varOmega,{\mathcal{F}},P)$) and evaluate its quality by the functional $\rho: {\mathcal{Z}}\to {\mathbbm{R}}$, called a *risk measure*. This leads to the problem $$\label{risk-opt}
\min_{x\in X} F(x) = \rho[H(x)].$$ Particularly simple and useful is the *mean–semideviation measure* [@OR:1999; @OR:2001], which is an example of a coherent measure of risk [@ADEH:1999] (see also [@follmer2011stochastic; @shapiro2009lectures] and the references therein). It has the following form: $$\label{msd}
\rho[Z] =  {\mathbbm{E}}[Z]  + \varkappa\, \bigg({\mathbbm{E}}\Big[\Big(\max\big(0,Z - {\mathbbm{E}}[Z]\big)\Big)^p\Big]\bigg)^{1/p},\quad p \ge 1.$$ For $p=1$, problem with the risk measure can be written in the form with two functions $$\begin{aligned}
f_1(x,u_2) &= {\mathbbm{E}}\Big[H(x) +  \varkappa\, \max\big(0,H(x) - u_2\big)\Big],\\
f_2(x) &= {\mathbbm{E}}[H(x)].\end{aligned}$$ For $p=2$, problem with the risk measure can be cast (with a minor regularization) in the form with three functions $$\begin{aligned}
f_1(x,u_2) &= {\mathbbm{E}}\big[H(x)\big] + \varkappa\,(\varepsilon+ u_2)^{1/2}, \quad \varepsilon>0,\\
f_2(x, u_3) &= {\mathbbm{E}}\Big[\Big(\max\big(0,H(x) - u_3\big)\Big)^2\Big],\\
f_3(x) &= {\mathbbm{E}}[H(x)].\end{aligned}$$ We added $\varepsilon>0$ into the definition of $f_1(\cdot,\cdot)$ to ensure local Lipschitz property for all $u_2\ge 0$. Frequently, $H(\cdot)$ is non-convex and non-differentiable in modern machine learning models. ]{}

\[e:SVI\] [We have a random mapping $H:{\mathbbm{R}}^n\times \Omega\to {\mathbbm{R}}^n$ on some probability space $(\varOmega,{\mathcal{F}},P)$ and a closed convex set $X$. The problem is to find ${x}\in X$ such that $$\label{SVI}
\big\langle \mathbb{E}[H(x)], \xi - {x} \big\rangle  \le 0, \quad \text{for all}\quad \xi \in X.$$ The reader is referred to the recent publications [@iusem2017extragradient] and [@koshal2013regularized] for a discussion of the challenges associated with this problem and its applications to stochastic equilibria. We may reformulate problem as by defining $f_1: {\mathbbm{R}}^n \times {\mathbbm{R}}^n\to {\mathbbm{R}}$ as $$f_1(x,u) = \max_{y \in X}\, \left\{ \langle u, y - x \rangle - \frac{r}{2} \| y - x\|^2 \right\},\quad r>0,\label{gad_function}$$ and $f_2(x) = \mathbb{E}[H(x)]$. In this case, we have access to the gradient of $f_1$, but the value and the Jacobian of $f_2$ must be estimated. We do not require $H(\cdot)$ to be monotone or differentiable. ]{}

The research on stochastic subgradient methods for nonsmooth and nonconvex functions started in the late 1970’s: see Nurminski [@Nurminski1979numerical] for weakly convex functions and a general methodology for studying convergence of non-monotonic methods, Gupal [@gupal1979stochastic] for convolution smoothing (mollification) of Lipschitz functions and resulting finite-difference methods, and Norkin [@norkin:phd] and [@mikhalevich1987nonconvex Ch. 3 and 7] for unconstrained problems with “generalized-differentiable” functions.

Recently, by an approach via differential inclusions, Duchi and Ruan [@duchi2018stochastic] studied proximal methods for sum-composite problems with weakly convex functions, Davis *et al.* [@davis2018stochastic] studied the subgradient method for locally Lipschitz Whitney $\mathcal{C}^1$-stratifiable functions with constraints, and Majewski *et al.* [@majewski2018analysis] studied several methods for subdifferentially regular Lipschitz functions.

The research on composition optimization problems started from penalty functions for stochastic constraints and composite regression models in [@Ermoliev71] and [@ermoliev1976methods Ch. V.4]. An established approach was to use two-level stochastic recursive algorithms with two stepsize sequences in different time scales: a slower one for updating the main decision variable $x$, and a faster one for tracking the value of the inner function(s). References [@wang2017stochastic; @WaLiFa17] provide a detailed account of these techniques and existing results. In [@yang2019multi] these ideas were extended to multilevel problems of form , albeit with multiple time scales and under continuous differentiability assumptions.

A Central Limit Theorem for problem has been established in [@dentcheva2017statistical]. Large deviation bounds for the empirical optimal value were derived in [@ermoliev2013sample].

The first single time-scale method for a two-level version ($M=2$) of problem with continuously differentiable functions has been recently proposed in [@ghadimi2018single]. It has the complexity of ${\cal O}(1/\epsilon^2)$ to obtain an $\varepsilon$-solution of the problem, the same as methods for one-level unconstrained stochastic optimization. However, the construction of the method and its analysis depend on the Lipschitz continuity of the gradients of the functions involved, and its parameters depend on the corresponding Lipschitz constants.

To the best of our knowledge, there has been no research on stochastic subgradient methods for composition problems of form where the functions involved may be neither convex nor smooth.

Our main objective is propose a single time-scale method for solving the multiple composition problem , and to establish its convergence with probability one on a broad class of problems, in which the functions $f_m(\cdot,\cdot)$, $m=1,\dots,M-1$, are assumed to be locally Lipschitz and admit a chain rule, while $f_M(\cdot)$ may be only differentiable in a generalized sense (to be defined in section \[s:2\]). The class of such functions is broader than the class of locally Lipschitz semismooth functions [@mifflin1977semismooth]. The method’s few parameters may be set to arbitrary positive constants. The main idea is to lift the problem to a higher dimensional space and to devise a single time scale scheme for not only estimating the solution, but also the values of all functions nested in , and the generalized gradient featuring in the optimality condition.

Our approach uses the differential inclusion method (see [@ljung1977analysis; @kushner2012stochastic; @kushner2003stochastic] and the references therein). Extension to differential inclusions was proposed in [@benaim2005stochastic; @benaim2006stochastic] and further developed, among others, in [@borkar2009stochastic; @duchi2018stochastic; @davis2018stochastic; @majewski2018analysis]. In our analysis, we use a specially tailored nonsmooth Lyapunov function for the method, which generalizes the idea of [@ruszczynski1986method; @ruszczynski87] to the composite setting.

The second contribution is the error analysis after finitely many iterations of the method. We prove that a non-optimality measure for problem , which corresponds to the squared norm of a gradient in the unconstrained one-level case, decreases at the rate $1/\sqrt{N}$, where $N$ is the number of iterations of the method. This matches the best rate estimates for general unconstrained one-level problems [@GhaLan12], and the statistical estimate of [@dentcheva2017statistical] for plug-in estimates of composition risk functionals; it outperforms the estimate of the method of [@yang2019multi].

The paper is organized as follows. In §\[s:2\], we recall the main facts on generalized differentiation and the chain rule on a path. In §\[s:3\], we describe our stochastic subgradient method for problem . In §\[s:4\] we introduce relevant multifunctions and prove the boundedness of the sequences generated by the method. §\[s:analysis\] contains the proof of its convergence for nonconvex and nonsmooth functions. Finally, in §\[s:5\], we provide solution quality guarantees after finitely many iterations with a constant stepsize, in the case when the functions have Lipschitz continuous derivatives.

Generalized Subdifferentials of Composite Functions {#s:2}
===================================================

We consider problem for locally Lipschitz continuous functions $f_m(\cdot)$ satisfying additional conditions of generalized differentiability, subdifferential regularity, or Whitney $\mathcal{C}^1$-stratification. Recall that $f:{\mathbbm{R}}^n\to {\mathbbm{R}}^m$ is *locally Lipschitz*, if for every $x_0\in {\mathbbm{R}}^n$ a constant $L$ and an open set $U$ containing $x_0$ exist, such that $\|f(x) - f(y)\| \le L\|x -y\|$ for all $x,y\in U$. A *Clarke generalized Jacobian* of $f(\cdot)$ at $x$ is defined as follows [@clarke1975generalized]: $$\partial\! f(x) = {\text{conv}}\big\{ \lim_{\substack{{y\to x}\\{y\in {\mathcal{D}}(f)}}} f'(y)\big\},$$ where ${\mathcal{D}}(f)$ is the set of points $y$ at which the usual Jacobian $f'(y)$ exists. To simplify further notation, $\partial\! f(x)$ of a function $f:{\mathbbm{R}}^n\to {\mathbbm{R}}^m$ is always understood as a set of $m\times n$ matrices (also for $m=1$). A locally Lipschitz function $f:{\mathbbm{R}}^n\to{\mathbbm{R}}$ is *subdifferentally regular* if for any $x\in {\mathbbm{R}}^n$ and $d\in {\mathbbm{R}}^n$ the directional derivative exists and satisfies the equation: $$\lim_{\tau\downarrow 0} \frac{f(x+\tau d)-f(x)}{\tau} = \max_{g \in \partial\!f(x)} gd.$$ *Whitney $\mathcal{C}^1$-stratification* is a partition of a graph of $f(\cdot)$ into finitely many subsets (strata), such that within each of them the function is continuously differentiable, and a special compatibility condition for the normals on the common parts of the closures of the strata is satisfied (see, [@bolte2007clarke] and the references therein).

We also recall the following definition.

\[d:Norkin\] A function $f:{\mathbbm{R}}^n\to{\mathbbm{R}}$ is *differentiable in a generalized sense at a point* $x\in {\mathbbm{R}}^n$, if an open set $U\subset {\mathbbm{R}}^n$ containing $x$, and a nonempty, convex, compact valued, and upper semicontinuous multifunction $G_f: U \rightrightarrows {\mathbbm{R}}^n$ exist, such that for all $y\in U$ and all $g \in G_f(y)$ the following equation is true: $$f(y) = f(x) +  g(y-x) + o(x,y,g),$$ with $$\lim_{y\to x} \sup_{g\in G(y)} \frac{o(x,y,g)}{\|y-x\|}=0.$$ The set $G_f(y)$ is the *generalized derivative* of $f$ at $y$. If a function is differentiable in a generalized sense at every $x \in {\mathbbm{R}}^n$ with the same generalized derivative $G_f:{\mathbbm{R}}^n\rightrightarrows {\mathbbm{R}}^n$, we call it simply *differentiable in a generalized sense*. For functions with values in ${\mathbbm{R}}^n$, generalized differentiability is understood component-wise.

The class of such functions is contained in the set of locally Lipschitz functions, and contains all semismooth locally Lipschitz functions. The generalized derivative $G_f(\cdot)$ is not uniquely defined in Definition \[d:Norkin\], which is essential for us, but the Clarke Jacobian $\partial\!f(x)$ is an inclusion-minimal generalized derivative. The class of such functions is closed with respect to composition and expectation, which allows for easy generation of stochastic subgradients in our case. For full exposition, see [@mikhalevich1987nonconvex Ch. 1 and 6].

An essential step in the analysis of stochastic recursive algorithms by the differential inclusion method is the *chain rule on a path* (see [@davis2019stochastic] and the references therein). For an absolutely continuous function $p:[0,\infty)\to{\mathbbm{R}}^n$ we denote by ${\overset{\hbox{\,\tiny${}_\bullet$}}{p}}(\cdot)$ its weak derivative: a measurable function such that $$p(t) = p(0) + \int_0^t {\overset{\hbox{\,\tiny${}_\bullet$}}{p}}(s)\;ds,\quad \forall \;  t \ge 0.$$

\[d:chain\] A locally Lipschitz continuous function $f:{\mathbbm{R}}^n \to {\mathbbm{R}}^m$ admits a chain rule on absolutely continuous paths, if $$\label{chain-path}
f(p(T))- f(p(0)) = \int_0^T   g(p(t)) \, {\overset{\hbox{\,\tiny${}_\bullet$}}{p}}(t)  \;dt,$$ for any absolutely continuous path $p:[0,\infty)\to {\mathbbm{R}}^n$, all selections $g(\cdot) \in \partial\!f(\cdot)$, and all $T>0$. It admits a chain rule on continuously differentiable paths, if is true for all continuously differentiable $p:[0,\infty)\to {\mathbbm{R}}^n$, all selections $g(\cdot) \in \partial\!f(\cdot)$, and all $T>0$.

Formula is true for convex functions [@brezis1971monotonicity] and, as recently demonstrated in [@drusvyatskiy2015curves], for subdifferentiably regular locally Lipschitz functions and Whitney $\mathcal{C}^1$-stratifiable locally Lipschitz functions. In [@Rusz2019] we proved that generalized differentiable functions admit the chain rule on generalized differentiable paths.

To formulate optimality conditions for our problem, and construct and analyze our method, we need to introduce several relevant multifunctions. For a point $x\in {\mathbbm{R}}^n$ we consider the generalized Jacobians $\partial\! f_m(\cdot)$, $m=1,\dots,M-1$, and we recursively define the sets and vectors: $$\label{GF}
\begin{aligned}
G_M(x) &= \partial\! f_M(x),\quad v_M = f_M(x);\\
G_m(x) &= {\rm conv} \big\{ z\in {\mathbbm{R}}^{n}: z = g_x + g_u J,\ (g_x,g_u) \in \partial\! f_m(x, v_{m+1}), \ J \in G_{m+1}(x) \big\}, \\
&\quad v_m = f_m(x,v_{m+1}),\quad m= M-1,\dots,1.
\end{aligned}$$ By [@mikhalevich1987nonconvex Thm. 1.6], each set $G_m(x)$ is a generalized Jacobian of the function $$F_m(x) = f_m\Big(x,f_{m+1}\big(x, \cdots f_{M-1}\big(x,f_M(x)\big)\cdots \big)\Big), \quad m=1,\dots,M,$$ at $x$. We also have $\partial F_m(x) \subseteq G_m(x)$ [@clarke1981generalized]. We call a point $x^*\in X$ *stationary* for problem , if $$\label{stationary}
0 \in  G_1(x^*) + N_X(x^*).$$ The set of stationary points of problem is denoted by $X^*$.

The single time-scale method with filtering {#s:3}
===========================================

The method generates $M+1$ random sequences: approximate solutions $\{x^k\}_{k\ge 0}\subset {\mathbbm{R}}^n$, path-averaged generalized subgradient estimates $\{z^k\}_{k\ge 0}\subset {\mathbbm{R}}^n$, and path-averaged inner functions estimates $\{u_m^k\}_{k\ge 0}\subset {\mathbbm{R}}^{d_{m}}$, $m=1,\dots,M$, all defined on a certain probability space $(\Omega,{\mathcal{F}},P)$. We let ${\mathcal{F}}_k$ to be the $\sigma$-algebra generated by $\{x^0,\dots,x^k,z^0,\dots,z^k,u^0,\dots,u^k\}$, with each $u^j=(u_1^j,\dots,u_{M}^j)$.

The method starts from $x^0 \in X$, $z^0 \in {\mathbbm{R}}^n$, $u^0 \in {\mathbbm{R}}^{d_1}\times\dots{\mathbbm{R}}^{d_M}$, and uses parameters $a>0$, $b>0$, and $\rho>0$. At each iteration $k=0,1,2,\dots$, we compute $$\label{QP}
y^k = {\mathop{\rm argmin}}_{y \in X}\  \left\{\langle z^k, y-x^k \rangle + \frac{\rho}{2} \|y-x^k\|^2\right\},$$ and, with an ${\mathcal{F}}_k$-measurable stepsize $\tau_k \in \big(0,\min(1,1/a,1/b)\big]$, we set $$\label{def_xk}
x^{k+1} = x^k + \tau_k (y^k-x^k).$$ Then, we obtain statistical estimates: $$\begin{aligned}
& \tilde{J}_m^{k+1}=\begin{bmatrix}\tilde{J}_{mx}^{k+1} & \tilde{J}_{mu}^{k+1}\end{bmatrix}
\text{ of }{J}_m^{k+1}=\begin{bmatrix} {J}_{mx}^{k+1} & {J}_{mu}^{k+1}\end{bmatrix}\in \partial\!f_m(x^{k+1},u_{m+1}^k),\quad m=1,\dots,M-1;\\
& \tilde{J}_M^{k+1}
\text{ of }{J}_M^{k+1}\in \partial\!f_M(x^{k+1});\\
& \tilde{h}_m^{k+1} \text{ of } f_m(x^{k+1},u_{m+1}^k),\quad  m=1,\dots,M-1;\\
& \tilde{h}_M^{k+1} \text{ of } f_M(x^{k+1}).\end{aligned}$$ We use the stochastic subgradients $\tilde{J}_m^{k+1}$ to construct a biased estimate of the composite function: $$\label{g-fold}
\tilde{g}_M^{k+1} = \tilde{J}_M^{k+1}, \quad
\tilde{g}_{m}^{k+1} = \tilde{J}_{mx}^{k+1} +  \tilde{J}_{mu}^{k+1} \tilde{g}_{m+1}^{k+1}, \quad m=M-1,\dots,1.$$ Finally, we update the path averages by backward recursion as follows: $$\begin{aligned}
z^{k+1} &= z^k + a\tau_k \Big( \big[\tilde{g}_{1}^{k+1}\big]^T  - z^k\Big),  \label{def_zk}\\
u_M^{k+1} &= u_M^k +  \tilde{J}_M^{\,k+1} (x^{k+1}-x^k) + b \tau_k \big(\tilde{h}_M^{k+1}-u_M^k\big),  \label{def_ukM}\\
u_m^{k+1} &= u_m^k +  \tilde{J}_m^{\,k+1} \begin{bmatrix}  x^{k+1}-x^k \\u_{m+1}^{k+1}-u_{m+1}^k \end{bmatrix}
+ b \tau_k \big(\tilde{h}_m^{k+1}-u_m^k\big),\quad  m=M-1,\dots, 1.\label{def_ukm}\end{aligned}$$ In fact, we do not need the sequence $\{u_1^k\}$ for the operation of the method, but we include it for uniformity of notation; it will also provide an estimate of the function value at the optimal solution.

We will analyze convergence of the algorithm – under the following conditions:

The set $X$ is convex and compact;

The functions $f_m(\cdot,\cdot)$, $m=1,\dots,M-1$, are Lipschitz continuous and admit the chain rule for every path $(x(\cdot),u_{m+1}(\cdot))$ with a continuously differentiable $x(\cdot)$ and absolutely continuous $u_{m+1}(\cdot)$. Moreover, for every $x\in X$ and $u_{m+1}\in {\mathbbm{R}}^{d_{m+1}}$ the $u$-part of the generalized Jacobian, $\partial_u f_m(x,u_{m+1})$, is single-valued.

The function $f_M(\cdot)$ is Lipschitz continuous and admits the chain rule for every continuously differentiable path $x(\cdot)$.

A constant $C$ exists, such that $\|g \| \le C( 1 + \|v\|)$ for all $g \in \partial f_m(x,v)$, all $x\in X$, and all $v\in {\mathbbm{R}}^{d_{m+1}}$, $m=1,\dots,M-1$.

The set $\{ F_1(x): x\in X^*\}$ does not contain an interval of nonzero length;

$\tau_k \in \big(0,\min(1,1/a,1/b)\big]$ for all $k$, $\lim_{k\to\infty}\tau_k= 0$, $\sum_{k=0}^\infty \tau_k = \infty$, $\sum_{k=0}^\infty {\mathbbm{E}}[\tau_k^2] < \infty$;

For all $k$,

[iii]{}

$\tilde{h}_m^{k+1} = f_m(x^{k+1},u_{m+1}^{k}) + e_m^{k+1} + \delta_m^{k+1}$, $m=2,\dots,M-1$, and\
$\tilde{h}_M^{k+1} = f_M(x^{k+1}) + e_M^{k+1} + \delta_M^{k+1}$, with $\lim_{k\to \infty} \delta_m^{k+1}=0$ and\
${\mathbbm{E}}\big\{e_m^{k+1}\big|{\mathcal{F}}_k\big\} = 0$, ${\mathbbm{E}}\big\{\|e_m^{k+1}\|^2|{\mathcal{F}}_k\big\}\le \sigma_e^2$, $m=2,\dots,M$;

$\tilde{J}_m^{\,k+1} = J_m^{k+1} + E_m^{k+1} + \Delta_m^{k+1}$,with\
$J_m^{k+1}\in \partial f_m(x^{k+1},u_{m+1}^k)$, $m=1,\dots,M-1$, $J_M^{k+1}\in \partial f_M(x^{k+1})$,\
${\mathbbm{E}}\big\{E_m^{k+1}\big|{\mathcal{F}}_k\big\} = 0$, ${\mathbbm{E}}\big\{\|E_m^{k+1}\|^2|{\mathcal{F}}_k\big\}\le \sigma_E^2$, $\lim_{k\to \infty}  \Delta_m^{k+1}=0$;

For $m=2,\dots,M-1$, the error $ E_{mu}^{k+1} $ is conditionally independent of $E_\ell^{k+1}$ and $e_\ell^{k+1}$, $\ell=m+1,\dots,M$, given ${\mathcal{F}}_k$.

\[r:chain\] Assumption (A2) is satisfied if $f_m(\cdot, \cdot)$, $m=1,\dots,M-1$, are continuously differentiable, or, more generally, if $$f_m(x,u_{m+1}) = \varphi_m(\psi_m(x),u_{m+1}),\quad m= 1,\dots, M-1,$$ where $\varphi_m(\cdot,\cdot)$ is continuously differentiable, and $\psi_m(\cdot)$ is differentiable in a generalized sense. Indeed, by virtue of [@Rusz2019 Thm. 1], if a path $x(t)$ is continuously differentiable, the function $\psi_m(\cdot)$ admits the chain rule on this path. This implies that the path $\psi_m(x(t))$ is absolutely continuous. Consequently, $\varphi(\psi_m(x(t)),u_{m+1}(t))$ admits the chain rule.

Basic Properties {#s:4}
================

As the calculation of generalized subgradients of a composition requires the knowledge of the values $v_m$ of the inner functions, we will also need a more general multifunction $\varGamma:{\mathbbm{R}}^n\times{\mathbbm{R}}^n\times{\mathbbm{R}}^{d_{1}}\times \dots \times {\mathbbm{R}}^{d_{M}} \rightrightarrows {\mathbbm{R}}^n\times{\mathbbm{R}}^{d_{1}}\times \dots \times {\mathbbm{R}}^{d_{M}}$: $$\label{Gamma-new-2}
\begin{aligned}
\varGamma(x,z,u) &= \Big\{ (g,v_1,\dots,v_M): \exists D_m\in \partial f_m(x,u_{m+1}), \;\exists D_M\in \partial f_M(x),\  v_M = D_M,\\
& v_m = D_{mx}\big(\bar{y}(x,z)-x\big) + D_{mu} v_{m+1}+ b\big(  f_m(x,u_{m+1})- u_m\big), \; m=1,\dots,M-1,\\
& g_M = D_M,\;g_m = D_{mx} + D_{mu} g_{m+1},\;m=1,\dots,M-1,\; g = a (g_1^T-z) \Big\}.
\end{aligned}$$ The multifunction $\varGamma(\cdot)$ is convex and compact valued, due to assumptions (A2) and (A3).

These multifunctions allow us to write the iterations of $\{u^k_m\}$ and $\{z^k\}$ in a more compact way, and to establish the boundedness of these sequences. In our analysis we consider the $x$- and $u$-components of $\varGamma(\cdot)$ separately, by writing $g \in \varGamma_x(x,z,u)$ and $v_m\in \varGamma_{mu}(x,z,u)$, but it is essential to keep in mind that they all derive from the same multifunction.

\[l:uk-bounded\] If conditions [(A1)]{} and [(A6)–(A8)]{} are satisfied, then for all $k$ the following relations are true: $$\label{u-abstract}
u_m^{k+1} \in  u_m^k + \tau_k \varGamma_{mu}(x^{k+1},z^k,u^k) + \tau_k \theta_m^{k+1} + \tau_k \alpha_m^{k+1},$$ where, for some constant $C^{\theta}_m$, $$\label{theta-abstract}
{\mathbbm{E}}\big[\theta_m^{k+1}\,\big|\, {\mathcal{F}}_k\big]=0,  \quad {\mathbbm{E}}\big[\|\theta_m^{k+1}\|^2\,\big|\,{\mathcal{F}}_k\big]\le C^{\theta}_m, \quad  k=0,1,\dots$$ and $$\label{alpha-abstract}
\lim_{k\to \infty} \alpha_m^{k+1} = 0.\quad \text{a.s.}.$$ Moreover, the sequence $\{u^k\}$ is bounded a.s..

For $m=M$, formula and assumption (A7) yield: $$\label{ukM}
u_M^{k+1}  =  u_M^k + \tau_k v_M^k + \tau_k \theta_M^{k+1} + \tau_k \alpha_M^{k+1},$$ with $$\begin{gathered}
v_M^k =  J_M^{k+1} \big(\bar{y}(x^k,z^k)-x^k\big) + b \big( f_M(x^{k+1})- u_M^k\big),\\
\theta_M^{k+1} =  E_M^{k+1}\big(\bar{y}(x^k,z^k)-x^k\big) + b e_M^{k+1},\\
\alpha_M^{k+1} = \Delta_M^{k+1} \big(\bar{y}(x^k,z^k)-x^k\big) + b \delta_M^{k+1}.\end{gathered}$$ Due to (A1) and (A7), relations – are true for coordinate $m=M$.

To verify the boundedness of $\{u_M^k\}$, we define the quantities $$\label{tildeu}
\tilde{u}_M^k = u_M^k +  \sum_{j=k}^\infty \tau_j \theta_M^{j+1}.$$ Owing to (A6) and , by virtue of the martingale convergence theorem, the series in the formula above is convergent a.s., and thus $\tilde{u}_M^k - u_M^k \to 0$ a.s., when $k\to \infty$. We can now use to establish the following recursive relation: $$\label{tildeu-simple}
\tilde{u}_M^{k+1} = (1-b\tau_k) \tilde{u}_M^k + \tau_k J_M^{k+1} \big(\bar{y}(x^k,z^k)-x^k\big)
+ b \tau_k f_M(x^{k+1})+ \tau_k \alpha_M^{k+1} \\
 + b\tau_k (\tilde{u}_M^k -u_M^k).$$ By (A1), the sequences $\{J_M^k\}$ and $\{f_M(x^k)\}$ are bounded, and thus $$\limsup_{k\to \infty} \big\| {u}_M^k \big\| = \limsup_{k\to \infty} \big\| \tilde{u}_M^k \big\|
\le \limsup_{k\to \infty} \big\| J_M^{k+1} \big(\bar{y}(x^k,z^k)-x^k\big)
+ b  f_M(x^{k+1}) \big\| < \infty.$$ Furthermore, it follows from and for $M$ that a constant $C_u$ exists, such that ${\mathbbm{E}}\big[ \|u_m^k\|^2\| \big] \le C_u$.

We now proceed by induction. Suppose the relations - are true for $m+1$, the sequence $\{u^k_{m+1}\}$ is bounded a.s., and ${\mathbbm{E}}\big[ \|u_{m+1}^k\|^2\| \big]$ is bounded as well. We shall verify these properties for $m$. From for $m$ and for $m+1$ we obtain $$u_m^{k+1} = u_m^k + \tau_k \tilde{J}_m^{\,k+1} \begin{bmatrix}  \bar{y}(x^k,z^k)-x^k \\v_{m+1}^{k} + \theta_{m+1}^{k+1} + \alpha_{m+1}^{k+1}\end{bmatrix}
+ b \tau_k \big(\tilde{h}_m^{k+1}-u_m^k\big).$$ This can be rewritten as follows: $$u_m^{k+1} =  u_m^k +  \tau_k v^k_m + \tau_k \theta_M^{k+1} + \tau_k \alpha_M^{k+1},$$ with $$\begin{gathered}
v_m^k =  J_m^{k+1} \begin{bmatrix} \bar{y}(x^k,z^k)-x^k \\ v^k_{m+1} \end{bmatrix} + b\big(  f_m(x^{k+1},u^k_{m+1}) - u^k_m\big), \label{uk-ind}\\
\theta_m^{k+1} = E_m^{k+1}\begin{bmatrix} \bar{y}(x^k,z^k)-x^k \\ v^k_{m+1} \end{bmatrix}
+ \big(J_{mu}^{k+1} + E_{mu}^{k+1}\big) \theta_{m+1}^{k+1}
 + E_{mu}^{k+1} \alpha_{m+1}^{k+1} + b e_m^{k+1},
\label{theta-ind}\\
\alpha_m^{k+1} = J_{mu}^{k+1} \alpha_{m+1}^{k+1}
+ \Delta_m^{k+1}\begin{bmatrix} \bar{y}(x^k,z^k)-x^k \\ v^k_{m+1} +\alpha_{m+1}^{k+1} \end{bmatrix}  + b \delta_m^{k+1}. \label{alpha-ind}\end{gathered}$$ As the sequences $\{x^k\}$ and $\{u_{m+1}^k\}$ are bounded a.s., the Jacobians $J_m^{k+1}$ are bounded as well. Moreover, by assumption (A5), ${\mathbbm{E}}\big[\| J_m^{k+1}\|^2\big]$ is bounded for all $k$.

To verify –, we only need to analyze the effect of various product terms in the formulae above. The product $J_{mu}^{k+1} v^k_{m+1}$ in is square integrable due to (A4). The same is true for $f_m(x^{k+1},u^k_{m+1})$. The conditional expectation ${\mathbbm{E}}\big[\theta_m^{k+1}\,\big|\,{\mathcal{F}}_k\big] = 0 $, thanks to assumption (A8), because $\theta^{k+1}_{m+1}$ depends only on observations of quantities associated with the functions $f_\ell(\cdot,\cdot)$, $\ell=m+1,\dots,M$. Furthermore, ${\mathbbm{E}}\big[\|\theta_m^{k+1}\|^2\,\big|\,{\mathcal{F}}_k\big]$ is bounded, due to (A4), (A7), and (A8). Consequently, we can define a sequence $\{\tilde{u}_m^k\}$ is a way analogous to , and establish for it a recursive relation of the form , with $M$ replaced by $m$. In the same way as above, we obtain the boundedness of $\{{u}_m^k\}$ with probability 1, and its square integrability. By induction, the assertion is true for all $m$.

We now pass to the analysis of the sequence $\{z^k\}$. Carrying out for exact subgradients, we would obtain $$\label{g-fold-exact}
{g}_M^{k+1} = {J}_M^{k+1}, \quad {g}_{m}^{k+1} = J_{mx}^{k+1}+J_{mu}^{k+1}{g}_{m+1}^{k+1}, \quad m=M-1,\dots,1.$$ Evidently, $a \big(\big[g_1^{k+1}\big]^T-z^k\big) \in \varGamma_{x}\big(x^{k+1},z^k,u^k\big)$.

\[l:zk-recursion\] If conditions [(A1)]{} and [(A6)–(A8)]{} are satisfied, then for all $k=0,1,\dots$ and all $m=1,\dots,M$, $$\label{gradient-error}
\tilde{g}_m^{k+1} = g_m^{k+1} + \varkappa_m^{k+1} + \beta_m^{k+1}, \quad m=1,\dots,M,$$ where, for some constant $C_\varkappa$, $$\label{gradient-error-2}
{\mathbbm{E}}\big[\varkappa_m^{k+1} \,\big|\, {\mathcal{F}}_k\big] = 0, \quad {\mathbbm{E}}\big[\|\varkappa_m^{k+1}\|^2 \,\big|\, {\mathcal{F}}_k\big] \le C_\varkappa, \quad \lim_{k\to \infty} \beta_m^{k+1} = 0.$$

Formulae – are true for $m=M$ directly by (A7)(ii). Supposing that they are true for $m+1$, we verify them for $m$. We have $$\begin{gathered}
\label{gradient-error-3}
\tilde{g}_m^{k+1} - g_m^{k+1} =
 \tilde{J}_{mx}^{k+1} - J_{mx}^{k+1}
 +\tilde{J}_{mu}^{k+1}\tilde{g}_{m+1}^{k+1}
-J_{mu}^{k+1} {g}_{m+1}^{k+1}\\
 =
 \tilde{J}_{mx}^{k+1} - J_{mx}^{k+1}
 +\big(\tilde{J}_{mu}^{k+1} -J_{mu}^{k+1} \big)\tilde{g}_{m+1}^{k+1}
+ J_{mu}^{k+1} \big( \tilde{g}_{m+1}^{k+1} -{g}_{m+1}^{k+1} \big).\end{gathered}$$ The first term in , $\tilde{J}_{mx}^{k+1} - {J}_{mx}^{k+1}$, admits decomposition of form – directly by (A7)(ii). The second term, by (A7)(ii), can be represented as follows: $$\big(\tilde{J}_{mu}^{k+1} -{J}_{mu}^{k+1}\big) \tilde{g}_{m+1}^{k+1} =
\big({E}_{mu}^{k+1} + \Delta_{mu}^{k+1}\big) \big({g}_{m+1}^{k+1} + \varkappa_{m+1}^{k+1} + \beta_{m+1}^{k+1}\big).$$ Owing to (A8), $${\mathbbm{E}}\big\{  {E}_{mu}^{k+1} \varkappa_{m+1}^{k+1}\,\big|\, {\mathcal{F}}_k\big\} =0,\quad
{\mathbbm{E}}\big\{  \big\|{E}_{mu}^{k+1}  \varkappa_{m+1}^{k+1}\big\|^2 \,\big|\, {\mathcal{F}}_k\big\} \le C_\varkappa\sigma_E^2.$$ Together with the boundedness of $\{x^k\}$ and $\{u^k\}$, this implies that the second term admits decomposition od form –. The third term in , by , can be represented as follows: $${J}_{mu}^{k+1} \big( \tilde{g}_{m+1}^{k+1} -{g}_{m+1}^{k+1} \big)
 =   {J}_{mu}^{k+1} \big( \varkappa_{m+1}^{k+1} + \beta_{m+1}^{k+1}\big).$$ Together with the boundedness of $\{J_m^k\}$, this implies that the third term admits decomposition od form – as well. Therefore, – are true for all $m$.

We can now establish the boundedness of the sequence $\{z^k\}$.

\[l:zk-bounded\] If conditions [(A1)]{} and [(A6)–(A8)]{} are satisfied, then with probability 1 the sequence $\{z^k\}$ is bounded.

We proceed as in the proof of lemma \[l:uk-bounded\]. We define the quantities $$\tilde{z}^k  = z^k +  a\sum_{j=k}^\infty \tau_j \big[\varkappa_1^{j+1}\big]^T, \quad k=0,1,\dots.$$ Due to Lemma \[l:zk-recursion\], by virtue of the martingale convergence theorem, the series in the formula above is convergent a.s., and thus $\tilde{z}^k - z^k \to 0$ a.s., when $k\to \infty$. We can now use to establish the following recursive relation: $$\tilde{z}^{k+1} = (1-a\tau_k) \tilde{z}^k
+  a\tau_k \big[g_1^{k+1}\big]^T +  a\tau_k(\big[\beta_1^{k+1}\big]^T+\tilde{z}^k -z^k).$$ Therefore, $$\limsup_{k\to \infty} \big\| {z}^k \big\| = \limsup_{k\to \infty} \big\| \tilde{z}^k \big\| \le \limsup_{k\to \infty} \big\| g_1^k \big\| < \infty\quad \text{a.s.},$$ as claimed.

Convergence analysis {#s:analysis}
====================

We start from a useful property of the gap function [$\eta:X\times {\mathbbm{R}}^n\to(-\infty,0]$]{}, $$\label{gap}
\eta(x,z) = \min_{y\in X} \left\{\langle z, y-x \rangle + \frac{\rho}{2} \|y-x\|^2\right\}.$$ We denote the minimizer in by $\bar{y}(x,z)$. Since it is a projection of $x-z/\rho$ on $X$, $$\label{opt-eta}
 \langle z ,\bar{y}(x,z)-x \rangle  + \rho \| \bar{y}(x,z)-x \|^2 \le 0.$$ Moreover, a point $x^*\in X^*$ if and only if $z^*\in G_1(x^*)$ exists such that $\eta(x^*,z^*)=0$.

We can now state the main result of the paper.

\[t:convergence\] If the assumptions [(A1)–(A8)]{} are satisfied, then with probability 1 every accumulation point $\hat{x}$ of the sequence $\{x^k\}$ is stationary, $\lim_{k\to\infty} (u_m^k-F_m(x^k))=0$, for $m=1,\dots,M$, and the sequence $\{F_1(x^k)\}$ is convergent.

We consider a specific trajectory of the method and divide the proof into three standard steps.

*Step 1: The Limiting Dynamical System.* We denote by $p^k=(x^k,z^k,u^k)$, $k=0,1,2,\dots$, a realization of the sequence generated by the algorithm. We introduce the accumulated stepsizes $t_k = \sum_{j=0}^{k-1}\tau_j$, $k=0,1,2 \dots$, and we construct the interpolated trajectory $$P_0(t) = p^k + \frac{t-t_{k}}{\tau_k}(p^{k+1}-p^k),\quad t_{k}\le t \le t_{k+1},\quad k=0,1,2,\dots.$$ For an increasing sequence of positive numbers $\{s_k\}$ diverging to $\infty$, we define shifted trajectories $P_k(t) =P_0(t+s_k)$. The sequence $\{p^k\}$ is bounded by Lemmas \[l:uk-bounded\] and \[l:zk-bounded\].

By [@majewski2018analysis Thm. 3.2], for any infinite set ${\mathcal{K}}$ of positive integers, there exist an infinite subset ${\mathcal{K}}_1 \subset {\mathcal{K}}$ and an absolutely continuous function $P_\infty:
[0,+\infty) \to X\times {\mathbbm{R}}^n\times {\mathbbm{R}}^m$ such that for any $T > 0$ $$\lim_{\substack{{k\to\infty}\\{k\in {\mathcal{K}}_1}}}\sup_{t\in[0,T]} \big\| P_k(t)-P_\infty(t)\big\|= 0,$$ and $P_\infty(\cdot)=\big(X_\infty(\cdot),Z_\infty(\cdot),U_\infty(\cdot)\big)$ is a solution of the system of differential equations and inclusions corresponding to , with , and : $$\begin{gathered}
{\overset{\hbox{\,\tiny${}_\bullet$}}{x}}(t) = \bar{y}\big(x(t),z(t)\big)-x(t), \label{dx}\\
\big({\overset{\hbox{\,\tiny${}_\bullet$}}{z}}(t), {\overset{\hbox{\,\tiny${}_\bullet$}}{u}}(t)\big) \in \varGamma(x(t),z(t),u(t)). \label{dz}
$$ Moreover, for any $t\ge 0$, the triple $\big(X_\infty(t),Z_\infty(t),U_\infty(t)\big)$ is an accumulation point of the sequence $\{(x^k,z^k,u^k)\}$.

In order to analyze the equilibrium points of the system –, we first study the dynamics of the functions $\varPhi_m(t) =  f_m(X(t),U_{m+1}(t))$, $m =1,\dots,M-1$, and $\varPhi_M(t) = f_M(X(t))$. It follows from that the path $X(\cdot)$ is continuously differentiable. By virtue of assumption (A3), for any $J_M(t) \in \partial f_M(X(t))$, $$\label{f-incr-M}
 {\overset{\hbox{\,\tiny${}_\bullet$}}{\varPhi}}_M(t)
=   J_{M}(t) {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t) .$$ Assumption (A2) means that for any $J_m(t) \in \partial f_m(X(t),U_{m+1}(t))$, $$\label{f-incr-m-2}
 {\overset{\hbox{\,\tiny${}_\bullet$}}{\varPhi}}_m(t)
=   J_{m}(t) \begin{bmatrix} {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t)  \\  {\overset{\hbox{\,\tiny${}_\bullet$}}{U}}_{m+1}(t) \end{bmatrix}, \quad m=1,\dots,M-1.$$ We need to understand the dynamics of $U_m(\cdot)$. From and we deduce that $$\label{UM}
{\overset{\hbox{\,\tiny${}_\bullet$}}{U}}_M(t) = \hat{J}_M(t) {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t) + b [\varPhi_M(t) - U_M(t)],$$ with some $\hat{J}_M(t) \in \partial f_M(X(t))$, and $${\overset{\hbox{\,\tiny${}_\bullet$}}{U}}_m(t) = \hat{J}_m(t) \begin{bmatrix} {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t) \\ {\overset{\hbox{\,\tiny${}_\bullet$}}{U}}_{m+1}(t)\end{bmatrix}  + b[ \varPhi_m(t) - U_m(t)], \quad m=1,\dots,M-1,$$ with some $\hat{J}_m(t) \in \partial f_m(X(t),U_{m+1}(t))$. Therefore, using $J_m(\cdot) = \hat{J}_m(\cdot)$, for $m=1,\dots,M$, in , we obtain $$\label{Um}
{\overset{\hbox{\,\tiny${}_\bullet$}}{U}}_m(t) =  {\overset{\hbox{\,\tiny${}_\bullet$}}{\varPhi}}_m(t) + b[ \varPhi_m(t) - U_m(t)].$$ We can verify by induction that the solution of – has the form: $$\label{Phi-sol}
   {\overset{\hbox{\,\tiny${}_\bullet$}}{\varPhi}}_m(t) = \hat{g}_m(t){\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t)
   + b \sum_{\ell=m}^{M-1} \prod_{q=m}^\ell \hat{J}_{q u}(t)\big[\varPhi_{\ell +1}(t)-U_{\ell +1}(t)\big],
   \quad m=1,\dots,M,$$ with $\hat{g}_m(t)$ defined by the recursive procedure: $$\hat{g}_M(t) =  \hat{J}_{Mx}(t), \quad \hat{g}_m(t) =  \hat{J}_{mx}(t) + \hat{J}_{mu}(t) \hat{g}_{m+1}(t),\quad m=M-1,\dots,1,$$ and $\hat{J}_{mx}$ and $\hat{J}_{mu}$ denoting the $x$-part and the $u$-part of $\hat{J}_m$, respectively. These observations will help us study the stability of the system.

*Step 2: Descent Along a Path.* We use the Lyapunov function $$\label{Lyapunov}
W(x,z,u) = a f_1(x,u_2) -  \eta(x,z) + \sum_{m=2}^{M-1} \gamma_m\, \big\|f_m(x,u_{m+1})-u_m\big\| + \gamma_M \big\|f_M(x)-u_M\big\|,$$ with the coefficients $\gamma_m >0$ to be specified later.

Directly from for $m=1$ we obtain $$\begin{gathered}
\label{f-incr2}
  f_1(X(T),U_2(T)) - f_1(X(0),U_2(0)) \\
 = \int_0^T  \hat{g}_1(t) {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t)  \;dt
  + b  \sum_{\ell=1}^{M-1} \int_0^T \prod_{q=1}^\ell \hat{J}_{q u}(t)\big[\varPhi_{\ell +1}(t)-U_{\ell +1}(t)\big]\;dt.\end{gathered}$$

We now estimate the change of $\eta(X(\cdot),Z(\cdot))$ from 0 to $T$. Since $\bar{y}(x,z)$ is unique, the function $\eta(\cdot,\cdot)$ is continuously differentiable. Therefore, the chain formula holds for it as well: $$\begin{gathered}
\eta(X(T),Z(T)) - \eta(X(0),Z(0)) \\
= \int_0^T  \big\langle \nabla_x \eta(X(t),Z(t)), {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t) \big\rangle \;dt +
\int_0^T  \big\langle \nabla_z \eta(X(t),Z(t)), {\overset{\hbox{\,\tiny${}_\bullet$}}{Z}}(t) \big\rangle \;dt.\end{gathered}$$ From and , we obtain $${\overset{\hbox{\,\tiny${}_\bullet$}}{Z}}(t) = a\big(\hat{g}_1^T(t) - Z(t)\big),$$ with the same $\hat{g}_1(\cdot)$ as in for $m=1$ and in . Substituting $\nabla_x \eta(x,z) = -z+\rho(x- \bar{y}(x,z))$, $\nabla_z \eta(x,z) =  \bar{y}(x,z)-x$, and using , we obtain $$\begin{aligned}
\lefteqn{\eta(X(T),Z(T)) - \eta(X(0),Z(0))}\quad  \\
&= \int_0^T  \big\langle -Z(t)+\rho(X(t)-\bar y(X(t),Z(t)))\, , \,\bar{y}(X(t),Z(t)) - X(t) \big\rangle \;dt \\
&{\quad } +
a \int_0^T  \big\langle \bar y(X(t),Z(t))-X(t)\,,\,\hat{g}^T_1(t) - Z(t) \big\rangle \;dt\\
&\ge\; a \int_0^T  \big\langle \bar y(X(t),Z(t))-X(t)\,,\, \hat{g}_1^T(t) - Z(t) \big\rangle \;dt\\
&\ge\; a \int_0^T  \hat{g}_1(t)\big( \bar y(X(t),Z(t))-X(t)\big) \;dt
 + a\rho \int_0^T  \big\| \bar y(X(t),Z(t))-X(t)\big\|^2 \;dt.
   \end{aligned}$$ With a view at , we conclude that $$\label{eta-incr}
 \eta(X(T),Z(T)) - \eta(X(0),Z(0))
  \ge \;  a \int_0^T  \hat{g}_1(t) {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t) \;dt
 + a\rho \int_0^T  \big\| {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t)\big\|^2 \;dt.$$ We now estimate the increment of $\big\| \varPhi_m(\cdot)-U_m(\cdot)\big\|$ from 0 to $T$. As $\| \cdot\| $ is convex and $\varPhi_m(\cdot)$ and $U_m(\cdot)$ are paths, the chain rule applies as well: for any $\lambda_m(t) \in \partial \| \varPhi_m(t)- U_m(t)\|$ we have $$\big\|\varPhi_m(T)-U_m(T)\big\| -  \big\|\varPhi_m(0)-U_m(0)\big\|
= \int_0^T \big\langle \lambda_m(t),  {\overset{\hbox{\,\tiny${}_\bullet$}}{\varPhi}}_m(t) - {\overset{\hbox{\,\tiny${}_\bullet$}}{U}}_m(t) \big\rangle\; dt.$$ By , ${\overset{\hbox{\,\tiny${}_\bullet$}}{\varPhi}}_m(t) - {\overset{\hbox{\,\tiny${}_\bullet$}}{U}}_m(t) = b\big[ U_m(t)-\varPhi_m(t)\big]$ for almost all $t$. Furthermore, $$\lambda_m(t)=\frac{\varPhi_m(t)-U_m(t)}{\|\varPhi_m(t)-U_m(t)\|},\quad \text{if}\quad  \varPhi_m(t)\ne U_m(t).$$ Therefore $$\label{norm-incr}
 \big\|\varPhi_m(T)-U_m(T)\big\| -  \big\|\varPhi_m(0)-U_m(0)\big\| = - b \int_0^T \big\| \varPhi_m(t)-U_m(t)\big\| \; dt.$$

We can now combine , , and to estimate the change of the function : $$\begin{gathered}
W\big(X(T),Z(T),U(T)\big) - W\big(X(0),Z(0),U(0)\big)\quad \\
\quad \le   a  b  \sum_{m=1}^{M-1} \int_0^T \prod_{q=1}^m \hat{J}_{q u}(t)\big[\varPhi_{m +1}(t)-U_{m +1}(t)\big]\;dt\\
 - a\rho \int_0^T  \big\| {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t)\big\|^2 \;dt
- b\sum_{m=2}^{M} \gamma_m  \int_0^T \big\|\varPhi_{m}(t)-U_{m}(t)\big\|\;dt.\end{gathered}$$ Because the paths $X(t)$ and $U(\cdot)$ are bounded a.s. and the functions $f_m$ are locally Lipschitz, a (random) constant $L$ exists, such that $\big \| J_{m u}(t)\big\| \le L$ for $m=1,\dots,M-1$. The last estimate entails: $$\begin{gathered}
\label{W-incr}
W\big(X(T),Z(T),U(T)\big) - W\big(X(0),Z(0),U(0)\big) \\
 \le - a\rho \int_0^T  \big\| {\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(t)\big\|^2 \;dt - b \sum_{m=2}^{M} (\gamma_m-aL^{m-1})  \int_0^T \|\varPhi_m(t)-U_m(t)\|\;dt.\end{gathered}$$ By choosing $\gamma_m > aL^{m-1}$ for all $m$, we ensure that $W(\cdot)$ has the descent property to be used in our stability analysis at Step 3. The fact that $L$ (and thus $\gamma_m$) may be different for different paths is irrelevant, because our analysis is path-wise.

*Step 3: Analysis of the Limit Points.* Define the set $$\begin{gathered}
\quad {\mathcal{S}}= \big\{ (x,z,u)\in X^* \times {\mathbbm{R}}^n\times {\mathbbm{R}}^m: \eta(x,z)=0,\\
 u_m=f_m(x,u_{m+1}), \ m=1,\dots,M-1,\ u_M=f_M(x)\big\}.\quad\end{gathered}$$ Suppose $(\bar{x},\bar{z},\bar{u})$ is an accumulation point of the sequence $\{(x^k,z^k,u^k)\}$. If $\eta(\bar{x},\bar{z}) <0$ or $\bar{u}_m\ne f_m(\bar{x},\bar{u}_{m+1})$, then every solution $(X(t),Z(t),U(t))$ of the system –, starting from $(\bar{x},\bar{z},\bar{u})$ has $\|{\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(0)\| > 0$ or $\|\varPhi_m(0)-U_m(0)\| >0$. Using and arguing as in [@duchi2018stochastic Thm. 3.20] or [@majewski2018analysis Thm. 3.5], we obtain a contradiction. Therefore, we must have $\eta(\bar{x},\bar{z})=0$ and $\bar{u}_m= f_m(\bar{x},\bar{u}_{m+1})$, $m=1,\dots,M-1$, and $\bar{u}_M= f_M(\bar{x})$. Suppose $\bar{x}\not \in X^*$. Then $$\label{non-opt}
{\mathop{\rm dist}}\big(0,  G_1(\bar{x}) + N_X(\bar{x})\big) >0.$$ Suppose the system – starts from $(\bar{x},\bar{z},\bar{u})$ and $X(t)=\bar{x}$ for all $t\ge 0$. From and , in view of the equations $\bar{y}(\bar{x},\bar{z})=\bar{x}$ and $\bar{u}_M=f_M(\bar{x})$, we obtain $U_M(t)=  f_M(\bar{x})$ for all $t \ge 0$. Proceeding by backward induction, we conclude that $U_m(t)=  f_m(\bar{x},\bar{u}_{m+1})$ for all $t \ge 0$ and all $m=1,\dots,M-1$. The inclusion , in view of , simplifies $${\overset{\hbox{\,\tiny${}_\bullet$}}{z}}(t) \in a \big( G_1(\bar{x}) - z(t)\big).$$ For the convex Lyapunov function $V(z) = {\mathop{\rm dist}}\big(z,G_1(\bar{x})\big)$, we apply the classical chain formula [@brezis1971monotonicity] on the path $Z(\cdot)$: $$V((Z(T)) - V(Z(0)) = \int_0^T \big\langle \partial V(Z(t)), {\overset{\hbox{\,\tiny${}_\bullet$}}{Z}}(t)\big\rangle \;dt.$$ For $Z(t)\notin G_1(\bar{x})$, we have $$\partial V(Z(t)) = \frac{Z(t) - {{\rm Proj}}_{G_1(\bar{x})}(Z(t))}{ \| Z(t) - {{\rm Proj}}_{G_1(\bar{x})}(Z(t))\|}$$ and ${\overset{\hbox{\,\tiny${}_\bullet$}}{Z}}(t) = a (d(t) - Z(t))$ with some $d(t)\in G_1(\bar{x})$. Therefore, $$\big\langle \partial V(Z(t)), {\overset{\hbox{\,\tiny${}_\bullet$}}{Z}}(t)\big\rangle \le  - a \| Z(t) - {{\rm Proj}}_{G_1(\bar{x})}(Z(t))\| = -a V(Z(t)).$$ It follows that $$V((Z(T)) - V(Z(0)) \le - a \int_0^T  V(Z(t)) \;dt,$$ and thus $$\label{Z-conv}
\lim_{t\to\infty} {\mathop{\rm dist}}\big(Z(t),G_1(\bar{x})\big) = 0.$$ It follows from – that $T>0$ exists, such that $ -Z(T) \not \in N_X(\bar{x})$, which yields ${\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(T)\ne 0$. Consequently, the path $X(t)$ starting from $\bar{x}$ cannot be constant (our supposition made right after cannot be true). But if is not constant, then again $T>0$ exists, such that ${\overset{\hbox{\,\tiny${}_\bullet$}}{X}}(T)\ne 0$. By Step 1, the triple $(X(T),Z(T),U(T))$ would have to be an accumulation point of the sequence $\{(x^k,z^k,u^k)\}$, a case already excluded. We conclude that every accumulation point $(\bar{x},\bar{z},\bar{u})$ of the sequence $\{(x^k,z^k,u^k)\}$ is in ${\mathcal{S}}$. The convergence of the sequence $\big\{W(x^k,z^k,u^k)\big\}$ then follows in the same way as [@duchi2018stochastic Thm. 3.20] or [@majewski2018analysis Thm. 3.5]. As $\eta(x^k,z^k)\to 0$, the convergence of $\{f_1(x^k,u_2^k)\}$ follows as well. Since $f_m(x^k)-u_m^k \to 0$, $m=2,\dots,M$, the sequence $\{F_1(x^k)\}$ is convergent as well.

Adapting the proof of Lemma \[l:zk-bounded\], we obtain convergence of path-averaged stochastic subgradients.

\[c:conv-zk\] If the sequence $\{x^k\}$ is convergent to a single point $\bar{x}$, then every accumulation point of $\{z^k\}$ is an element of the generalized gradient $G_{F_1}(\bar{x})$ satisfying the optimality condition at $\bar{x}$.

In fact, if we introduced path-averaging of the vectors $\tilde{g}_{m}^{k}$ used in , in a way similar to : $$z_m^{k+1} = z_m^k + a\tau_k \Big( \big[\tilde{g}_{m}^{k+1}\big]^T  - z_m^k\Big),  \quad m=1,\dots,M,$$ we could extend Corollary \[c:conv-zk\] to the convergence of each $\{z_m^k\}$ to an element of the generalized gradient $G_{F_m}(\bar{x})$.

Fixed Stepsize Performance {#s:5}
==========================

We are now interested in performance guarantees after finitely many steps with a constant stepsize. We make stronger assumptions about the functions $f_m$ in problem and about the estimation bias in (A6):

(A9)

:   The functions $f_m(\cdot,\cdot)$, $m=1,\dots,M-1$, and $f_M(\cdot)$ are continuously differentiable with Lipschitz continuous derivatives;

(A10)

:   A constant $C_\delta>0$ exists, such that $\max\big(\|\delta_m^{k+1}\|,\|\Delta_m^{k+1}\|\big) \le C_\delta \tau_k$ with probability 1, for all $k$ and $m=1,\dots,M$;

(A11)

:   A constant $L$ exists, such that $\| {J}_{mu}^{k+1}\| \le L$ with probability 1, for all $k=0,1,2\dots$.

In (A10) and (A11) we need $C_\delta$ and $L$ to be the same for all paths, because we need one Lyapunov function for all paths. Assumptions (A9)–(A11) imply (A2)–(A4). We also do not need assumptions (A5) and (A6), since we are interested in finitely many iterations.

Our reasoning is similar to the line of argument in [@ghadimi2018single Sec. 3], albeit with the multi-level nested structure and with the new Lyapunov function .

By using assumptions (A7)–(A9), and adapting to discrete time our analysis from Step 2 of the proof of Theorem \[t:convergence\], we estimate the three terms of the difference $$\begin{gathered}
{\mathbbm{E}}\big[W(x^{k+1},z^{k+1},u^{k+1})\,\big|\, {\mathcal{F}}_k\big] - W(x^k,z^k,u^k)\\
=
a \big( {\mathbbm{E}}\big[f_1(x^{k+1},u_2^{k+1})\,\big|\, {\mathcal{F}}_k\big] - f_1(x^k,u_2^k)\big)
 - \big( {\mathbbm{E}}\big[\eta(x^{k+1},z^{k+1})\,\big|\, {\mathcal{F}}_k\big] - \eta(x^k,z^k) \big)\\
 +
\sum_{m=2}^{M-1} \gamma_m \big( {\mathbbm{E}}\big[ \|f_m(x^{k+1},u_{m+1}^{k+1})-u_m^{k+1} \| \,\big|\, {\mathcal{F}}_k\big] - \|f_m(x^{k},u_{m+1}^{k})-u_m^{k} \|\big) \\
+ \gamma_M \big( {\mathbbm{E}}\big[ \|f_M(x^{k+1})-u_M^{k+1} \| \,\big|\, {\mathcal{F}}_k\big] - \|f_M(x^{k})-u_M^{k} \|\big).\end{gathered}$$ In our considerations below, $C$ is a sufficiently large deterministic constant dependent on various constants featuring in the problem definition, such as Lipschitz constants, variances of errors, etc.

Denote $d^k=y^k-x^k$. By a discrete-time version of the argument that lead to , the decrease of the first part of the Lyapunov function can be estimated as follows: $$\label{f-incr-dis}
\begin{aligned}
\lefteqn{ {\mathbbm{E}}\big[ f_1(x^{k+1},u_2^{k+1})\,\big|\,{\mathcal{F}}_k\big]  - f_1(x^{k},u_2^{k}) }\qquad\\
 &\le \tau_k \hat{g}_1^{k+1} d^k
  + b \tau_k \sum_{m=1}^{M-1} \prod_{q=1}^m \hat{J}_{q u}^{k+1}\big(\varPhi_{m +1}^k-U_{m +1}^k\big)
  + C \tau_k^2 \\
  &\le \tau_k \hat{g}_1^{k+1} d^k
  + b \tau_k \sum_{m=2}^{M} L^{m-1} \big\|\varPhi_{m}^k-U_{m}^k\big\|
  + C \tau_k^2 .
\end{aligned}$$ Due to (A6) and the continuous differentiability of the function $\eta(\cdot,\cdot)$, the change of the second part of the Lyapunov function can be esitmated as follows (the inequalities are derived from , as in ): $$\label{eta-incr-dis}
 {\mathbbm{E}}\big[ \eta(x^{k+1},z^{k+1}) - \eta(x^{k},z^{k}))
  \ge \; a \tau_k    \hat{g}_1^{k+1} d^k   + a\rho \tau_k  \big\| d^k \big\|^2 - C\tau_k^2.$$ The remaining part of the Lyapunov function requires more attention, because it is a sum of nonsmooth functions. We first analyze the term associated with $f_M(\cdot)$ to expose our ideas. Expanding the square, we obtain $$\begin{gathered}
\label{H-first}
\big\|f_M(x^{k+1})-u_M^{k+1}\big\|^2 \\
= \big\|f_M(x^{k})-u_M^{k}\big\|^2 + 2 \big\langle f_M(x^{k})-u_M^{k} ,  f_M(x^{k+1})-u_M^{k+1} - f_M(x^{k}) + u_M^{k} \big\rangle \\
 + \big\|f_M(x^{k+1})-u_M^{k+1} - f_M(x^{k}) + u_M^{k}\big\|^2.\end{gathered}$$ Due to assumption (A9), $$f_M(x^{k+1})- f_M(x^k) = \tau_k J_M^{k+1} d^k + r_M^{k+1},$$ with ${\mathbbm{E}}\big[\|r_M^{k+1}\| \,\big|\, {\mathcal{F}}_k\big] \le C \tau_k^2$. By , $$\begin{aligned}
u_M^{k+1} - u_M^k &= \tau_k \tilde{J}^{\,k+1}d^k + b \tau_k \big(\tilde{h}_M^{k+1}-u_M^k\big) \\
&= \tau_k {J}^{k+1}d^k + \tau_k E_M^{k+1}d^k + \tau_k \Delta_M^{k+1} d^k + b \tau_k (f_M(x^{k+1}) - u^k) + b \tau_k e_M^{k+1} + b\tau_k \delta_M^{k+1}.\end{aligned}$$ Therefore $$\begin{aligned}
\lefteqn{f_M(x^{k+1})- f_M(x^k) - u_M^{k+1} + u_M^k}\\
&= -b\tau_k (f_M(x^{k+1}) - u_M^k)- \tau_k E_M^{k+1}d^k - \tau_k \Delta^{k+1} d^k - b \tau_k e_h^{k+1} - b\tau_k \delta_h^{k+1}+  r_M^{k+1} \\
&= -b\tau_k (f_M(x^{k}) - u_M^{k})  - \tau_k E_M^{k+1}d^k  - b \tau_k e_h^{k+1}+  r^{k+1},\end{aligned}$$ where $\|r^{k+1}\| \le C \tau_k^2$. Then $$\big\| f_M(x^{k+1})- f_M(x^k) - u_M^{k+1} + u_M^k\big\|^2 \le 2b \tau_k^2 \big\| f_M(x^{k}) - u_M^{k}\big\|^2 + 2 C \tau_k^2.$$ We can now substitute the last two expressions into the scalar product in and take the conditional expectation of both sides with respect to ${\mathcal{F}}_k$: $${\mathbbm{E}}\big[ \|f_M(x^{k+1})-u_M^{k+1}\|^2 \,\big|\, {\mathcal{F}}_k\big]  \le (1- b \tau_k +2 b^2\tau_k^2)\|f_M(x^{k})-u_M^{k}\|^2   + C \tau_k^2.$$ The Jensen inequality yields, $$\label{H-last}
\begin{aligned}
\lefteqn{{\mathbbm{E}}\big[ \|f_M(x^{k+1})-u_M^{k+1}\| \,\big|\, {\mathcal{F}}_k\big]
 \le  \Big( {\mathbbm{E}}\big[ \|f_M(x^{k+1})-u_M^{k+1}\|^2 \,\big|\, {\mathcal{F}}_k\big]\Big)^{\frac{1}{2}}  } \qquad\\
&\le \big( (1- b \tau_k+2 b^2\tau_k^2)\|f_M(x^{k})-u_M^{k}\|^2   + C \tau_k^2 \big)^{\frac{1}{2}} \\
&\le \sqrt{1- b \tau_k+2 b^2\tau_k^2}\; \|f_M(x^{k})-u_M^{k}\|   + \sqrt{C} \tau_k\\
&\le
\Big( 1- \frac{1}{2} b \tau_k + b^2\tau_k^2\Big)\|f_M(x^{k})-u_M^{k}\|   + \sqrt{C} \tau_k.
\end{aligned}$$ For $1 \le m < M$ the analysis follows the same steps. The expansion is the same: $$\begin{gathered}
\big\|f_m(x^{k+1},u_{m+1}^{k+1})-u_M^{k+1}\big\|^2 =  \big\|f_m(x^{k},u_{m+1}^k)-u_m^{k}\big\|^2 \\
 +
2 \big\langle f_m(x^{k},u_{m+1}^k)-u_m^{k} ,  f_m(x^{k+1},u_{m+1}^{k+1})-u_m^{k+1} - f_m(x^{k},u_{m+1}^k) + u_M^{k} \big\rangle \\
 + \big\|f_m(x^{k+1},u_{m+1}^{k+1})-u_m^{k+1} - f_m(x^{k},u_{m+1}^k) + u_M^{k}\big\|^2.\end{gathered}$$ The only difference is the presence of the arguments $u_{m+1}^{k}$ and $u_{m+1}^{k+1}$. Due to assumption (A9), $$f_m(x^{k+1},u_{m+1}^{k+1})- f_m(x^{k},u_{m+1}^k) = \tau_k J_{mx}^{k+1} d^k + J_{mu}^{k+1}(u_{m+1}^{k+1}-u_{m+1}^{k}) + r_m^{k+1}.$$ The additional term, $J_{mu}^{k+1}(u_{m+1}^{k+1}-u_{m+1}^{k})$ is compensated (with stochastic errors) by the corresponding term in the update of $u_m$ by : $$u_m^{k+1} - u_m^k = \tau_k \tilde{J}_{mx}^{\,k+1} d^k + \tilde{J}_{mu}^{\,k+1} (u_{m+1}^{k+1}-u_{m+1}^{k})  + b \tau_k \big(\tilde{h}_m^{k+1}-u_m^k\big).$$ The remaining steps are identical and lead to the estimate corresponding to : $$\label{Hm-last}
{\mathbbm{E}}\big[ \|f_m(x^{k+1},u_{m+1}^{k+1})-u_M^{k+1}\| \,\big|\, {\mathcal{F}}_k\big]
 \le
\Big( 1- \frac{1}{2} b \tau_k + b^2\tau_k^2\Big)\|f_m(x^{k},u_{m+1}^k)-u_m^{k}\|   + \sqrt{C} \tau_k.$$ Integrating , , , and , we obtain the following inequality, corresponding to (we write $\varPhi_m^k$ for $f_m(x^k,u_{m+1}^k)$ and $\varPhi_M^k$ for $f_M(x^k)$): $$\begin{gathered}
\label{W-discrete}
{\mathbbm{E}}\big[ W(x^{k+1},z^{k+1},u^{k+1})\,\big|\,{\mathcal{F}}_k\big]  - W(x^{k},z^{k},u^{k}) \le -a\rho \tau_k\| y^k-x^k\|^2  \\
- b \tau_k \sum_{m=2}^{M} \Big( \frac{\gamma_m}{2} - b \gamma_m\tau_k - aL^{m-1}\Big) \|\varPhi_m^{k}-u_m^{k}\|
 + C \tau_k^2 + \sqrt{C}\sum_{m=2}^M \gamma_m\tau_k,\end{gathered}$$ where $C$ is a constant independent of $\gamma_m$.

Suppose we carry out $N$ iterations with constant stepsize $\tau$. In , we can always assume that $b\tau < \frac{1}{4}$ and $\frac{\gamma_m}{4} > aL^{m-1} + c$, with some $c>0$. Setting $\bar{\gamma} =\sum_{m=2}^M \gamma_m$, we obtain $$\begin{gathered}
\label{W-discrete-2}
{\mathbbm{E}}\big[ W(x^{k+1},z^{k+1},u^{k+1})\,\big|\,{\mathcal{F}}_k\big]  - W(x^{k},z^{k},u^{k})\\
 \le -a\rho \tau\| y^k-x^k\|^2
- b c \tau \sum_{m=2}^{M}\|\varPhi_m^{k}-u_m^{k}\|
 + C \tau^2 + \sqrt{C}\bar{\gamma}\tau.\end{gathered}$$ Taking the expected value of both sides and summing from 0 to $N-1$, we get $$\begin{gathered}
{\mathbbm{E}}\big[ W(x^{N},z^{N},u^{N})\big]   - W(x^{0},z^{0},u^{0}) \\
\le
- \pi \tau  \sum_{k=0}^{N-1} {\mathbbm{E}}\bigg[ \| y^k-x^k\|^2 + \sum_{m=2}^{M}\|\varPhi_m^{k}-u_m^{k}\| \bigg]  + NC \tau^2 + N\sqrt{C}\bar{\gamma}\tau,\end{gathered}$$ where $\pi=  \min(a\rho,bc)>0$. By construction, $W(x^N,z^N,u^N) \ge a\,\min_{x\in X} F(x)$, and thus $$\begin{gathered}
 \frac{1}{N} \sum_{k=0}^{N-1} {\mathbbm{E}}\bigg[ \| y^k-x^k\|^2 + \sum_{m=2}^{M}\|\varPhi_m^{k}-u_m^{k}\| \bigg] \\
 \le \frac{1}{N\pi\tau} \left(W(x^{0},z^{0},u^{0}) - a\,\min_{x\in X} F(x)\right)  +  \frac{C \tau^2}{\pi}
 + \frac{\sqrt{C}\bar{\gamma}\tau}{\pi}.\end{gathered}$$ Setting $\tau = \alpha/\sqrt{N}$ we obtain $$\begin{gathered}
 \frac{1}{N} \sum_{k=0}^{N-1} {\mathbbm{E}}\bigg[ \| y^k-x^k\|^2 + \sum_{m=2}^{M}\|\varPhi_m^{k}-u_m^{k}\| \bigg] \\
\le \frac{1}{\pi\alpha\sqrt{N}} \left(W(x^{0},z^{0},u^{0}) - a\,\min_{x\in X} F(x)\right)  + \frac{C\alpha^2}{\pi N}
 + \frac{ \sqrt{C} \bar{\gamma} \alpha}{\pi\sqrt{N}}.\end{gathered}$$ Therefore, at a random iteration number $R$, drawn from the uniform distribution in $\{0,1,\break\dots, N-1\}$, independently from other random quantities in the method, we have $$\label{rate}
{\mathbbm{E}}\bigg[ \| y^R-x^R\|^2 + \sum_{m=2}^{M-1}\|f_m(x^{R},u_{m+1}^R)-u_m^{R}\| + \|f_M(x^{R})-u_M^{R}\|\bigg] \le \frac{\text{Const}}{\sqrt{N}}.$$ If we consider the quantity on the left hand side of as a measure of non-optimality of the triple $(x^R,z^R,u^R)$, we can deduce that $\mathcal{O}(1/\varepsilon^2)$ iterations are needed to achieve an error of size $\varepsilon$. The reference [@ghadimi2018single Sec. 2] provides a thorough discussion of similar optimality measures used for constrained stochastic optimization. In fact, compared to the measure used in [@ghadimi2018single], we see from that the “tracking error” $\sum_{m=2}^{M-1}\|f_m(x,u_{m+1})-u_m\| + \|f_M(x)-u_M\|$ decreases, on average, an order of magnitude faster than the “stationarity error” $\| \bar{y}(x,z)- x\|$. This reinforces our point that multiple time scales are not needed in this setting.

[10]{}

P. Artzner, F. Delbaen, J.-M. Eber, and D. Heath. Coherent measures of risk. , 9:203–228, 1999.

M. Bena[ï]{}m, J. Hofbauer, and S. Sorin. Stochastic approximations and differential inclusions. , 44(1):328–348, 2005.

M. Bena[ï]{}m, J. Hofbauer, and S. Sorin. Stochastic approximations and differential inclusions, part ii: Applications. , 31(4):673–695, 2006.

J. Bolte, A. Daniilidis, A. Lewis, and M. Shiota. Clarke subgradients of stratifiable functions. , 18(2):556–572, 2007.

V. S. Borkar. . Springer, New York, 2009.

H. Br[é]{}zis. Monotonicity methods in [H]{}ilbert spaces and some applications to nonlinear partial differential equations. In [*Contributions to Nonlinear Functional Analysis*]{}, pages 101–156. Elsevier, 1971.

F. H. Clarke. Generalized gradients and applications. , 205:247–262, 1975.

F. H. Clarke. Generalized gradients of [L]{}ipschitz functionals. , 40(1):52–67, 1981.

D. Davis and D. Drusvyatskiy. Stochastic model-based minimization of weakly convex functions. , 29(1):207–239, 2019.

D. Davis, D. Drusvyatskiy, S. Kakade, and J. D. Lee. Stochastic subgradient method converges on tame functions. , pages 1–36, 2018.

D. Dentcheva, S. Penev, and A. Ruszczy[ń]{}ski. Statistical estimation of composite risk functionals and risk optimization problems. , 69(4):737–760, 2017.

D. Drusvyatskiy, A. D. Ioffe, and A. S. Lewis. Curves of descent. , 53(1):114–138, 2015.

J. C. Duchi and F. Ruan. Stochastic methods for composite and weakly convex optimization problems. , 28(4):3229–3259, 2018.

Yu. M. Ermoliev. A general stochastic programming problem. , 1(4):106–112, 1971.

Yu. M. Ermoliev. . Nauka, Moscow, 1976.

Yu. M. Ermoliev and V. I. Norkin. Sample average approximation method for compound stochastic optimization problems. , 23(4):2231–2263, 2013.

H. F[ö]{}llmer and A. Schied. . Walter de Gruyter, 2011.

S. Ghadimi and G. Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic programming. , 23(4):2341–2368, 2013.

S. Ghadimi, A. Ruszczy[ń]{}ski, and M. Wang. A single time-scale stochastic approximation method for nested stochastic optimization. , 2020. Accepted for publication (arXiv preprint 1812.01094).

A. M. Gupal. . Naukova Dumka, Kiev, 1979.

A. N. Iusem, A. Jofr[é]{}, R. I. Oliveira, and P. Thompson. Extragradient method with variance reduction for stochastic variational inequalities. , 27(2):686–724, 2017.

J. Koshal, A. Nedic, and U. V. Shanbhag. Regularized iterative stochastic approximation methods for stochastic variational inequality problems. , 58(3):594–609, 2013.

H. Kushner and G. G. Yin. . Springer, New York, 2003.

H. J. Kushner and D. S. Clark. . Springer, New York, 1978.

L. Ljung. Analysis of recursive stochastic algorithms. , 22(4):551–575, 1977.

S. Majewski, B. Miasojedow, and E. Moulines. Analysis of nonsmooth stochastic approximation: the differential inclusion approach. , 2018.

R. Mifflin. Semismooth and semiconvex functions in constrained optimization. , 15(6):959–972, 1977.

V. S. Mikhalevich, A. M. Gupal, and V. I. Norkin. . Nauka, Moscow, 1987.

V. I. Norkin. Generalized-differentiable functions. , 16(1):10–12, 1980.

V. I. Norkin. . PhD thesis, Institute of Cybernetics of the Ukrainian Academy of Sciences, Kiev, 1983.

E. A. Nurminski. . Naukova Dumka, Kiev, 1979.

W. Ogryczak and A. Ruszczyński. From stochastic dominance to mean–risk models: semideviations as risk measures. , 116:33–50, 1999.

W. Ogryczak and A. Ruszczyński. On consistency of stochastic dominance and mean–semideviation models. , 89:217–232, 2001.

A. Ruszczy[ń]{}ski. A method of feasible directions for solving nonsmooth stochastic programming problems. In F. Archetti, G. Di Pillo, and M. Lucertini, editors, [ *Stochastic Programming*]{}, pages 258–271. Springer, 1986.

A. Ruszczy[ń]{}ski. A linearization method for nonsmooth stochastic programming problems. , 12(1):32–49, 1987.

A. Ruszczyński. Convergence of a stochastic subgradient method with averaging for nonsmooth nonconvex constrained optimization. , 2020 (in press). arXiv preprint arXiv:1912.07580.

A. Shapiro, D. Dentcheva, and A. Ruszczy[ń]{}ski. . SIAM, Philadelphia, 2009.

M. Wang, E. X. Fang, and B. Liu. Stochastic compositional gradient descent: Algorithms for minimizing compositions of expected-value functions. , 161(1-2):419–449, 2017.

M. Wang, J. Liu, and E. X. Fang. Accelerating stochastic composition optimization. , 18:1–23, 2017.

S Yang, M. Wang, and E. X. Fang. Multilevel stochastic gradient methods for nested composition optimization. , 29(1):616–659, 2019.

[^1]: Rutgers University, Department of Management Science and Information Systems, Piscataway, NJ 08854, USA; email: [rusz[@]{}rutgers.edu]{}
