Huber Penalty Approach to Problems with Linear ConstraintsA. Nedić and T. Tatarenko

Huber Loss-Based Penalty Approach to Problems with Linear ConstraintsSubmitted to the editors on ......

Angelia NedićThe School of Electrical, Computer and Energy Engineering at Arizona State University, USA
(Angelia.Nedich@asu.edu). A. Nedić gratefully acknowledges support of this work by the
Office of Naval Research grant no. N00014-12-1-0998.
Tatiana TatarenkoDepartment of Control Theory and Robotics, TU Darmstadt, Germany
(tatiana.tatarenko@tu-darmstadt.de). T. Tatarenko gratefully acknowledges support of this work by the German Research Foundation
(Deutsche Forschungsgemeinschaft, DFG) within the SPP 1984 “Hybrid and
multimodal energy systems: System theoretical methods for the transformation
and operation of complex networks”.

We consider an optimization problem with many linear inequalities constraints.
To deal with a large number of constraints, we provide a
penalty reformulation of the problem, where the penalty is a variant of the one-sided Huber loss function.
We study the infeasibility properties of the solutions of penalized problems for nonconvex and convex objective functions, as the penalty parameters vary with time.
Then, we propose a random incremental penalty method for solving convex problem,
and investigate its convergence properties for convex and strongly convex objective functions. We show that the iterates of the method converge to a solution of the original problem almost surely and in expectation for suitable choices of the penalty parameters and the stepsize. Also, we establish convergence rate of the method for appropriately defined weighted averages of the iterates for the expected function values.
We establish  $O(\ln^{1/2+\e} k/{\sqrt k})$-convergence rate when the objective function is convex, and  $O(\ln^{\e} k/k)$-convergence rate when the objective function is strongly convex with $\e$ being any small positive number.
To the best of our knowledge, these are the first results on the convergence rate for the penalty-based incremental subgradient method with time-varying penalty parameters.

Convex minimization, linear constraints, inexact penalty, incremental methods

90C25, 90C06, 65K05

§ INTRODUCTION

In  this  paper,  we  study  the  problem  of minimizing  a
convex function $f:\R^n\to\R$ over a convex and closed set $X$
that is the intersection of finitely many sets $X_i$, represented by linear inequalities, $i=1,\ldots,m,$
where $m\ge 2$ is large,

min   f(x),   s.t.    x∈X = ∩_i=1^m X_i.

Optimization problems of the form (<ref>) arise in many areas of research,
such as digital filter settings in communication systems [1], energy consumption in Smart Grids [7], convex  relaxations of  various  combinatorial  optimization  problems in machine learning applications [19, 39].

Our interest is in case when $m$ is large, which prohibits us from using
projected gradient and
augmented Lagrangian methods [2, 38],
that require either computation of the (Euclidean) projection or
an estimation of the gradient for the sum of many functions, at each iteration.
To reduce the complexity, one may consider a method that operates on a  single set $X_i$
from the constraint set collection $\{X_1,\ldots,X_m\}$ at each iteration.
Algorithms using random constraint sampling for general convex optimization problems (<ref>)  have been first considered  in [22] and were extended in [36] to a broader class of randomization over the sets of constraints. The convergence rate analysis in [22, 36]
demonstrates that the expected optimality error (in terms of function values) diminishes to zero with the rate of $O(1/\sqrt k)$
which is the optimal convergence rate for merely convex functions.
The work [26] further extends the results from [22, 36] to non asymptotic analysis of problems with potentially infinitely many constraints. On the other hand, a primal-dual method updating a random coordinate of the dual variable pro iteration was proposed for generalized linear programming with equality constraints in [32].

A possible reformulation of the problem (<ref>) is through the use of the indicator functions of
the constraint sets, resulting in the following unconstrained problem

min_x∈^n∑_i=1^m {1/mf(x) + χ_i(x)},
where $\chi_i(\cdot):\R^n\to\R\cup\{+\infty\}$ is the indicator function of the set $X_i$
(taking value $0$ at the points
$x\in X_i$ and, otherwise, taking value $+\infty$).
The advantage of this reformulation is that the objective function is the sum of convex functions and incremental methods can be employed that compute only a (sub)-gradient of one of the component functions at each iteration.
The traditional incremental methods do not have memory, and their origin can be traced back to work of Kibardin [15]. They have been studied for
smooth least-square problems [3, 18], for training the neural networks [10, 11], for smooth convex problems [31, 35] and
for non-smooth convex problems [23, 9, 12, 14, 37]
(see [5] for a more comprehensive survey of these methods).
However, no rate of convergence to the exact solution has been obtained for such procedures.
Reformulation (<ref>) has been considered in [16] as a departure point
toward an exact penalty reformulation using the set-distance functions.
This exact penalty formulation has been motivated by a simple exact penalty model proposed in [4] (using only the set-distance functions)
and a more general penalty model considered in [5].
In [16], a lower bound on the penalty parameter
has been identified guaranteeing that the optimal solutions of the penalized problem are also optimal solutions of the original problem (<ref>). However, this bound depends on a so-called  regularity constant for the constraint set, which might be difficult to estimate. Moreover, the proposed approaches in [16] do not utilize
incremental processing, but rather primal-dual approaches where a full (sub)-gradient of the penalized function
is used.
On a broader scale, our work is related to random methods for solving linear feasibility problems [21, 33], and their extensions to solving convex  inequality systems [20].

In contrast to the penalized formulation in (<ref>) and the
works mentioned above, this paper deals with a penalized reformulation of the problem (<ref>), where the penalized problems vary with time.
This is done by varying penalty parameters so as to gradually decrease the infeasibility of the iterates. In this way,  we can guarantee convergence of the  single time scale
procedure incremental procedure to an exact solution of the original problem (<ref>). Our choice of the penalty functions is a variant of the one-sided Huber losses [17], which have Lipschitz continuous gradients.
In the work [34], existence of the fixed penalty choices for this penalized reformulation has been shown under which the fast incremental algorithms can be applied to achieve convergence to a feasible point in a predefined neighborhood of the optimal solution of the original problem, with a linear convergence rate. However, to guarantee this convergence, some problem specific parameters need to be known, which are difficult to estimate in practice. The recent paper [25] deals with the penalty parameters which vary with time. However, that work considers a strongly convex objective function $f(\cdot)$ exclusively.
In this present paper, we show that some of the results in [25] extend to nonconvex objective functions $f(\cdot)$ and to merely convex functions, as given in Section <ref> and
Section <ref>, respectively.
The random incremental penalty-based gradient method has been proposed in [25] and analyzed for a strongly convex objective function $f(\cdot)$ with Lipschitz continuous gradients. In this paper, however,
we consider a subgradient variant of the method in Section <ref> for both convex and strongly convex $f(\cdot)$. Therein, we analyze its almost sure  convergence and show its convergence rate of the order $O(\ln^{1/2+\e} k/\sqrt{k})$ for merely convex objective function, where  $\epsilon>0$ is arbitrarily small. Moreover, for strongly convex objective function, we establish convergence rate
of the order $O(\ln^{\e} k /{k})$ for arbitrarily small $\epsilon>0$, which improves the convergence rate of $O(1/\sqrt{k})$ provided in [25]. Note that these rates possess better dependence on the logarithmic term than the rates obtained in [8]. In that work the authors focus on an optimization problem over a convex set with a special structure and present a  penalty-based method with a double loop structure, where at each iteration the corresponding penalized problem has
to be solved up to some accuracy, and prove the convergence rates of the order $O(\ln k /{\sqrt k})$ and $O(\ln k /{k})$ for the purely convex and strongly convex cases respectively.

The outline of the paper is as follows. In Section <ref>,
we provide the penalty based formulation of the original problem (<ref>)
and some basic properties of the chosen penalty functions.
In Section <ref>, we investigate the relations for the solutions of the penalized problems, as we vary the penalty parameters, for continuous nonconvex objective function $f(\cdot)$ and for convex objective function $f(\cdot)$.
In Section <ref>, we propose a random incremental penalty method and show that its iterates converge almost surely and in expectation to a solution of the original problem, under suitable assumptions on the penalty parameters and the stepsize. In Section 4, we also provide convergence rate estimates using appropriate weighted averages of the iterates. To the best of our knowledge, these are the first results on the convergence rate for the penalty-based incremental subgradient method with time-varying penalty parameters.
In Section <ref>, we conclude the paper.

§ PROBLEM FORMULATION AND ITS PENALTY-BASED REFORMULATION

We consider the following optimization problem:

min    f(x),   s.t.    a_i,x- b_i≤0,   i∈[m],
where $[m]=\{1,\ldots,m\}$ and
the vectors $a_i\in\R^n$, $i\in[m]$, are nonzero.
We will assume that the problem is feasible.
Throughout the paper, we use $X_i$ to denote the set of points satisfying the $i$-th inequality constraint, i.e.,
$X_i=\{x\in\R^n\mid  \la a_i,x\ra - b_i\le 0\} \quad\hbox{for all  }i\in[m],$
and $X$ to denote the nonempty intersection of the sets $X_i$, $i\in[m]$, i.e., $X=\cap_{i\in[m]}X_i$.
Associated with problem (<ref>), we
consider a penalized problem

\begin{equation}\label{eq:penfun}
F_{\g\d}(x) = f(x) + \frac{\gamma}{m} \sum_{i=1}^{m} h_\d\left(x; a_i, b_i\right).\end{equation}

Here, $\gamma>0$ and $\d\ge0$ are penalty parameters.
The vectors $a_i$ and scalars $b_i$ are the same as those characterizing the constraints in problem (<ref>).
For a given nonzero vector $a\in\R^n$ and $b\in\R$, the penalty function
$h_\d(\cdot;a,b)$ is given by

\begin{align}\label{eq:hfun}
h_\delta(x; a,b) = \begin{cases}
\frac{\la a,x\ra -b}{\|a\|} &\text{ if }  \la a,x\ra - b>\de,\\
\frac{(\la a,x\ra -b + \de)^2}{4\de\|a\|} &\text{ if }  -\de\le\la a,x\ra - b\le\de,\\
0 &\text{ if }  \la a,x\ra - b<-\de,
\end{cases}
\end{align}

(see Figure <ref> for an illustration).

(22,34)$\delta =1$
(22,31.4)$\delta =\frac12$
(22,28.5)$\delta =\frac14$

Penalty functions $h_\delta(x;1,1)$ for the constraint $x-1\le 0$, $x\in\R$,
with $\d\in\left\{\frac{1}{4},\frac{1}{2},1\right\}$.

The penalty function $h_\d(\cdot; a,b)$ is a variant of the one-sided Huber loss functions.
Originally, such functions have been introduced in applications of robust regression models
to make them less sensitive to outliers in data in comparison with the squared error loss [17].
In contrast, here,
we use this type of penalty function to smoothen the distance-based penalties (the average distance to the sets $X_i$)
proposed in [5].
Furthermore, an appropriate choice of the parameter $\delta\ge0$ allows us to increase the penalty values as compared to the penalties based on the squared distances to the sets $X_i$, which have a small penalty values around an optimum lying close to the boundary of the constraint set [30].

For any $\d\ge0$, the function $h_\d(\cdot;a,b)$ satisfies the following relations:

\begin{align}\label{eq:hfunineq}
h_\delta(x; a,b) \ge 0\qquad\hbox{for all }x\in\R^n,
\end{align}

\begin{align}\label{eq:hfunineq1}
h_\delta(x; a,b) \le \frac{\delta}{4\|a\|}\qquad\hbox{for all $x$ satisfying }
\la a,x\ra\le b,
\end{align}

\begin{align}\label{eq:hfunineq2}
h_\delta(x; a,b) > \frac{\delta}{4\|a\|}\qquad\hbox{for all $x$ satisfying }
\la a,x\ra> b.
\end{align}

For a vector $a\ne0$, the function $h_\delta(\cdot; a,b)$ is a composition of
a scalar function

\begin{align}\label{eq:sfun}
p_\d(s)= \begin{cases}
s &\text{ if } \quad s>\d,\\
\frac{(s + \de)^2}{4\d} &\text{ if } \quad -\de\le s \le\de,\\
0 &\text{ if } \quad s<-\de,
\end{cases}
\end{align}

with a linear function  $x\mapsto \la a,x\ra-b$, scaled by $\frac{1}{\|a\|}$, i.e.,

h_δ(x; a,b) =1/ap_(̣a,x-b)    for all x∈^n.The function $p_\d(\cdot)$ is convex on $\R$ for any $\delta\ge0$,
implying that the function $h_\d(\cdot;a,b)$ is convex on $\R^n$ for any $\d\ge0$.

Furthermore, the function
$p_\d(\cdot)$ is twice differentiable for any $\d>0$,
with the first and second derivatives given by

\begin{align}\label{eq:pderiv}
p'_\d(s)= \begin{cases}
1 &\text{ if } \quad s>\d,\\
\frac{(s + \de)}{2\d} &\text{ if } \quad -\de\le s \le\de,\\
0 &\text{ if } \quad s<-\de,
\end{cases}
\end{align}

\begin{align*}
p''_\d(s) = \begin{cases}
\frac{1}{2\d} &\text{ if } \quad -\d\le s\le\d,\\
0 &\text{ if } \quad s<-\de \quad \text{or}\quad s>\d.
\end{cases}
\end{align*}

Thus, for $\d>0$, the function $p_\d(\cdot)$ has Lipschitz continuous derivatives with a constant $\frac{1}{2\d}$.
Hence, the function $h_\delta(\cdot; a,b)$ is differentiable for any $\d>0$, and
its gradient is given by

∇h_δ(x; a,b) =1/a p'_(̣a,x-b) a    for all x∈^n.Moreover,  $\nabla h_\delta(\cdot; a,b)$ is Lipschitz continuous with a constant $\frac{\|a\|}{2\d}$, i.e.,
for all $x,y\in\R^n$,

∇h_δ(x; a,b) -∇h_δ(y; a,b) ≤a/2 x-y.
In view of the definition of the penalty function $F_{\g\d}$ in (<ref>)
and relation (<ref>),
we can see that the magnitude of the “slope" of the penalty function is controlled by the parameter $\g>0$,
while the ratio of the parameters $\g$ and $\d$ is controlling the ``curvature" of the penalty function.
Finally, in the following lemma, we provide some additional properties of the gradients $∇h_(·;a,b)$.
\begin{lem}\label{lem:pderiv}
Consider the function $h_\d(\cdot; a,b)$ as given in~\eqref{eq:hfun}.
Then, we have $\|\nabla h_\d(x;a,b)\|\le 1$ for all $x\in\R^n$. Additionally,
if $\d_1 \ge \d_2$, then
\[\max_{x\in\R^n}\|\nabla h_{\d_1}(x;a,b)- \nabla h_{\d_2}(x;a,b)\|
\le \frac{\d_1-\d_2}{2\d_1}.\]
\end{lem}
\begin{proof} Can be found in [34]. \end{proof}

Throughout the rest of the article,
we let $Π_Y [x]$ denote the (Euclidean) projection of a point $x$ on a convex closed set $Y$,
$(x, Y) = x - Π_Y [x].$
Also, the smallest norm of the vectors $a_i$, $i∈[m]$, is denoted by $_min$, i.e.,

\section{Time-varying Penalty Parameters}\label{sec-var-par}
Consider positive scalar sequences ${δ_k}$ and ${_k}$
and, for each $k$, let  $F_k$ denote the penalty function $F__̣k_k(x)$, i.e.,
\begin{equation}\label{eq-penalized-function}
F_k(x) = f(x) + \frac{\g_k}{m} \sum_{i=1}^{m} h_k\left(x; a_i, b_i\right),
\end{equation}
where we use $h_k(·; a_i, b_i)$ to denote $h__̣k(·; a_i, b_i)$ for each~$i$.
We let $X_k^*$ denote the optimal set for the penalized problem $min_x∈^n F_k(x)$ and $X^*$  denote the optimal set for the original problem~\eqref{eq:problem}. We next investigate the properties of the optimal solutions of the penalized problems for the case of a general function $f(·)$  and the case when $f(·)$ is convex.
\subsection{General Function $f$}\label{subsec-var-par-a}
When $f(·)$ is continuous and has bounded lower-level sets,
the optimal set  $X^*$ for the original problem is nonempty and compact. Also, for each $k$,
the optimal set $X_k^*$ of the penalized problem is also nonempty and compact.

