# An inexact restoration-nonsmooth algorithm with variable accuracy for
stochastic nonsmooth convex optimization problems in machine learning and
stochastic linear complementarity problems

Nataša Krejić111Department of Mathematics and Informatics, Faculty of
Sciences, University of Novi Sad, Trg Dositeja Obradovića 4, 21000 Novi Sad,
Serbia. e-mail<EMAIL_ADDRESS>Nataša Krklec Jerinkić222Department of
Mathematics and Informatics, Faculty of Sciences, University of Novi Sad, Trg
Dositeja Obradovića 4, 21000 Novi Sad, Serbia. e-mail:
<EMAIL_ADDRESS>Tijana Ostojić 333Department of Fundamental
Sciences, Faculty of Technical Sciences, University of Novi Sad, Trg Dositeja
Obradovića 6, 21000 Novi Sad, Serbia. e-mail<EMAIL_ADDRESS>444Corresponding author

(November 2, 2022)

###### Abstract

We study unconstrained optimization problems with nonsmooth and convex
objective function in the form of a mathematical expectation. The proposed
method approximates the expected objective function with a sample average
function using Inexact Restoration-based adapted sample sizes. The sample size
is chosen in an adaptive manner based on Inexact Restoration. The algorithm
uses line search and assumes descent directions with respect to the current
approximate function. We prove the a.s. convergence under standard
assumptions. Numerical results for two types of problems, machine learning
loss function for training classifiers and stochastic linear complementarity
problems, prove the efficiency of the proposed scheme.

Key words: Nonsmooth optimization, Subgradient, Inexact Restoration, Sample
Average Approximation, Variable sample size.

##  1 Introduction

Let us observe an unconstrained optimization problem which objective function
takes the form of a mathematical expectation

$\min_{x}f(x)=E\left[F(x,\xi)\right],$ (1.1)

where $F:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ is
continuous and convex function with respect to $x$, bounded from below,
$\xi:\Omega\rightarrow\mathbb{R}^{m}$ is random vector and
$(\Omega,\mathcal{F},P)$ is probability space. Convexity implies that $F$ is
locally Lipschitz, [1]. No additional smoothness assumption is imposed. A
number of important problems can be stated in the form (1.1) - starting from
data analytics with huge data sets which require working with subsamples or
online training with the permanently increasing data sets [2], to simulations
of natural and industrial processes with number of random parameters [3, 4, 5,
6].

The objective function in (1.1) can rarely be computed exactly and might be
nonsmooth. Thus, the main issues that arise in iterative methods for solving
(1.1) are the approximation of the objective function and the choice of search
directions. The most common approximation of the mathematical expectation is
the Sample Average Approximation (SAA). For a given independent and
identically distributed, i.i.d., sample $\\{\xi_{1},\ldots,\xi_{N}\\}$ of the
size $N,$ the SAA approximate objective function is defined as

$f_{N}(x)=\frac{1}{N}\sum_{i=1}^{N}f_{i}(x),$ (1.2)

where $f_{i}(x)=F(x,\xi_{i})$. The sample vectors $\xi_{1},\ldots,\xi_{N}$ are
assumed to be i.i.d. and the sample size $N$ determines the accuracy of the
approximation (1.2), [7]. Naturally, larger $N$ implies higher accuracy of the
approximate function $f_{N},$ but makes any optimization algorithm more costly
as the cost of computing $f_{N},$ as well as search directions, increases with
$N.$ There is a vast literature dealing with variable sample size methods for
SAA approximations, [8, 9, 10, 11, 12], which range from simple heuristics to
complex schemes, all of them with the idea of using cheaper, lower accuracy
approximations of the objective function whenever possible, in order to save
the computational effort.

The second issue one needs to address is the choice of search directions. In
the case of smooth problems we can choose between relatively slow but cheap
first order methods or more elaborate and more costly second order methods,
depending on a particular problem structure, needed accuracy etc. In the case
of nonsmooth problems the gradient is generally replaced by a subgradient or
more elaborate schemes like gradient sampling, [13, 14], bundle methods, [15],
proximal methods, [16], and so on. A number of recent papers deals with second
order search directions [17, 18, 19].

The method presented in this paper addresses both issues by using an adaptive
variable accuracy and descent directions with respect to the current
approximate functions. The sample size is governed by Inexact Restoration (IR)
framework introduced by Martinez and Pilota [20] and consists of two phases:
the restoration and the optimality phase. The main idea of IR is to treat the
phases, restoration and optimality, in a modular way and then to use a merit
function, which combines feasibility and optimality and enforces progress
towards a feasible optimal point. As IR is constrained optimization tool, the
problem (1.1) is reformulated into a constrained problem as follows

$\min f_{N}(x),\text{ s.t. }f_{N}(x)=f(x),$ (1.3)

where $f$ and $f_{N}$ are defined in (1.1) and (1.2), respectively.

Notice that (1.3) is equivalent to (1.1) if the constraint is satisfied.
However if we consider methods that are not strictly feasible, i.e., not all
iterations satisfy the constraint, then we can treat $N$ as an additional
variable in the constraint. That is precisely what we will do in the IR
approach - in each iteration $k$ of the method we will determine a suitable
$N_{k}.$ There are numerous studies that have confirmed the benefits of using
the IR approach in the varying accuracy approximations framework, [21, 22].
The key advantage of this approach is the fact that feasibility and optimality
are kept in balance through merit function. Therefore, the accuracy of the
approximate objective function depends on the progress towards optimality in
each iteration. Obviously, the accuracy is adaptive, endogenous to the
algorithm and there is no need for additional parameters or heuristics in the
sample size determination. Furthermore, the sequence of sample sizes is very
often nonmonotone, increasing the accuracy (and the computational cost)
whenever we approach the solution to ensure good quality of the approximate
solution, and decreasing the accuracy (and the costs) when the current iterate
is far away from the solution. The approach has been used for variable
accuracy approximations for the first time in [22] for the problem of finite
sum minimization coupled with line search descent direction method, based on
results from [23]. It is extended to trust region framework and constrained
problems, [20, 21, 24]. An approach for solving problems with variable
accuracy in both objective function and constraints is analyzed in [25].

The step size is a challenging issue in stochastic analysis and it was a
subject of research in many papers, [26, 27, 28, 29, 30, 31]. Line search
methods, which are an important tool in deterministic optimization, are not
easily extended to the stochastic case due to the mutual dependence of step
size and search direction, which are both random variables in the stochastic
framework. An important study on this topic is given in [30] where the
approximations of the objective function and its gradient are assumed to be
good enough with a fixed high probability. Under these settings, the
complexity analysis in terms of expected number of iterations to reach near-
optimal solution is provided. In [32] a second order direction is considered
but an additional sampling is used in Armijo-like condition to overcome the
bias issue. The approach presented here differs in several aspects. First of
all, we consider the approximate objective of the form (1.2) and prove that
the algorithm introduced here yields $N\to\infty.$ In other words we approach
the objective function almost surely under some standard conditions. This
property of the algorithm is a direct consequence of IR strategy. Furthermore,
the conditional expectation of the relevant SAA estimator is equal to the
objective function under our settings (for details see the final paragraph of
Section 2 and the proof of Lemma 3.1), and the step size is not directly
involved. Another important difference lies in the fact that the objective
function and its approximations are not differentiable, and thus the step size
analysis is more complicated even in the strongly convex case.

Our contributions are the following. We define Inexact Restoration - Nonsmooth
(IR-NS) algorithm for nonsmooth optimization with variable accuracy and prove
a.s. convergence of the algorithm under the set of standard assumptions. By
using Inexact Restoration for sample size selection we generalize the results
from [33]. More precisely, since IR-NS pushes the SAA error to zero, in the
case of finite sum problems where the objective function is given by (1.2)
with the finite full sample size $N$, the true objective function is reached
eventually and the convergence results from [33] hold. IR-NS algorithm also
covers wider class of problems than finite sums, including infinite sums. Our
experiments confirm the intuitive reasoning that working with variable,
adaptive sample size is more effective than working with predefined or full
sample size as in [33]. To emphasize this fact we present experiments with the
same search direction as in [33] \- the nonsmooth BFGS descent direction, and
demonstrate the advantages of variable sample size approach proposed in IR-NS.
In general, an arbitrary descent direction in the sense of Assumption 2 stated
below is applicable. From theoretical point of view, the complexity of order
$\varepsilon^{-2}$ is proved, which also applies to the method from [33]. The
obtained complexity is in line with the results from [34] where the complexity
of IR is analyzed. The result in [34] is obtained for smooth constrained
problems and is of the form $\varepsilon_{feas}^{-1}+\varepsilon_{opt}^{-2},$
with $\varepsilon_{feas}$ being the constant for feasibility and
$\varepsilon_{opt}$ coincides with the $\varepsilon$ that we consider here.
Notice that the problems considered in [34] are smooth and deterministic. The
complexity results obtained in [24] are not comparable to the complexity
results for IR-NS as the methods analyzed in [24] are specialized for smooth
problems and problems with regularization. It is important to notice that the
choice of sample size we propose here introduces stochastic iterative sequence
which might seem as an unnecessary complication if one is dealing with finite
sum problems. However we will show that the complexity remains the same and
asymptotically we get a.s. convergence, so the stochastic nature does not
alter the expected theoretical results. On the other hand, the intrinsic
nature of the sample size variation, based on the progress of the iterative
process, yields significant computational cost savings as demonstrated in the
numerical results.

The paper is organized as follows. The algorithm and some preliminaries are
given in Section 2, while Section 3 contains convergence analysis. Numerical
results are presented in Section 4. Some conclusions are drawn in Section 5.

##  2 Proposed algorithm

The following assumption summarizes the properties of the problem (1.1).

######  Assumption A 1.

Assume that $f_{i}(x)=F(x,\xi_{i}),$ $i=1,2,\ldots,$ are continuous, convex
and bounded from below with a constant $C$ for all $\xi_{i}.$

Notice that Assumption A1 implies that $f$ is convex and continuous function
as well as $f_{N}.$ Following the standard line search method, we assume that
a descent direction can be provided for any given function $f_{N}$.

######  Assumption A 2.

For any given $N$, $x$ and $B$ such that $mI\preceq B(x)\preceq MI,$ for some
positive and bounded constants $m\leq M$ we can compute a direction
$p_{N}\in\mathbb{R}^{n}$ such that

$p_{N}(x)=-B(x)\bar{g}_{N}(x)\quad\mbox{and}\sup_{g\in\partial
f_{N}(x)}g^{T}p_{N}(x)\leq-\frac{m}{2}\|\bar{g}_{N}(x)\|^{2},$

where $\bar{g}_{N}(x)\in\partial f_{N}(x).$

Let us briefly discuss the plausibility of the above assumption. One
possibility to generate such direction is presented in [33] where $B$ is the
BFGS matrix. If an oracle for calculating $\sup_{g\in\partial
f_{N}(x)}g^{T}p_{N}(x)$ is available, then we can take the subgradient descent
direction. Another approach would be to use gradient subsampling techniques
[35]. For directions that satisfy Assumption 2 the following result holds,
[33]. We provide the proof for the sake of completeness.

######  Lemma 2.1.

Let Assumptions A1 and A2 hold. Then there exists $\tau_{N}(x)>0$ and
$\gamma\in(0,1)$ such that the subgradient Armijo condition

$f_{N}(x+\alpha p_{N}(x))\leq f_{N}(x)-\gamma\alpha\|p_{N}(x)\|^{2}.$

holds for all $\alpha\in[0,\tau_{N}(x)].$

###### Proof.

Let us fix an arbitrary $N$ and an arbitrary $x\in\mathbb{R}^{n}$. If
$\bar{g}_{N}(x)=0$ the statement is obviously true. In the case
$\bar{g}_{N}(x)\neq 0$ we can define $\delta(\alpha):=f_{N}(x+\alpha
p_{N}(x))$, where $p_{N}(x)$ is a descent direction satisfying Assumption 2.
For such $p_{N}(x)$ there holds

$\delta^{\prime}(0)=\sup_{g\in\partial f_{N}(x)}g^{T}p_{N}(x)<0.$

Consider

$l(\alpha):=f_{N}(x)+\alpha\eta\sup_{g\in\partial f_{N}(x)}g^{T}p_{N}(x),$

for some $\eta\in(0,1)$. Given that $\sup_{g\in\partial
f_{N}(x)}g^{T}p_{N}(x)<0,$ $f_{N}$ is bounded from below and convex by
Assumption A1, there exists an unique intersection of the functions $\delta$
and $l$ on the interval $\alpha\in(0,\infty)$. Let us denote this intersection
by $\tau_{N}(x)$. Then, for all $\alpha\in[0,\tau_{N}(x)]$ there holds

$f_{N}(x+\alpha p_{N}(x))\leq f_{N}(x)+\alpha\eta\sup_{g\in\partial
f_{N}(x)}g^{T}p_{N}(x).$

Furthermore, Assumption 2 implies

$f_{N}(x+\alpha p_{N}(x))\leq
f_{N}(x)-\alpha\eta\frac{m}{2}\|\bar{g}_{N}(x)\|^{2}\leq
f_{N}(x)-\alpha\eta\frac{m}{2M^{2}}\|p_{N}(x)\|^{2}$

and the statement holds for $\gamma=\eta m/(2M^{2})$. ∎

The problem we are solving is defined by (1.3). Clearly the feasibility
condition $f_{N}(x)=f(x)$ can not be enforced in the general case of expected
value as in that case we should have $N\to\infty.$ Furthermore, neither the
deviation from feasible condition $|f(x)-f_{N}(x)|$ can be computed. Thus we
introduce an approximate infeasibility measure as a function $h(N)$ for
arbitrary integer $N.$ Assume that $h:\mathbb{N}\to\mathbb{R_{+}}\cup\\{0\\}$
is monotonically decreasing function such that
$\lim_{N\rightarrow\infty}h(N)=0.$ In other words, $h(N)$ is a proxy for
$|f(x)-f_{N}(x)|$. If we are solving a finite sum problem, i.e. if
$f(x)=f_{N_{\max}}(x)$ for a fixed $N_{\max}$ then for arbitrary $N\leq
N_{\max}$ we can define $h(N)=(N_{max}-N)/N_{max}.$ For the case of unbounded
$N$ one possible simple choice is $h(N)=N^{-1}.$ The merit function for IR is
defined in the usual way

$\Phi(x,N,\theta):=\theta f_{N}(x)+(1-\theta)h(N),$

where $\theta\in(0,1)$ is the penalty parameter used to give different weights
to the objective function and the measure of infeasibility and $N$ is an
integer that defines the level of accuracy in the approximate function
$f_{N}.$

At each iteration $k$ we have the accuracy parameter as an integer $N_{k},$
the solution estimate $x_{k}$, the penalty parameter $\theta_{k}$ and the
approximate objective function $f_{N_{k}}.$ The algorithm is as follows.

Algorithm 1: IR-NS (Inexact Restoration - Nonsmooth)

  * S0

Given
$x_{0}\in\mathbb{R}^{n},N_{0}\in\mathbb{N},\theta_{0},r\in(0,1),\beta,\gamma,\overline{\gamma}>0$.
Set $k=0$.

  * S1

Restoration phase. Find $\tilde{N}_{k+1}\geq N_{k}$ such that

$h(\tilde{N}_{k+1})\leq rh(N_{k}),$
$f_{\tilde{N}_{k+1}}(x_{k})-f_{N_{k}}(x_{k})\leq\beta h(N_{k}).$ (2.1)

  * S2

If
$\Phi(x_{k},\tilde{N}_{k+1},\theta_{k})-\Phi(x_{k},N_{k},\theta_{k})\leq\frac{1-r}{2}\left(h(\tilde{N}_{k+1})-h(N_{k})\right)$
set $\theta_{k+1}=\theta_{k}.$ Else

$\theta_{k+1}:=\frac{(1+r)(h(N_{k})-h(\tilde{N}_{k+1}))}{2\left[f_{\tilde{N}_{k+1}}(x_{k})-f_{N_{k}}(x_{k})+h(N_{k})-h(\tilde{N}_{k+1})\right]}.$

  * S3

Optimization Phase. Choose $N_{k+1}\leq\tilde{N}_{k+1},$
$p_{N_{k+1}}\in\mathbb{R}^{n}$ and $\alpha_{k}\in(0,1]$ such that

$f_{N_{k+1}}(x_{k}+\alpha_{k}p_{N_{k+1}}(x_{k}))-f_{\tilde{N}_{k+1}}(x_{k})\leq-\gamma\alpha_{k}||p_{N_{k+1}}(x_{k})||^{2},$
(2.2) $h(N_{k+1})\leq
h(\tilde{N}_{k+1})+\bar{\gamma}\alpha_{k}^{2}||p_{N_{k+1}}(x_{k})||^{2},$
(2.3)
$\Phi(x_{k}+\alpha_{k}p_{N_{k+1}}(x_{k}),N_{k+1},\theta_{k+1})-\Phi(x_{k},N_{k},\theta_{k+1})\leq\frac{1-r}{2}\bigl{(}h(\tilde{N}_{k+1})-h(N_{k})\bigr{)}.$
(2.4)

  * S4