As a solution $x_k^*∈X_k^*$ need not be feasible, we next estimate the distance of any $x_k^*∈X_k^*$
from the feasible set $X$.
To simplify the notation, we denote by $H_k(·)$ the average of the penalty functions $h_k(·;a_i,b_i)$, $i∈[m]$,
\begin{equation}\label{eq-hk}
H_k(x)= \frac{1}{m}\sum_{i=1}^{m}h_k(x;a_i,b_i)\qquad\hbox{for all }x\in\R^n.
\end{equation}
Thus, the penalty function in~\eqref{eq-penalized-function}
is written as
\begin{equation}\label{eq-penfun-hk}
F_k(x) = f(x) + \g_kH_k(x)\qquad\hbox{for all }x\in\R^n.
\end{equation}
Regarding the average penalty $H_k(·)$, we have an upper bound on its values at feasible points $x∈X$ via relation~\eqref{eq:hfunineq1}.
Specifically, setting $=̣_̣k$ in relation~\eqref{eq:hfunineq1}, we have $h_k(x; a_i,b_i) ≤_̣k/4a_i≤_̣k/4_min$  for all $x∈X$ and $i∈[m]$,
implying that
\begin{equation}\label{eq-hk-atfeasx}
H_k(x)\le \frac{\d_k}{4\a_{\min}}\qquad\hbox{for all $x\in X$ and all }k\ge1.
\end{equation}

To lower bound the value $H_k(x)$ at any $x∈^n$, we use
the following result regarding the intersection of linear sets.
\begin{lem}[Hoffman lemma~[13]]\label{lem-hoffman}
Given a collection of sets $X_i=\{x\in\R^n\mid \la a_i,x\ra -b\le 0\}$, $i\in[m]$,
with a nonempty intersection $X=\cap_{i=1}^m X_i$,
there exists a scalar
$\beta=\beta(a_1,\ldots,a_m)>0$ such that
$\beta \sum_{i=1}^{m} \dist(x, X_i) \ge \dist(x,X)$ for all  $x\in\R^n$.
\end{lem}
The following lemma and its corollary provide some additional
properties of the penalty function $h_(̣·;a,b)$ that we will use later on.
The proof can be found in Lemma~1 in~[34].