Set $p_{k}=p_{N_{k+1}}(x_{k})$, $x_{k+1}=x_{k}+\alpha_{k}p_{k}$, $k:=k+1$ and
go to S1.

Let us briefly discuss the key points of IR-NS algorithm. In Step S1 the
feasibility is improved, i.e. a new sample size candidate $\tilde{N}_{k+1}$ is
chosen. Additionally, the value $f_{\tilde{N}_{k+1}}(x_{k})$ might increase
with respect to $f_{N_{k}}(x_{k})$ by at most $\beta h(N_{k}).$ Thus,
optimality can deteriorate with respect to the previous iteration but the
deterioration is controlled by the function $h,$ i.e., it depends on the
accuracy of the objective function. So, for smaller $N_{k}$ \- which means
looser approximation of the true objective function, the deterioration of
optimality can be relatively large, as we assume that we are still far away
from solution. Parameter $\beta$ can be arbitrary large, but finite. In some
applications (ex. finite sums) one can prove that such $\beta$ exists under
standard conditions. However, in general, since we do not impose
differentiability of the objective function nor any other special property,
the following assumption is needed.

######  Assumption A 3.

Suppose that there exists $\beta$ such that (2.1) holds for each $k.$

The penalty parameter is updated in such way that it ensures a decrease of the
merit function as stated in Lemma 2.2. Moreover, it can also be shown that the
sequence of $\theta_{k}$ is non-increasing and bounded away from zero which
prevents the optimality part to vanish from the merit function. The proof of
Lemma 2.2 is fundamentally the same as in [22, Lemma 2.1] and thus we omit it
here.

######  Lemma 2.2.

[22] Let Assumptions A1\- A3 hold. Then the sequence $\\{\theta_{k}\\}$
generated by Algorithm IR-NS is positive and non-increasing, the inequality

$\Phi(x_{k},\tilde{N}_{k+1},\theta_{k+1})-\Phi(x_{k},N_{k},\theta_{k+1})\leq\frac{1-r}{2}\left(h(\tilde{N}_{k+1})-h(N_{k})\right)$

holds and there exists $\theta^{*}>0$ such that
$\lim_{k\rightarrow\infty}\theta_{k}=\theta^{*}.$

In Step S3 we chose the sample size to be used in the subsequent iteration.
Notice that one possible choice is $N_{k+1}=\tilde{N}_{k+1}$ since (2.2)-(2.3)
are satisfied due to Lemma 2.1 and, as we will prove in Lemma 2.3, there
exists $\alpha_{k}$ which satisfies inequality (2.4) in that case as well. On
the other hand, in order to decrease the overall costs, we try to decrease the
sample size if it still provides the decrease in the merit function (2.4). The
resulting sample size $N_{k+1}$ can be larger, equal or smaller than $N_{k}$.
Our numerical study shows that allowing the decrease of a sample size is
beneficial in terms of overall function evaluations. In practical
implementations, we estimate the sample size lower bound $N^{trial}_{k+1}$
derived from (2.4) and let
$N_{k+1}\in\\{N^{trial}_{k+1},\lceil(N^{trial}_{k+1}+\tilde{N}_{k+1})/2\rceil,\tilde{N}_{k+1}\\}$.
We use the backtracking technique for finding $\alpha_{k}$, but at each
backtracking step we try all three candidate values for $N_{k+1}$. This is
just one possible approach and the optimal strategy remains an open question,
probably problem-dependent.

######  Lemma 2.3.

Let Assumptions A1\- A3 hold. Then, there exists $\gamma>0$ such that Step 3
of Algorithm IR-NS is well-defined.

###### Proof.

The algorithm is well defined if there exists a choice of
$N_{k+1}\leq\tilde{N}_{k+1}$ and a descent direction $p_{k}$ such that (2.2) -
(2.4) hold for some $\alpha_{k}>0$ and a suitable $\gamma>0$ for each $k.$ Let
us take $N_{k+1}=\tilde{N}_{k+1}$ and retain the same sample so that
$f_{N_{k+1}}=f_{\tilde{N}_{k+1}}$. In that case Lemma 2.1 implies the
existence of $\tau_{k}:=\tau_{N_{k+1}}(x_{k})>0$ such that the inequality
(2.2) holds for all $\alpha\in[0,\tau_{k}]$. Since (2.3) is trivially
satisfied for this choice of $N_{k+1}$, it remains to prove the existence of
$\alpha_{k}\in[0,\tau_{k}]$ such that (2.4) holds. By (2.2), (2.3) and Lemma
2.2, for all $\alpha\in[0,\tau_{k}],$

$\displaystyle\Phi(x_{k}+\alpha
p_{k},N_{k+1},\theta_{k+1})-\Phi(x_{k},N_{k},\theta_{k+1})$
$\displaystyle=\Phi(x_{k}+\alpha
p_{k},N_{k+1},\theta_{k+1})-\Phi(x_{k},\tilde{N}_{k+1},\theta_{k+1})+\Phi(x_{k},\tilde{N}_{k+1},\theta_{k+1})-\Phi(x_{k},N_{k},\theta_{k+1})$
$\displaystyle\leq\Phi(x_{k}+\alpha
p_{k},N_{k+1},\theta_{k+1})-\Phi(x_{k},\tilde{N}_{k+1},\theta_{k+1})+\frac{1-r}{2}\left(h(\tilde{N}_{k+1})-h(N_{k})\right)$
$\displaystyle=\theta_{k+1}\left(f_{N_{k+1}}(x_{k}+\alpha
p_{k})-f_{\tilde{N}_{k+1}}(x_{k})\right)+\frac{1-r}{2}\left(h(\tilde{N}_{k+1})-h(N_{k})\right)$
$\displaystyle\leq-\theta_{k+1}\gamma\alpha||p_{k}||^{2}+\frac{1-r}{2}\left(h(\tilde{N}_{k+1})-h(N_{k})\right)\leq\frac{1-r}{2}\left(h(\tilde{N}_{k+1})-h(N_{k})\right).$

Therefore, (2.4) holds for all $\alpha\in[0,\tau_{k}]$. ∎

Notice that in the above Lemma 2.3 we proved only that the algorithm is well
defined, i.e., we can always take $N_{k+1}=\tilde{N}_{k+1}$ and the $(k+1)$th
iteration is well defined. However, other possibilities for $N_{k+1}$ exists
and we discuss some of them in Section 4. Since the sample size sequence is
not monotonically increasing in general, it is not obvious that $N_{k}$ tends
to infinity. Nevertheless, using essentially the same proof as in [22, Theorem
2.1], we conclude that infeasibility measure tends to zero yielding the result
of $\lim_{k\to\infty}N_{k}=\infty$. Specially, for the finite sum problem we
conclude that the full sample is reached after a finite number of iterations.
The proof of Theorem 2.1 in [22] contains an important relation stated below

$\sum_{k=0}^{\infty}h(N_{k})\leq C_{1}<\infty,$ (2.5)

where $C_{1}>0$ is a constant, that we will use in further convergence
analysis presented in the next Section.

Let us now provide more insights regarding the stochastic concept of the
proposed algorithm. IR-NS yields stochastic sequence of iterates $x_{k}$. The
stochastic nature comes from the sequence of random variables $N_{k}$ that
determine the samples to be used for the SAA functions. Assume that we are at
iteration $k$ and $x_{k}$ is known. Denote by $\mathcal{F}_{k}$ the
$\sigma$-algebra generated by $x_{0},...,x_{k}$, i.e., by random variables
that determine $f_{\tilde{N}_{j}},j=1,...,k$ and $f_{N_{j}},j=0,...,k.$ Since
the samples are assumed to be i.i.d., we have conditionally unbiased
estimators. More precisely, at the beginning of Step S1 of the algorithm a new
sample size $\tilde{N}_{k+1}$ is chosen and a random sample is generated to
obtain $f_{\tilde{N}_{k+1}}$. Thus, since $x_{k}$ is
$\mathcal{F}_{k}$-measurable (i.e., known at that point of the algorithmic
procedure), there holds

$E\left[f_{\tilde{N}_{k+1}}(x_{k})|\mathcal{F}_{k}\right]=f(x_{k}),$ (2.6)

where $E\left[\cdot|\mathcal{F}_{k}\right]$ denotes the conditional
expectation with respect to $\mathcal{F}_{k}$ [35]. Also
$E\left[f_{N_{k+1}}(x_{k})|\mathcal{F}_{k}\right]=f(x_{k})$. However,
$E\left[f_{N_{k+1}}(x_{k+1})|\mathcal{F}_{k}\right]$ is not equal to
$f(x_{k+1})$ in general because $x_{k+1}$ is dependent on $N_{k+1}$. More
precisely, the second round of stochastic influence within iteration $k$ comes
at the Step S3 where we choose $N_{k+1}$ which may yield totally different
sample for $f_{N_{k+1}}$ with respect to $f_{\tilde{N}_{k+1}}$ in general
(each trial sample size may yield different sample). Moreover, the direction
$p_{N_{k+1}}(x_{k+1})$ and the step size $\alpha_{k}$ directly depend on the
generated samples and thus we lose the martingale property. This is a common
situation in stochastic line search (see [14] for instance). In Step S4, we
set the next iterate and return to Step S1, repeating the procedure.

##  3 Convergence analysis

The convergence analysis is performed under the set of standard assumptions
for stochastic problems stated below. We analyze conditions needed for a.s.
convergence of IR-NS and provide complexity result at the end of this section.
The two assumptions stated in this Section are needed to ensure that the
Uniform Law of Large Numbers (ULLN) holds.

######  Assumption A 4.

The objective function $f$ has bounded level sets.

This assumption holds if the objective function is strongly convex for
example, and we have the following result.

######  Lemma 3.1.

Let Assumptions A1-A4 hold. Suppose that there exists a constant $C_{0}$ such
that $F(x_{0},\xi)\leq C_{0}$ for any $\xi$. Then $f(x_{k})\leq C_{2}$ holds
for all $k$, i.e., $\\{x_{k}\\}_{k\in\mathbb{N}}\subseteq D,$ where

$D=\\{x\in\mathbb{R}^{n}\;|\;f(x)\leq C_{2}\\}$

and $C_{2}=C_{0}+2\beta C_{1}$.

###### Proof.

The set $D$ is compact by Assumption A4. Using inequalities (2.1)-(2.2), for
all $k$ we obtain

$f_{N_{k+1}}(x_{k+1})\leq
f_{\tilde{N}_{k+1}}(x_{k})-\gamma\alpha_{k}||p_{N_{k+1}}(x_{k})||^{2}\leq
f_{N_{k}}(x_{k})+\beta h(N_{k}).$

Furthermore, using the induction argument and (2.5) we get

$f_{N_{k+1}}(x_{k+1})\leq f_{N_{0}}(x_{0})+\beta\sum_{j=0}^{k}h(N_{j})\leq
f_{N_{0}}(x_{0})+\beta C_{1},$

for all $k=0,1,...$. Obviously, the assumption of uniformly bounded $F$ at the
initial point $x_{0}$ implies that $f_{N_{0}}(x_{0})\leq C_{0}$ and we obtain

$f_{N_{k}}(x_{k})\leq C_{0}+\beta C_{1},$ (3.1)

for all $k=1,2,...$. Finally, by (2.6) and inequalities (2.1) and (3.1) we get

$f(x_{k})=E\left[f_{\tilde{N}_{k+1}}(x_{k})|{\mathcal{F}}_{k}\right]\leq
E\left[f_{N_{k}}(x_{k})+\beta h(N_{k})|{\mathcal{F}}_{k}\right]\leq
C_{0}+2\beta C_{1}:=C_{2},$

which completes the proof. ∎

######  Assumption A 5.

The function $F$ is dominated by an integrable function on a bounded open set
$\tilde{D}^{0}$ such that $D\subset\tilde{D}^{0}.$

Under the stated assumptions the ULLN [7] implies that
$\lim_{N\rightarrow\infty}\sup_{x\in D}|f_{N}(x)-f(x)|=0$ a.s. Notice that
this equality holds trivially if the sample is finite and the full sample is
eventually achieved and retained. Denote by
$X^{*}=\\{x\in\mathbb{R}^{n}:f(x)=\inf_{y}f(y):=f^{*}\\}$ the set of solutions
for problem (1.1). Define

$t_{k}:=\max_{x,y\in\tilde{D}}\\{|f(x)-f_{N_{k+1}}(x)|+|f(y)-f_{\tilde{N}_{k+1}}(y)|\\},$
(3.2)

where $\tilde{D}$ is a compact enlargement of $D,$ i.e., $\tilde{D}$ is the
closure of an open set $\tilde{D}^{0}\supset D.$ Therefore, both $D$ and
$\tilde{D}$ are compact sets and $D\mathchar 13608\relax\tilde{D}.$ Notice
that ULLN and the fact $h(N_{k})\to 0$ imply that $t_{k}\to 0$ a.s. if
$N_{k}\to\infty.$ Let us analyse the convergence depending on properties of
the step size sequence $\\{\alpha_{k}\\}$ and the error sequence
$\\{t_{k}\\}$.

######  Theorem 3.1.

Let Assumptions A1-A5 hold and $\\{x_{k}\\}$ be a sequence generated by
Algorithm IR-NS. If $\alpha_{k}\geq\overline{\alpha}>0$ for all
$k\in\mathbb{N}$ then there exists an accumulation point $x^{*}$ of
$\\{x_{k}\\}$ which is a solution of problem (1.1) a.s.

###### Proof.

Denote $\bar{g}_{k}=\bar{g}_{N_{k}}(x_{k}).$ Then Assumption A2 and (2.2)
imply

$f_{N_{k+1}}(x_{k+1})\leq
f_{\tilde{N}_{k+1}}(x_{k})-\gamma\alpha_{k}||p_{k}||^{2}\leq
f_{\tilde{N}_{k+1}}(x_{k})-\eta\alpha_{k}||\overline{g}_{k}||^{2},$

where $\eta=\gamma m^{2}.$ Furthermore,

$\displaystyle f(x_{k+1})$ $\displaystyle\leq
f_{\tilde{N}_{k+1}}(x_{k})-\eta\alpha_{k}||\overline{g}_{k}||^{2}+f(x_{k+1})-f_{N_{k+1}}(x_{k+1})$
$\displaystyle\leq
f(x_{k})-\eta\alpha_{k}||\overline{g}_{k}||^{2}+|f(x_{k+1})-f_{N_{k+1}}(x_{k+1})|+|f_{\tilde{N}_{k+1}}(x_{k})-f(x_{k})|.$

From the definition of $t_{k}$ (3.2), we obtain

$f(x_{k+1})\leq f(x_{k})-\eta\bar{\alpha}||\overline{g}_{k}||^{2}+t_{k}.$
(3.3)

We will show that $\liminf_{k\rightarrow\infty}||\overline{g}_{k}||^{2}=0.$
Assume the contrary, i.e., that $||\overline{g}_{k}||^{2}\geq\varrho>0$ for
some $\varrho>0$ and all $k.$ Then
$\eta\overline{\alpha}||\overline{g}_{k}||^{2}\geq\eta\overline{\alpha}\varrho>0.$
Since $t_{k}\to 0$ a.s., there exists $\overline{k}$ such that for all
$k\geq\overline{k}$ there holds
$t_{k}\leq\frac{1}{2}\eta\overline{\alpha}||\overline{g}_{k}||^{2}\text{
a.s.}$ and thus (3.3) implies $f(x_{k+1})\leq
f(x_{k})-\eta\overline{\alpha}/2$ a.s. Equivalently, for all $s\in\mathbb{N}$
we have

$f(x_{\overline{k}+s})\leq
f(x_{\overline{k}})-\frac{s}{2}\eta\overline{\alpha}\varrho\quad\mbox{a.s.}$
(3.4)

Letting $s\to\infty$ yields a contradiction with the Assumption A1 which
implies that $f$ is bounded from bellow. Therefore, we conclude that there
there exists $K\subseteq\mathbb{N}$ such that $\lim_{k\in
K}\overline{g}_{k}=0$ a.s. Since $\\{x_{k}\\}\subset D$ and $D$ is compact
there follows that there exist $K_{1}\subseteq K$ and $x^{*}\in D$ such that
$x^{*}=\lim_{k\in K_{1}}x_{k}.$ Now, using the fact that
$\overline{g}_{k}\in\partial f_{N_{k+1}}(x_{k})$, for all $x\in\mathbb{R}^{n}$
we have $f_{N_{k+1}}(x)\geq f_{N_{k+1}}(x_{k})+\overline{g}^{T}_{k}(x-x_{k}).$
Thus, for arbitrary $x\in\tilde{D}$ we have