\begin{lem}\label{lem:penalty}
Given a nonzero vector $a\in\R^n$ and a scalar $b\in\R$,
consider the penalty function $h_\d(\cdot;a,b)$ defined in~\eqref{eq:hfun} with $\d\ge0$.
Let $Y=\{x\in\R^n\mid \la a,x\ra-b\le0\}.$
Then, we have for $\d=0$, $h_0(x;a,b)=\dist(x,Y)$ for all $x\in\R^n$,
and for any $0\le \d \le \d'$,
$h_\d(x;a,b) \le h_{\d'}(x;a,b)$ for all $x\in\R^n$.
\end{lem}

The following corollary is an immediate consequence of Lemma~\ref{lem:penalty},
which shows that any
feasible point $x̂∈X$ can be used to construct non-empty level sets of
$F_(·)$ and $f(·)$.
The proof can be found in Corollary~2 in~[34].
\begin{cor}\label{cor:lset}
Let $\g>0$ and $\d\ge 0$ be arbitrary, and let
$\hat x$ be a feasible point for the original problem~\eqref{eq:problem}.
Then, for the scalar $t_{\g\d}(\hat x)$ defined by
$t_{\g\d}(\hat x)= f(\hat x)+\g\d/(4\a_{\min}),$
the level set
$\{x\in\R^n\mid F_{\g\d}(x)\le t_{\g\d}(\hat x)\}$ is nonempty and
$\{x\in\R^n\mid F_{\g\d}(x)\le t_{\g\d}(\hat x)\}\subseteq
\{x\in\R^n\mid f(x)\le t_{\g\d}(\hat x)\}$.
Moreover, the solution set $X^*_{\g\d}$ of
the penalized problem~\eqref{eq:pen-problem} is contained in the level set
$\{x\in\R^n\mid f(x)\le t_{\g\d}(\hat x)\}$.
\end{cor}
Now, a lower bound on the value $H_k(x)$ at any $x∈^n$ follows from
Hoffman Lemma~\ref{lem-hoffman} and Lemma~\ref{lem:penalty}, as follows.
By Lemma~\ref{lem:penalty} with $=̣0$, $'̣=_̣k$, $a=a_i$, $b=b_i$, and $Y=X_i$, we obtain $h_k(x;a_i,b_i)≥h_0(x;a_i,b_i)=(x,X_i)  for all  x∈^n$ for all $i∈[m]$.
Therefore, for all $x∈^n$,
$H_k(x)=1/m∑_i=1^m h_k(x; a_i,b_i)
≥1/m∑_i=1^m (x,X_i)$.
By using Hoffman Lemma~\ref{lem-hoffman}, we obtain
\begin{equation}\label{eq-hk-atanyx}
H_k(x)\ge \frac{1}{m\b} \dist(x,X)\quad\hbox{for all } x\in\R^n.
\end{equation}

The following result provides an upper bound on the distance of $x_k^*$ from the feasible set $X$,
and shows that this distance goes to 0 if $_k$ tends to infinity and $_̣k$ tends to 0.

\begin{prop}\label{prop-distancetofeas}
Let $f(\cdot)$ be continuous with bounded lower-level sets,
and let $\g_k>0$, $\d_k>0$, and $\g_k\d_k\le c$ for all $k$ and for some $c>0$.
Then, for arbitrary $\hat x\in X$ and all $k$,
\[\dist(x_k^*,X) \le \frac{m\b}{\g_k} (f(\hat x) - f(x_k^*))+\frac{m\b\d_k}{4\a_{\min}},\]
where %$\a_{\min}=\min_{i\in[m]}\|a_i\|$ and
$\beta$ is the Hoffman constant from Lemma~\ref{lem-hoffman}. In particular,
$\lim_{k\to\infty}\dist(x_k^*,X)=0$ as  $\g_k\to\infty$  and  $\d_k\to0$,
with the convergence rate of the order $O(\g_k^{-1} +\d_k)$.
\end{prop}
We refer the reader to Appendix~\ref{app:distancetofeas} for the proof.

\subsection{Convex Function $f$ not Necessarily Differentiable}\label{subsec-var-par-b}
Here, we consider the case when $f(·)$ is convex on $^n$. Since the domain of $f(·)$ is $^n$,
the function $f$ is continuous. The subdifferential set $∂f(x)$ is nonempty at every $x∈^n$
since $f(·)$ is convex and $(f)=^n$.
The function $f(·)$ is strongly convex with a constant $μ>0$ if and only if
for all $u,v∈^n$ and all subgradients $∇̃f(x)∈∂f(x)$, we have
\begin{equation}\label{eq-strconvex}
f(v)+\la \tilde\nabla f(v),u-v\ra +\frac{\mu}{2}\|u-v\|^2\le f(u).
\end{equation}
Moreover, when $f(·)$ is strongly convex with a constant $μ$, we also have for all $x,y∈^n$, and all subgradients $∇̃f(x)∈∂f(x)$
and $∇̃f(y)∈∂f(y)$,
\begin{equation}\label{eq-subgrad-strong}
\mu\|x-y\|^2\le \la \tilde \nabla f(x)-\tilde \nabla f(y),x-y\ra.
\end{equation}
If $f(·)$ is just convex relations~\eqref{eq-strconvex} and~\eqref{eq-subgrad-strong} hold with $=0$.
To capture both cases when $f$ is strongly convex and when $f$ is just convex,
we will sometimes abuse the definition of strong convexity in~\eqref{eq-strconvex} by allowing
the possibility that $μ=0$.

When $f(·)$ is strongly convex with $μ>0$, then so is every penalty function $F_k(·)$ in~\eqref{eq-penalized-function} with the same $μ>0$.
In this case,
the original problem~\eqref{eq:problem} and each penalized problem $min_x∈^nF_k(x)$, $k≥1$, have unique solutions, respectively, denoted by $x^*∈X$ and $x_k^*∈^n$, respectively.
Moreover, under mild conditions on the penalty parameters $_̣k$ and $_k$, the optimal set sequence ${X_k^*}$ is uniformly bounded in the case of continuous function with bounded lower-level sets, as seen in the following lemma. Its proof relies on Corollary~\ref{cor:lset} and is provided in Appendix~\ref{app:solsbded}.
\begin{lem}\label{lem-solsbded}
Let $f(\cdot)$ be continuous and have bounded lower-level sets.
Let $\g_k>0$,
$\d_k>0$, and $\g_k\d_k\le c$ for all $k\ge 1$ and for some $c>0$.
Then, the optimal set $X_k^*$ for the penalized problem $\min_{x\in \R^n} F_k(x)$
is nonempty compact set and
$X_k^*\subseteq\left\{x\in\R^n\mid f(x)\le t_c(\hat x)\right\}$ for all~$k$, where
$t_c(\hat x)=f(\hat x)+c/(4\a_{\min})$, with $\hat x\in X$.
In particular, the set sequence $\{X_k^*\}$ is uniformly bounded.
\end{lem}
\begin{rem}\label{remsgd-bound}
When the conditions of Lemma~\ref{lem-solsbded} are satisfied, the optimal solutions of the penalized problems
$\min_{x\in\R^n} F_k(x)$, for all $k\ge1$,
are uniformly bounded, i.e., there exists $D>0$ such that
$\|x_k^*\|\le D$ for all $x_k^*\in X_k^*$ and all $k\ge1$.
Therefore, the projections $\Pi_X[x_k^*]$ of these optimal solutions on the feasible set $X$
are also uniformly bounded, i.e., there exists $R>0$ such that
$\|\Pi_X[x_k^*]\|\le R$ for all $x_k^*\in X_k^*$ and all $k\ge1$.
Hence, the subgradients $\tilde \nabla f(x)\in\partial f(x)$, for all $x$ with $\|x\|\le R$,
are bounded, i.e., \begin{equation}\label{eq:sgdbound}
L=\max_{\|x\|\le R}\{\|\tilde \nabla f(x)\|\,\mid \tilde \nabla f(x)\in\partial f(x)\}<\infty.
\end{equation}
\end{rem}

We next provide a set of conditions on parameters $_̣k$ and $_k$ ensuring
that the sequence ${x_k^*}$ converges to $x^*$ as $k→∞$ when $f(·)$ is strongly convex.
When $f(·)$ is just convex, we obtain a special bound on $(x_k^*,X)$ for any solution $x_k^*∈X_k^*$ to the penalized problem.
This bound yields an improved convergence rate of $(x_k^*,X)→0$ compared to
that for a general function $f(·)$ provided in Proposition~\ref{prop-distancetofeas}.

\begin{prop}\label{prop-sols}
Let $f(\cdot)$ be strongly convex with $\mu\ge0$. If $\mu=0$, assume that $f(\cdot)$
has bounded lower-level sets.
Let $\g_k>0$, $\d_k>0$, and $\g_k\d_k\le c$ for all $k$ and some $c>0$.
Then, for all $k$, we have\vspace{-0.25cm}
\begin{align*}
&\frac{\mu}{2}\|x^* -x_k^*\|^2 + \frac{\mu}{2}\|x^*-\Pi_X[x_k^*]\|^2  +\left(\frac{\mu}{2}+\frac{\g_k}{m\b} -L\right)\dist(x_k^*,X)
\le \frac{\g_k\d_k}{4\a_{\min}},
\end{align*}
where $L$ is given by~\eqref{eq:sgdbound}.% and $\a_{\min}=\min_{ i\in[m]}\|a_i\|$.
\end{prop}
The proof is in Appendix~\ref{app:sols}.

Proposition~\ref{prop-sols} indicates that, when $μ>0$,
by letting $_k→∞$,
we will have $_k/4m≥L$ for all large enough $k$,
\begin{align}\label{eq-disttofeas-convex}
\dist(x_k^*,X)&\le \frac{\g_k\d_k}{4\a_{\min} \left(\frac{\mu}{2}+\frac{\g_k}{4m\b}-L\right)}\approx O(\d_k).
\end{align}
Thus, if $_̣k→0$,  the distance of $x_k^*$ to the feasible set $X$ will go to 0 asymptitically
at the rate of  $O(_̣k)$, independent of $_k$.
The preceding relation holds with $μ=0$ if $f(·)$ is merely convex.
The convergence rate order is better than $O(_k^-1)+O(_̣k)$ obtained for a general $f(·)$
in Proposition~\ref{prop-distancetofeas}.

When $f(·)$ is strongly convex with $μ>0$,
Proposition~\ref{prop-sols} shows that $x^* -x_k^*^2≤_k_̣k/2μ_min$ for large enough $k$.
If additionally $_k_̣k→0$, then the solutions $x_k^*$ of the penalized problems converge to
the optimal solution $x^*$ of the original problem, with the rate in the order of $O(_k_̣k)$.

The common requirement for $(x_k^*,X)→0$, as $k→∞$, in Proposition~\ref{prop-distancetofeas} and Proposition~\ref{prop-sols} is that $_̣k→0$.
The main difference between these propositions is in the requirement for the penalty parameter $_k$.
to ensure $(x_k^*,X)→0$, the penalty $_k$ has to increase to $+∞$
for a general function~$f(·)$ (Proposition~\ref{prop-distancetofeas}).
In contrast,
to ensure that $(x_k^*,X)→0$ for a convex function $f(·)$, one can choose a fixed penalty value $_k=$, for all $k$, with $γ$ large enough so that $/4m> L$, as seen from relation~\eqref{eq-disttofeas-convex}.
However, determining such a value of $$ is challenging as it is hard to obtain upper estimates for the subgradient norm bound $L$ and the Hoffman constant $$̱. Estimating $L$ requires knowing a region that contains the projections of the solutions to the penalized problems on the feasible set (see (<ref>)). Determining an upper estimate of the Hoffman constant $\b$ is also a difficult problem, which has recently been addressed in~[27] via computational approaches.

Based on Proposition~\ref{prop-sols}, one can construct a two-loop iterative approach to compute
the optimal point $x^*$ of the original problem {in the case of strongly convex $f$}. The outer loop is on the index $k$ where the penalty values $_k$ and $_̣k$ are set.
For any given $k$, the inner loop of iterations compute the optimal point $x_k^*$ for the penalized problem $min_x∈^nF_k(x)$.
This naive two-loop approach is quite inefficient. Later in Section~\ref{sec-algo}, we propose a more efficient single-loop  algorithm,
where at each iteration $k$, we
adjust the parameters $_k$ and $_̣k$, and use a (stochastic) gradient of the penalty function $F_k(·)$ for the update.

\section{Random Incremental Penalty Algorithm}\label{sec-algo}
Assuming that the function $f(·)$ is convex over $^n$,
we consider an algorithm that takes one gradient step for minimizing $F_k(·)$ at iteration $k$,
as opposed to determining $x^*_k$ for each function $F_k(·)$.
To deal with the large number of component functions $h_k(·;a_i,b_i)$ involved in
$F_k(·)$, we consider a random incremental subgradient algorithm using only
one randomly chosen constraint (indexed by $i_k$) to estimate
a subgradient $∇̃F_k(x_k)∈∂F_k(x_k)$ at iteration $k$, when $x_k$ is available.
This estimation is employed to construct $x_k+1$ via {\textit {random penalty}} corresponding to subgradient sampling,
as opposed to determining a full subgradient of $F_k(x_k)$.

We apply random incremental update to the penalty function~\eqref{eq-penalized-function}
represented in the following form:
\[F_k(x) = \frac{1}{m} \sum_{i=1}^m\left( f(x) + \g_k  h_k\left(x; a_i, b_i\right) \right),\]
and the random incremental penalty method is:
for $k≥1$,
\begin{equation}\label{eq:gradmet}
x_{k+1} = x_k-s_k [\tilde \nabla f(x_k) + \g_k\nabla h_k(x_k;a_{i_k},b_{i_k})],
\end{equation}
where $s_k>0$ is a stepsize, $∇̃f(x_k)$ is a subgradient of $f(·)$ at $x=x_k$,
%($\tilde\nabla f(x_k)\in\partial f(x_k)$), and
and the index $i_k∈{1,…,m}$ is chosen uniformly at random at every iteration $k$.
\an{The algorithm is initiated with a random initial point $x_1\in\R^n$, for which we assume that

Note that $∇̃f(x_k) + _k∇h_k(x_k;a_i_k,b_i_k)$  is an unbiased estimate of
a subgradient $∇̃F_k(x_k) ∈∂F_k(x_k)$,
since by the uniform distribution of $i_k$ we have
\begin{align}\label{eq-nobias}
\ex{\tilde \nabla f(x_k) + \g_k\nabla h_k(x;a_{i_k},b_{i_k})|\EuScript F_k}=\tilde \nabla F_k(x_k),
\end{align}
where $F_k$ is the $σ$-algebra generated by the random variables
${i_j,  1 ≤j≤k-1}$ and the random initial iterate $x_1$, which is equivalent to
the $σ$-algebra $F_k = σ({x_t}_t=1^k)$ for all $k≥1$.

Unlike the standard random incremental method,
the random incremental penalty method in~\eqref{eq:gradmet} selects one
random component from the time-varying function $F_k(·)$.
This makes the analysis of the method more challenging since the iterates ${x_k}$ need not be feasible for the original problem,
which poses difficulties, and the penalty parameters $_k$ and $_̣k$ have to be carefully tuned
to ensure the convergence of the iterates.

\subsection{Preliminary results}

We first establish a basic result for the iterates $x_k$ of the random method~\eqref{eq:gradmet}.
\begin{lem}\label{lem-basic-iter}
Let $f(\cdot)$ be strongly convex with $\mu\ge 0$, and
let $\g_k>0$, $\d_k>0$, and $s_k>0$ for all $k$.
Then, the iterates $x_k$ of the method~\eqref{eq:gradmet} surely satisfy
for all $y\in X$ and $k\ge1$,
\begin{align*}
&\|x_{k+1}- y\|^2 \le(1-\mu s_k)\|x_k - y\|^2
&+\frac{s_k\g_k \d_k}{2\a_{\min}} - 2s_k\g_k\dist(x_k,X_{i_k})
+s_k^2(\|\tilde\nabla f(x_k)\| + \g_k)^2.
\end{align*}
\end{lem}

{\it Proof}:
By the definition of $x_k+1$ in~\eqref{eq:gradmet}, we surely have
for any $y∈X$ and all $k≥1$,
\begin{align*}
\|x_{k+1}- y\|^2 =\|x_k - y\|^2 &-2s_k\la g_k(x_k),x_k-y\ra +s_k^2\| g_k(x_k)\|^2,
\end{align*}
with $g_k(x_k)=∇̃f(x_k) + _k∇h_k(x_k;a_i_k,b_i_k)$.
By the convexity of $h_k(·;a_i_k,b_i_k)$, we have
for all $y∈X$ and all $k≥1$,
\begin{align*}
&\|x_{k+1}- y\|^2 \le\|x_k - y\|^2
-2s_k\la \tilde\nabla f(x_k), x_k - y\ra \cr
&+2s_k \g_k(h_k(y;a_{i_k},b_{i_k})-h_k(x_k;a_{i_k},b_{i_k}))
+s_k^2\| g_k(x_k)\|^2.
\end{align*}
Since $y$ is feasible, we have $h_k(y;a_i_k,b_i_k)≤_̣k / 4_min$, with $_min=min_i∈[m]a_i$,
which follows by relation~\eqref{eq:hfunineq1} (where $=̣_̣k$ and $x=y$).
By the monotonicity of $h_(̣·;a_i,b_i)$ with respect to $$̣, we have that
$\dist(x,X_i)=h_0(x;a_{i},b_{i})\le h_{\d_k}(x;a_i,b_i)$ for all $x$ and $i$
(see Lemma <ref> where $\d=0$, $\d'=\d_k$).
Hence, it follows that for all $y\in X$ and all $k\ge1$,

\begin{align*}
&\|x_{k+1}- y\|^2 \le\|x_k - y\|^2
-2s_k\la \tilde\nabla f(x_k), x_k - y\ra\cr
&+\frac{s_k\g_k \d_k}{2\a_{\min}} - 2s_k\g_k\dist(x_k,X_{i_k})
+s_k^2\| g_k(x_k)\|^2.
\end{align*}

By the strong convexity relation (where $u=y$, $v=x_k$), it follows that

\[-2s_k\la \tilde\nabla f(x_k), x_k - y\ra\le 2s_k(f(y)-f(x_k))-\mu s_k\|x_k-y\|^2.\]

By combining the preceding two relations,
we obtain that, surely, for all $y\in X$ and all $k\ge1$,

\begin{align*}
&\|x_{k+1}- y\|^2 \le(1-\mu s_k)\|x_k - y\|^2
&+\frac{s_k\g_k \d_k}{2\a_{\min}} - 2s_k\g_k\dist(x_k,X_{i_k})
+s_k^2\| g_k(x_k)\|^2.
\end{align*}

To estimate $\| g_k(x_k)\|$, we write

\[\|g_k(x_k)\|\le \|\tilde \nabla f(x_k)\| + \g_k\|\nabla h_k(x_k;a_{i_k},b_{i_k})\|
\le \|\tilde \nabla f(x_k)\| + \g_k,\]

where the last inequality uses the fact that
$\|\nabla h_k(x;a_{i},b_{i})\|\le 1$ for all $x$ and all $i$ (see Lemma <ref>).
The stated relation follows from the preceding two inequalities.

Lemma <ref> is important for both convergence and convergence rate analysis.
We now introduce additional assumptions and refine Lemma <ref>.
Specifically, we assume that the objective function is convex and has bounded level sets, and the subgradient norms
$\|\tilde \nabla f(x)\|$ grow at most linearly with $\|x\|$.

The function $f(\cdot)$ is convex and has bounded level sets.

Note that if $f(\cdot)$ is strongly convex with $\mu>0$, then Assumption <ref> is satisfied.
The assumption is also satisfied if $f(\cdot)$ is convex and coercive, i.e., $\lim_{\|x\|\to\infty}f(x)=+\infty$.
We make the following assumption regarding the subgradients of $f(\cdot)$.

There exist scalars $M_1,M_2>0$
such that $\|\tilde \nabla f(x)\|\le M_1\|x\|+M_2$
for all subgradients $\tilde \nabla f(x)\in\partial f(x)$ and for $x\in\R^n$.

The following lemma will be important in establishing the convergence properties of the method.

Let Assumption <ref> and Assumption <ref> hold.
Then, it surely holds
for all $x^*\in X^*$  and $k\ge1$,

\begin{align*}
&\|x_{k+1}- x^*\|^2 \le(1+4s_k^2M_1^2)\|x_k - x^*\|^2 \cr
&\ \
+2s_k\left((1-\rho) B +\rho( M_1B+M_2)\right)\dist(x_k,X)\cr
&\ \
+2s_k\rho(f^*- f(\Pi_X[x_k]))
+2s_k\rho M_1\|x_k-x^*\|\,\dist(x_k,X)\cr
& \ \ +\frac{s_k\g_k \d_k}{2\a_{\min}}- 2s_k\g_k\dist(x_k,X_{i_k})
+4s_k^2(M_1^2 B^2 +M_2^2+\g_k^2),\qquad
\end{align*}

where $\rho\in[0,1]$ is arbitrary scalar and $B$ is a norm bound for the optimal solutions $x^*\in X^*$ and the subgradients $\tilde\nabla f(x^*)$ at any $x^*\in X^*$.

By Lemma <ref>, where we omit the term $-\mu s_k$ in the coefficient of
$\|x_k - y\|^2$, we surely have
for all $y\in X$ and $k\ge1$,

\begin{align}\label{eq-jedan}
&\|x_{k+1}- y\|^2 \le \|x_k - y\|^2
&+\frac{s_k\g_k \d_k}{2\a_{\min}} - 2s_k\g_k\dist(x_k,X_{i_k})
+s_k^2(\|\tilde\nabla f(x_k)\| + \g_k)^2.
\end{align}

Under our assumption on the subgradient-norm growth (Assumption <ref>)
we have that

\[\|\tilde\nabla f(x_k)\| + \g_k
\le M_1\|x_k\|+M_2 + \g_k
\le M_1(\|x_k-y\|+\|y\|)+M_2 + \g_k\]

where the last inequality is obtained by using $\|x_k\|\le \|x_k-y\|+\|y\|$.
Next, using $(a+b+c+d)^2\le 4 (a^2 + b^2+c^2+d^2)$, which is valid for any scalars $a,b,c,$ and $d$,
we have

\[(\|\tilde\nabla f(x_k)\| + \g_k)^2\le 4 M_1^2(\|x_k-y\|^2+\|y\|^2)+4M_2^2+4\g_k^2.\]

Hence, by substituting the preceding estimate back in relation (<ref>),
we surely have for all $y\in X$ and all $k\ge1$,

\begin{align}\label{eq-rel-perk}
&\|x_{k+1}- y\|^2 \le(1+4s_k^2 M_1^2)\|x_k - y\|^2
+2s_k(f(y)- f(x_k)) \cr
& +\frac{s_k\g_k \d_k}{2\a_{\min}}- 2s_k\g_k\dist(x_k,X_{i_k})
+4s_k^2(M_1^2 \|y\|^2 +M_2^2+\g_k^2).\qquad
\end{align}

By our assumption that the function $f(\cdot)$ has bounded level sets (Assumption <ref>),
the problem $\min_{x\in X}f(x)$ has a nonempty compact convex solution set $X^*$, i.e.,  $\|x^*\|\le B$
for some $B>0$ and for all $x^*\in X^*$.
Let $x^*\in X^*$ be an arbitrary solution,
and let $y=x^*$ in (<ref>). Thus, we  surely have for all $x^*\in X^*$  and $k\ge1$,

\begin{align}\label{eq-rel-perk-last}
&\|x_{k+1}- x^*\|^2 \le(1+4s_k^2M_1^2)\|x_k - x^*\|^2
+2s_k(f^*- f(x_k)) \cr
& +\frac{s_k\g_k \d_k}{2\a_{\min}}- 2s_k\g_k\dist(x_k,X_{i_k})
+4s_k^2(M_1^2 B^2 +M_2^2+\g_k^2),\qquad
\end{align}

where $f^*$ is the optimal value of the original problem.

The remaining part of the proof relies on using two different ways to upper bound the value $f^*- f(x_k)$ in (<ref>) by using the convexity of $f(\cdot)$. One way is to
write  for any $x^*\in X^*$,

\begin{eqnarray}\label{eq:lowerbound}
f^*- f(x_k)
&=& f(x^*)- f(x_k)\cr
&\le &\la\tilde\nabla f(x^*),x^*-x_k\ra \cr
&=&\la\tilde\nabla f(x^*),x^*-\Pi_X[x_k]\ra +\la\tilde\nabla f(x^*),\Pi_X[x_k]-x_k\ra \cr
&\le& \|\tilde\nabla f(x^*)\|\, \|\Pi_X[x_k]-x_k\|,
\end{eqnarray}

where the last inequality is obtained by using $\la\tilde\nabla f(x^*),x^*-\Pi_X[x_k]\ra\le0$,
which holds by the optimality of $x^*$ (since $\Pi_X[x_k]$ is feasible),
and by applying the Cauchy-Schwarz inequality to estimate the other inner product term.
Since $X^*$ is bounded by Assumption <ref>, we may assume without loss of generality that $B$ is large enough so that
$\|\tilde\nabla f(x^*)\|\le B$ and $\|x^*\|\le B$ for all subgradients $\tilde\nabla f(x^*)$ and all $x^*\in X^*$, so that we surely have

\begin{equation}\label{eq-oneway}
f^*- f(x_k)\le B \|\Pi_X[x_k]-x_k\|
\end{equation}

Another way is to write

\begin{eqnarray*}
f^*- f(x_k) &=&
f^*-f(\Pi_X[x_k])+f(\Pi_X[x_k]) -f(x_k)\cr
&\le &f^*-f(\Pi_X[x_k])+
\la\tilde\nabla f(\Pi_X[x_k]),\Pi_X[x_k]-x_k\ra\cr
&\le & f^*-f(\Pi_X[x_k])+
\|\tilde\nabla f(\Pi_X[x_k])\|\,\|\Pi_X[x_k]-x_k\|.
\end{eqnarray*}

Using the assumption that the subgradient-norm growth is at most linear
(Assumption <ref>),
we have that

\[\|\tilde\nabla f(\Pi_X[x_k])\|\le M_1\|\Pi_X[x_k]\| +M_2\le M_1(\|\Pi_X[x_k]-x^*\| +B)+M_2,\]

where we use the fact that $\|x^*\|\le B$
for any $x^*\in X^*$.
By combining the preceding two relations and using
the fact that $\|\Pi_X[x_k]-x^*\|\le \|x_k-x^*\|$,
we obtain

\begin{equation}\label{eq-anotherway}
f^*- f(x_k) \le  f^*-f(\Pi_X[x_k]) +
(M_1\|x_k-x^*\| +M_1B+M_2 )\|\Pi_X[x_k]-x_k\|
\end{equation}

Multiplying the inequality (<ref>) with $1-\rho$ and the inequality (<ref>) with $\rho$, for some $\rho\in[0,1]$, and combining these
with inequality (<ref>), we obtain surely
for all $x^*\in X^*$  and $k\ge1$,

\begin{align*}
&\|x_{k+1}- x^*\|^2 \le(1+4s_k^2M_1^2)\|x_k - x^*\|^2 +2s_k(1-\rho)
B \|\Pi_X[x_k]-x_k\|\cr
&\ \
+2s_k\rho(f^*- f(\Pi_X[x_k]))
+2s_k\rho\left(M_1\|x_k-x^*\| +M_1B+M_2\right)\,\|\Pi_X[x_k]-x_k\|\cr
& \ \ +\frac{s_k\g_k \d_k}{2\a_{\min}}- 2s_k\g_k\dist(x_k,X_{i_k})
+4s_k^2(M_1^2 B^2 +M_2^2+\g_k^2),\qquad
\end{align*}

which gives the stated relation by noting that
$\|\Pi_X[x_k]-x_k\|=\dist(x_k,X)$, and by grouping the terms accordingly.

§.§ Almost Sure and in-Expectation Convergence

In this section, we establish almost sure convergence of the random incremental penalty method (<ref>).
In the forthcoming discussion we will often use
almost surely.
Our convergence analysis makes use of
the following result on semi-supermartingale convergence,
which is due to Robbins and Siegmund [29] (it can also be found in [28], Lemma 11, page 50).

Let $\{v_k\}$, $\{u_k\}$, $\{\a_k\}$, and $\{\b_k\}$ be random nonnegative scalar sequences such that
$\sum_{k=0}^\infty \a_k<\infty$ and $\sum_{k=0}^\infty \b_k<\infty$ ,

\[\mathbb{E}\left[v_{k+1}\mid \mathcal{F}_k \right]\le(1+\a_k) v_k -u_k+\b_k\quad \hbox{for all }k\geq 0\ \as,\]

$\mathcal{F}_k=\{v_\ell,u_\ell,\a_\ell,\b_\ell;\, 0\le \ell\le k\}$.
Then, we have  that $\sum_{k=0}^{\infty}u_k<\infty$ and $\lim_{k\to\infty}v_k=v$ for a random variable $v\geq 0$.

Using Lemma <ref>, we establish the following result which will be used in the convergence analysis of the method.

Consider a minimization problem $\min_{x \in Z} \phi(z)$, where $\phi:\mathbb{R}^n\to\mathbb{R}$ is a
continuous function and $Z\subseteq\mathbb{R}^n$ is a closed convex set. Assume that
the solution set $Z^*$ of the problem is nonempty. Let
$\{z_k\}\subset \mathbb{R}^n$ be a random sequence,
$\{a_k\}$ and $\{c_k\}$ be random nonnegative scalar sequences, and $\{b_k\}$ and  $\{b'_k\}$ be deterministic nonnegative scalar sequences
such that  for all $z^*\in  Z^*$ and for all $k\ge1,$

\[\mathbb{E}[\|z_{k+1}-z^*\|^2\mid\mathcal{F}_k] \le (1+a_k)\|z_k-z^*\|^2
- b_k\left (\phi(\Pi_Z[z_k]) -\phi^*\right) +c_k,\]

\[\sum_{k=1}^\infty b'_k\dist(z_k,Z_k)<\infty\qquad\as\]

where $\mathcal{F}_k=\{z_1,\ldots,z_k\}$ for all $k$,
$\phi^*=\min_{x \in Z} \phi(z)$, while the scalar sequences satisfy $\sum_{k=1}^\infty a_k<\infty$ and
$\sum_{k=1}^\infty c_k<\infty$ , and
$\sum_{k=1}^\infty b_k=\infty$ and $\sum_{k=1}^\infty b'_k=\infty$.
Then, the sequence $\{z_k\}$ converges to some (random)  optimal solution $z^*\in X^*$ almost surely. Moreover, if the solution set $Z^*$ is bounded, then there exists a scalar $M_0$ such that $\|z_k\|\le M_0$ for all $k\ge0$ and convergence to $z^*$ is also in expectation.

We note that the conditions of Lemma <ref> are satisfied with $v_k=\|z_k-z^*\|^2$, for every $z^*\in Z^*$, and $u_k=b_k\left (\phi(\Pi_Z[z_k]) -\phi^*\right)$, and by this lemma
we obtain the following statements:

\begin{equation}\label{eq:xtconv}
\hbox{$\{\|z_k - z^*\|^2\}$ converges \as\ for each $z^*\in Z^*$},
\end{equation}

\begin{equation}\label{eq:sumf}
\sum_{k=1}^\infty b_k \left(\phi(\Pi_Z[z_k]) - \phi^* \right)<\infty \qquad \as.\end{equation}

Since $\sum_{k=1}^\infty b_k=\infty$, it follows from (<ref>) that ,
$\liminf_{k\to\infty} \phi(\Pi_Z[z_k])=\phi^*$.
The conditions $\sum_{k=1}^\infty b'_k\dist(z_k,Z_k)<\infty$ , and $\sum_{k=1}^\infty b'_k=\infty$
imply that
$\liminf_{k\to\infty} \|z_k-\Pi_Z[z_k]\|=0.$

Let $\{z_{k_\ell}\}$ be a subsequence of $\{z_k\}$ such that ,

\begin{equation}\label{eq-lim1}
\lim_{\ell\to\infty} \phi(\Pi_Z[z_{k_\ell}])
=\liminf_{k\to\infty} \phi(\Pi_Z[z_k])=\phi^*,\end{equation}

\begin{equation}\label{eq-lim2}
\lim_{\ell\to\infty} \|z_{k_\ell}-\Pi_Z[z_{k_\ell}]\| =\liminf_{k\to\infty} \|z_k-\Pi_Z[z_k]\|=0.
\end{equation}

Now, relation (<ref>) implies that the sequence $\{z_k\}$ is bounded , so
without loss of generality, we can assume that $\{z_{k_\ell}\}$ is converging  to some
random point $\tilde z$ (for otherwise, we
can in turn select  convergent subsequence of  $\{z_{k_\ell}\}$). Since the projection mapping $z\mapsto \Pi_Z[z]$ is continuous, it follows that

\[\lim_{\ell\to\infty}\Pi_Z[z_{k_\ell}]
=\Pi_Z[\tilde z]\qquad\as\]

The preceding relation and relation (<ref>)
imply that

\[0=\lim_{\ell\to\infty}\|z_{k_\ell}-\Pi_Z[z_{k_\ell}]\|=
\|\tilde z -\Pi_Z[\tilde z]\|\qquad\as\]

Therefore, we have that $\tilde z\in Z$ Moreover, by the continuity of $\phi(\cdot)$,

\[\lim_{\ell\to\infty} \phi(z_{k_\ell})=\phi(\tilde z)\qquad\as,\]

which by relation (<ref>) and the fact $\tilde z\in Z$  implies that $\tilde z\in Z^*$  By letting $z^*=\tilde z$ in (<ref>)
we obtain that $\{z_k\}$ converges to $\tilde z$

When the set $Z^*$ is bounded, the convergence point $\tilde z$ of the iterates $\{z_k\}$ is bounded by a (deterministic) scalar almost surely, implying that $\{\|z_k\|\}$ is   bounded by some deterministic scalar $M_0$. Finally, using  Lebesgue's dominated convergence theorem, we conclude that $\{z_k\}$ converges to $\tilde z$ also in expectation.

Having Lemma <ref>,
Lemma <ref>, and Lemma <ref> in place, we  show next the convergence of the method
under some conditions on the stepsize $s_k$ and the penalty parameters $\g_k,\d_k$, as given in the following assumption.

Let $\g_k>0$, $\d_k>0$, and $s_k>0$ for all $k$, and assume that
$\lim_{k\to\infty}\g_k=\infty$, $\sum_{k=1}^\infty s_k=\infty$,
$\sum_{k=1}^\infty s_k\g_k \d_k<\infty$, $\sum_{k=1}^\infty s_k^2\g_k^2<\infty$.

The intuition behind the conditions in Assumption <ref> is as follows. The condition $\lim_{k\to\infty}\g_k=\infty$ ensures that the penalty function pushes the iterations into the feasible set as time runs, whereas $\sum_{k=1}^\infty s_k=\infty$ allows the algorithm to make sufficient progress toward an optimal solution of the original problem. The last two conditions, $\sum_{k=1}^\infty s_k\g_k \d_k<\infty$ and $\sum_{k=1}^\infty s_k^2\g_k^2<\infty$, keep the perturbations caused by the penalty function under control.

In the following proposition, we establish almost sure iterates' convergence
and their boundedness properties. For this purpose we define the sigma-algebra $\EuScript F_k$ relevant to the random method (<ref>) as follows:

\[\EuScript F_k=\{x_1,\ldots,x_k\}\qquad\hbox{for all }k\ge 1.\]

Let Assumptions <ref>–<ref> hold. Then, the iterates of the random method (<ref>) converge to a (random)  optimal solution  and in expectation. Moreover, there exists a scalar $M>0$ such that
$\|x_k\|\le M$ for all $k\ge 1$

Proof: We use Lemma <ref>, according to which we surely have for all $x^*\in X^*$, all $\rho\in[0,1]$, and all $k\ge 1$,

\begin{align*}
&\|x_{k+1}- x^*\|^2 \le(1+4s_k^2M_1^2)\|x_k - x^*\|^2 \cr
&\ \
+2s_k\left((1-\rho) B +\rho( M_1B+M_2)\right)\dist(x_k,X)\cr
&\ \
+2s_k\rho(f^*- f(\Pi_X[x_k]))
+2s_k\rho M_1\|x_k-x^*\|\,\dist(x_k,X)\cr
& \ \ +\frac{s_k\g_k \d_k}{2\a_{\min}}- 2s_k\g_k\dist(x_k,X_{i_k})
+4s_k^2(M_1^2 B^2 +M_2^2+\g_k^2),\qquad
\end{align*}

We will take the conditional expectation in the preceding relation with respect to $\EuScript F_k=\{x_1,\ldots,x_k\}$. In doing so, we note that, given $x_k$, the index
$i_k$ is uniformly distributed over $\{1,\ldots,m\}$, for all $k\ge1$, implying that

\[\ex{\dist(x_k,X_{i_k})\mid \EuScript F_k} =\frac{1}{m}\sum_{i=1}^m\dist(x_k,X_i).\]

By using Hoffman's Lemma (Lemma <ref>), we obtain

\[\ex{\dist(x_k,X_{i_k})\mid \EuScript F_k}\ge \frac{1}{\b m}\,\dist(x_k,X).\]

Therefore, it follows that  for all $x^*\in X^*$, all $\rho\in[0,1]$, and all $k\ge 1$,

\begin{align}\label{eq-keys}
\ex{\|x_{k+1}- x^*\|^2 \mid \EuScript F_k}\le(1+4s_k^2M_1^2)\|x_k - x^*\|^2 \cr
&\ \
+2s_k\left((1-\rho) B +\rho( M_1B+M_2)\right)\dist(x_k,X)\cr
&\ \
+2s_k\rho(f^*- f(\Pi_X[x_k]))
+2s_k\rho M_1\|x_k-x^*\|\,\dist(x_k,X)\cr
& \ \ +\frac{s_k\g_k \d_k}{2\a_{\min}}- \frac{2s_k\g_k}{\b m}\,\dist(x_k,X)
+4s_k^2(M_1^2 B^2 +M_2^2+\g_k^2).\qquad
\end{align}

The proof proceeds in two major steps both using relation (<ref>): (Step 1) we show that
$\sum_{k=1}^\infty s_k\dist(x_k,X)<\infty$ , and (Step 2) we show the almost sure convergence of the method by means of Lemma <ref>.

Step 1:
In relation (<ref>), we let $x^*\in X^*$ be arbitrary but fixed, and we set
$\rho=0$. In the resulting relation, we group the common coefficients with
$\dist(x_k,X)$ and obtain  for all $k\ge 1$,

\begin{align*}
\ex{\|x_{k+1}- x^*\|^2 \mid \EuScript F_k}\le(1+4s_k^2M_1^2)\|x_k - x^*\|^2
+\frac{s_k\g_k \d_k}{2\a_{\min}}\cr
&\ \ - 2s_k\left(\frac{\g_k}{\b m} - B\right)\,\dist(x_k,X)
+4s_k^2(M_1^2 B^2 +M_2^2+\g_k^2).\qquad
\end{align*}

By Assumption <ref>, we have that $\lim_{k\to\infty}\g_k=\infty$, implying that
for some sufficiently large $k_1\ge1$, we have $\frac{\g_k}{m\b} - B\ge B$. Therefore,
for all $k\ge k_1$ ,

\begin{align*}
&\ex{\|x_{k+1}- x^*\|^2 \mid \EuScript F_k}\le(1+4s_k^2M_1^2)\|x_k - x^*\|^2
+\frac{s_k\g_k \d_k}{2\a_{\min}}\cr
&\ \  - 2s_kB\dist(x_k,X)
+4s_k^2(M_1^2B^2 +M_2^2+\g_k^2).
\end{align*}

Under Assumption <ref>,
the preceding relation
satisfies all the conditions of Lemma <ref>, starting at time $k_1$,
with $v_k=\|x_k - x^*\|^2$,
$u_k=2s_kB\dist(x_k,X)$, and
$\b_k=\frac{s_k\g_k \d_k}{2\a_{\min}} +4s_k^2(M_1^2B^2 +M_2^2+\g_k^2)$.
Thus, by Lemma <ref>,

\begin{equation}\label{eq-step1}
\sum_{k=1}^\infty s_k\dist(x_k,X)<\infty,\qquad\as
\end{equation}

Step 2:
We now use relation (<ref>) with $\rho=1$
and obtain  that  for all $x^*\in X^*$ and all $k\ge 1$,

\begin{align*}
&\ex{\|x_{k+1}- x^*\|^2 \mid \EuScript F_k}\le(1+4s_k^2M_1^2)\|x_k - x^*\|^2 \cr
&\ \ +2s_k(f^*- f(\Pi_X[x_k]))
+2s_k M_1\|x_k-x^*\|\,\dist(x_k,X)\cr
& \ \ +\frac{s_k\g_k \d_k}{2\a_{\min}}- 2s_k\left(\frac{\g_k}{\b m}-(M_1B+M_2)\right)\,\dist(x_k,X)
+4s_k^2(M_1^2 B^2 +M_2^2+\g_k^2).\qquad
\end{align*}

By Assumption <ref>, we have that $\lim_{k\to\infty}\g_k=\infty$. Thus,
for some sufficiently large $k_2\ge1$, we have $\frac{\g_k}{m\b} - (M_1B+M_2)\ge 0$ and $M_1^2 B^2 +M_2^2\le \g_k$ for all $k\ge k_2$. Using this and
noting that $2\|x_k-x^*\|\le 1+\|x_k-x^*\|^2$, we find that  for all $x^*\in X^*$ and all $k\ge k_2$,

\begin{align}\label{eq-keys1}
\ex{\|x_{k+1}- x^*\|^2 \mid \EuScript F_k}&\le\left(1+4s_k^2M_1^2+s_k M_1\,\dist(x_k,X)\right)
\|x_k - x^*\|^2 \cr
&\ \ +2s_k(f^*- f(\Pi_X[x_k]))
+s_k M_1\,\dist(x_k,X) \cr
&\ \ +\frac{s_k\g_k \d_k}{2\a_{\min}}
\end{align}

By Assumption <ref> we have
$\sum_{k=1}^\infty s_k=\infty$.
By Assumption <ref>
and the fact that $\sum_{k=1}^\infty s_k\dist(x_k,X)<\infty$  (see  (<ref>)),
it follows that relation (<ref>)
satisfies the conditions of Lemma <ref>
with $z_k=x_k$, $f(\cdot)=\phi(\cdot)$, and

\[a_k=4s_k^2M_1^2+s_k M_1\,\dist(x_k,X),
\qquad
b_k=2s_k,\qquad b'_k=s_k,\]

\[c_k=s_k M_1\,\dist(x_k,X) +\frac{s_k\g_k \d_k}{2\a_{\min}}+8s_k^2\g_k^2.\]

Hence, the results follow from Lemma <ref>.

Next we provide some sufficient conditions for the parameters to satisfy Assumption <ref> and formulate their properties which we will use in the forthcoming convergence rate analysis.

Let us consider the following setting for the parameters:
\[\mbox{$s_k = O\left(\frac{1}{k^c\ln^{(1+3g)/2}(k+1)}\right)$, $\gamma_k = O\left(\ln^g (k+1)\right)$, $\delta_k =  O\left(\frac{1}{k^d}\right)$ with $c,g,d>0$}.\]

Then, if $c\in[1/2,1)$ and $d> 1/2$, Assumption <ref> holds. Moreover, given this choice of the parameters,
$\sum_{k=1}^{t} s_k \ge \frac{(t+1)^{1-c}-1}{(1-c)\ln^{(1+3g)/2}(t+1)}$,
$\sum_{k=1}^{t} s^2_k = O(1)$,  $\sum_{k=1}^{t} s_k\gamma_k\delta_k = O(1)$,  $\sum_{k=1}^{t} s^2_k\gamma^2_k = O(1)$ as $t\to\infty$.

Obviously, given the settings, $\lim_{k\to\infty}\gamma_k = \lim_{k\to\infty} \ln^g (k+1) = \infty$.

\begin{align*}
\sum_{k=1}^{t} s_k = \sum_{k=1}^{t} \frac{1}{k^c\ln^{(1+3g)/2}(k+1)} &\ge \frac{1}{\ln^{(1+3g)/2}(t+1)}\int_{1}^{t+1}\frac{dx}{x^c}\cr& = \frac{(t+1)^{1-c}-1}{(1-c)\ln^{(1+3g)/2}(t+1)}.
\end{align*}

On the other hand,

\begin{align*}
\sum_{k=1}^{t} s_k = \sum_{k=1}^{t} \frac{1}{k^c\ln^{(1+3g)/2}(k+1)} &\le \frac{1}{\ln^{(1+3g)/2}(2)}\left(1+\int_{1}^{t}\frac{dx}{x^c}\right)\cr& = \frac{t^{1-c}-c}{(1-c)\ln^{(1+3g)/2}(2)}.
\end{align*}

Thus, $\sum_{k=1}^{t} s_k = O\left(\frac{t^{1-c}}{\ln^{(1+3g)/2}(2)}\right)$, which implies also the condition $\sum_{k=1}^{\infty} s_k = \infty$. Next,

\begin{align*}
\sum_{k=1}^{t} s^2_k\gamma^2_k = \sum_{k=1}^{t} \frac{1}{k^{2c}\ln^{1+g}(k+1)} &\le \frac{1}{\ln^{1+g}2}+\int_{1}^{t}\frac{dx}{x\ln^{1+g}(x+1)} \cr&\le \frac{1}{\ln^{1+g}2}+\frac{g}{\ln^g 2}.
\end{align*}

Hence, $\sum_{k=1}^{t} s^2_k\le\sum_{k=1}^{t} s^2_k\gamma^2_k = O(1)$ and the corresponding series converge.

\begin{align*}
\sum_{k=1}^{t} s_k\gamma_k\delta_k = \sum_{k=1}^{t} \frac{1}{k^{c+d}\ln^{(1+g)/2}(k+1)} &\le \frac{1}{\ln^{(1+g)/2}(2)}\left(1+\int_{1}^{t}\frac{dx}{x^{c+d}}\right)\cr& = \frac{t^{1-c-d}-c-d}{(1-c-d)\ln^{(1+3g)/2}(2)}.
\end{align*}

Thus, taking into account that $c+d>1$, we conclude that $\sum_{k=1}^{t} s_k\gamma_k\delta_k = O(1)$ and $\sum_{k=1}^{\infty} s_k\gamma_k\delta_k<\infty$.

§.§ Convergence Rate Results

In this section,
we establish convergence rate results for the iterates of the method (<ref>), under
the following condition:
there is a scalar $M>0$ such that

\begin{align}\label{eq:gradBound}
\|\tilde\grad f(x_k)\|\le M,\qquad
\|\tilde\grad f(\Pi_X[x_k])\|\le M
\qquad\hbox{for all  $k\ge1$}\ \ \as
\end{align}

This condition is satisfied under Assumptions <ref>–<ref>. Specifically, as seen in Proposition <ref>, the norms of iterates are  bounded by a constant.
As a consequence, the subgradients $\tilde\grad f(x_k)$
as well as the subgradients at the projection points $\tilde\grad f(\Pi_X[x_k])$ are also  bounded by a constant.

In what follows, we allow for the strong convexity constant to take value 0 ($\m=0$) in order to provide a unified treatment of both cases when $f(\cdot)$ is strongly convex and when it is just convex.

Using relation (<ref>), we prove the following relations which will be useful in the further analysis.

Let $f(\cdot)$ be strongly convex with $\mu\ge 0$ and let the condition in (<ref>) hold.

Then, the following relations are  valid for the iterates $x_k$ of the method (<ref>)
for all $k\ge1$,

$f(\Pi_X[x_k])-f(x_k)\le M \dist(x_k,X)-\frac{\mu}{2}\dist^2(x_k,X),$

$f(\Pi_X[x_k])-f(x_k)\ge -M \dist(x_k,X)+\frac{\mu}{2}\dist^2(x_k,X).$

Here, $\mu\ge0$ is the strong convexity constant for $f(\cdot)$ and $M>0$ is subgradient norm bound along the iterates and their projections on $X$ from (<ref>).

Let $k\ge 1$ be arbitrary.
By using the convexity of $f(\cdot)$ (see (<ref>)),
we estimate $f(\Pi_X[x_k])-f(x_k)$ as follows:

\begin{align*}
&\le \la \tilde \nabla f(\Pi_X[x_k]), \Pi_X[x_k]- x_k\ra-\frac{\mu}{2}\|\Pi_X[x_k]-x_k\|^2\cr
&\le M \|\Pi_X[x_k]-  x_k\|-\frac{\mu}{2}\|\Pi_X[x_k]-x_k\|^2,
\end{align*}

where the last inequality follows by the Cauchy-Schwarz inequality and (<ref>).
By  using $\|\Pi_X[x_k]-x_k\|=\dist(x_k,X)$ we obtain the relation in part (a).
The relation in part (b) is obtained similarly by using (<ref>), where
the roles of $x_k$ and $\Pi_X[x_k]$ are exchanged.

The next lemma provides an iterate relation that is crucial in our convergence rate analysis of the method.
In the lemma, we use a scalar  $\eta\in[0,1]$ in order to estimate $f(y)-f(x_k)$ by a convex combination
of $f(y)-f(x_k)$ and $f(y)-f(\Pi_X[x_k])$ plus some additional terms.
We will use $\eta=0$ or $\eta=1$,
depending on which quantities we want to estimate.
The lemma is a refinement of Lemma <ref>,
which possible due to the condition (<ref>) and Lemma <ref>.

Let $f(\cdot)$ be strongly convex with $\mu\ge 0$ and let $\g_k>0$, $\d_k>0$, and $s_k>0$ for all $k$. Also, let the condition in (<ref>) hold. Then, the iterate sequence $\{x_k\}$ of the method (<ref>)  satisfies
for all $y\in X$ and $k\ge1$,

\begin{align*}
\ex{\|x_{k+1}- y\|^2\mid \EuScript F_k}
&\le(1-\mu s_k)\|x_k - y\|^2
+2s_k(1-\eta)(f(y) - f(x_k))\cr
&\quad+2s_k\eta (f(y) - f(\Pi_X[x_k]))\cr
&\quad +\frac{s_k\g_k \d_k}{2\a_{\min}}
- 2s_k\left(\frac{\g_k}{m\b} -\eta M\right)\dist(x_k,X)\cr
&\quad - \eta\mu s_k\dist^2(x_k,X)
\end{align*}

where $\eta\in[0,1]$ is arbitrary,

$\beta>0$ is the Hoffman constant from Lemma <ref>, and
$M>0$ is subgradient norm bound from (<ref>).

By Lemma <ref>, we surely have for all
$y\in X$ and all $k\ge1$,

\begin{align*}
&\|x_{k+1}- y\|^2 \le(1-\mu s_k)\|x_k - y\|^2
&+\frac{s_k\g_k \d_k}{2\a_{\min}} - 2s_k\g_k\dist(x_k,X_{i_k})
+s_k^2(\|\tilde\nabla f(x_k)\| + \g_k)^2.
\end{align*}

By the condition in (<ref>)
the subgradients of $f(x_k)$ are bounded , so that $\|\tilde\nabla f(x_k)\| \le M$ Hence, we have  for all $y\in X$ and all $k\ge1$,

\begin{align}\label{eq-rel-perk1}
&\|x_{k+1}- y\|^2 \le(1-\mu s_k)\|x_k - y\|^2
+2s_k(f(y)- f(x_k)) \cr
& +\frac{s_k\g_k \d_k}{2\a_{\min}}- 2s_k\g_k\dist(x_k,X_{i_k})
\end{align}

Taking the conditional expectation with respect to $\EuScript F_k$ in relation (<ref>),
we obtain  that   for all $y\in X$ and all $k\ge1$,

\begin{align*}
&\ex{\|x_{k+1}- y\|^2\mid \EuScript F_k} \le
(1-\mu s_k) \|x_k - y\|^2 +2s_k(f(y) - f(x_k))\cr
& \ \ +\frac{s_k\g_k \d_k}{2\a_{\min}} - 2s_k\g_k \ex{\dist(x_k,X_{i_k})\mid \EuScript F_k}
\end{align*}

Given $x_k$, the index $i_k$ is uniformly distributed over $\{1,\ldots,m\}$,  so it follows that
$\ex{\dist(x_k,X_{i_k})\mid \EuScript F_k} =\frac{1}{m}\sum_{i=1}^m\dist(x_k,X_i),$
and by using Hoffman's Lemma (Lemma <ref>), we obtain

\[\ex{\dist(x_k,X_{i_k})\mid \EuScript F_k}\ge \frac{1}{\b m}\dist(x_k,X).\]

Therefore, we have   for all $y\in X$ and all $k\ge1$,

\begin{align}\label{eq-prima}
& \ex{\|x_{k+1}- y\|^2\mid \EuScript F_k} \le
(1-\mu s_k) \|x_k - y\|^2 +2s_k(f(y) - f(x_k))\cr
& \ \ +\frac{s_k\g_k \d_k}{2\a_{\min}}
- \frac{2s_k\g_k}{m\b} \dist(x_k,X)+2s_k^2(M^2+\g_k^2).
\end{align}

Let $\eta\in[0,1]$ and lets write

\[f(y) - f(x_k)=(1-\eta)(f(y) - f(x_k))+\eta(f(y) - f(x_k)).\]

Further, we write

\[\eta(f(y) - f(x_k))=\eta(f(y) - f(\Pi_X[x_k] +f(\Pi_X[x_k]) - f(x_k)).\]

By Lemma <ref>(a), we have

\[f(\Pi_X[x_k]) - f(x_k)\le M \dist(x_k,X)-\frac{\mu}{2}\dist^2(x_k,X).\]

Combining the preceding relations yields

\begin{align*}
f(y) - f(x_k)&\le(1-\eta)(f(y) - f(x_k)) + \eta(f(y) - f(\Pi_X[x_k]) \cr
&+ \eta M \dist(x_k,X)- \frac{\eta\mu}{2}\dist^2(x_k,X).
\end{align*}

By substituting the preceding estimate back in (<ref>), we obtain the stated relation.

§.§.§ Convergence rate in merely convex case

Here, we provide convergence rate result for the method in the case when $f(\cdot)$ is merely convex, i.e.,
$\mu=0$ in relation (<ref>).
To obtain a convergence rate estimate for the method, we will
consider the weighted averages and the truncated weighted averages of the iterates. We define the weighted averages with respect to a generic sequence $\{\nu_k\}$  of positive scalars.
Given a positive sequence $\{\nu_k\}$, we define the $\nu$-weighted averages, as follows:

\begin{align}\label{eq-waver}
\xtav &=S_t^{-1}\sum_{k=1}^t \nu_k x_k,\quad S_t =\sum_{k=1}^t \nu_k,\qquad\hbox{for all }t\ge 1.
\end{align}

Let us focus on $\xtav$ defined in (<ref>) as follows.
First, let us note that the truncated $\nu$-weighted average $\xttav$ of the iterates $x_k$ for $k=\tau,\ldots,t$, is defined by

\begin{align}\label{eq-truncwaver}
x_{\tau, t}^{\textrm av}=(S_t-S_\tau)^{-1}\sum_{k=\tau}^t \nu_k x_k\qquad\hbox{for all }t\ge \tau\ge 1.
\end{align}

The $\nu$-weighted average $\xtav$ can be related to
the average $x_\tau^{\textrm av}$, with $\tau\le t$, and the truncated $\nu$-weighted average $\xttav$. Specifically, for any $t\ge \tau$, we have

\[\xtav=\frac{\sum_{k=1}^t \nu_k x_k}{\sum_{k=1}^t \nu_k}
=\frac{\sum_{k=1}^\tau \nu_k x_k +\sum_{k=\tau}^t \nu_k x_k}{\sum_{k=1}^t \nu_k}.\]

By using the definitions of $S_t$, $x_\tau^{\textrm av}$, and $\xttav$ i
n (<ref>)–(<ref>), it follows that

+ (S_t-S_τ)/S_t    for all t≥τ≥1.Thus, $\xtav$ is a convex combination of $x_{\tau}^{\textrm av}$ and
the truncated-weighted average $\xttav$.

The following lemma provides some preliminary estimates regarding $\xtav$.

Let $f(\cdot)$ be convex, and
let the condition in (<ref>) hold. Assume that the problem $\min_{x\in X}f(x)$ has an optimal solution. Also,
assume that $s_k>0$, $\g_k>0$, $\d_k>0$, and
$\g_{k+1}\ge \g_k$ for all $k$, and $\lim_{k\to\infty}\g_k=\infty$.
Let $\{x_k\}$ be generated by the method (<ref>). Consider
the $s$-weighted averages $\{\xtav\}$ of the iterates obtained by using
$\nu_k=s_k$ in (<ref>) for all $k\ge1$.
Then, the following estimates are valid:
for all $t\ge1$,

\[\ex{f(\xtav)}-f^*\le \frac{\ex{\|x_1 - x^*\|^2}}{2S_t}
+\frac{\sum_{k=1}^t\b_k}{2S_t}\ \hbox{for  all $x^*\in X^*$},\]

where $\b_k=\frac{s_k\g_k \d_k}{2\a_{\min}}+2s_k^2(M^2+\g_k^2)$
for all $k$, with $M$ being the constant from condition (<ref>).
Moreover, there exists $k_1\ge 1$ such that for all $t\ge k_1$,

\begin{align*}
\ex{\dist(\xtav,X)}&\le
\frac{S_{k_1}}{S_t}\ex{\dist(x_{k_1}^{\textrm av},X)}
+\frac{m\b\sum_{k=k_1}^t \g_k^{-1}\b_k}{S_t}.\end{align*}

In Lemma <ref>,
we let $\mu=0$ and $\eta=0$, and thus obtain
for all $y\in X$ and $k\ge1$,

\begin{align*}
&\ex{\|x_{k+1}- y\|^2\mid \EuScript F_k} \le\|x_k - y\|^2 +2s_k(f(y) - f(x_k))\cr
& +\frac{s_k\g_k \d_k}{2\a_{\min}}- \frac{2s_k\g_k}{m\b} \dist(x_k,X)+2s_k^2(M^2+\g_k^2).
\end{align*}

By taking the total expectation
in the preceding relation, we obtain for all $y\in X$ and $k\ge1$,

\begin{align}\label{eq-777}
&\ex{\|x_{k+1}- y\|^2} \le \ex{\|x_k - y\|^2} +2s_k\ex{f(y) - f(x_k)}
- \frac{2s_k\g_k}{m\b} \ex{\dist(x_k,X)}+\b_k,\qquad
\end{align}

$\b_k=\frac{s_k\g_k \d_k}{2\a_{\min}}+2s_k^2(M^2+\g_k^2)$.
By summing these relations over $k=1,\ldots, t$ and
by dropping the distance-related terms, after re-arranging the remaining terms,
we obtain for all $y\in X$ and all $t\ge 1,$

\begin{align*}
&2 \sum_{k=1}^t s_k\left(\ex{f(x_k)}- f(y)\right)
+\ex{\|x_{t+1}- y\|^2} \le\ex{\|x_1 - y\|^2}+\sum_{k=1}^t\b_k.
\end{align*}

Next, we divide the preceding relation with $2 S_t$, where $S_t=\sum_{k=1}^t s_k$, and we use the definition of $\xtav$  (i.e., (<ref>) with $\nu_k=s_k$).
This and the convexity of $f(\cdot)$ yield for $y=x^*$, where $x^*\in X^*$ is arbitrary,
and for all $t\ge1$,

\[\ex{f(\xtav)}-f^*\le \frac{\ex{\|x_1 - x^*\|^2}}{2S_t}

To show the estimate for $\ex{\dist(\xtav,X)}$, we use
relation (<ref>) with $y=\Pi_X[x_k]\in X$, and
obtain for all $k\ge1$,

\begin{align*}
&\ex{\|x_{k+1}- \Pi_X[x_k]\|^2} \le\ex{\|x_k - \Pi_X[x_k]\|^2} \cr
&+2s_k\ex{f(\Pi_X[x_k]) - f(x_k)}
- \frac{2s_k\g_k}{m\b} \ex{\dist(x_k,X)}+\b_k.
\end{align*}

Since $\ex{\dist(x_{k+1},X)}\le \ex{\|x_{k+1}- \Pi_X[x_k]\|}$
and $\|x_k - \Pi_X[x_k]\|=\dist(x_k,X)$ it follows that
for all $k\ge1$,

\begin{align*}
\ex{\dist^2(x_{k+1},X)} \le&\ex{\dist^2(x_k,X)} +2s_k\ex{f(\Pi_X[x_k]) - f(x_k)}\cr
&- \frac{2s_k\g_k}{m\b}\ex{\dist(x_k,X)}+\b_k.
\end{align*}

By Lemma <ref>(a), where $\mu=0$,
we have that

\[f(\Pi_X[x_k]) - f(x_k)\le M\dist(x_k,X),\]

implying that  for all $k\ge1$,

\begin{align*}
\ex{\dist^2(x_{k+1},X)} \le\ex{\dist^2(x_k,X)} -2s_k\left(\frac{\g_k}{m\b}-M\right)\ex{\dist(x_k,X)}+\b_k.
\end{align*}

Since $\g_k\to+\infty$, there exists a large enough $k_1$ such that
$\g_k/(m\b)-M\ge \g_k/(2m\b)$, implying that for all $k\ge k_1$,

\begin{align*}
\ex{\dist^2(x_{k+1},X)} \le\ex{\dist^2(x_k,X)} -\frac{s_k\g_k}{m\b}\ex{\dist(x_k,X)}+\b_k.
\end{align*}

Dividing the preceding relation with $\g_k$ and using the assumption that
$\g_{k+1}\ge \g_k$, we obtain for all $k\ge k_1$,

\begin{align*}
\g_{k+1}^{-1}\ex{\dist^2(x_{k+1},X)} &\le\g_k^{-1}\ex{\dist^2(x_k,X)} -\frac{s_k}{m\b}\ex{\dist(x_k,X)}+\g_k^{-1}\b_k.
\end{align*}

By summing these relations over $k=k_1,\ldots,t,$
after re-arranging the terms, we can see that for all $t\ge k_1$, we have

\begin{align*}
\sum_{k=k_1}^t \frac{s_k}{m\b}\ex{\dist(x_k,X)} &\le\g_{k_1}^{-1}\ex{\dist^2(x_{k_1},X)}
\end{align*}

Since $X$ is convex, the distance function $\dist(\cdot,X)$ is convex, and
by dividing the preceding relation with $\sum_{k=k_1}^t s_k$, we obtain
for all $t\ge k_1$,

\begin{align}\label{eq-88}
\frac{1}{m\b}\ex{\dist(\xtktav,X)}
\le\frac{\g_{k_1}^{-1}\ex{\dist^2(x_{k_1},X)}}{S_t-S_{k_1}}
+\frac{\sum_{k=k_1}^t \g_k^{-1}\b_k}
\end{align}

where we used the notation $S_t=\sum_{k=1}^t s_k$.

By expressing $\xtav$ as a convex combination of $x_{k_1}^{\textrm av}$ and $\xtktav$, i.e., by using (<ref>) with $\nu_k=s_k$ and $\tau=k_1$,

\[\xtav=\frac{S_\tau x_{\tau}^{\textrm av}
+ (S_t-S_\tau)\xttav}{S_t}\qquad\hbox{for all }t\ge k_1\ge 1.\]

Then, by using the convexity of the distance function $\dist(\cdot,X)$,
we find that for all $t\ge k_1$,

\begin{align*}
\ex{\dist(\xtav,X)}
\le
\frac{S_{k_1}}{S_t}\ex{\dist(x_{k_1}^{\textrm av},X)} \ +\frac{S_t-S_{k_1}}{S_t}\ex{\dist(\xtktav,X)}.\end{align*}

The preceding relation and (<ref>) yield for all $t\ge k_1$,

\begin{align*}
\ex{\dist(\xtav,X)}\le&
\frac{S_{k_1}}{S_t}\ex{\dist(x_{k_1}^{\textrm av},X)}\cr &+\frac{m\b\g_{k_1}^{-1}}{S_t}\ex{\dist^2(x_{k_1},X)}
+\frac{m\b\sum_{k=k_1}^t \g_k^{-1}\b_k}{S_t}.\end{align*}

The following proposition provides the convergence rate for the expected function values estimated at $\xtav$.

Let $f(\cdot)$ be convex, and
let the condition in (<ref>) hold. Assume that the problem $\min_{x\in X}f(x)$ has an optimal solution. Also, assume that
\[\mbox{$s_k = O\left(\frac{1}{k^c\ln^{(1+3g)/2}(k+1)}\right)$, $\gamma_k = O\left(\ln^g (k+1)\right)$, $\delta_k =  O\left(\frac{1}{k^d}\right)$ }\]

with $g>0$, $c\in[1/2,1)$ and $d> 1/2$.
Let $\{x_k\}$ be generated by the method (<ref>). Consider
the $s$-weighted averages $\{\xtav\}$ of the iterates obtained by using
$\nu_k=s_k$ in (<ref>) for all $k\ge1$.
Then, the following estimate is valid, as $t\to\infty$,

\[|\ex{f(\xtav)}-f^*| = O\left(\frac{\ln^{1/2+3g/2} t}{t^{1-c}}\right).\]

Let us notice that the setting for the parameters $s_k$, $\gamma_k$, and $\delta_k$ above is the same as one in Lemma <ref>.
Thus, by combining the result provided by that lemma with the result for function values
from Lemma <ref>, namely,

\[\ex{f(\xtav)}-f^*\le \frac{\ex{\|x_1 - x^*\|^2}}{2S_t}
+\frac{\sum_{k=1}^t\b_k}{2S_t}\quad \hbox{for  all $x^*\in X^*$},\]

where $\b_k=\frac{s_k\g_k \d_k}{2\a_{\min}}+2s_k^2(M^2+\g_k^2)$, we conclude that there exists a scalar $C_1>0$ such that for sufficiently large $t$

\begin{align}\label{eq:upper}
\ex{f(\xtav)}-f^*\le \frac{C_1 \ln^{1/2+3g/2} t}{t^{1-c}}.
\end{align}

The constant $C_1$ depends on $M^2$ (where $M$ is from the condition (<ref>)), the smallest norm $\a_{\min}$ of the vectors $a_i,i\in[m]$, and the squared distance between the initial point and the corresponding solution $\ex{\|x_1 - x^*\|^2}$.

On the other hand, using the result on expected distances of Lemma <ref> and the inequality $\ex{f(\xtav)}-f^*
\ge -\|\tilde \nabla f(x^*)\|\,\ex{\dist(\xtav,X)}$ (see (<ref>)), we conclude that there exists finite $k_1$ such that for all $t>k_1$,

\begin{align*}&\ex{f(\xtav)}-f^*
\ge -\|\tilde \nabla f(x^*)\|\,\ex{\dist(\xtav,X)} \cr
& \ge -\|\tilde \nabla f(x^*)\|\left(\frac{S_{k_1}}{S_t}\ex{\dist(x_{k_1}^{\textrm av},X)}+\frac{m\b\g_{k_1}^{-1}}{S_t}\ex{\dist^2(x_{k_1},X)}
+\frac{m\b\sum_{k=k_1}^t \g_k^{-1}\b_k}{S_t}\right),
\end{align*}

where $\tilde \nabla f(x^*)$ is a subgradient of $f(\cdot)$ at some optimal point $x^*\in X^*$.
By using Lemma <ref> again, we conclude the existence of some constant $C_2>0$ such that

\begin{align}\label{eq:lower}
\ex{f(\xtav)}-&f^*\ge -\frac{C_2\ln^{1/2+3g/2} t}{t^{1-c}}.
\end{align}

The constant $C_2$ in its turn depends on the index $k_1$, the upper bound $M$ of the subgradient norms, the number of the constraints $m$, the value $\a_{\min}$, Hoffman's constant $\beta$, and distance between $x_{k_1}$ and the feasible set $X$.
Combining (<ref>) and (<ref>), we obtain the stated result.

By optimizing the parameters $c$ and $d$ in Proposition <ref>,
the convergence rate of the order $O\left(\frac{\ln^{1/2+\epsilon} t}{t^{1/2}}\right)$ for any small positive $\varepsilon$ is obtained for $c=1/2$ and $g = 2\epsilon/3$.
The constant involved in $O$-notation in Proposition <ref> can be explicitly derived using Lemma <ref>, and it linearly depends on $M^2$ (with $M$ from (<ref>)), the Hoffman constant $\beta$, the number of the constraints $m$, the smallest inverse norm $\a^{-1}_{\min}$ of the vectors $a_i,i\in[m]$, and the expected squared distance from the initial point $x_1$ to the solution set $X^*$.

§.§.§ Convergence rate in strongly convex case

Here, we analyze the convergence rate of the method (<ref>) for strongly convex
$f(\cdot)$ with $\mu>0$ (see (<ref>)). We note that in this case the original
problem of minimizing $f(x)$ over $x\in X$ has a unique solution (as $X$ is assumed to be nonempty).
To establish the convergence rate of the method, we will consider the $s^{-1}$-
weighted averages $\xtav$ defined by (<ref>) with $\nu_k=s_k^{-1}$.
To differentiate these averages from those used in the preceding section, we define

\begin{align}\label{eq-waver-new}
\bxtav &=\bar {S}_t^{-1}\sum_{k=1}^t s_k^{-1} x_k\quad, \bar{S}_t=\sum_{k=1}^t s_k^{-1},\qquad\hbox{for all }t\ge 1.
\end{align}

The $s^{-1}$-weighted averages of the form (<ref>) have been proposed in [24] to analyze
the convergence rate of a stochastic gradient method in terms of the expected function values.
As seen in [24], using the $s^{-1}$-weighted averages for a stochastic gradient method applied to
minimizing merely convex function $f(\cdot)$ yields
the convergence rate of $O(1/\sqrt{k})$, where $k$ is the number of iterations.
However, in this case, an additional assumption that the iterates are bounded is needed, which is the main reason why we did not consider the $s^{-1}$-weighted averages in the preceding section.

The $s^{-1}$-weighted averages of the form (<ref>) have the convergence rate of $O(1/k)$
for a stochastic subgradient method as applied to a strongly convex function $f(\cdot)$, as shown in [24]. Unlike the setting in [24], here in addition to time-varying functions,  we have to deal with the infeasibility of the iterates.

Using the weighted average $\bxtav$ given in (<ref>), we define
the truncated weighted average $\bxttav$ for $k=\tau,\ldots,t$, as follows

\begin{align}\label{eq-truncwaver-new}
\bxttav=(S_t-S_\tau)^{-1}\sum_{k=\tau}^t \nu_k x_k\qquad\hbox{for all }t\ge \tau\ge 1.
\end{align}
with $\nu_k=s_k^{-1}$ for all $k$.
Note that the $s^{-1}$-weighted average $\bxtav$ and the truncated
weighted average $\bar{x}_\tau^{\textrm av}$, with $\tau\le t$, satisfy the relation

+ (S̅_t-S̅_τ)/S̅_t    for all t≥τ≥1.A preliminary result regarding the function value at $\bxtav$  is given in the following lemma.

Let $f(\cdot)$ be strongly convex with $\mu>0$, and let the condition in (<ref>)) hold.
Assume that $s_k=\frac{2}{\mu k}$.

Also, assume that $\g_k>0$, and
Let $\{x_k\}$ be the iterate sequence  generated by method (<ref>).
Then, the following estimates are valid for the $s^{-1}$-weighted averages $\bxtav$:
for all $t\ge1$,

\[\ex{f(\bxtav)}-f^*\le \frac{s_0^{-2}\ex{\|x_1 - x^*\|^2}}{2\bar{S}_t}
+\frac{\sum_{k=1}^t c_k}{2\bar{S}_t},\]

where $s_0=2\mu^{-1}$,
$c_k=\frac{\g_k \d_ks_k^{-1}}{2\a_{\min}}+2(M^2+\g_k^2)$
for all $k$, with $M$ being the constant from (<ref>).
Moreover, there exists $k_1\ge 1$ such that for all $t\ge k_1$,

\begin{align*}
\ex{\dist(\bxtav,X)}\le
\frac{\bar{S}_{k_1}}{\bar{S}_t}\ex{\dist(\bar{x}_{k_1}^{\textrm av},X)}
&+\frac{m\b\sum_{k=k_1}^t \g_k^{-1}c_k}{\bar{S}_t}.\end{align*}

In Lemma <ref>, where $\mu>0$, we let $\eta=0$, and thus obtain
for all $y\in X$ and $k\ge1$,

\begin{align*}
\ex{\|x_{k+1}- y\|^2\mid \EuScript F_k}
&\le(1-\mu s_k)\|x_k - y\|^2
+2s_k(f(y) - f(x_k))+\frac{s_k\g_k \d_k}{2\a_{\min}}\cr
& \quad- \frac{2s_k\g_k}{m\b} \dist(x_k,X)+2s_k^2(M^2+\g_k^2).
\end{align*}

By taking the total expectation and dividing with $s_k^2$, we obtain for all $y\in X$ and $k\ge1$,

\begin{align*}
&s_k^{-2}\ex{\|x_{k+1}- y\|^2}
\le s_k^{-2}(1-\mu s_k)\ex{\|x_k - y\|^2} \cr
&+2s_k^{-1}\ex{f(y) - f(x_k)}
- \frac{2s_k^{-1}\g_k}{m\b} \ex{\dist(x_k,X)}
\end{align*}

\[c_k=\frac{\g_k \d_ks_k^{-1}}{2\a_{\min}}+2(M^2+\g_k^2).\]

Note that for the given stepsize $s_k=\frac{2}{\mu k}$,
we have for sufficiently large $k$

\[ s_k^{-2}(1-\mu s_k)=\frac{\mu^2k(k-2)}{4}
\le \frac{\mu^2 (k-1)^{2}}{4}=s_{k-1}^2.\]

Therefore, we obtain for all $y\in X$ and sufficiently large $k\ge1$,

\begin{align}\label{eq-77ag}
&s_k^{-2}\ex{\|x_{k+1}- y\|^2}
\le s_{k-1}^{-2}\ex{\|x_k - y\|^2} \cr
&+2s_k^{-1}\ex{f(y) - f(x_k)}
- \frac{2s_k^{-1}\g_k}{m\b} \ex{\dist(x_k,X)}
\end{align}

By summing the relations in (<ref>) over $k=1,\ldots, t$ and omitting the terms with $\ex{\dist(x_k,X)}$,
after re-arranging the terms,
we find that for all $y\in X$ and all $t\ge 1,$

\begin{align*}
2 \sum_{k=1}^t s_k^{-1}\left(\ex{f(x_k)}- f(y)\right)
+s_t^{-2}\ex{\|x_{t+1}- y\|^2}
\le s_0^{-2}\ex{\|x_1 - y\|^2}+\sum_{k=1}^t c_k,
\end{align*}

where we define $s_0=2\mu^{-1}$.
By dividing with $2\bar S_t$, where $\bar{S}_t=\sum_{k=1}^t s_k^{-1}$, using the definition of
$\bxtav$ (see (<ref>)) and the convexity of $f(\cdot)$,
we obtain for $y=x^*$ and all $t\ge1$,

\[\ex{f(\bxtav)}-f^*\le \frac{s_0^{-2}\ex{\|x_1 - x^*\|^2}}{2\bar{S}_t}
+\frac{\sum_{k=1}^t c_k}{2\bar{S}_t}.\]

Next, we establish the estimate for $\ex{\dist(\bxtav,X)}$. We use
relation (<ref>) with $y=\Pi_X[x_k]\in X$. Noting that
$\ex{\dist(x_{k+1},X)}\le \ex{\|x_{k+1}- \Pi_x[x_k]\|^2}$
and $\|x_k - \Pi_X[x_k]\|=\dist(x_k,X)$ for all $k$, we obtain for all $k\ge1$,

\begin{align*}
&\le s_{k-1}^{-2}\ex{\dist^2(x_k,X)} +2s_k^{-1}\ex{f( \Pi_X[x_k]) - f(x_k)}\cr
&\quad- \frac{2s_k^{-1}\g_k}{m\b} \ex{\dist(x_k,X)}
\end{align*}

By the convexity of $f(\cdot)$ and
the definition of the subgradient bound $M$ in (<ref>),
we have that

\[f(\Pi_X[x_k]) - f(x_k)\le M\dist(x_k,X)\]

(see Lemma <ref>(a) where the term $\dist^2(x_k,X)$ is omitted). Therefore, for all $k\ge1$,

\begin{align*}
\le s_{k-1}^{-2}\ex{\dist^2(x_k,X)}
&- 2s_k^{-1} \left(\frac{\g_k}{m\b}-M\right)\ex{\dist(x_k,X)}\cr
\end{align*}

As $\g_k\to+\infty$, there is a large enough $k_1$ so that
$\g_k/(m\b)-M\ge \g_k/(2m\b)$, implying that for all $k\ge k_1$,

\begin{align*}
\le s_{k-1}^{-2}\ex{\dist^2(x_k,X)} -&\frac{s_k^{-1}\g_k}{m\b}\ex{\dist(x_k,X)}\cr
\end{align*}

Now, we divide the preceding relation with $\g_k$. By using the assumption that
$\g_{k+1}\ge \g_k$, we obtain for all $k\ge k_1$,

\begin{align*}
\g_{k+1}^{-1}s_k^{-2} \ex{\dist^2(x_{k+1},X)} \le\g_k^{-1}s_{k-1}^{-2}\ex{\dist^2(x_k,X)}
\end{align*}

By summing these relations over $k=k_1,\ldots,t,$
after re-arranging the terms, we can see that for all $t\ge k_1$,

\begin{align*}
\sum_{k=k_1}^t \frac{s_k^{-1}}{m\b}\ex{\dist(x_k,X)} &\le\g_{k_1}^{-1}s_{k_1-1}^{-2}\ex{\dist^2(x_{k_1},X)}
\end{align*}

The distance function $\dist(\cdot,X)$ is convex since $X$ is convex. Hence,
upon dividing the preceding relation with $\sum_{k=k_1}^t s_k^{-1}$, and by using the convexity of $\dist(\cdot,X)$
and the definition of the truncated $s^{-1}$-weighted averages (see (<ref>), and (<ref>) with $\nu_k=s_k^{-1}$), we obtain
for all $t\ge k_1$,

\begin{align}\label{eq-88ag}
\frac{1}{m\b}\ex{\dist(\bxtktav,X)}
\le\frac{\g_{k_1}^{-1}s_{k_1-1}^{-2}\ex{\dist^2(x_{k_1},X)}}{\bar{S}_t-\bar{S}_{k_1}}
+\frac{\sum_{k=k_1}^t \g_k^{-1}c_k}
\end{align}

where we used the notation $\bar{S}_t=\sum_{k=1}^t s_k^{-1}$ (see (<ref>)).

By expressing $\bxtav$ as a convex combination of $\bar{x}_{k_1}^{\textrm av}$ and $\bxtktav$, i.e., by using (<ref>) with $\nu_k=s_k^{-1}$ and $\tau=k_1$, we have

\[\bxtav=\frac{\bar{S}_{k_1} \bar{x}_{{k_1}}^{\textrm av}
+ (\bar{S}_t-\bar{S}_{k_1})\bar{x}_{{k_1},t}^{\textrm av} }{\bar{S}_t}\qquad\hbox{for all }t\ge k_1\ge 1.\]

Then, by using the convexity of the distance function $\dist(\cdot,X)$,
we find that for all $t\ge k_1$,

\begin{align*}
\ex{\dist(\bxtav,X)}\le
\frac{\bar{S}_{k_1}}{\bar{S}_t}\ex{\dist(\bar{x}_{k_1}^{\textrm av},X)}
+\frac{\bar{S}_t - \bar{S}_{k_1}}{\bar{S}_t}\ex{\dist(\bxtktav,X)}.\end{align*}

The preceding relation and (<ref>) yield for all $t\ge k_1$,

\begin{align*}
\ex{\dist(\bxtav,X)}\le
\frac{\bar{S}_{k_1}}{\bar{S}_t}\ex{\dist(\bar{x}_{k_1}^{\textrm av},X)}
&+\frac{m\b\sum_{k=k_1}^t \g_k^{-1}c_k}{\bar{S}_t}.\end{align*}

The next proposition provides the convergence rate of the procedure (<ref>) in terms of the $s^{-1}$-
weighted averages $\bxtav$ for a strongly convex function $f$.

Let $f(\cdot)$ be strongly convex with $\mu>0$, and
let Assumption <ref> hold.
Assume $s_k=\frac{2}{\mu k}$, $\d_k = \frac{1}{k^d}$ with $d>1$, and
$\g_{k} = \ln^g k$ with $g>0$ for all $k$.
Let $\{x_k\}$ be the iterate sequence  generated by method (<ref>).
Then, the following estimate is valid for the $s^{-1}$-weighted averages $\bxtav$, $t\to\infty$,

\[|\ex{f(\bxtav)}-f^*| = O\left(\frac{\ln^{2g} t}{t}\right). \]

Given the choice $s_k =\frac{2}{\mu k}$,

\begin{align}\label{eq:barS_t}
\bar{S}_t = \sum_{k=1}^{t}\frac{\mu k}{2} = \frac{\mu}{2} \frac{(1+t)t}{2}.
\end{align}

Next, taking into account that $\g_k = \ln^g k$ and $\d_k = \frac{1}{k^d}$, $d>1$, we obtain

\begin{align}\label{eq:c_t}
\sum_{k=1}^t s^{-1}_k\g_k\d_k &=\sum_{k=1}^t \frac{\mu\ln^g k}{2k^{d-1}} \cr
&\le \frac{\mu\ln^g t}{2} [\int_{1}^t\frac{dx}{x^{d-1}} + 1] =\begin{cases}
\frac{\mu\ln^g t}{2}\left(\frac{t^{2-d}-1}{2-d}+1\right), \, \mbox{if $d\ne 2$},\\
\frac{\mu\ln^g t}{2}\left(\ln t+1\right), \, \mbox{if $d= 2$}.
\end{cases}
\end{align}

By combining the inequalities in (<ref>) and (<ref>) with the estimate for the function values
from Lemma <ref>, namely,

\[\ex{f(\bxtav)}-f^*\le \frac{s_0^{-2}\ex{\|x_1 - x^*\|^2}}{2\bar{S}_t}
+\frac{\sum_{k=1}^t c_k}{2\bar{S}_t},\]

where $s_0=2\mu^{-1}$,
$c_k=\frac{\g_k \d_ks_k^{-1}}{2\a_{\min}}+2(M^2+\g_k^2)$, we conclude that there exists the constant $D_1>0$ such that

\begin{align}\label{eq:upper1}
\ex{f(\bxtav)}-f^*\le \frac{D_1 \ln^{2g} t}{t}.
\end{align}

Note that $D_1$ depends on the squared upper bound $M^2$ of the subgradients defined in relation (<ref>), the value $\a_{\min}$, and the expected squared distance between the initial point and the corresponding solution $\ex{\|x_1 - x^*\|^2}$.

To determine a lower bound for $\ex{f(\bxtav)}-f^*$, we use the estimate for the expected distance from Lemma <ref> and $\ex{f(\bxtav)}-f^*
\ge -\|\tilde \nabla f(x^*)\|\,\ex{\dist(\bxtav,X)}$ (see (<ref>)), taking into account that the solution set $X^*$ is singleton, i.e., $X^*=\{x^*\}$. Without loss of generality, we may assume that the subgradient $\tilde\nabla f(x^*)$  satisfies
$\|\tilde\nabla f(x^*)\|\le M$, where  $M$ is defined in (<ref>). In this way, we can conclude that there exists finite $k_1$ such that for all $t\ge k_1$,

\begin{align*}\ex{f(\bxtav)}-f^*
&\ge -M\, \an{\ex{\dist(\bxtav,X)}}\cr
&\ge -M\Big(\frac{\bar{S}_{k_1}}{\bar{S}_t}\ex{\dist(\bar{x}_{k_1}^{\textrm av},X)}
&\quad+\frac{m\b\sum_{k=k_1}^t \g_k^{-1}c_k}{\bar{S}_t}\Big).
\end{align*}

By using the estimation analogous to one in (<ref>) and (<ref>), we conclude existence of some constant $D_2>0$ such that

\begin{align}\label{eq:lower1}
\ex{f(\bxtav)}-f^*\ge -\frac{D_2\ln^g t}{t}.
\end{align}

The constant $D_2$ in its turn depends on the index $k_1$, the upper bound $M$ of the subgradient norms, number of the constraints $m$, $\a_{\min}$, the Hoffman constant $\beta$, and distance between $x_{k_1}$ and the feasible set $X$.
Hence, combining (<ref>) and (<ref>), we conclude the result.

The constant in the $O$-notation can be obtained using Lemma <ref>. In particular, it can be seen that, it depends linearly on the  strong-convexity constant $\mu$ and its inverse counterpart $\mu^{-1}$, the squared upper bound of the gradients over a region containing the iterates and their projections on the constraint set $X$ (see (<ref>)),  the
Hoffman constant $\beta$, the number of the constraints $m$,
the inverse of the smallest norm $\a^{-1}_{\min}$ of the vectors $a_i,i\in[m]$, and the expected distance $\ex{\|x_1-x^*\|^2}$ between the initial point $x_1$ and the optimal solution $x^*$. Given the choice of the parameters, we conclude that the rate achieved by the proposed procedure in the case of strongly convex $f$ is $O\left(\frac{\ln^{\e} t}{t}\right)$ for any small positive $\e$.

As seen in Propositions <ref> and <ref>, we obtain the convergence rate that is within a logarithmic factor matching the best convergence rates of $O(1/\sqrt{t})$ and $O(1/t)$ known for a stochastic subgradient method applied to a convex and strongly convex function respectively. The extra logarithmic factor is incurred due to the fact that the method (<ref>) uses stochastic subgradients of time-varying functions $F_k(\cdot)$ yielding the sequence of infeasible iterates.

As seen from 	the proofs of these results, the extra
logarithmic factor comes from the penalty parameters,
which control the infeasibility of the iterates. However, it is worth noting that the factors $\ln^{1/2+\e} t$ and $\ln^{\e} t$ for any small positive $\e$ in the merely convex and strongly convex cases respectively improve the previously obtained results on convergence rates of incremental procedures with many constraints [8].

§ SIMULATION RESULTS

(6,18)90Relative Error

(40,-1)time, $k$

Strongly convex function $f$ with the optimum inside the feasible set $X$.

(6,18)90Relative Error

(40,-1)time, $k$

Strongly convex function $f$ with the optimum outside the feasible set $X$.

(6,18)90Relative Error

(40,-1)time, $k$

Merely convex function $f$ with the optimum inside the feasible set $X$.

(6,18)90Relative Error

(40,-1)time, $k$

Merely convex function $f$ with the optimum outside the feasible set $X$.

(6,18)90Relative Error

(40,-1)time, $k$

Full gradient method for the case of optimum outside the feasible set $X$.

In this section we present some simulation results obtained for the proposed procedure (<ref>).
To implement the algorithm for a strongly convex problem, we set up the objective function as follows:

\[f(x) = \|x-x^0\|^2,\]

where $f:\R^{10}\to\R$ and $x^0\in\R^{10}$ is some randomly generated vector. The vectors $a_i\in\R^{10}$ and the scalars $b_i$, $i\in\{1,\ldots,m\}$, defining the constraint set $X$, are generated according to a normal distribution. We distinguish between two settings: 1) $x^0\in\R^{10}$, as well as $a_i\in\R^{10}$ and $b_i\in\R$, are generated in such a way that $x^0\in X$; 2) $x^0\in\R^{10}$, as well as $a_i\in\R^{10}$ and $b_i\in\R$, are generated in such a way that $x^0\notin X$. The simulation results in these cases, given $m=100, 500, 1000$, are demonstrated by Figures <ref> and <ref> respectively. In these simulations, we estimate the relative error in terms of the actual iterates, namely $\frac{\|x_k - x^*\|}{\|x^*\|}$. The parameters $s_t$, $\g_t$, and $\d_t$ are chosen according to the conditions in Proposition <ref>, namely $s_t = \frac{1}{k^{0.99}}$, $\g_t = c\ln k$, $s_t = \frac{1}{k^{2}}$, where the parameter $c$ is tuned.

For the case of non-strongly convex optimization, we choose the following objective function:

\[f(x) = \|x-x^0\|_{1},\]

where $\|\cdot\|_1$ denotes the $l_1$-norm. Analogously to the strongly convex case, $f:\R^{10}\to\R$ and $x^0\in\R^{10}$ is some randomly generated vector. The vectors $a_i\in\R^{10}$ and the scalars $b_i$, $i\in\{1,\ldots,m\}$, defining the constraint set $X$, are generated according to a normal distribution. As before, we distinguish between two settings: 1) $x^0\in\R^{10}$, as well as $a_i\in\R^{10}$ and $b_i\in\R$, are generated in such a way that $x^0\in X$; 2) $x^0\in\R^{10}$, as well as $a_i\in\R^{10}$ and $b_i\in\R$, are generated in such a way that $x^0\notin X$. The simulation results in these cases, given $m=100, 500, 1000$, are demonstrated by Figures <ref> and <ref> respectively. Here, the parameters $s_t$, $\g_t$, and $\d_t$ are chosen according to the conditions in Proposition <ref>, namely $s_t = \frac{S}{k^{0.5}}$, $\g_t = c\ln k$, $s_t = \frac{1}{k^{2}}$, where the parameters $S$ and $c$ are tuned.

For the completeness of simulations, we implemented the non-incremental version of the procedure (<ref>), where at each iteration the full gradient of the penalized objective function is used. The corresponding results are presented by Figure <ref>.

§ CONCLUSION

This work deals with penalty reformulation of the optimization problems subjected to linear constraints, where the penalty is a variant of the one-sided Huber loss function. The infeasibility properties of the solutions of penalized problems for nonconvex and convex objective functions under time-varying penalty parameters are analyzed.
A random incremental penalty method for solving convex problem is proposed. This method is proven to converge to a solution of the original problem almost surely and in expectation for suitable choices of the penalty parameters and the stepsize. Moreover,  $O(\ln^{1/2+\e} k/{\sqrt k})$-convergence rate when the objective function is convex, and  $O(\ln^{\e} k/k)$-convergence rate when the objective function is strongly convex are established, given any positive $\e$.

Some interesting questions for the future work include applicability of
accelerated incremental algorithms for the proposed penalty
reformulation in the case of both strongly and non-strongly
convex optimization as well as extension of the presented penalty approach to other structures of constraints.

[1]
J. W. Adams.
$\mbox{FIR}$ digital filters with least-squares stopbands subject to
peak-gain constraints.
IEEE Transactions on Circuits and Systems, 38(4):376–388,

[2]
D. P. Bertsekas.
Constrained Optimization and Lagrange Multiplier Methods
(Optimization and Neural Computation Series).
Athena Scientific, 1 edition, 1996.

[3]
D. P. Bertsekas.
A hybrid incremental gradient method for least squares.
SIAM Journal on Optimization, 7:913–926, 1997.

[4]
D. P. Bertsekas.
Incremental proximal methods for large scale convex optimization.
Mathematical Programming, 129(2):163–195, 2011.

[5]
D. P. Bertsekas.
Incremental gradient, subgradient, and proximal methods for convex
optimization: A survey.
available on arxiv at https://arxiv.org/abs/1507.01030, 2015.

[6]
P. Billingsley.
Probability and Measure.
John Wiley & Sons, Inc., New York, NY, USA, 3rd edition, 1995.

[7]
G. Dorini, P. Pinson, and H. Madsen.
Chance-constrained optimization of demand response to price signals.
IEEE Transactions on Smart Grid, 4(4):2072–2080, 2013.

[8]
O. Fercoq, A. Alacaoglu, I. Necoara, and V. Cevher.
Almost surely constrained convex optimization.
In ICML, 2019.

[9]
M. Gaudioso, G. Giallombardo, and G. Miglionico.
An incremental method for solving convex finite min-max problems.
Mathematics of Operations Research, 31:173–187, 2006.

[10]
L. Grippo.
A class of unconstrained minimization methods for neural network
Optimization Methods and Software, 4:135–150, 1994.

[11]
L. Grippo.
Convergent on-line algorithms for supervised learning in neural
IEEE Transactions on Neural Networks, 11:1284–1299, 2000.

[12]
E. S. Helou and A. R. De Pierro.
Incremental subgradients for constrained convex optimization, a
unified framework and new methods.
SIAM Journal on Optimization, 20:1547–1572, 2009.

[13]
A. J. Hoffman.
On approximate solutions of systems of linear inequalities.
Journal of Research of the National Bureau of Standards,
49(4):263–265, 1952.

[14]
B. Johansson, M. Rabi, and M. Johansson.
A randomized incremental subgradient method for distributed
optimization in networked systems.
SIAM Journal on Optimization, 20:1157–1170, 2009.

[15]
V. M. Kibardin.
Decomposition into functions in the minimization problem.
Automation and Remote Control, 40:1311–1323, 1980.

[16]
A. Kundu, F. Bach, and C. Bhattacharyya.
Convex optimization over intersection of simple sets: improved
convergence rate guarantees via an exact penalty approach.
In Proceedings of the 21st International Conference on
Artificial Intelligence and Statistics, volume 84 of Proceedings of
Machine Learning Research, pages 958–967. PMLR, 2018.

[17]
W. Li and J. Swetits.
The linear $\ell$-1 estimator and the Huber M-estimator.
SIAM Journal on Optimization, 8(2):457–475, 1998.

[18]
Z. Q. Luo.
On the convergence of the lms algorithm with adaptive learning rate
for linear feedforward networks.
Neural Computation, 3:226–245, 1991.

[19]
C. Mathieu and W. Schudy.
Correlation clustering with noisy input.
In SODA, pages 712–728. SIAM, 2010.

[20]
I. Necoara, P. Richtárik, and A. Patrascu.
Randomized projection methods for convex feasibility: Conditioning
and convergence rates.
SIAM Journal on Optimization, 29(4):2814–2852, 2019.

[21]
A. Nedić.
Random projection algorithms for convex set intersection problems.
In Proceedings of the 49th IEEE Conference on Decision and
Control, pages 7655–7660, 2010.

[22]
A. Nedić.
Random algorithms for convex minimization problems.
Mathematical Programming, 129(2):225–253, 2011.

[23]
A. Nedić and D. P. Bertsekas.
Incremental subgradient methods for nondifferentiable optimization.
SIAM Journal on Optimization, 12:109–138, 2001.

[24]
A. Nedić and S. Lee.
On stochastic subgradient mirror-descent algorithm with weighted
SIAM Journal on Optimization, 24(1):84–107, 2014.

[25]
A. Nedić and T. Tatarenko.
Convergence rate of a penalty method for strongly convex problems
with linear constraints.
In Proceedings of the 59th IEEE Conference on Decision and
Control, Jeju Island, Republic of Korea, Dec. 14–18, 2020, pages 372–377,

[26]
A. Patrascu and I. Necoara.
Nonasymptotic convergence of stochastic proximal point methods for
constrained convex optimization.
J. Mach. Learn. Res., 18(1):7204–7245, jan 2017.

[27]
J. Pena, J. Vera, and L. F. Zuluaga.
New characterizations of hoffman constants for systems of linear
Mathematical Programming, 187:79–109, 2021.

[28]
B. T. Polyak.
Introduction to optimization.
Optimization Software, Inc., New York, 1987.

[29]
H. Robbins and D. Siegmund.
A convergence theorem for nonnegative almost supermartingales and
some applications.
In Optimizing Methods in Statistics, pages 233–257. Academic
Press, New York, 1971.

[30]
W. Siedlecki and J. Sklansky.
Constrained genetic optimization via dynamic reward-penalty balancing
and its use in pattern recognition.
In Proceedings of the Third International Conference on Genetic
Algorithms, pages 141–150, San Francisco, CA, USA, 1989. Morgan Kaufmann
Publishers Inc.

[31]
M. V. Solodov.
Incremental gradient algorithms with stepsizes bounded away from
Comput. Opt. Appl., 11:28–35, 1998.

[32]
C. Song, Cheuk Yin Lin, Stephen J. Wright, and Jelena Diakonikolas.
Coordinate linear variance reduction for generalized linear
programming, 2021.

[33]
T. Strohmer and R. Vershynin.
A randomized kaczmarz algorithm with exponential convergence.
Journal of Fourier Analysis and Applications, 15(2), 2008.
article number 262.

[34]
T. Tatarenko and A. Nedić.
A smooth inexact penalty reformulation of convex problems with linear
SIAM Journal on Optimization, 31(3):2141–2170, 2021.

[35]
P. Tseng.
An incremental gradient(-projection) method with momentum term and
adaptive stepsize rule.
SIAM Journal on Optimization, 8:506–531, 1998.

[36]
M. Wang and D. P. Bertsekas.
Incremental constraint projection methods for variational
Mathematical Programming, 150(2):321–363, 2015.

[37]
S. J. Wright, R. D. Nowak, and M. A. T. Figueiredo.
Sparse reconstruction by separable approximation.
In Proceedings of the IEEE International Conference on
Acoustics, Speech and Signal Processing (ICASSP), pages 3373–3376, 2008.

[38]
Y. Xu.
Primal-dual stochastic gradient method for convex programs with many
functional constraints.
SIAM J. Optim., 30:1664–1692, 2020.

[39]
M. Zaslavskiy, F. Bach, and J. P. Vert.
A path following algorithm for the graph matching problem.
IEEE Transactions on Pattern Analysis and Machine Intelligence,
31(12):2227–2242, 2009.

§ PROOF OF PROPOSITION <REF>

Let $\hat x\in X$ be an arbitrary feasible point, and $k\ge1$ be arbitrary but fixed.
By the optimality of $x_k^*$, we have

\[F_k(x_k^*)\le F_k(\hat x).\]

By using $F_k(x)=f(x)+\g_k H_k(x)$ (see (<ref>)),
after re-arranging the terms, we obtain

\[\g_kH_k(x_k^*)\le f(\hat x) - f(x_k^*)+\g_k H_k(\hat x).\]

By relation (<ref>) with $x=\hat x\in X$ and
relation (<ref>) with $x=x^*_k$, we obtain

\[\frac{\g_k}{m\b}\dist(x_k^*,X) \le f(\hat x) - f(x_k^*)+\frac{\g_k\d_k}{4\a_{\min}}.\]

By multiplying the preceding relation with $m\b/\g_k$, we have

\begin{equation}\label{eq-disto}
\dist(x_k^*,X) \le \frac{m\b}{\g_k} (f(
\hat x) - f(x_k^*))+\frac{m\b\d_k}{4\a_{\min}},\end{equation}

which is the stated relation.
By Lemma <ref>, the optimal sets $X_k^*$ are uniformly bounded, i.e., there exist some $D>0$, such that
$\|x_k^*\|\le D$ for all $x_k^*\in X_k^*$ and for all $k\ge1$. Hence, the values $f(x_k^*)$ are also uniformly bounded
implying that
$\frac{m\b}{\g_k} (f(\hat x) - f(x_k^*))\to0$ and $\g_k\to\infty$. The rest of the statement follows immediately from relation (<ref>).

§ PROOF OF LEMMA <REF>

The set $X_k^*$ is closed by the continuity of $F_k(\cdot)$ for all $k$. It remains to show that $X_k^*$ is bounded and nonempty for all $k$. Let $k\ge 1$ be arbitrary.
By Corollary <ref>, with $\d=\d_k$,
$\g=\g_k$, and $\hat x\in X$, we have

\begin{align}\label{eq-bounded}
\{x\in\R^n\mid F_k(x)\le t_{k}(\hat x)\}\subseteq
\{x\in\R^n\mid f(x)\le t_{k}(\hat x)\},\end{align}

with $t_k(\hat x)= f(\hat x)+\g_k\d_k/(4\a_{\min}).$
Since $\g_k\d_k\le c$, it follows that

\[t_{k}(\hat x)\le t_c(\hat x)\triangleq
f(\hat x)+\frac{c}{4\a_{\min}},\]

implying that

\[\{x\in\R^n\mid F_k(x)\le t_{k}(\hat x)\}\subseteq
\{x\in\R^n\mid f(x)\le  t_c(\hat x)\}.\]

The set $\{x\in\R^n\mid f(x)\le t_{c}(\hat x)\}$
is bounded since $f$ has bounded lower-level sets, implying by the preceding relation that the set $\{x\in\R^n\mid F_k(x)\le t_{k}(\hat x)\}$ is bounded. Since
$X_k^*\subseteq \{x\in\R^n\mid F_k(x)\le t_{k}(\hat x)\}$, the set $X_k^*$ is bounded.

To show that $X_k^*\ne\emptyset$, we observe that the penalized problem
$\min_{x\in \R^n} F_k(x)$ and its related constrained problem

\[\min_{F_k(x)\le t_{k}(\hat x)} F_k(x)\]

have the same optimal values and the same optimal sets.
The optimal set of the problem above is nonempty by the Weierstrass Theorem, implying that $X_k^*$ is nonempty.

To show that $\{X_k^*\}$ is uniformly bounded, we use the fact that
$X_k^*\subseteq\{x\in\R^n\mid F_k(x)\le t_{k}(\hat x)\}$ for all $k$, which by (<ref>) implies that
$X_k^*\subseteq\{x\in\R^n\mid f(x)\le t_{c}(\hat x)\}$ for all $k$.
The boundedness of the set $\{x\in\R^n\mid f(x)\le t_{c}(\hat x)\}$ implies that $\{X_k^*\}$ is uniformly bounded.

§ PROOF OF PROPOSITION <REF>

In case $\mu=0$, the conditions of Lemma <ref> are satisfied.
If $f(\cdot)$ is strongly convex with $\mu>0$, then it has bounded lower-level sets.

In this case, the conditions of Lemma <ref> are also satisfied.
Thus, by Lemma <ref> the optimal solutions $x_k^*$ of the penalized problems are uniformly bounded, i.e., there exists $D>0$ such that
$\|x_k^*\|\le D$ for any optimal solution $x_k^*\in X_k^*$ of the penalized problem $\min_{x\in\R^n} F_k(x)$
and for all $k$. Since
$\{x_k^*\mid x_k^*\in X_k^*, k\ge 1\}$
is bounded, the set $\{\Pi_X[x_k^*]\mid x_k\in X_k^*,\ k\ge1\}$ is also bounded, i.e.,
there is $R>0$ such that $\Pi_X[x_k^*]\le R$ for all $x_k^*\in X_k^*$ and all $k$.
Thus, the constant $L$ in (<ref>) is finite.

Since each $F_k(\cdot)$ is strongly convex with a constant $\mu\ge0$,
by the optimality of $x_k^*$ we have

\[\frac{\mu}{2}\|x^* -x_k^*\|^2\le F_k(x_k^*)- F_k(x)\quad\hbox{for any $x\in\R^n$}.\]

Using $x=x^*$, where $x^*\in X^*$  and the expression for $F_k$ in (<ref>), we obtain for all $k$,

\begin{align}\label{eq-11}
\frac{\mu}{2}&\|x^* -x_k^*\|^2\le f(x^*)-f(x_k^*) + \g_kH_k(x^*) - \g_kH_k(x_k^*).
\end{align}

Adding and subtracting $f(\Pi_X[x_k^*])$ yields

\begin{align}\label{eq-12}
f(x^*)-f(x_k^*)&=f(x^*)-f(\Pi_X[x_k^*]) + f(\Pi_X[x_k^*]) - f(x_k^*)\cr
&\le -\frac{\mu}{2}\|x^*-  \Pi_X[x_k^*]\|^2
\end{align}

where the inequality is obtained using
$\frac{\mu}{2}\|\Pi_X[x_k^*]-x^*\|^2 + \la \tilde \nabla f(x^*), \Pi_X[x_k^*]-x^*\ra +f(x^*)\le f(p_k)$
for a subgradient $\tilde \nabla f(x^*)\in\partial f(x^*)$ (see (<ref>)),
and the fact that $\la \tilde \nabla f(x^*), \Pi_X[x_k^*]-x^*\ra\ge0$,
which holds since $\Pi_X[x_k^*]$ is feasible
and $x^*$ is the optimal point of $f(\cdot)$ over $X$.
By the (strong) convexity relation (<ref>), we also have

\begin{align*}
f(\Pi_X[x_k^*])-f(x_k^*)&\le -\frac{\mu}{2}\|\Pi_X[x_k^*]-x_k^*\|^2 \cr
&\ + \la \tilde\nabla f(\Pi_X[x_k^*]), \Pi_X[x_k^*] - x_k^*\ra.\end{align*}

Since the subgradients of $f(\cdot)$ at points $\Pi_X[x^*_k]$ are uniformly bounded (see (<ref>)),
it follows that

\[f(\Pi_X[x_k^*])-f(x_k^*)\le -\frac{\mu}{2}\|\Pi_X[x_k^*]-x_k^*\|^2 + L\|\Pi_X[x_k^*] - x_k^*\|.\]

Combining the preceding inequality with (<ref>) and (<ref>), after
re-arranging nonpositive terms, we obtain

\begin{align*}
\frac{\mu}{2}&\|x^* -x_k^*\|^2+\frac{\mu}{2}\|x^*-\Pi_X[x_k^*]\|^2 +\frac{\mu}{2}\|\Pi_X[x_k^*] - x_k^*\|^2\cr
&\le L \|\Pi_X[x_k^*] - x_k^*\|+ \g_k H_k(x^*)-\g_k H_k(x_k^*).
\end{align*}

Since $x^*$ is feasible, we apply relation (<ref>) with $x=x^*$.
Also, we use (<ref>) with $x=x_k^*$, and obtain

\begin{align*}
\frac{\mu}{2}&\|x^* -x_k^*\|^2+\frac{\mu}{2}\|x^*-\Pi_X[x_k^*]\|^2 +\frac{\mu}{2}\|\Pi_X[x_k^*] - x_k^*\|^2\cr
&\le L \|\Pi_X[x_k^*] - x_k^*\| + \frac{\g_k\d_k}{4\a_{\min}} -\frac{\g_k}{m\b} \dist(x_k^*,X).
\end{align*}

$\|\Pi_X[x_k^*] -x_k^*\|=\dist(x_k^*,X)$ and grouping the terms accordingly, we arrive at the stated relation.

[1]
J. W. Adams.
$\mbox{FIR}$ digital filters with least-squares stopbands subject to
peak-gain constraints.
IEEE Transactions on Circuits and Systems, 38(4):376–388,

[2]
D. P. Bertsekas.
Constrained Optimization and Lagrange Multiplier Methods
(Optimization and Neural Computation Series).
Athena Scientific, 1 edition, 1996.

[3]
D. P. Bertsekas.
A hybrid incremental gradient method for least squares.
SIAM Journal on Optimization, 7:913–926, 1997.

[4]
D. P. Bertsekas.
Incremental proximal methods for large scale convex optimization.
Mathematical Programming, 129(2):163–195, 2011.

[5]
D. P. Bertsekas.
Incremental gradient, subgradient, and proximal methods for convex
optimization: A survey.
available on arxiv at https://arxiv.org/abs/1507.01030, 2015.

[6]
P. Billingsley.
Probability and Measure.
John Wiley & Sons, Inc., New York, NY, USA, 3rd edition, 1995.

[7]
G. Dorini, P. Pinson, and H. Madsen.
Chance-constrained optimization of demand response to price signals.
IEEE Transactions on Smart Grid, 4(4):2072–2080, 2013.

[8]
O. Fercoq, A. Alacaoglu, I. Necoara, and V. Cevher.
Almost surely constrained convex optimization.
In ICML, 2019.

[9]
M. Gaudioso, G. Giallombardo, and G. Miglionico.
An incremental method for solving convex finite min-max problems.
Mathematics of Operations Research, 31:173–187, 2006.

[10]
L. Grippo.
A class of unconstrained minimization methods for neural network
Optimization Methods and Software, 4:135–150, 1994.

[11]
L. Grippo.
Convergent on-line algorithms for supervised learning in neural
IEEE Transactions on Neural Networks, 11:1284–1299, 2000.

[12]
E. S. Helou and A. R. De Pierro.
Incremental subgradients for constrained convex optimization, a
unified framework and new methods.
SIAM Journal on Optimization, 20:1547–1572, 2009.

[13]
A. J. Hoffman.
On approximate solutions of systems of linear inequalities.
Journal of Research of the National Bureau of Standards,
49(4):263–265, 1952.

[14]
B. Johansson, M. Rabi, and M. Johansson.
A randomized incremental subgradient method for distributed
optimization in networked systems.
SIAM Journal on Optimization, 20:1157–1170, 2009.

[15]
V. M. Kibardin.
Decomposition into functions in the minimization problem.
Automation and Remote Control, 40:1311–1323, 1980.

[16]
A. Kundu, F. Bach, and C. Bhattacharyya.
Convex optimization over intersection of simple sets: improved
convergence rate guarantees via an exact penalty approach.
In Proceedings of the 21st International Conference on
Artificial Intelligence and Statistics, volume 84 of Proceedings of
Machine Learning Research, pages 958–967. PMLR, 2018.

[17]
W. Li and J. Swetits.
The linear $\ell$-1 estimator and the Huber M-estimator.
SIAM Journal on Optimization, 8(2):457–475, 1998.

[18]
Z. Q. Luo.
On the convergence of the lms algorithm with adaptive learning rate
for linear feedforward networks.
Neural Computation, 3:226–245, 1991.

[19]
C. Mathieu and W. Schudy.
Correlation clustering with noisy input.
In SODA, pages 712–728. SIAM, 2010.

[20]
I. Necoara, P. Richtárik, and A. Patrascu.
Randomized projection methods for convex feasibility: Conditioning
and convergence rates.
SIAM Journal on Optimization, 29(4):2814–2852, 2019.

[21]
A. Nedić.
Random projection algorithms for convex set intersection problems.
In Proceedings of the 49th IEEE Conference on Decision and
Control, pages 7655–7660, 2010.

[22]
A. Nedić.
Random algorithms for convex minimization problems.
Mathematical Programming, 129(2):225–253, 2011.

[23]
A. Nedić and D. P. Bertsekas.
Incremental subgradient methods for nondifferentiable optimization.
SIAM Journal on Optimization, 12:109–138, 2001.

[24]
A. Nedić and S. Lee.
On stochastic subgradient mirror-descent algorithm with weighted
SIAM Journal on Optimization, 24(1):84–107, 2014.

[25]
A. Nedić and T. Tatarenko.
Convergence rate of a penalty method for strongly convex problems
with linear constraints.
In Proceedings of the 59th IEEE Conference on Decision and
Control, Jeju Island, Republic of Korea, Dec. 14–18, 2020, pages 372–377,

[26]
A. Patrascu and I. Necoara.
Nonasymptotic convergence of stochastic proximal point methods for
constrained convex optimization.
J. Mach. Learn. Res., 18(1):7204–7245, jan 2017.

[27]
J. Pena, J. Vera, and L. F. Zuluaga.
New characterizations of hoffman constants for systems of linear
Mathematical Programming, 187:79–109, 2021.

[28]
B. T. Polyak.
Introduction to optimization.
Optimization Software, Inc., New York, 1987.

[29]
H. Robbins and D. Siegmund.
A convergence theorem for nonnegative almost supermartingales and
some applications.
In Optimizing Methods in Statistics, pages 233–257. Academic
Press, New York, 1971.

[30]
W. Siedlecki and J. Sklansky.
Constrained genetic optimization via dynamic reward-penalty balancing
and its use in pattern recognition.
In Proceedings of the Third International Conference on Genetic
Algorithms, pages 141–150, San Francisco, CA, USA, 1989. Morgan Kaufmann
Publishers Inc.

[31]
M. V. Solodov.
Incremental gradient algorithms with stepsizes bounded away from
Comput. Opt. Appl., 11:28–35, 1998.

[32]
C. Song, Cheuk Yin Lin, Stephen J. Wright, and Jelena Diakonikolas.
Coordinate linear variance reduction for generalized linear
programming, 2021.

[33]
T. Strohmer and R. Vershynin.
A randomized kaczmarz algorithm with exponential convergence.
Journal of Fourier Analysis and Applications, 15(2), 2008.
article number 262.

[34]
T. Tatarenko and A. Nedić.
A smooth inexact penalty reformulation of convex problems with linear
SIAM Journal on Optimization, 31(3):2141–2170, 2021.

[35]
P. Tseng.
An incremental gradient(-projection) method with momentum term and
adaptive stepsize rule.
SIAM Journal on Optimization, 8:506–531, 1998.

[36]
M. Wang and D. P. Bertsekas.
Incremental constraint projection methods for variational
Mathematical Programming, 150(2):321–363, 2015.

[37]
S. J. Wright, R. D. Nowak, and M. A. T. Figueiredo.
Sparse reconstruction by separable approximation.
In Proceedings of the IEEE International Conference on
Acoustics, Speech and Signal Processing (ICASSP), pages 3373–3376, 2008.

[38]
Y. Xu.
Primal-dual stochastic gradient method for convex programs with many
functional constraints.
SIAM J. Optim., 30:1664–1692, 2020.

[39]
M. Zaslavskiy, F. Bach, and J. P. Vert.
A path following algorithm for the graph matching problem.
IEEE Transactions on Pattern Analysis and Machine Intelligence,
31(12):2227–2242, 2009.