$\displaystyle f(x)$ $\displaystyle\geq
f_{N_{k+1}}(x_{k})+\overline{g}^{T}_{k}(x-x_{k})+f(x)-f_{N_{k+1}}(x)$
$\displaystyle=f(x_{k})+\overline{g}^{T}_{k}(x-x_{k})-\left(f_{N_{k+1}}(x)-f(x)+f(x_{k})-f_{N_{k+1}}(x_{k})\right)$
$\displaystyle\geq
f(x_{k})+\overline{g}^{T}_{k}(x-x_{k})-\left(|f(x)-f_{N_{k+1}}(x)|+|f(x_{k})-f_{N_{k+1}}(x_{k})|\right).$
(3.5)

Therefore, $f(x)\geq f(x_{k})-||\overline{g}_{k}||||x-x_{k}||-2t_{k}$. Taking
the limit over $K_{1}$ and using the fact that $||x-x_{k}||$ is bounded, we
obtain that for every $x\in\tilde{D}$ there holds

$f(x)\geq f(x^{*}),\mbox{ a.s.}$ (3.6)

Recall that $x^{*}\in D$ and $\tilde{D}$ is a compact enlargement of $D$ so
$x^{*}$ cannot be on the boundary of $\tilde{D}$ and there exists $\epsilon>0$
such that $\mathcal{B}(x^{*},\epsilon)\subset\tilde{D}$ and we conclude that
$x^{*}$ is a local minimizer of $f$ a.s. Since $f$ is assumed to be convex, we
conclude that $x^{*}\in X^{*}$ a.s. ∎

We can also prove that every strictly strong accumulation point [36] is a
solution a.s. A point $x^{*}$ is called strictly strong accumulation point of
the sequence $\\{x_{k}\\}_{k\in\mathbb{N}}$ if there exists a subsequence
$K\subseteq\mathbb{N}$ and a constant $b\in\mathbb{N}$ such that
$\lim_{k_{i}\in K}x_{k_{i}}=x^{*}$ and $k_{i+1}-k_{i}\leq b$ for any two
consecutive elements $k_{i},k_{i+1}\in K.$ According to the available
literature, [7, 37], and up to the best of our knowledge, stronger statement
in a.s. sense is not possible without some additional assumptions on the rate
of increase of $N_{k}.$

######  Theorem 3.2.

Assume that the conditions of Theorem 3.1 hold. Then every strictly strong
accumulation point of the sequence $\\{x_{k}\\}$ is a solution of problem
(1.1) a.s.

###### Proof.

Let $x^{*}$ be an arbitrary strictly strong accumulation point of the sequence
$\\{x_{k}\\}$, i.e., $x^{*}=\lim_{i\to\infty}x_{k_{i}}$ and
$s_{i}:=k_{i+1}-k_{i}\leq b$ for every $i\in\mathbb{N}.$ Since (3.3) holds for
each $k\in\mathbb{N},$ we obtain

$f(x_{k_{i+1}})\leq
f(x_{k_{i}})-\eta\overline{\alpha}\sum_{j=0}^{s_{i}-1}||\overline{g}_{k_{i}+j}||^{2}+\sum_{j=0}^{s_{i}-1}t_{k_{i}+j}\leq
f(x_{k_{i}})-\eta\overline{\alpha}||\overline{g}_{k_{i}}||^{2}+\omega_{i},$

where $\omega_{i}=\sum_{j=0}^{b-1}t_{k_{i}+j}.$ Notice that
$\omega_{i}\rightarrow 0,i\rightarrow\infty$ a.s. We want to show that

$\liminf_{i\to\infty}||\overline{g}_{k_{i}}||^{2}=0\mbox{ a.s.}$ (3.7)

Assume the contrary, i.e., for all $i\in\mathbb{N}$ there holds
$||\overline{g}_{k_{i}}||^{2}\geq\varrho>0$ for some $\varrho>0.$ Then,
$\eta\overline{\alpha}||\overline{g}_{k_{i}}||^{2}\geq\eta\overline{\alpha}\varrho>0$
for all $i\in\mathbb{N}.$ Therefore, there exists $\overline{i}$ such that for
all $i\geq\overline{i}$ there holds
$\omega_{i}\leq\frac{1}{2}\eta\overline{\alpha}\varrho$ a.s. and thus
$f(x_{k_{i+1}})\leq f(x_{k_{i}})-\frac{1}{2}\eta\overline{\alpha}\varrho$ a.s.
Letting $i\to\infty$ in the last inequality we obtain

$f(x^{*})\leq f(x^{*})-\frac{1}{2}\eta\overline{\alpha}\varrho<f(x^{*}),$

which is contradiction. So, (3.7) holds and repeating the steps (3.4)-(3.6)
from the proof of Theorem 3.1, we obtain the result, i.e. $x^{*}\in X^{*}$
a.s. ∎

Next, we show that the convergence result as in Theorem 3.1 can be obtained
under weaker assumptions on the step size sequence, but assuming that the
sample size $N_{k}$ is eventually increased fast enough such that
$\sum_{k=0}^{\infty}t_{k}<\infty$. For instance, if the sample is cumulative,
the log bound given in Proposition 3.5 of [10] holds and
$\sum_{k=0}^{\infty}t_{k}<\infty$ is true if $N_{k}\geq e^{k}$. Therefore, one
can switch to exponential growth after a certain number of iterations of IR-NS
algorithm, taking advantage of cheap iterations in early stages and
theoretically proved convergence for fast increase of the sample size sequence
in the later stages of algorithm. The switching point is an interesting
problem itself, but beyond the scope of this paper.

######  Theorem 3.3.

Let Assumptions A1-A5 hold and $\\{x_{k}\\}$ be a sequence generated by
Algorithm IR-NS. If $\sum_{k=0}^{\infty}\alpha_{k}=\infty$ and
$\sum_{k=0}^{\infty}t_{k}<\infty$ then there exists an accumulation point
$x^{*}$ of $\\{x_{k}\\}$ which is a solution of problem (1.1).

###### Proof.

Following the steps of the proof of Theorem 3.1 we obtain $f(x_{k+1})\leq
f(x_{k})-\eta\alpha_{k}||\overline{g}_{k}||^{2}+t_{k}$ for every $k$ and thus

$f(x_{k+1})\leq
f(x_{0})-\eta\sum_{i=0}^{k}\alpha_{i}||\overline{g}_{i}||^{2}+\sum_{i=0}^{k}t_{i}.$

The function $f$ is bounded from below and $\sum_{k=0}^{\infty}t_{k}<\infty,$
so we conclude

$\sum_{k=0}^{\infty}\alpha_{k}||\overline{g}_{k}||^{2}<\infty.$ (3.8)

Furthermore, the assumption $\sum_{k=0}^{\infty}\alpha_{k}=\infty$ implies the
existence of a subset $K_{1}$ such that $\lim_{k\in K_{1}}\overline{g}_{k}=0$.
Indeed, if we assume the contrary, i.e., that there exists $\varepsilon>0$
such that $\|\overline{g}_{k}\|\geq\varepsilon>0$ for $k$ large enough, then
we obtain

$\sum_{k=0}^{\infty}\alpha_{k}||\overline{g}_{k}||^{2}\geq\sum_{k=0}^{\infty}\alpha_{k}\varepsilon^{2}=\varepsilon^{2}\sum_{k=0}^{\infty}\alpha_{k}=\infty,$

which is in contradiction with (3.8). Since the whole sequence
$\\{x_{k}\\}_{k\in\mathbb{N}}$ is bounded due to Lemma 3.1, there exist
$K_{2}\subseteq K_{1}$ and $x^{*}\in D$ such that $\lim_{k\in
K_{2}}x_{k}=x^{*}$. Now, repeating the proof of Theorem 3.1 \- the part after
(3.4), we conclude that $x^{*}\in X^{*}.$ ∎

The following result is based on considerations in [9] and [38] and
essentially yields worst-case complexity analysis with respect to the expected
objective function value.

######  Theorem 3.4.

Let Assumptions A1-A5 hold, $\varepsilon>0$ and $\\{x_{k}\\}$ be a sequence
generated by Algorithm IR-NS. Furthermore, assume that
$\alpha_{k}\geq\overline{\alpha}>0$ for all $k\in\mathbb{N}$ and
$\sum_{k=0}^{\infty}t_{k}\leq\overline{t}<\infty.$ Then, after at most

$\overline{k}=\big{\lceil}\frac{R^{2}(\overline{t}+f(x_{0})-f^{*})}{\eta\overline{\alpha}}\varepsilon^{-2}\big{\rceil}$

iterations, we have $E\left[f(x_{\overline{k}})-f^{*}\right]\leq\varepsilon,$
where $R$ is the diameter of $D.$

###### Proof.

First, notice that (3.8) holds and since $\alpha_{k}\geq\overline{\alpha}$ we
obtain $\lim_{k\to\infty}||\overline{g}_{k}||^{2}=0.$ Take arbitrary
$\varepsilon>0$ and define $\varepsilon_{1}=\varepsilon/R.$ Since
$\overline{g}_{k}$ tends to zero, there exists $\overline{k}$ such that
$||\overline{g}_{\overline{k}}||\leq\varepsilon_{1}.$ Let $\overline{k}$ be
the first such iteration. Then for $k=0,1,\ldots,\overline{k}-1$ we have
$||\overline{g}_{k}||>\varepsilon_{1}.$ Moreover, from (3.3) we get
$t_{k}+f(x_{k})-f(x_{k+1})\geq\eta\overline{\alpha}\varepsilon_{1}^{2}$ for
$k=0,1,\ldots,\overline{k}-1$ and by summing up both sides of this inequality
and using $\sum_{k=0}^{\infty}t_{k}\leq\overline{t}<\infty$ we obtain

$\eta\overline{\alpha}\varepsilon_{1}^{2}\overline{k}\leq\overline{t}+f(x_{0})-f(x_{\overline{k}})\leq\overline{t}+f(x_{0})-f^{*},$

i.e.,
$\overline{k}\leq(\overline{t}+f(x_{0})-f^{*})/(\varepsilon_{1}^{2}\eta\overline{\alpha})=\varepsilon^{-2}(R^{2}(\overline{t}+f(x_{0})-f^{*}))/(\eta\overline{\alpha}).$
Since $f_{N_{k+1}}$ is convex and $\overline{g}_{k}\in\partial
f_{N_{k+1}}(x_{k})$ there holds $f_{N_{\overline{k}+1}}(x^{*})\geq
f_{N_{\overline{k}+1}}(x_{\overline{k}})+\overline{g}^{T}_{\overline{k}}(x^{*}-x_{\overline{k}})$,
i.e.,

$f_{N_{\overline{k}+1}}(x_{\overline{k}})-f_{N_{\overline{k}+1}}(x^{*})\leq\overline{g}^{T}_{\overline{k}}(x_{\overline{k}}-x^{*})\leq||\overline{g}_{\overline{k}}||||x^{*}-x_{\overline{k}}||\leq\varepsilon_{1}R=\varepsilon.$

Denote by ${\mathcal{F}}_{\bar{k}}$ the $\sigma$-algebra generated by
$x_{0},...,x_{\bar{k}}$. Since the sample is assumed to be i.i.d. and the
approximate functions are computed as sample average, we obtain

$E\left[f(x_{\overline{k}})-f(x^{*})\right]=E\left[E\left[f_{N_{\overline{k}+1}}(x_{\overline{k}})-f_{N_{\overline{k}+1}}(x^{*})|{\mathcal{F}}_{\bar{k}}\right]\right]\leq\varepsilon.$

∎

Let us conclude this section by considering finite sum case which falls into
the IR-NS framework. Recall that $h(N_{k})\to 0.$ So, in the case of finite
sum we have $N_{k}=N_{\max}$ for all $k\geq k_{0}$ where $k_{0}$ is random,
but finite. Moreover, $t_{k}$ becomes zero eventually, so the summability of
$t_{k}$ holds. Furthermore, (3.3) reveals that $f(x_{k+1})\leq f(x_{k})$ for
all $k\geq k_{0}$ and thus the iterations remain in the level set
$\mathcal{L}=\\{x|f(x)\leq f(x_{k_{0}})\\}.$ If the level set is compact then
the Assumption A4 is obviously satisfied. Finally, notice that Assumption A5
is needed only to ensure that $t_{k}$ tends to zero a.s. which is obviously
true in the finite sum case. Also, notice that in the strongly convex finite
sum case there exists $C$ such that all $f_{i}$ functions are bounded from
bellow by $C.$ Therefore the following result holds.

######  Corollary 3.1.

Let Assumptions A2-A3 hold and assume $\sum_{k}\alpha_{k}=\infty.$ If
$f=f_{N_{max}}$ and $f_{i},i=1,...,N_{max}$ are continuous and strongly
convex, then there exists an accumulation point $x^{*}$ of $\\{x_{k}\\}$ which
is a solution of problem (1.1). Moreover, if
$\alpha_{k}\geq\overline{\alpha}>0$ for all $k\in\mathbb{N},$ then the worst-
case complexity is of order ${\mathcal{O}}(\varepsilon^{-2})$.

##  4 Numerical experiments

In this section, we test IR-NS variable sample size scheme on two classes of
nonsmooth convex problems: 1) Finite Sums (FS), i.e., bounded sample size with
real-world data, and 2) Expected Residual Minimization (ERM) reformulation of
Stochastic Linear Complementarity Problems (SLCP) with unbounded sample size
and simulated data. The first class belongs to the machine learning framework
and considers $L_{2}$-regularized binary hinge loss functions (see [33] and
the references therein) for binary classification. The considered data sets
are given in Table 1 and the problem is of the form

$\min_{x\in\mathbb{R}^{n}}f(x):=\frac{\lambda}{2}||x||^{2}+\frac{1}{N_{max}}\sum_{i=1}^{N_{max}}\max(0,1-z_{i}x^{T}w_{i}),$

where $\lambda=10^{-5}$ is a regularization constant, $w_{i}\in\mathbb{R}^{n}$
are the input features, $z_{i}\in\\{\pm 1\\}$ the corresponding labels,
$N_{max}$ is the size of relevant data set (testing or training).

| Data set | $N$ | $n$ | $N_{train}$ | $N_{test}$ | $Max_{FEV}$
---|---|---|---|---|---|---
1 | SPLICE [39] | 3175 | 60 | 2540 | 635 | $10^{6}$
2 | MUSHROOMS [40] | 8124 | 112 | 6500 | 1624 | $10^{6}$
3 | ADULT9 [39] | 32561 | 123 | 26049 | 6512 | $10^{7}$
4 | MNIST(binary) [41] | 70000 | 784 | 60000 | 10000 | $10^{7}$
Table 1: Properties of the data sets used in the experiments.

SLCP consists of finding a vector $x\in\mathbb{R}^{n}$ such that

$x\geq 0,M(\xi)x+q(\xi)\geq 0,x^{T}(M(\xi)x+q(\xi))=0,\xi\in\Omega,$

where $\Omega$ is the underlying sample space, $M(\xi)\in\mathbb{R}^{n,n}$ is
a random matrix and $q(\xi)\in\mathbb{R}^{n}$ is a random vector. ERM
reformulation (see [42] for example) is defined as follows

$\min f(x)=E\left[||\tilde{F}(x,\xi)||^{2}\right],\quad\mbox{s. t. }\quad
x\geq 0,\ $

where $\tilde{F}(x,\xi):\mathbb{R}^{n}\times\Omega\rightarrow\mathbb{R}^{n}$,
$\tilde{F}(x,\xi)=\phi(x,M(\xi)x+q(\xi))$ and
$\phi:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is the NCP function defined as
$\phi(a,b)=\min\\{a,b\\}$.

The SAA approximate objective function (1.2) is defined as

$f_{N_{k}}(x)=\frac{1}{N_{k}}\sum_{j=1}^{N_{k}}f_{j}(x)$

with
$f_{j}(x)=\|\tilde{F}(x,\xi_{j})\|^{2}=\sum_{l=1}^{n}\left(\min\\{x_{l},[M(\xi_{j})x]_{l}+[q(\xi_{j})]_{l}\\}\right)^{2}.$

Since numerical results for deterministic (full sample) problem provided in
[33] reveal the advantages of BFGS-type methods in nonsmooth optimization, we
chose to use the method proposed therein for finding a descent direction
satisfying Assumption A2. The functions in consecutive iterations differ in
general, and $y_{k}$ needed for BFGS update is the difference of subgradients
of different SAA functions, a safeguard is needed to ensure that the resulting
matrices are uniformly positive definite. Thus we start with the identity
matrix and skip the BFGS update if
$y_{k}(x_{k+1}-x_{k})<10^{-4}\|y_{k}\|^{2}$. Both types of tested problem, FS
and ERM allow us to calculate $\sup_{g\in\partial f_{N}(x)}p^{T}g$ which is
crucial for finding the descent BFGS direction. We denote the proposed
algorithm by IRBFGS to emphasize the fact that the BFGS directions are used.

The parameters of IRBGFS algorithm are
$\theta_{0}=0.9,\;r=0.95,\overline{\gamma}=1$ and $\gamma=10^{-4}.$ The
function $h$ is defined as $h(N_{k})=(N-N_{k})/N$ for FS and
$h(N_{k})=1/N_{k}$ for ERM problem. Thus, we have
$\tilde{N}_{k+1}=\min\\{N,\lceil N-r(N-N_{k})\rceil\\}$ for bounded and
$\tilde{N}_{k+1}=\lceil N_{k}/r\rceil$ for unbounded sample case.
$N_{0}=\lceil 0.1N\rceil$ for FS, while for ERM problems we take $N_{0}=1000$.
Step S3 is performed as already stated: we estimate the sample size lower
bound $N^{trial}_{k+1}$ derived from (2.4) and let
$N_{k+1}\in\\{N^{trial}_{k+1},\lceil(N^{trial}_{k+1}+\tilde{N}_{k+1})/2\rceil,\tilde{N}_{k+1}\\}$.
The backtracking technique for finding $\alpha_{k}=0.5^{j}$ is used, but at
each backtracking step we try all three candidate values for $N_{k+1}$. We use
cumulative samples, although other approaches are feasible as well. The value
$N^{trial}_{k+1}$ is calculated as follows: for FS

$N_{k+1}^{trial}:=N_{k}+\frac{1-r}{2}\cdot\frac{\tilde{N}_{k+1}-N_{k}}{1-\theta_{k+1}}-\hat{\theta}_{k+1}\left(\gamma\alpha||p_{k-1}||^{2}-f_{\tilde{N}_{k+1}}(x_{k})+f_{N_{k}}(x_{k})\right),$

where $\hat{\theta}_{k+1}=N\cdot\frac{\theta_{k+1}}{1-\theta_{k+1}}$; for ERM

$N_{k+1}^{trial}:=\frac{1-\theta_{k+1}}{\frac{1-r}{2}\cdot\frac{N_{k}-\tilde{N}_{k+1}}{\tilde{N}_{k+1}N_{k}}+\frac{1-\theta_{k+1}}{N_{k}}+\theta_{k+1}\left(\gamma\alpha||p_{k-1}||^{2}-f_{\tilde{N}_{k+1}}(x_{k})+f_{N_{k}}(x_{k})\right)}.$

The motivation for these choices comes from condition (2.4) from Step S3. The
merit function at new point should be decreased for at least
$\frac{1-r}{2}(h(\tilde{N}_{k+1})-h(N_{k}))$. Therefore, approximating
$||p_{k}||$ with $||p_{k-1}||$ and using (2.2) and (2.3) from Step S3, we
obtain the lower bound $N_{k+1}^{trial}$ for $N_{k+1}$. If this value falls
below $N_{0}$, we simply take $N_{k+1}^{trial}=N_{0}$.

Our numerical study has two goals:

  * 1)

to investigate if the variable sample size approach is beneficial in terms of
overall optimization costs;

  * 2)

to investigate if the potential decrease of the sample size coming from S3 is
beneficial.

This is why we compare the proposed IRBFGS method to: 1) FBFGS which takes the
full sample (when applicable) at each iteration, i.e., in FS problems
$N_{k}=N_{max}$ for each $k$; 2) HBFGS which takes $N_{k+1}=\tilde{N}_{k+1}$
for each $k$. The criterion for comparison is the number of scalar products
denoted by FEV. We report the average values of 10 independent runs. The
algorithms are stopped when the maximum number of scalar products, $Max_{FEV}$
is reached. In the FS case, we track the value of the (full sample) objective
function, while in the ERM case we track the Euclidean difference between
$x_{k}$ and the solution $x^{*}$ since the objective function is not
computable while the solution is known in advance.

Figure 1: FS Problem. Testing loss versus function evaluations. Figure 2: ERM
Problem. The error $\|x_{k}-x^{*}\|$ versus function evaluations

Figure 3: IRBFGS sample size versus HBFGS sample size sequence: FS Problem -
SPLICE data set (left) and ERM Problem (right).

Fig. 1 shows the results on FS problems with uniform random $x_{0}$. Since
training and testing errors behave similarly, we report only the testing
error. The $y$-axes are in logarithmic scale. The plots demonstrate the
computational savings obtained by IRBFGS in almost all cases. In fact, both
subsampled method, IRBFGS and HBFGS use smaller FEV to obtain the solutions of
the same quality as the full BFGS - FBFGS. Comparing IRBFGS and HBFGS, one can
see that IRBFGS is more efficient and occasional decrease of $N_{k}$ in Step
S3 is beneficial in terms of computational effort measured by FEV. Typical
behavior of the sample size sequence is plotted in Fig. 3 (left).

ERM problems are formed as in [42, 43, 44] where the first order methods were
tested. Here we proceed with the nonsmooth BFGS direction. We report the
results for problem with $n=100$ and volatility measure $\sigma=10$.
$Max_{FEV}$ is set to $10^{5}$ and the average ending sample size is 4714 for
IRBFGS and 3110 for HBFGS. The results and typical behavior of the sample size
sequence are presented in Fig. 2 and 3 (right), respectively. As we can see,
IRBFGS algorithm significantly outperforms the heuristic scheme HBFGS.

##  5 Conclusions

We proposed a framework for minimization of nonsmooth convex function in the
form of mathematical expectation. The general algorithm is defined within
Inexact Restoration approach, using a suitable approximate function computed
as the sample average approximation in each iteration. The sample size is
determined adaptively, taking into account the progress toward the stationary
point and thus balancing the computational cost and accuracy in endogenous way
without heuristic elements. The Armijo line search rule, adapted to the
nonsmooth function, is used for step sizes. Algorithm is defined with a
general descent direction for nonsmooth function, assuming that a suitable
oracle for direction computation is available. It is proved, using the
standard IR methodology, that the sample size tends to infinity or attains the
fixed maximal value. Therefore, the method generates the approximate solution
of desired accuracy but with lower computational costs. The theoretical
analysis reveals a.s. convergence towards stationary points under the set of
standard assumptions. The numerical experiments are based on the BFGS
direction adapted to the nonsmooth environment [33]. The oracle for computing
the direction is taken from literature for the hinge loss problems and
Expected Residual Minimization of Stochastic Linear Complementarity Problem.
The obtained numerical results are in line with the theoretical considerations
and confirm the efficiency of the algorithm.

Acknowledgement. We are grateful to the editor and the anonymous referees for
their comments that helped us to improve the paper.
Funding The work of Krejić and Krklec Jerinkić is supported by Provincial
Secretariat for Higher Education and Scientific Research of Vojvodina, grant
no. 142-451-2593/2021-01/2. The work of Ostojić is supported by the Ministry
of Education, Science and Technological Development, Republic of Serbia.
Availability statement The datasets analysed during the current study are
available in the MNIST database of handwritten digits [41], LIBSVM Data:
Classification (Binary Class) [39] and UCI Machine Learning Repository [40].

## Declarations

Conflict of interest The authors declare no competing interests.

## References

  * [1] A. Bagirov, N.Karmitsa, M. Mäkelä, Introduction to Nonsmooth Optimization, Springer, (2014), https://doi.org/10.1007/978-3-319-08114-4.
  * [2] V. Cevher, S. Becker, M. Schmidt, Convex Optimization for Big Data: Scalable, randomized, and parallel algorithms for big data analytics, IEEE Signal Process. Mag. 31(5) (2014) 32-43, DOI: 10.1109/MSP.2014.2329397.
  * [3] P. J. Carrington, J. Scott, S. Wasserman, eds., Models and Methods in Social Network Analysis, Structural Analysis in the Social Sciences, Cambridge University Press (2005), https://doi.org/10.1017/CBO9780511811395.
  * [4] K. Marti, Stochastic optimization methods, Springer, Heidelberg, third ed. (2015), Applications in engineering and operations research, https://doi.org/10.1007/978-3-662-46214-0.
  * [5] L. Martinez, R. Andrade, E.G. Birgin, J.M. Martinez, Packmol: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem. 30 (2009) 2157-2164, https://doi.org/10.1002/jcc.21224.
  * [6] D. Vicari, A. Okada, G. Ragozini, C. Weihs, eds., Analysis and Modeling of Complex Data in Behavioral and Social Sciences, Springer, Cham (2014), https://doi.org/10.1007/978-3-319-06692-9.
  * [7] A. Shapiro, D. Dentcheva, A. Ruszczynski, Lectures on stochastic programming: modeling and theory, Society for Industrial and Applied Mathematics (2021), https://doi.org/10.1137/1.9781611976595.
  * [8] F. Bastin, C. Cirillo, P.L. Toint, An adaptive Monte Carlo algorithm for computing mixed logit estimators, Comput. Manag. Sci. 3(1) (2006) 55-79, https://doi.org/10.1007/s10287-005-0044-y.
  * [9] S. Bellavia, N. Krejić, N. Krklec Jerinkić, Subsampled Inexact Newton methods for minimizing large sums of convex function, IMA J. Numer. Anal. 40(4) (2018) 2309-2341, https://doi.org/10.1093/imanum/drz027.
  * [10] T. Homem-de-Mello, Variable-Sample Methods for Stochastic Optimization, ACM Trans. Model. Comput. Simul. 13(2) (2003) 108–133, https://doi.org/10.1145/858481.858483.
  * [11] N. Krejić, N. Krklec, Line search methods with variable sample size for unconstrained optimization, J. Comput. Appl. Math. 245 (2013) 213-231, https://doi.org/10.1016/j.cam.2012.12.020.
  * [12] N. Krklec Jerinkić, A. Rožnjik, Penalty variable sample size method for solving optimization problems with equality constraints in a form of mathematical expectation, Numer. Algorithms 83(2) (2020) 701-718, https://doi.org/10.1007/s11075-019-00699-6.
  * [13] F.E. Curtis, X. Que,  An adaptive gradient sampling algorithm for nonsmooth optimization, Optim. Methods Softw. 28(6) (2013) 1302-1324, DOI: 10.1080/10556788.2012.714781.
  * [14] K.C. Kiwiel, Convergence of the gradient sampling algorithm for nonsmooth nonconvex optimization, SIAM J. Optim. 18(2) (2007) 379-388, https://doi.org/10.1137/050639673.
  * [15] L. Lukšan, J. Vlček, A bundle-Newton method for nonsmooth unconstrained minimization, Math. Program. 83(1) (1998) 373 – 391, https://doi.org/10.1007/BF02680566.
  * [16] C. Lemarechal, C. Sagastizabal, Variable metric bundle methods: From conceptual to implementable forms, Math. Program. 76(3) (1997) 393 – 410, https://doi.org/10.1007/BF02614390.
  * [17] A. Asl, M.L. Overton, Analysis of limited-memory BFGS on a class of nonsmooth convex functions, IMA J. Numer. Anal. 41(1) (2021) 1-27, https://doi.org/10.1093/imanum/drz052.
  * [18] A. Asl, M.L. Overton, Analysis of the gradient method with an Armijo–Wolfe line search on a class of nonsmooth convex functions, Optim. Methods Softw. 35(2) (2020) 223-242, https://doi.org/10.1080/10556788.2019.1673388.
  * [19] A. Jalilzadeh, A. Nedić, U. V. Shanbhag, F. Yousefian, A Variable Sample-Size Stochastic Quasi-Newton Method for Smooth and Nonsmooth Stochastic Convex Optimization, Proc. IEEE Conf. Decis. Control, Miami Beach, FL, (2018) 4097-4102, doi: 10.1109/CDC.2018.8619209.
  * [20] J.M. Martinez, E.A. Pilotta, Inexact restoration algorithms for constrained optimization, J. Optim. Theory Appl. 104 (2000) 135-163, https://doi.org/10.1023/A:1004632923654.
  * [21] S. Bellavia, N. Krejić, B. Morini, Inexact restoration with subsampled trust-region methods for finite-sum minimization, Comput. Optim. Appl. 76 (2020) 701–736, https://doi.org/10.1007/s10589-020-00196-w
  * [22] N. Krejić, J.M. Martinez, Inexact Restoration approach for minimization with inexact evaluation of the objective function, Math. Comput. 85 (2016) 1775-1791, https://doi.org/10.1090/mcom/3025.
  * [23] A. Fischer, A. Friedlander,  A new line search inexact restoration approach for nonlinear programming, Comput. Optim. Appl. 46(2) (2010) 333-346, https://doi.org/10.1007/s10589-009-9267-0
  * [24] E. Birgin, N. Krejić, J.M. Martinez, Iteration and evaluation complexity on the minimization of functions whose computation is intrinsically inexact, Math. Comput. 89 (2020) 253-278, https://doi.org/10.1090/mcom/3445.
  * [25] L.F. Bueno, J.M. Martinez, Inexact Restoration for Minimization with Inexact Evaluation both of the Objective Function and the Constraints, arXiv preprint arXiv:2201.01162, (2022).
  * [26] J.C. Duchi, E. Hazan, Y. Singer, Adaptive subgradient methods for online learning and stochastic optimization, J. Mach. Learn. Res. 12 (2011) 2121-2159.
  * [27] A.P. George, W.B. Powell, Adaptive stepsizes for recursive estimation with applications in approximate dynamic programming, Mach. Learn. 65 (2006) 167-198, https://doi.org/10.1007/s10994-006-8365-9.
  * [28] P. Hennig, Fast probabilistic optimization from noisy gradients, Proc. 30th ICML (2013) 62-70.
  * [29] D.P. Kingma, J. Ba, Adam: A method for stochastic optimization, Pro. ICLR (2015).
  * [30] C. Paquette, K. Scheinberg, A stochastic line search method with expected complexity analysis, SIAM J. Optim. 30(1) (2020) 349-376, https://doi.org/10.1137/18M1216250.
  * [31] T. Schaul, S. Zhang, Y. LeCun, No more pesky learning rates, PMLR (2013) 343-351.
  * [32] D. di Serafino, N. Krejić, N. Krklec Jerinkić, M. Viola, LSOS: Line-search Second-Order Stochastic optimization methods for nonconvex finite sums, Math. Comput. (to appear), arXiv:2007.15966v2 (2021).
  * [33] J. Yu, S. Vishwanathan, S. Guenter, N. Schraudolph, A Quasi-Newton Approach to Nonsmooth Convex Optimization Problems in Machine Learning, J. Mach. Learn. Res. 11 (2010) 1145-1200.
  * [34] L.F. Bueno, J.M. Martinez, On the complexity of an inexact restoration method for constrained optimization, SIAM J. Optim. 30(1) (2020) 80-101, https://doi.org/10.1137/18M1216146
  * [35] J. Burke, A. Lewis, M. Overton,  Approximating subdifferentials by random sampling of gradients, Math. Oper. Res. 27(3) (2002) 567-584, https://doi.org/10.1287/moor.27.3.567.317.
  * [36] D. Yan, H. Mukai, Optimization Algorithm with Probabilistic Estimation, J. Optim. Theory Appl. 64, 79(2) (1993) 345-371, https://doi.org/10.1007/BF00940585.
  * [37] Y. Wardi, Stochastic algorithms with Armijo stepsizes for minimization of functions, J. Optim. Theory Appl. 64 (1990) 399–417, https://doi.org/10.1007/BF00939456.
  * [38] G.N. Grapiglia, E.W. Sachs, On the worst-case evaluation complexity of nonmonotone line search algorithms, Comput. Optim. Appl. 68(3) (2017) 555-577, DOI: 10.1007/s10589-017-9928-3.
  * [39] LIBSVM Data: Classification (Binary Class), https://www.csie.ntu.edu.tw/$\sim$cjlin/libsvmtools/datasets/binary.html
  * [40] M. Lichman, UCI machine learning repository (2013), https://archive.ics. uci.edu/ml/index.php
  * [41] Y. LeCun, C. Cortes, C.J.C. Burges, The MNIST database of handwritten digits (1998), http://yann.lecun.com/exdb/mnist/
  * [42] N. Krejić, N. Krklec Jerinkić, S. Rapajić, Barzilai-Borwein method with variable sample size for stochastic linear complementarity problems, Optimization 65(2) (2016) 479-499, https://doi.org/10.1080/02331934.2015.1062008.
  * [43] X. Chen, C. Zhang and M. Fukushima, Robust solution of monotone stochastic linear complementarity problems, Math. Program. 117(1) (2009) 51-80, https://doi.org/10.1007/s10107-007-0163-z
  * [44] X. Li, H. Liu, X. Sun, Feasible smooth method based on Barzilai-Borwein method for stochastic linear complementarity problem, Numer. Algorithms 57 (2011) 207-215, https://doi.org/10.1007/s11075-010-9424-7