\section{Introduction.}
Outstanding growth in the use of computers and sensors has
attracted interest in myriad scientific and engineering fields
for solving difficult optimization problems involving functions available only through a zeroth-order oracle (i.e., the functions are {\it black boxes}~\cite{AuHa2017}).
Derivative-free optimization (DFO~\cite{AuHa2017,CoScVibook,LMW2019AN}) addresses such situations where closed-form expressions and derivatives are not available.
This paper focuses on unconstrained stochastic DFO, wherein the objective function values are accessible only through a blackbox oracle corrupted by stochastic noise. We consider the problem
\begin{equation}\label{probl1}
	\underset{\bm{x}\in\mathbb{R}^n}{\min} f(\bm{x}), \qquad\mbox{with}\quad f(\bm{x})=\mathbb{E}_{\ubar{\theta}}\left[f_{\ubar{\theta}}(\bm{x})\right],
\end{equation}
where the values of the continuously differentiable function $f:\mathbb{R}^n\to\R$ are available only via $f_{\ubar{\theta}}$, a stochastically noisy version of $f$, and where $\ubar{\theta}$ is a random variable whose distribution is possibly unknown.
Stochastic gradient descent (SGD~\cite{RoMo1951}) is arguably the most used algorithm to solve~\eqref{probl1}. However, SGD is very sensitive to the choice of its learning rate; and not only does the method use step directions that are not necessarily descent ones~\cite{paquette2018stochastic}, but also it exhibits slow convergence~\cite{chen2018stochastic}. Similar remarks have motivated recent work in stochastic DFO, leading to many algorithms such as model-based methods~\cite{blanchet2016convergence,chen2018stochastic,shashaani2018astro}, those using a direct-search approach~\cite{anagnostidis2021direct,audet2019stomads,Ch2012,dzahini2020expected,dzahini2020constrained,rinaldi2022weak}, and  others~\cite{berahas2019global,RBSWlong21,miaoSchNeur2021,paquette2018stochastic} that employ stochastic estimation of the gradient $\nabla f(\bm{x})$. Because of the unavailability of true function values and the resulting need to use stochastic estimates and/or models, the performance of these methods is highly dependent on how accurate those stochastic quantities are. For example,
STORM~\cite{blanchet2016convergence,chen2018stochastic} is a stochastic trust-region algorithm using so-called $\beta$-probabilistically $\varepsilon_f$-accurate estimates of unknown function values as well as $\alpha$-probabilistically $\kappa$-fully linear models of the objective function. For convergence purposes, such estimates and models whose accuracies are dynamically controlled by means of a trust-region radius $\delta_k$ need to be sufficiently accurate probabilistically. More precisely, at a given iteration~$k$, the aforementioned estimates and models require, respectively, at least
$\mathcal{O}\left(\delta_k^{-4}/(1-\sqrt{\beta})\right)$
and
$\mathcal{O}\left((n+1)\max\accolade{\delta_k^{-2},\delta_k^{-4}}/(1-\alpha^{1/(n+1)})\right)$
function evaluations in a simple stochastic noise framework where no noisy gradient values are available; these quantities grow rapidly as $\alpha,\beta\in (0,1)$ or the dimension $n$ gets larger or when $\delta_k$ tends to zero. On the other hand, ASTRO-DF~\cite{shashaani2018astro} is a class of stochastic DFO trust-region algorithms where random estimates and polynomial interpolation models are constructed adaptively by using a Monte Carlo sampling whose extent is determined by continuously balancing and monitoring measures of sampling error and model bias. Recalling that the construction of {\it full space} quadratic interpolation models requires $q(n):=(n+1)(n+2)/2$ ($\ell(n):=n+1$ in the linear case) points~\cite{CoScVibook} and that the accuracy of Monte Carlo-based estimates crucially depends on the sampling size, as suggested by the strong law of large numbers~\cite[Theorem~2.1.8]{tao2012topicsInRMT}, ASTRO-DF also does not have low per-iteration costs in terms of function evaluations.

In light of these observations,
a key question naturally arises regarding model-based stochastic DFO methods: is there a way to improve the ability of these methods to handle large-scale problems.
To answer a question similar to the one above in a context where the objective function is deterministic,  Cartis and  Roberts  \cite{CRsubspace2021} recently introduced RSDFO, a general framework of scalable subspace methods for model-based DFO. To achieve scalability, RSDFO approximates the deterministic objective only in subspaces using so-called $\Qk$-fully linear models, with~$\Qk\in \mathbb{R}^{n\times p}$ denoting a matrix whose entries are randomly selected and where $p\leq n$ (ideally, $p\ll n$) is a user-determined parameter. The RSDFO framework was then specialized to deterministic nonlinear least-squares problems, and high-probability worst-case complexity bounds of the methods were derived under mild assumptions. As highlighted in~\cite{CRsubspace2021}, another model-based subspace DFO method with similarities to RSDFO but admitting no convergence analysis is the {\it moving ridge function} approach~\cite{GrossParks2021}, where an interpolation model is built in an {\it active space} that is determined by using existing objective function evaluations. Neumaier  et al.~\cite{NeFeSchLei2011} also proposed the VXQR method with no convergence analysis, where line searches are performed along directions selected in a subspace determined by previous iterates. More recently, a direct-search method based on probabilistic descent in {\it reduced spaces} was introduced in~\cite{RobRoy2022RedSpace} for the optimization of deterministic objective functions, where the polling directions are constructed in random subspaces; complexity bounds were also derived, making use of probabilistic properties related to both the latter directions and subspaces.

Of the cited methods, the only ones that achieve scalability using random subspace strategies are developed for deterministic objective functions; to our knowledge, no such method exists for stochastic DFO. We address this gap by introducing a stochastic trust-region algorithm in random subspaces (STARS), in which scalability arises from constructing and then minimizing stochastic models that approximate the objective function in low-dimensional random subspaces. By defining these random subspaces through user-provided random matrices $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}\in\mathbb{R}^{n\times p}$ with $p\ll n$, only $q(p)$ and $\ell(p)$ points are required for the construction of quadratic and linear models, respectively, which are cheap models in terms of function evaluations since $p$ can always be chosen independently of $n$ in STARS.  While both the convergence and expected complexity analyses of STARS are inspired by those of STORM~\cite{blanchet2016convergence,chen2018stochastic}, dealing with the additional difficulties introduced by the randomness stemming from the entries of $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ is not a trivial task. As an example, no prior work exists showing how so-called $\beta_m$-{\it probabilistically} $(\kappa_{ef},\kappa_{eg};\Qk)$-{\it fully linear models} in random subspaces, crucial for the present analysis, can be made available. Moreover, the theoretical analyses of STORM utilizing random steps $\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}_k}\in\mathbb{R}^n$ do not straightforwardly hold when the latter are replaced with  $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}_k}\in\mathbb{R}^n$ used by STARS. Our main results trivially imply those related to STORM in the full-space case where $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}=\bm{I_n}\in\mathbb{R}^{n\times n}$ for all $k\in\mathbb{N}$, with $\bm{I_n}$
denoting the identity matrix.

One of the key contributions of the present work is that it extends the analysis of a STORM-like framework~\cite{audet2019stomads,BerDioKunRoy2022LvnbrgMrqdt,blanchet2016convergence,chen2018stochastic,curtis2020adaptive,dzahini2020expected,dzahini2020constrained,paquette2018stochastic} to settings where additional randomness stems from internal mechanics (here, the entries of the random matrices~$\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$) of an algorithm for stochastic DFO.
To our knowledge, STARS is the first DFO algorithm for stochastic objectives that achieves scalability using a random subspace strategy; STARS both extends STORM techniques to settings of scalable subspace methods and extends RSDFO techniques to stochastic objective functions.
Another key contribution of this work is the results of Theorem~\ref{subspaceModel1} together with Corollary~\ref{CorQkFully}, which are, to the best of our knowledge, the first to rigorously provide a detailed strategy for the construction of $\beta_m$-probabilistically $(\kappa_{ef},\kappa_{eg};\Qk)$-fully linear models
for stochastic objective functions in random subspaces.
Moreover, relying on a framework introduced in~\cite{blanchet2016convergence} using results derived from martingale theory, the analysis of STARS demonstrates that while using random subspace models, the expected complexity of the algorithm is surprisingly similar to that of STORM~\cite{blanchet2016convergence}, stochastic direct-search~\cite{dzahini2020expected}, and line-search-based~\cite{paquette2018stochastic} methods up to constants.

The manuscript is organized as follows. Section~\ref{Section2} introduces the general framework of STARS and explains how it results in a stochastic process. Section~\ref{Section3} discusses various strategies for the selection of random subspaces, demonstrates how probabilistic estimates and models can be constructed in these subspaces, and provides conditions related to these random quantities that are necessary for the convergence of STARS. Sections~\ref{Section4} and~\ref{Section5},  respectively, present the convergence and expected complexity analyses of STARS. Section~\ref{Section6} presents numerical results, followed by a discussion and suggestions for future work.

\section{Random subspace trust-region and resulting stochastic process.}\label{Section2}
$ $
This section presents the general framework of STARS and explains how the proposed method results in a stochastic process.
\subsection{The random subspace stochastic trust-region method.}
Unlike the stochastic trust-region framework STORM~\cite{blanchet2016convergence,chen2018stochastic}, which builds full space random models $\ubar{m}_k(\bm{x})$, $\bm{x}\in\mathbb{R}^n$, to approximate $f(\bm{x})$, STARS operates as follows. On iteration $k$, given a current iterate $\bm{{x}_k}$, let $\mathcal{Y}_k\subseteq \mathbb{R}^n$ be the affine space randomly chosen through the range of a matrix $\Qk\in\mathbb{R}^{n\times p}$
whose entries are randomly selected; that is,
\[\mathcal{Y}_k = \accolade{\bm{{x}_k}+\Qk\bm{{\bm{s}}}:\bm{{\bm{s}}}\in\R^p}.\]
Then, given a trust-region radius $\delta_k>0$, a subspace model $\mhatk(\bm{s}), \bm{s}\in\mathbb{R}^p$ is built only on $\mathcal{Y}_k$ using realizations of the stochastically noisy function $f_{\ubar{\theta}}$; the model  serves as an approximation of $f(\bm{x})$ in $\mathcal{B}(\bm{{x}_k},\delta_k;\Qk):=\accolade{\bm{{x}_k}+\Qk\bm{s} \in \mathcal{Y}_k:\normii{\bm{s}}\leq\delta_k}$. For concreteness, here we employ a quadratic subspace model given by
\begin{equation*}
f(\bm{{x}_k}+\Qk\bm{s})\approx\mhatk(\bm{\bm{s}}):= f_k+\bm{\hat{g}_k}^\top\bm{s}+\bm{s}^\top\bm{\hat{H}}_k\bm{s},
\end{equation*}
where $\bm{\hat{g}_k}\in \R^p$ and $\bm{\hat{H}}_k\in \Sym^{p}$ are the low-dimensional model gradient and Hessian, respectively. A {\it tentative step} $\Qk\bm{s_k}\in\mathbb{R}^n$ is produced by using the solution $\bm{s_k}\in\mathbb{R}^p$ obtained by approximately minimizing $\mhatk$ inside the trust  region; that  is, $\bm{s_k}\approx{\argmin}\left\{ \mhatk(\bm{s}): \bm{s}\in\mathbb{R}^p, \, \normii{\bm{{\bm{s}}}}\leq \delta_k\right\}$. Inspired by~\cite{CRsubspace2021,chen2018stochastic}, the {\it trial step} $\bm{s_k}$ has to provide a sufficient decrease in $\mhatk$ by satisfying the following standard fraction of the Cauchy decrease condition.

\begin{assumption}
For every iteration $k$, a trial step $\bm{s_k}\in\mathbb{R}^p$ is computed so that\footnote{Throughout the manuscript, the matrix norm $\norme{\cdot}$ is supposed to be {\it consistent} with the Euclidean norm; that is, $\normii{\bm{A}\bm{x}}\leq \norme{\bm{A}}\normii{\bm{x}}$, which holds for both Frobenius and spectral matrix norms.}
\begin{equation}\label{CauchyDecrease}
\mhatk(\bm{0})-\mhatk(\bm{\bm{s_k}})\geq \frac{\kappa_{fcd}}{2}\normii{\bm{\hat{g}_k}}\min\accolade{\delta_k,\frac{\normii{\bm{\hat{g}_k}}}{\max\left(\norme{\bm{\hat{H}}_k},1\right)}},
\end{equation}
for some constant $\kappa_{fcd}\in (0,1]$.
\end{assumption}
Estimates $f_k^0$ and $f_k^s$ of $f(\bm{{x}_k})$ and $f(\bm{{x}_k}+\Qk\bm{s_k})$, respectively, are constructed by using evaluations of the noisy function $f_{\ubar{\theta}}$. The possible change in $f$ given by $\bm{s_k}$ is measured by comparing the estimates $f_k^0$ and $f_k^s$ through the value of a ratio $\rho_k$. An iteration is called successful if a decrease in the estimates is deemed sufficient, in which case the current solution is updated by $\bm{{x}_k}+\Qk\bm{s_k}$ and the trust-region radius is not decreased. Otherwise, the latter is decreased while the former is not updated. Because of inaccurate estimates, successful iterations qualified as {\it false} could lead to an increase in $f$. While the present algorithmic framework is not able to identify such iterations, as was the case in~\cite{audet2019stomads,chen2018stochastic,dzahini2020expected,dzahini2020constrained,paquette2018stochastic}, it will  be possible to prove later in the proof of Theorem~\ref{liminfTheorem} combined with~\eqref{randomWalk} that {\it true} iterations occur sufficiently often for convergence of the proposed method to hold.

\subsection{Stochastic process generated by the algorithm.}

The random variables considered here are all defined on the same probability space $(\Omega,\mathcal{F}, \mathbb{P})$, where $\Omega$, referred to as the {\it sample space}, is a nonempty set whose elements $\omega$ and subsets are called {\it sample points} and {\it events}, respectively. $\mathcal{F}$ is a collection of such events, which is called a {\it $\sigma$-algebra}. $\mathbb{P}$ is a finite measure defined on the measurable space $(\Omega,\mathcal{F})$, satisfying $\prob{\Omega}=1$ and referred to as the {\it probability measure}. An increasing subsequence $\accolade{\mathcal{S}_k}_k$ of $\sigma$-algebras of $\mathcal{F}$ will be called a {\it filtration}. Let $\mathscr{B}(\mathbb{R}^n)$ be the $\sigma$-algebra generated by the open sets of $\mathbb{R}^n$, also known as the Borel $\sigma$-algebra of $\mathbb{R}^n$. A random variable $\ubar{z}$ is a measurable map defined on $(\Omega,\mathcal{F}, \mathbb{P})$ into the measurable space $(\mathbb{R}^n, \mathscr{B}(\mathbb{R}^n))$, where measurability means that $\accolade{\ubar{z}\in A}:=\accolade{\omega\in\Omega:\ubar{z}(\omega)\in A} =:{\ubar{z}}^{-1}(A)\in\mathcal{F}\ \forall A\in \mathscr{B}(\mathbb{R}^n)$~\cite{bhattacharya2007basic}. For the remainder of the manuscript, vectors will be written in lowercase boldface (e.g., $\bm{x}\in\mathbb{R}^n, n\geq 2$) while matrices will be written in uppercase boldface (e.g., $\Qk\in\mathbb{R}^{n\times p}$), and underlined letters (e.g., $\ubar{z}, \ubar{\bm{x}}, \bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$) will be used to denote random quantities; $z=\ubar{z}(\omega)$ will denote a realization of $\ubar{z}$.

\begin{algorithm}[htb]
\caption{STARS}
\label{algoStoScalTR}
\textbf{[0] Initialization}\\
\hspace*{10mm}Choose constants $\gamma>1$, $\eta_1\in(0,1)$, $\eta_2>0$, $j_{\max}\in\mathbb{N}$, $\varepsilon\in(0,1)$, $\beta\in(0,1)$, \\
\hspace*{10mm}$c_1\geq 1$, initial trust-region radius $\delta_0>0$, and maximum trust-region radius\\
\hspace*{10mm}$\delta_{\max}=\gamma^{j_{\max}}\delta_0$, starting point $\bm{{\bm{x}}_0}\in\mathbb{R}^n$, and dimension $p\in\intbracket{1,n}$.\\
\hspace*{10mm}Set the iteration counter $k \gets 0$.\\
\textbf{[1] Construction of subspace model}\\
\hspace*{10mm}Generate $\Qk$: a realization of a\footnotemark {\em WAM}($1-\varepsilon,\beta$), using a distribution $\mathbb{Q}_{k}$.\\
\hspace*{10mm}Build model $\mhatk:\mathbb{R}^p\to\R$ that is $(\kef,\keg;\Qk)$-fully linear
in  $\mathcal{B}(\bm{{x}_k},c_1\delta_k;\Qk)$.\\
\textbf{[2] Step calculation}\\
\hspace*{10mm}Compute $\bm{s_k}\approx {\argmin}\left\{ \mhatk(\bm{s}) : \; \bm{s}\in\mathbb{R}^p:\normii{\bm{s}}\leq \delta_k\right \}$ satisfying~\eqref{CauchyDecrease}.\\
\textbf{[3] Estimate computation}\\
\hspace*{10mm}Obtain estimates $f_k^0\approx f(\bm{{x}_k})$ and $f_k^s\approx f(\bm{{x}_k}+\Qk\bm{s_k})$ satisfying~\eqref{accEstimates}.\\
\textbf{[4] Updates}\\
\hspace*{10mm}Compute $\rho_k=\frac{f_k^0-f_k^s}{\mhatk(\bm{0})-\mhatk(\bm{s_k})}.$\\
\hspace*{10mm}If $\rho_k\geq \eta_1$ and $\normii{\bm{\hat{g}_k}}\geq \eta_2\delta_k$ (\textbf{success}): \\
\hspace*{16mm}set $\bm{\bm{{x}_{k+1}}}=\bm{{x}_k}+\Qk\bm{s_k}$ and ${\delta_{k+1}}=\min\accolade{\gamma\delta_k,\delta_{\max}}$. \\
\hspace*{10mm}Otherwise (\textbf{failure}): set $\bm{\bm{{x}_{k+1}}}=\bm{{x}_k}$ and ${\delta_{k+1}}=\gamma^{-1}\delta_k$.\\
\hspace*{10mm}Update the iteration counter $k \gets k+1$, and go to~\textbf{[1]}.\\
\end{algorithm}
\footnotetext{See Definition~\ref{wellAlignedDefinition} for the meaning of {\em WAM}($1-\varepsilon,\beta$).}

The deterministic models $\mhatk$ and the function estimates $f_k^0$ and $f_k^s$ are constructed by using evaluations of the noisy function $f_{\ubar{\theta}}$ and the randomly selected entries of the deterministic matrix~$\Qk=\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}(\omega)$.  Hence, $\mhatk$, $f_k^0$, and $f_k^s$ can be considered as realizations of random models and estimates $\Mhatk$, ${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0$, and ${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s$, respectively. Each iteration of Algorithm~\ref{algoStoScalTR} is therefore influenced by the behavior of these random quantities; and consequently, the algorithm results in a stochastic process. The present work shows that under certain conditions on the sequences~$\accoladekinN{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}}$, $\accoladekinN{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0,{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s}$, and $\accoladekinN{\Mhatk}$, the resulting stochastic process has desirable convergence properties, conditioned on the past, where {\it past} means the {\it past history} of the algorithm (see~Remark~\ref{whyTheseSigmaAlgebras}). In particular, for the needs of convergence and expected complexity analyses presented  in Sections~\ref{Section4} and~\ref{Section5}, and as will be seen  in Section~\ref{Section3}, the random models~$\Mhatk$ and estimates~${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0,{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s$ are required to be probabilistically sufficiently accurate conditioned on the past, as was the case in~\cite{chen2018stochastic}. Moreover, $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ must guarantee a high quality of the selected subspace probabilistically and conditioned on the past, as will be detailed in Section~\ref{Section31}.

For all $k\in\mathbb{N}$, denote by $(\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k})_{ij},  (i,j)\in \intbracket{1,n}\times\intbracket{1,p}=:\mathbb{I}_{n,p}$ the entries of the random matrix $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$, where for $a\leq b$, $\intbracket{a,b}:=[a,b]\cap\mathbb{Z}$  throughout the manuscript.
To formalize conditioning on the past, we consider the three filtrations $\accolade{\mathcal{F}_{k-1}}_{k\in\mathbb{N}}$, $\accolade{\mathcal{F}_{k-1}^Q}_{k\in\mathbb{N}}$, and $\accolade{\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}_{k\in\mathbb{N}}$ defined respectively by
\begin{equation*}
\begin{array}{rcl}
\mathcal{F}_{k-1} &:=& \sigma\left({\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_{\ell}^0, {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_{\ell}^s, \hat{{{\stackunder[1.1pt]{$m$}{\rule{1.1ex}{.09ex}}}}}_{\ell},  (\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_{\ell}})_{ij} \mbox{ with}\  \ell \in \intbracket{0,k-1} \mbox{ \&}\ (i,j)\in \mathbb{I}_{n,p} \right)\\
\mathcal{F}_{k-1}^Q &:=& \sigma\left({\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_{\ell}^0, {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_{\ell}^s, \hat{{{\stackunder[1.1pt]{$m$}{\rule{1.1ex}{.09ex}}}}}_{\ell}, (\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_{\ell}})_{ij}, (\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k})_{ij} \mbox{ with}\  \ell \in \intbracket{0,k-1} \mbox{ \&}\ (i,j)\in \mathbb{I}_{n,p} \right)\\
\mathcal{F}^{\hat{m}\cdot Q}_{k-1} &:=& \sigma\left({\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_{\ell}^0, {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_{\ell}^s, \hat{{{\stackunder[1.1pt]{$m$}{\rule{1.1ex}{.09ex}}}}}_{\ell}, (\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_{\ell}})_{ij}, (\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k})_{ij}, \ \Mhatk \mbox{ with}\  \ell \in \intbracket{0,k-1} \mbox{ \&}\ (i,j)\in \mathbb{I}_{n,p}  \right).
\end{array}
\end{equation*}
Here, for {\it completeness}~\cite{billingsley1995ThirdEdition}, $\mathcal{F}_{-1}=\sigma(\bm{\bm{x}_0})$.
By construction, $\mathcal{F}_{k-1}\subset\mathcal{F}_{k-1}^Q\subset\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$, $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ is both~$\mathcal{F}_{k-1}^Q$ and~$\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$-measurable, and $\Mhatk$ is~$\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$-measurable, which imply that $\E{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}|\mathcal{F}_{k-1}^Q}=\E{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}=\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ and $\E{\Mhatk|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}=\Mhatk$.

\begin{remark}\label{whyTheseSigmaAlgebras}
At a given iteration~$k$, when generating~$\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$, because no model or estimates were generated since the start of the iteration, it obviously follows from the construction of Algorithm~\ref{algoStoScalTR} that its past history can be formalized by~$\mathcal{F}_{k-1}$. On the other hand, when constructing~$\Mhatk$, the past of the algorithm that includes~$\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$, given that the latter is already generated since the beginning of iteration~$k$, can therefore be formalized by~$\mathcal{F}_{k-1}^Q$. A similar observation  explains why~$\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$ can formalize the past history of Algorithm~\ref{algoStoScalTR} before the construction of~${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0$ and~${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s$. These observations will also play an important role in the formalization of sufficient accuracy of the random models and estimates in Definitions~\ref{defProbQkfully} and~\ref{defProbAccEstim} and the formalization of high quality of the random subspaces in Assumption~\ref{wellAlignedAssumption}.
\end{remark}

\section{Subspace selection and probabilistic models and estimates in random subspaces.}\label{Section3}
The random subspaces in which the models $\Mhatk$ are built are defined by specific random matrices. For efficiency of the proposed method and for convergence needs, not only must these subspaces be probabilistically rich enough, but also the accuracies of the models and estimates of unknown function values must hold with sufficiently high, but fixed, probabilities. This section discusses strategies for such random subspace selection, demonstrates by means of rigorous results how probabilistic estimates and subspace models can be constructed, and provides conditions that are necessary for the convergence of STARS.

\subsection{Subspace selection strategies.}\label{Section31}
To choose a subspace,
there must be enough gradient living in the subspace to ensure analogous reduction of $f$. In other words, the matrix $\Qk\in\mathbb{R}^{n\times p}$ determining the subspace needs to satisfy (probabilistically) $\normii{\Qk^\top\nabla f(\bm{{x}_k})}\geq \alpha_{Q}\normii{\nabla f(\bm{{x}_k})}$ at every iteration, for some constant $\alpha_{Q}\in (0,1)$ independent of $k$. This requirement motivates the following modified definition suggested by~\cite{CRsubspace2021,shao2022Thesis}, which basically says that at least some fraction of the gradient of~$f$ must be maintained after projecting it into the selected subspace.

\begin{definition}\label{wellAlignedDeterministic}
For $\alpha_{Q}\in (0,1)$, a matrix $\Q$ is $\alpha_{Q}$-well aligned if, for any vector $\bm{{\bm{v}}}\in\mathbb{R}^n$, $\normii{\Q^\top\bm{{\bm{v}}}}\geq \alpha_{Q}\normii{\bm{{\bm{v}}}}$.
\end{definition}

Recalling Remark~\ref{whyTheseSigmaAlgebras} and inspired by~\cite{CRsubspace2021}, we propose a probabilistic quantification of the quality of the subspace selection as follows.

\begin{definition}\label{wellAlignedDefinition}
For fixed $\alpha_Q,\beta_Q\in (0,1)$, a sequence $\accoladekinN{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}}$ of random matrices
is $(1-\beta_Q)$-probabilistically $\alpha_{Q}$-well aligned if, for any $\mathcal{F}_{k-1}$-measurable random vector $\bm{{\stackunder[1.0pt]{$v$}{\rule{0.8ex}{.075ex}}}}$ with realizations $\bm{v}\in \mathbb{R}^n$, the events
$\mathcal{A}_k:= \accolade{\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\bm{{\stackunder[1.0pt]{$v$}{\rule{0.8ex}{.075ex}}}}}\geq \alpha_{{Q}}\normii{\bm{{\stackunder[1.0pt]{$v$}{\rule{0.8ex}{.075ex}}}}}}$ satisfy
\begin{equation}\label{Akineq}
\prob{\mathcal{A}_k|\mathcal{F}_{k-1}}=\E{\mathds{1}_{\mathcal{A}_k}|\mathcal{F}_{k-1}}\geq 1-\beta_Q,
\end{equation}
where $\mathds{1}_{\mathcal{A}_k}$ denotes the indicator function of the event $\mathcal{A}_k$; that  is, $\mathds{1}_{\mathcal{A}_k}(\omega)=1$ if $\omega\in \mathcal{A}_k$ and $\mathds{1}_{\mathcal{A}_k}(\omega)=0$ otherwise. Random matrices satisfying~\eqref{Akineq} will be referred to as {WAM}($\alpha_Q,\beta_Q$).
\end{definition}
For convergence purposes  in Sections~\ref{Section4} and~\ref{Section5}, the following will be assumed throughout the manuscript.
\begin{assumption}\label{wellAlignedAssumption}
The sequence $\accoladekinN{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}}$ of matrices used by Algorithm~\ref{algoStoScalTR} is $(1-\beta_Q)$-probabilistically $\alpha_{Q}$-well aligned for some fixed $\alpha_Q\in (0,1)$ and $\beta_Q\in (0,1/2)$.
\end{assumption}

The next result from~\cite[Lemma~1]{LauBeaMaP2000} provides a slight generalization of an exponential inequality for chi-square distributions. It will be required later for the proof of Theorem~\ref{GaussianJLTmatrix}, inspired by~\cite[Theorem~2.1]{WoodrDavP2014}, providing a random {\it matrix ensemble} that will be shown in Corollary~\ref{JLCorollary1} to satisfy Assumption~\ref{wellAlignedAssumption}.
\begin{lemma}\label{BeatriceLaurent}
Let $({\stackunder[0.6pt]{$y$}{\rule{1.0ex}{.075ex}}}_1,{\stackunder[0.6pt]{$y$}{\rule{1.0ex}{.075ex}}}_2,\dots,{\stackunder[0.6pt]{$y$}{\rule{1.0ex}{.075ex}}}_p)$ be i.i.d.\ Gaussian random variables with mean zero and variance one. Let $a_1$, $a_2$, $\dots$, $a_p$ be nonnegative real numbers. Let $\ubar{z}$ be the random variable defined by $\ubar{z}=\sum_{i=1}^{p}a_i({\stackunder[0.6pt]{$y$}{\rule{1.0ex}{.075ex}}}_i^2-1)$. Then it holds that for any $t\geq 0$,
\[\prob{\ubar{z}\leq -2\normii{\bm{a}}\sqrt{t}}\leq e^{-t}.\]
\end{lemma}

\begin{theorem}\label{GaussianJLTmatrix}
Let $\varepsilon>0$, $\beta<1$, $p\geq 4\varepsilon^{-2}\log(1/\beta)$, and $\bm{{\stackunder[1.0pt]{$R$}{\rule{1.3ex}{.075ex}}}}=\left({{\stackunder[1.0pt]{$r$}{\rule{0.8ex}{.075ex}}}}_{ij}\right)_{\underset{1\leq j\leq n}{1\leq i\leq p}}\in \mathbb{R}^{p\times n}$ be a random matrix with i.i.d.\ standard Gaussian entries. Let $\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}=\left({\stackunder[1.0pt]{$s$}{\rule{0.8ex}{.075ex}}}_{ij}\right)_{\underset{1\leq j\leq n}{1\leq i\leq p}}\in \mathbb{R}^{p\times n}$ be a random matrix defined by $\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}=\frac{1}{\sqrt{p}}\bm{{\stackunder[1.0pt]{$R$}{\rule{1.3ex}{.075ex}}}}$; that is, ${\stackunder[1.0pt]{$s$}{\rule{0.8ex}{.075ex}}}_{ij}=\frac{1}{\sqrt{p}}{\stackunder[1.0pt]{$r$}{\rule{0.8ex}{.075ex}}}_{ij}\sim\mathcal{N}(0,1/p)$. Then for all $\bm{v}\in\mathbb{R}^n$,
$\mathbb{P}\left[\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}^2\geq (1-\varepsilon)\normii{\bm{v}}^2\right]\geq 1-\beta$.
\end{theorem}

\begin{proof}
The proof is inspired by that of~\cite[Lemma~2.12]{WoodrDavP2014}. Consider any deterministic vector\\ $\bm{v}=\left(v_1,v_2,\dots,v_n\right)^\top\in\mathbb{R}^n$,  and denote by $\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}}_{: j}\in\R^p, \, j \in\intbracket{1,n}$ the columns of the random matrix $\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}$. Then the $i$th component of the $p$-tuple $\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}=\sum_{j=1}^{n}v_j\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}}_{: j}$ is the Gaussian random variable $(\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v})_i=\sum_{j=1}^{n}v_j{{\stackunder[1.0pt]{$s$}{\rule{0.8ex}{.075ex}}}}_{i j}$ with mean zero and variance $\mathbb{V}\left[(\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v})_i\right]= \sum_{j=1}^{n}v_j^2\var{{{\stackunder[1.0pt]{$s$}{\rule{0.8ex}{.075ex}}}}_{i j}}=\frac{\normii{\bm{v}}^2}{p}$. It follows from the independence of $(\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v})_i, \, i\in\intbracket{1,p}$ that the random variable $\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}^2=\sum_{i=1}^p(\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v})_i^2$ is equal in distribution to $\frac{\normii{\bm{v}}^2}{p}{\ubar{w}}$, where ${\ubar{w}}$ is a $\chi^2_p$ random variable with $p$ degrees of freedom. Applying Lemma~\ref{BeatriceLaurent} with $t=\frac{\varepsilon^2p}{4}$ and $a_i=1,\ i\in\intbracket{1,p}$, yields $\prob{{\ubar{w}}-p\leq -\varepsilon p}\leq e^{-\frac{\varepsilon^2p}{4}}.$
Hence, $\mathbb{P}\left[\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}^2\geq (1-\varepsilon)\normii{\bm{v}}^2\right]=1-\mathbb{P}\left[{\ubar{w}}\leq(1-\varepsilon)p\right]\geq 1-e^{-\frac{\varepsilon^2p}{4}}$.
The proof is completed by noticing that $1-2e^{-\frac{\varepsilon^2p}{4}}\geq 1-\beta$ if $p\geq 4\varepsilon^{-2}\log(1/\beta)$.
\end{proof}

The next corollary shows that Assumption~\ref{wellAlignedAssumption} can be satisfied by using matrices resulting from Theorem~\ref{GaussianJLTmatrix}, and whose number~$p$ of columns does not depend on~$n$.
\begin{corollary}\label{JLCorollary1}
At a given iteration $k$, let the subspace selection matrix $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}=\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\in \mathbb{R}^{p\times n}$ be provided by Theorem~\ref{GaussianJLTmatrix}, with $\varepsilon=1-\alpha_Q$ and $\beta=\beta_Q$, for some $\alpha_Q\in (0,1)$ and $\beta_Q\in (0,1/2)$. Assume that $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ is independent of $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_0}, \bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_1},\dots, \bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_{k-1}}$ and ${\ubar{\theta}}$. Then Assumption~\ref{wellAlignedAssumption} holds; that is, for any $\bm{{\stackunder[1.0pt]{$v$}{\rule{0.8ex}{.075ex}}}}$ $\mathcal{F}_{k-1}$-measurable,
\begin{equation*}
\prob{\accolade{\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\bm{{\stackunder[1.0pt]{$v$}{\rule{0.8ex}{.075ex}}}}}\geq \alpha_{{Q}}\normii{\bm{{\stackunder[1.0pt]{$v$}{\rule{0.8ex}{.075ex}}}}}}|\mathcal{F}_{k-1}}\geq 1-\beta_Q.
\end{equation*}
\end{corollary}
\begin{proof} Since $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ is independent of $\mathcal{F}_{k-1}$ and $\bm{{\stackunder[1.0pt]{$v$}{\rule{0.8ex}{.075ex}}}}$ is $\mathcal{F}_{k-1}$-measurable, it suffices to show that\\ $\prob{\accolade{\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\bm{v}}\geq \alpha_{{Q}}\normii{\bm{v}}}}\geq 1-\beta_Q$
for any deterministic vector $\bm{v}\in\mathbb{R}^n$, which easily follows from Theorem~\ref{GaussianJLTmatrix} and the inclusion $\accolade{\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}^2\geq (1-\varepsilon)\normii{\bm{v}}^2} \subseteq \accolade{\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}\geq (1-\varepsilon)\normii{\bm{v}}}$ due to the inequality $1-\varepsilon\geq (1-\varepsilon)^2$.
\end{proof}

Inspired by~\cite{CRsubspace2021}, a technical assumption needed for convergence and expected complexity analyses later in Sections~\ref{Section4} and~\ref{Section5} is the following.
\begin{assumption}\label{QmaxAssumption}
Any realization $\Qk$ of the random matrix used by Algorithm~\ref{algoStoScalTR} satisfies $\norme{\Qk}\leq \Qmax$ for all $k\in\mathbb{N}$, for some constant $\Qmax >0$ independent of~$k$.
\end{assumption}
The next theorem, partially proved in the Appendix and formulated based on results from~\cite{KaNel2014SparseLidenstrauss}, provides another technique for constructing {\em WAM}($1-\varepsilon,\beta$) satisfying Assumption~\ref{wellAlignedAssumption} through so-called Johnson--Lindenstrauss (JL) transforms~\cite{JohnsonLind1984,KaNel2014SparseLidenstrauss}, and whose number~$p$ of columns does not depend on~$n$. Moreover, one can easily see that the resulting matrix ensemble satisfies Assumption~\ref{QmaxAssumption} with $\Qmax=\sqrt{n}$, unlike the one from Theorem~\ref{GaussianJLTmatrix} where $\norme{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}}\leq \Qmax$ holds with high probability with\footnote{The reader is referred to~\cite{knuth1976big} for details regarding the notation~$\Theta(\cdot)$.} $\Qmax=\Theta(\sqrt{n/p})$~\cite[Corollary~3.11]{bandeira2016sharp}.
\begin{theorem}\label{sHashingTheorem}
Let $\varepsilon>0$, $\beta<1/2$, and $\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}=\left({\stackunder[1.0pt]{$s$}{\rule{0.8ex}{.075ex}}}_{ij}\right)_{\underset{1\leq j\leq n}{1\leq i\leq p}}\in \mathbb{R}^{p\times n}$ be a real-valued random matrix defined by
${\stackunder[1.0pt]{$s$}{\rule{0.8ex}{.075ex}}}_{ij}=\frac{1}{\sqrt{r}}{\stackunder[1.0pt]{$\eta$}{\rule{0.8ex}{.075ex}}}_{ij}{\stackunder[1.0pt]{$\sigma$}{\rule{0.8ex}{.075ex}}}_{ij},$
where $r\in\mathbb{N}\setminus\accolade{0}$ and ${\stackunder[1.0pt]{$\sigma$}{\rule{0.8ex}{.075ex}}}_{ij}$ are independent Rademacher random variables. The ${\stackunder[1.0pt]{$\eta$}{\rule{0.8ex}{.075ex}}}_{ij}$ satisfying $\prob{\sum_{i=1}^p{\stackunder[1.0pt]{$\eta$}{\rule{0.8ex}{.075ex}}}_{ij}=r}=1$ for all $j\in\intbracket{1,n}$ are negatively correlated indicator random variables for the events $\accolade{{\stackunder[1.0pt]{$s$}{\rule{0.8ex}{.075ex}}}_{ij}\neq 0}$; that  is,  they satisfy
$\E{\underset{(i,j)\in T}{\prod}{\stackunder[1.0pt]{$\eta$}{\rule{0.8ex}{.075ex}}}_{ij}}\leq (r/p)^{\abs{T}}$
for all $T\subseteq\intbracket{1,p}\times\intbracket{1,n}$ with $\abs{T}\leq\ell$, where $\ell\geq \ln\left[\beta^{-1}(\ell+1)\ell/2\right]$. Then $\mathbb{P}\left[(1-\varepsilon)\normii{\bm{v}}\leq\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}\leq (1+\varepsilon)\normii{\bm{v}}\right] > 1-\beta$ for any $\bm{v}\in\mathbb{R}^n$, provided $r\geq 8e^4\sqrt{e}(\ell+1)/(2\varepsilon-\varepsilon^2)$ and $p=2r^2/(e\ell)$.
\end{theorem}

The next corollary shows that {\em WAM}($\alpha_Q,\beta_Q$) satisfying Assumptions~\ref{wellAlignedAssumption} and~\ref{QmaxAssumption} can be provided by specific JL transforms $\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}$ resulting from Theorem~\ref{sHashingTheorem}.
\begin{corollary
At a given iteration $k$, let the subspace selection matrix $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}=\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\in \mathbb{R}^{p\times n}$ be provided by Theorem~\ref{sHashingTheorem}, with $\varepsilon=1-\alpha_Q$ and $\beta=\beta_Q$. Assume that $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ is independent of $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_0}, \bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_1},\dots, \bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_{k-1}}$ and~${\ubar{\theta}}$. Then Assumption~\ref{wellAlignedAssumption} holds.
\end{corollary}
\begin{proof}
As was explained in the proof of Corollary~\ref{JLCorollary1}, it suffices to show that for all deterministic vectors $\bm{v}\in\mathbb{R}^n$, $\prob{\accolade{\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\bm{v}}\geq \alpha_{{Q}}\normii{\bm{v}}}}\geq 1-\beta_Q$. This immediately follows from the result of Theorem~\ref{sHashingTheorem}, together with the inclusion
$\accolade{(1-\varepsilon)\normii{\bm{v}}\leq\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}\leq (1+\varepsilon)\normii{\bm{v}}} \subseteq \accolade{\normii{\bm{{\stackunder[1.0pt]{$S$}{\rule{1.3ex}{.075ex}}}}\bm{v}}\geq (1-\varepsilon)\normii{\bm{v}}}.$
\end{proof}

\subsection{Probabilistic models and estimates in random subspaces.}
As mentioned at the beginning of Section~\ref{Section3}, the models used by Algorithm~\ref{algoStoScalTR} need to be sufficiently accurate. Motivated by~\cite{chen2018stochastic}, this sufficient accuracy is formalized for deterministic subspace models by the following measure of accuracy inspired by~\cite{CRsubspace2021}.
\begin{definition}\label{QkFullyLinearModels}
Assume $\nabla f$ is  Lipschitz continuous. Given $\Qk\in\mathbb{R}^{n\times p}$ and $\bm{{x}_k}\in\mathbb{R}^n$, a function $\mhatk:\mathbb{R}^p\to\R$ is a $(\kef,\keg; \Qk)$-fully linear model of $f$ in $\mathcal{B}(\bm{{x}_k},\delta_k;\Qk)$ for some $\delta_k>0$ if there exist constants $\kappa_{ef}, \kappa_{eg} >0$ independent of $k$ such that for all $\bm{s}\in\mathbb{R}^p$ with $\normii{\bm{s}}\leq \delta_k$,
\begin{equation}\label{Eq36}
\abs{f(\bm{{x}_k} +\Qk\bm{s})-\mhatk(\bm{s})}\leq \kappa_{ef}\delta_k^2\quad \mbox{and}\quad
\normii{\Qk^{\top}\nabla f(\bm{{x}_k} +\Qk\bm{s})-\nabla\mhatk(\bm{s})}\leq \kappa_{eg}\delta_k.
\end{equation}
\end{definition}
As pointed out in~\cite{CRsubspace2021}, the gradient condition in~\eqref{Eq36} is due to the fact that given any fixed~$\bm{x}$ and~$\Q$, if $\hat{f}(\bm{s}):=f(\bm{x}+\Q\bm{s})$, then $\nabla_{\bm{s}} \hat{f}(\bm{s})=\Q^{\top}\nabla_{\bm{x}} f(\bm{x}+\Q\bm{s})$. Moreover, denoting by $\bm{I}$ the identity matrix in full-dimensional subspaces where $p=n$, then $(\kef,\keg; \bm{I})$-fully linear models correspond to standard fully linear models~\cite[Definition~6.1]{CoScVibook}.

Recalling Remark~\ref{whyTheseSigmaAlgebras}, the next definition formalizing sufficient accuracy of probabilistic subspace models $\Mhatk$ is a stochastic variant of Definition~\ref{QkFullyLinearModels}, which generalizes~\cite[Definition~3.4]{chen2018stochastic} to low-dimensional subspaces. The existence of such models will be rigorously demonstrated  in Theorem~\ref{subspaceModel1} and Corollary~\ref{CorQkFully}.
\begin{definition}\label{defProbQkfully}
A sequence $\accolade{\Mhatk}$ of random models is $\beta_m$-probabilistically \\ $(\kef,\keg; \Qk)$-fully linear with respect to the random sequence $\accolade{\mathcal{B}(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k},{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k;\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k})}$ if the events
\[I_k^Q=\accolade{\mbox{Given}\ \Qk\ \mbox{and}\ \bm{{x}_k}, \Mhatk\ \mbox{is}\ (\kef,\keg;\Qk)\mbox{-fully linear for}\ f\ \mbox{in}\ \mathcal{B}(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k},{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k;\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k})}\]
satisfy the supermartingale-like condition
\begin{equation}\label{ikQProb}
\prob{I^Q_k|\mathcal{F}_{k-1}^Q}=\E{\mathds{1}_{I_k^Q}|\mathcal{F}_{k-1}^Q}\geq \beta_m.
\end{equation}
\end{definition}
Recalling that estimates of unknown function values also need to be sufficiently accurate, we introduce a formalization of sufficient accuracy inspired by~\cite{audet2019stomads,chen2018stochastic,dzahini2020expected,dzahini2020constrained,paquette2018stochastic} and motivated by~\cite{CRsubspace2021} as follows.
\begin{definition}
Given $\varepsilon_f>0$ and $\Qk\bm{s_k}\in\mathbb{R}^n, \, \bm{{x}_k}\in\mathbb{R}^n$, $f_k^0$ and $f_k^s$ are called $\varepsilon_f$-accurate estimates of $f(\bm{{x}_k})$ and $f(\bm{{x}_k}+\Qk\bm{s_k})$, respectively, for a given $\delta_k$ if
\begin{equation}\label{accEstimates}
\abs{f_k^0-f(\bm{{x}_k})}\leq\varepsilon_f\delta_k^2\quad \mbox{ and }\quad  \abs{f_k^s-f(\bm{{x}_k}+\Qk\bm{s_k})}\leq\varepsilon_f\delta_k^2.
\end{equation}
\end{definition}

Recalling Remark~\ref{whyTheseSigmaAlgebras},
we  extend  this definition to the following stochastic variant.
\begin{definition}\label{defProbAccEstim}
A sequence of random estimates $\accolade{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0,{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s}$ is said to be $\beta_f$-probabilistically $\varepsilon_f$-accurate with respect to the corresponding sequence $\accolade{\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k},{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k,\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}_k}}$ if the events
\begin{eqnarray*}
J_k^Q &=&\left\lbrace\mbox{Given}\ \Qk\ \mbox{and}\ \bm{{x}_k}, {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0 \mbox{ and } {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s \mbox{ are }\varepsilon_f\mbox{-accurate estimates of } f(\bm{{x}_k})  \mbox{ and}\right.\\
& &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \left.  f(\bm{{x}_k}+\Qk\bm{s_k}), \mbox{ respectively, for } {\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k \right\rbrace
\end{eqnarray*}
satisfy the supermartingale-like condition
$\prob{J^Q_k|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}\geq \beta_f.$
\end{definition}

Based on Definitions~\ref{wellAlignedDefinition}, \ref{defProbQkfully}, and~\ref{defProbAccEstim}, {\it true} iterations, in other words, those for which $\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}=1$, occur with a total probability of at least $\tilde{\beta}:=\beta_f\beta_m(1-\beta_Q)$ conditioned on $\mathcal{F}_{k-1}$. Noting that by construction $\mathds{1}_{\mathcal{A}_k}$ is both $\mathcal{F}_{k-1}^Q$-measurable and $\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$-measurable since $\mathcal{F}_{k-1}^Q\subset\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$, and that $\mathds{1}_{I_k^Q}$ is $\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$-measurable, then for $F:= \accolade{\mathcal{A}_k\cap I^Q_k\cap J_k^Q}$, we have
\begin{equation}\label{prodBetaIneq}
\begin{split}
\mathbb{P}(F|\mathcal{F}_{k-1})&=\E{\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}|\mathcal{F}_{k-1}}
=\E{\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\E{\mathds{1}_{J_k^Q}|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}} |\mathcal{F}_{k-1}}\\
&\geq \beta_f\E{\mathds{1}_{\mathcal{A}_k}\E{\mathds{1}_{I_k^Q}|\mathcal{F}_{k-1}^Q} |\mathcal{F}_{k-1}}\geq \beta_f\beta_m\E{\mathds{1}_{\mathcal{A}_k}|\mathcal{F}_{k-1}}\geq \tilde{\beta}.
\end{split}
\end{equation}
The second equality and the first two inequalities used the  {\it William's Tower Property} for conditional expectations~\cite[Theorem~34.4.]{billingsley1995ThirdEdition}, namely, the fact that if ${\ubar{z}}$ is integrable and the $\sigma$-algebras $\mathcal{G}_1$ and $\mathcal{G}_2$ satisfy $\mathcal{G}_1\subset\mathcal{G}_2$, then $\E{{\ubar{z}}|\mathcal{G}_1}=\E{\E{{\ubar{z}}|\mathcal{G}_2}|\mathcal{G}_1}$. Moreover, for a $\sigma$-algebra $\mathcal{G}$,  $\E{{\ubar{r}}{\ubar{u}}|\mathcal{G}}={\ubar{r}}\E{{\ubar{u}}|\mathcal{G}}$ if~${\ubar{r}}$ is $\mathcal{G}$-measurable and~${\ubar{r}}$ and~${\ubar{r}}{\ubar{u}}$ are integrable~\cite[Theorem~34.3.]{billingsley1995ThirdEdition}.

Note that even though the present algorithmic framework does not distinguish true iterations from {\it false} ones, it holds that $\prob{\underset{k\to\infty}{\limsup}\ \ubar{w}_k=\infty}=1$, where
\begin{equation}\label{randomWalk}
\ubar{w}_k:=\sum_{i=0}^k\left(2\mathds{1}_{\mathcal{A}_i}\mathds{1}_{I^Q_i}\mathds{1}_{J^Q_i}-1\right),
\end{equation}
as will be seen later in the proof of Theorem~\ref{liminfTheorem} provided $\tilde{\beta}>1/2$. This shows that true iterations occur sufficiently often for the convergence of Algorithm~\ref{algoStoScalTR} to hold (see, e.g.,~\cite[Theorem~3.6]{dzahini2020constrained} for a similar result).

\begin{remark
Calculations similar to those in~\eqref{prodBetaIneq} easily yield
\begin{eqnarray}\label{prodBetafBetam}
\prob{I^Q_k\cap J_k^Q|\mathcal{F}_{k-1}^Q}\geq \beta_f\beta_m, & & \prob{J_k^Q|\mathcal{F}_{k-1}^Q}\geq \beta_f,\\
\label{otherProbIneqs}
\prob{I^Q_k|\mathcal{F}_{k-1}}\geq\beta_m & \qquad \mbox{ and } \qquad & \prob{J^Q_k|\mathcal{F}_{k-1}}\geq\beta_f,
\end{eqnarray}
which will be useful in the proofs of Theorems~\ref{zerothOrderScal} and~\ref{liminfTheorem}, respectively.
\end{remark}

\subsection{Construction of probabilistic estimates and models in random subspaces.}
The construction of probabilistic estimates of Definition~\ref{defProbAccEstim} trivially follows the strategies described in~\cite[Section~2.3]{audet2019stomads}, \cite[Section~5]{chen2018stochastic}, and~\cite[Section~5.1]{dzahini2020constrained} and hence are not presented here again.

Next is stated an extension to subspaces of~\cite[Lemma~9.4]{AuHa2017}, which will be useful for the proof of Theorem~\ref{subspaceModel1}, one of the main results of the present work. Its proof is presented in the Appendix.

\begin{lemma}\label{FundTheorCalc}
Let $f$ be differentiable with a $L_g$-Lipschitz continuous gradient, and let $\bm{s},\bm{d}\in\mathbb{R}^p$. For a given matrix $\Qk$ and vector $\bm{{x}_k}$, let $\hat{f}(\bm{s})=f(\bm{{x}_k}+\Qk\bm{s})$, and assume that Assumption~\ref{QmaxAssumption} holds. Then
\[\abs{\hat{f}(\bm{s}+\bm{d})-\hat{f}(\bm{s})-\bm{d}^{\top}\nabla \hat{f}(\bm{s})}\leq \frac{1}{2}L_g\Qmax^2\normii{\bm{d}}^2.\]
\end{lemma}
The next result, which is a stochastic variant of~\cite[Theorem~9.5]{AuHa2017}, shows, together with Corollary~\ref{CorQkFully}, how $\beta_m$-probabilistically $(\kef,\keg; \Qk)$-fully linear models of Definition~\ref{defProbQkfully} can be made available.
\begin{theorem}\label{subspaceModel1}
Let the assumptions of Lemma~\ref{FundTheorCalc} hold.
Assume that there exists a finite constant $V_f>0$ such that
$\mathbb{V}_{\ubar{\theta}}[f_{\ubar{\theta}}(\bm{x})]\leq V_f$ for all $\bm{x}\in\mathbb{R}^n$. For all $i=0,1,\dots, p$, let ${\ubar{\theta}}^i_\ell, \ell=1,\dots,\pi_k$ be independent random samples of the independent random variables ${\ubar{\theta}}^i$ following the same distribution as ${\ubar{\theta}}$, and define ${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i):=\frac{1}{\pi_k}\sum_{\ell=1}^{\pi_k} f_{{\ubar{\theta}}^i_\ell}(\bm{{x}_k}+\Qk\bm{s_i})$, where the points  $\SSk:=\accolade{\szero,\sone,\dots,\spi}\subseteq\mathbb{R}^p\cap\mathcal{B}(\bm{0},c_1\delta_k)$ are affinely independent with $\szero=\bm{0}$ and $c_1\geq1$ a~constant. Let $\rho_k:=\underset{\bm{s_i} \in \SSk}{\max}\normii{\bm{s_i}-\szero}$ denote the approximate diameter of $\SSk$ and
define $\bm{\hat{L}_k}:=\frac{1}{\rho_k}\left[\sone-\szero\quad  \cdots\quad \spi-\szero \right]\in\mathbb{R}^{p\times p}$. Let $\bm{1}:=[1,\dots,1]^{\top}\in\R^{p+1}$, and define $\bm{Y}:=\left[\szero\ \, \sone\ \cdots\ \spi\right]\in\R^{p\times (p+1)}$. Consider the random vector  $\ubar{f}_k(\SSk,{\ubar{\theta}}):=\left[{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0), {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\sone,{\ubar{\theta}}^1),\dots, {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\spi,{\ubar{\theta}}^p)\right]^{\top}$ $\in\R^{p+1}$ and the random linear model $\Mhatk(\bm{s}):={\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}_0+\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}^{\top}\bm{s}$, where $({\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}_0,\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}})\in\R\times\mathbb{R}^p$ is the unique solution of $\left[\bm{1}\quad \bm{Y}^{\top}\right]\comb{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}_0}{\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}}=\ubar{f}_k(\SSk,{\ubar{\theta}})$.
Defining\footnote{While $\kef$ and $\keg$ seem to always depend on $k$, this is not the case since $\norme{\bm{\hat{L}_k}^{-1}}$ can be controlled by the geometry of the set $\SSk$, as will be seen later by means of Corollary~\ref{CorQkFully}.}
\begin{equation}\label{kefkeg}
\kef:= \frac{1}{2}\left(1+\frac{1}{2}\textcolor{black}{c_1^2\delta_{\max}}+2\textcolor{black}{c_1}\sqrt{p}\norme{\bm{\hat{L}_k}^{-1}}\right)L_g\Qmax^2,\quad \quad
\keg:= \left(1+\textcolor{black}{c_1}\sqrt{p}\norme{\bm{\hat{L}_k}^{-1}}\right)L_g\Qmax^2,
\end{equation}
consider the events
\begin{equation*}
E_f:= \accolade{\abs{\hat{f}(\bm{s})-\Mhatk(\bm{s})}\leq \kef\delta_k^2}, \qquad  E_g:=\accolade{\normii{\nabla \hat{f}(\bm{s})-\nabla\Mhatk(\bm{s})}\leq \keg\delta_k}.
\end{equation*}
Then
\begin{equation}\label{gradientBoundQlinear}
\mathbb{P}\left[E_f\cap E_g\right]\geq \beta_m \quad \mbox{for some}\ \beta_m\in(0,1)\ \mbox{ and all }\ \normii{\bm{s}}\leq \delta_k,
\end{equation}
provided
\begin{equation}\label{ppkChoice}
\pi_k\geq \frac{16V_f}{\textcolor{black}{c_1^2}L_g^2\Qmax^4\rho_k^2\min\{\delta_k^2,\delta_k^4\}\left(1-\beta_m^{1/(p+1)}\right)}.
\end{equation}
\end{theorem}
\begin{proof}
We build the random linear model $\Mhatk(\bm{s})={\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}_0+\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}^{\top}\bm{s}$ by seeking values for $({{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}_0},\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}})\in\R\times \R^{p}$ such that
$\Mhatk(\bm{s_i})={\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)$  for all $\bm{s_i}$ in the interpolation set $\SSk$,
which is equivalent to the $(p+1)\times(p+1)$ linear system of equations
\begin{equation}\label{systemTwo}
\left[\bm{1}\quad \bm{Y}^{\top}\right]\comb{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}_0}{\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}}=\ubar{f}_k(\SSk,{\ubar{\theta}}).
\end{equation}
Since the points in $\SSk$ are affinely independent, the matrix $\bm{\hat{L}_k}$ is invertible (see, e.g., \cite[Proposition~9.1]{AuHa2017}). Hence, system~\eqref{systemTwo} has a unique solution. Moreover, $\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}=\nabla\Mhatk(\bm{s})$ can also be computed by solving the $p\times p$ linear system
\begin{equation}\label{systemThree}
\rho_k \bm{\hat{L}_k}^{\top}\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}=\bm{{\ubar{\delta}}^{f_k(\SSk,{\ubar{\theta}})}}, \quad\mbox{where}\quad
\bm{{\ubar{\delta}}^{f_k(\SSk,{\ubar{\theta}})}}={\begin{bmatrix}
		{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\sone,{\ubar{\theta}}^1)-{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)}\\
		{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\stwo,{\ubar{\theta}}^2)-{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)}\\
		{\vdots}\\
		{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\spi,{\ubar{\theta}}^p)-{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)}
\end{bmatrix}}\in \mathbb{R}^p.
\end{equation}

To prove~\eqref{gradientBoundQlinear}, we consider the events
\[B_i:=\accolade{\abs{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)-\hat{f}(\bm{s_i})}\leq \frac{1}{2}\kappa_{eg}'\rho_k\min\{\delta_k,\delta_k^2\}},\quad i=0,1,\dots,p,\]
where $\kappa_{eg}':=\frac{1}{2}\textcolor{black}{c_1}L_g\Qmax^2$; and we assume that
\begin{equation}\label{BiProb}
\prob{B_i}\geq \beta_m^{1/(p+1)}\quad\mbox{for some }\ \beta_m\in (0,1)\  \mbox{and all}\ i=0,1,\dots,p.
\end{equation}
It follows from the independence of the random variables ${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i),\ i=0,1,\dots,p$, that the event $B:=\bigcap_{i=0}^pB_i$ satisfies $\prob{B}=\prod_{i=0}^{p}\prob{B_i}\geq \beta_m$. Then the remainder of the proof considers three  parts. The first two show respectively that $B\subseteq E_g$ and $B\subseteq E_f$, which imply that $B\subseteq E_f\cap E_g$ and hence $\prob{E_f\cap E_g}\geq \beta_m$. Part~3 provides the condition under which~\eqref{BiProb} holds.\\
$ $\\
\textbf{Part~1 ($B\subseteq E_g$)}. To demonstrate the latter inclusion, we will show that
\begin{equation}\label{FirstEvent}
B\subseteq E_1:= \accolade{\normii{\bm{\hat{L}_k}^{\top}\left[\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)\right]}\leq 2\sqrt{p}\kappa_{eg}'\delta_k}.
\end{equation}

Then, since the matrix norm is consistent with the Euclidean norm, the inequality \[\normii{\bm{\hat{L}_k}^{-\top}\bm{\hat{L}_k}\left(\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)\right)}\leq \norme{\bm{\hat{L}_k}^{-\top}}\normii{\bm{\hat{L}_k}\left(\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)\right)} \]
implies that $E_1\subseteq E_2:=\accolade{\normii{\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)}\leq 2\sqrt{p}\kappa_{eg}'\norme{\bm{\hat{L}_k}^{-1}}\delta_k}$, where we used the fact that $\norme{\bm{\hat{L}_k}^{-\top}}=\norme{\bm{\hat{L}_k}^{-1}}$. Using the $L_g$-Lipschitz continuity of $\nabla f$ and the fact that $\nabla \hat{f}(\bm{s})=\Qk^{\top}\nabla f(\bm{{x}_k}+\Qk\bm{s})$ with $\norme{\Qk}\leq \Qmax$, we get
\begin{eqnarray*}
\normii{\nabla \hat{f}(\bm{s})-\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}}&\leq&\normii{\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)}+\normii{\nabla \hat{f}(\szero)-\nabla \hat{f}(\bm{s})}\\
&\leq& 2\sqrt{p}\kappa_{eg}'\norme{\bm{\hat{L}_k}^{-1}}\delta_k + L_g\norme{\Qk^{\top}}\norme{\Qk}\normii{\szero-\bm{s}}\\
&\leq& \left(2\sqrt{p}\kappa_{eg}'\norme{\bm{\hat{L}_k}^{-1}}+L_g\Qmax^2\right)\delta_k=\left(1+\textcolor{black}{c_1}\sqrt{p}\norme{\bm{\hat{L}_k}^{-1}}\right)L_g\Qmax^2\delta_k,
\end{eqnarray*}
which implies that $E_2\subseteq E_g:=\accolade{\normii{\nabla \hat{f}(\bm{s})-\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}}\leq\left(1+\textcolor{black}{c_1}\sqrt{p}\norme{\bm{\hat{L}_k}^{-1}}\right)L_g\Qmax^2\delta_k}$.

To show~\eqref{FirstEvent}, we first notice, using~\eqref{systemThree}, that $\bm{\hat{L}_k}^{\top}\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}=\frac{1}{\rho_k}\bm{{\ubar{\delta}}^{f_k(\SSk,{\ubar{\theta}})}}$, and then
\begin{equation*
\bm{\hat{L}_k}^{\top}\nabla \hat{f}(\szero)=\frac{1}{\rho_k}\left[(\sone-\szero)^{\top}\nabla \hat{f}(\szero),\; \dots, \;  (\spi-\szero)^{\top}\nabla \hat{f}(\szero)\right]^{\top}\in\mathbb{R}^p.
\end{equation*}
Thus, the $i$th component of the vector $\bm{\hat{L}_k}^{\top}\left(\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)\right)$ is given by
\begin{equation*
\ubar{\Psi}_k(\bm{s_i},\szero):=\frac{1}{\rho_k}\left({\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)-{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)-(\bm{s_i}-\szero)^{\top}\nabla \hat{f}(\szero)\right).
\end{equation*}
We notice that
$
\abs{{\ubar{\Psi}}_k(\bm{s_i},\szero)}
\leq {{\ubar{\Psi}}_k^{1,i}} + {\psi_k^{2,i}},
$
where
\begin{eqnarray*}
{{\ubar{\Psi}}_k^{1,i}}&:=&\frac{1}{\rho_k}\abs{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)-\hat{f}(\bm{s_i})}+\frac{1}{\rho_k}\abs{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)-\hat{f}(\szero)}\\
\mbox{and }\quad {{{\psi}}_k^{2,i}}&:=&\frac{1}{\rho_k} \abs{\hat{f}(\bm{s_i})-\hat{f}(\szero)-(\bm{s_i}-\szero)^{\top}\nabla \hat{f}(\szero)}.
\end{eqnarray*}
It follows from Lemma~\ref{FundTheorCalc} that
\begin{equation}\label{psiktwoIneq}
{{{\psi}}_k^{2,i}}\leq\frac{1}{2\rho_k}L_g\Qmax^2\normii{\bm{s_i}-\szero}^2 \leq\frac{1}{2}L_g\Qmax^2\rho_k\leq \frac{1}{2}L_g\Qmax^2\textcolor{black}{c_1}\delta_k.
\end{equation}
Now assume that the event $B$ occurs. Then ${\Psi_k^{1,i}}\leq \kappa_{eg}'\delta_k$ for all $i=0,1,\dots,p$. It follows from~\eqref{psiktwoIneq} and the inequality $\abs{{\ubar{\Psi}}_k(\bm{s_i},\szero)}^2\leq 2\left(\left({\ubar{\Psi}}_k^{1,i}\right)^2+ \left({{\psi}}_k^{2,i}\right)^2\right)$ that

\begin{equation*
\begin{split}
\normii{\bm{\hat{L}_k}^{\top}\left(\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)\right)}^2&= \sum_{i=1}^p{\ubar{\Psi}}_k(\bm{s_i},\szero)^2\\
&\leq 2p\kappa_{eg}^{\prime\, 2}\delta_k^2+\frac{1}{2}pL_g^2\Qmax^4\textcolor{black}{c_1^2}\delta_k^2=pL_g^2\Qmax^4\textcolor{black}{c_1^2}\delta_k^2,
\end{split}
\end{equation*}
where the last equality follows from the definition of $\kappa_{eg}'$. This means that $B\subseteq E_1$.\\
$ $\\

\textbf{Part~2 ($B\subseteq E_f$)}. To show the latter inclusion, we recall that $\Mhatk(\szero)={\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)$ and $\Mhatk(\szero)-\Mhatk(\bm{s})=(\szero-\bm{s})^{\top}\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}$. Then the following holds:
\begin{eqnarray*}
\abs{\hat{f}(\bm{s})-\Mhatk(\bm{s})}&=&\abs{\hat{f}(\bm{s})-{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)+\Mhatk(\szero)-\Mhatk(\bm{s})}\\
&=&\abs{\hat{f}(\bm{s})-{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)+(\szero-\bm{s})^{\top}\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}}\\
&=&\left\lvert\hat{f}(\bm{s})-\hat{f}(\szero)-(\bm{s}-\szero)^{\top}\nabla\hat{f}(\szero)-\left({\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)-\hat{f}(\szero)\right) +(\bm{s}-\szero)^{\top}\left(\nabla\hat{f}(\szero)-\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}\right)\right\rvert\\
&\leq& \abs{\hat{f}(\bm{s})-\hat{f}(\szero)-(\bm{s}-\szero)^{\top}\nabla\hat{f}(\szero)}+ \abs{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)-\hat{f}(\szero)} + \normii{\bm{s}-\szero}\normii{\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla\hat{f}(\szero)}.
\end{eqnarray*}
As before, the first term on the right-hand side of the last inequality is bounded by using Lemma~\ref{FundTheorCalc} as follows:
\begin{equation}\label{lemmabound2}
\abs{\hat{f}(\bm{s})-\hat{f}(\szero)-(\bm{s}-\szero)^{\top}\nabla\hat{f}(\szero)}\leq \frac{1}{2}L_g\Qmax^2\normii{\bm{s}-\szero}^2\leq\frac{1}{2}L_g\Qmax^2\delta_k^2.
\end{equation}
To bound the last two terms, we note that the inclusions $B\subseteq B_0$ and $B\subseteq E_1\subseteq E_2$ yield  $\abs{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\szero,{\ubar{\theta}}^0)-\hat{f}(\szero)}$ $\leq \frac{1}{2}\kappa_{eg}'\rho_k\delta_k^2\ $ and $\ \normii{\bm{{\stackunder[1.0pt]{$a$}{\rule{0.8ex}{.075ex}}}}-\nabla \hat{f}(\szero)}\leq 2\sqrt{p}\kappa_{eg}'\norme{\bm{\hat{L}_k}^{-1}}\delta_k$, respectively. Thus, these inequalities combined with~\eqref{lemmabound2} and the inequalities $\rho_k\leq c_1\delta_{\max}$ and $\normii{\bm{s}-\szero}\leq \delta_k$ lead to

\begin{eqnarray*}
\abs{\hat{f}(\bm{s})-\Mhatk(\bm{s})}&\leq&\frac{1}{2}L_g\Qmax^2\delta_k^2+\frac{1}{2}\kappa_{eg}'\rho_k\delta_k^2+2\sqrt{p}\kappa_{eg}'\norme{\bm{\hat{L}_k}^{-1}}\delta_k^2\\
&=&\frac{1}{2}\left(1+\frac{1}{2}\textcolor{black}{c_1^2\delta_{\max}}+2\textcolor{black}{c_1}\sqrt{p}\norme{\bm{\hat{L}_k}^{-1}}\right)L_g\Qmax^2\delta_k^2,
\end{eqnarray*}
which shows that $B\subseteq E_f$.\\
$ $\\
\textbf{Part~3$\ $} To complete the proof, we show~\eqref{BiProb} using the Chebyshev inequality as follows:
\begin{eqnarray*}
\prob{B_i^c}&=&\prob{\abs{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)-\mathbb{E}\left[{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)\right]}> \frac{1}{2}\kappa_{eg}'\rho_k\min\{\delta_k,\delta_k^2\}}\\
&\leq& \frac{4 \mathbb{V}\left[{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)\right]}{\kappa_{eg}^{\prime\, 2}\rho_k^2\min\{\delta_k^2,\delta_k^4\}}\leq \frac{4 V_f}{\pi_k\kappa_{eg}^{\prime\, 2}\rho_k^2\min\{\delta_k^2,\delta_k^4\}}\leq 1-\beta_m^{1/(p+1)},
\end{eqnarray*}
provided~\eqref{ppkChoice} holds; that is, $\pi_k$ is chosen according to
\[\pi_k\geq \frac{4 V_f}{\kappa_{eg}^{\prime\, 2}\rho_k^2\min\{\delta_k^2,\delta_k^4\}\left(1-\beta_m^{1/(p+1)}\right)}= \frac{16V_f}{\textcolor{black}{c_1^2}L_g^2\Qmax^4\rho_k^2\min(\delta_k^2,\delta_k^4)\left(1-\beta_m^{1/(p+1)}\right)}.\]
\end{proof}
The next result shows that for a particular geometry of the interpolation set $\SSk$ of Theorem~\ref{subspaceModel1}, the resulting model gradient is a forward finite-difference stochastic gradient estimator of $\Qk^{\top}\nabla f(\bm{{x}_k})$ and its corresponding parameters $\kef$ and $\keg$ do not depend on $k$, as mentioned above.
\begin{corollary}\label{CorQkFully}
Under all the assumptions of Theorem~\ref{subspaceModel1}, assume further that $\bm{s_i}=h\bm{e_i}\in\mathbb{R}^p, i\in \intbracket{1,p}$, where $\bm{e_i}$ is the $i$th standard basis vector of $\mathbb{R}^p$ and $h=\min\{\hopt, \delta_k\}$ for some $\hopt>0$ sufficiently small. Let ${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0(\bm{{x}_k}):={\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(0,{\ubar{\theta}}^0)$ and ${\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{{x}_k}+h(\Qk)_{: i}):={\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{s_i},{\ubar{\theta}}^i)$ for all $i\in \intbracket{1,p}$, where $(\Qk)_{: i}$ denotes the $i$th column of $\Qk$. Then $\Qk^{\top}\nabla f(\bm{{x}_k})\approx \bm{\hat{g}_k}=\bm{\hat{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}}_k}(\omega)$, where $\bm{\hat{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}}_k}$ is defined by the  forward finite-difference scheme
$(\bm{\hat{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}}_k})_i=\frac{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s(\bm{{x}_k}+h(\Qk)_{:i})-{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0(\bm{{x}_k})}{h}$, $i\in \intbracket{1,p}$.
Moreover, the corresponding model parameters~$\kef$ and~$\keg$ given by~\eqref{kefkeg} are independent of $k$.
\end{corollary}

\begin{proof}
Note that $\SSk=\accolade{\bm{0},h\bm{e_1},h\bm{e_2},\dots,h\bm{e_p}}\subseteq\mathbb{R}^p\cap\B{0}{c_1\delta_k}$, $\rho_k=h$, and hence $\bm{\hat{L}_k}=\bm{I_p}$. The proof immediately follows from~\eqref{systemThree} by replacing $\ubar{\bm{a}}$ with $\bm{\hat{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}}_k}$. The nondependence of~$\kef$ and~$\keg$ on $k$ trivially follows from the fact that $\bm{\hat{L}_k}^{-1}=\bm{I_p}$.
\end{proof}

\section{Convergence analysis.}\label{Section4}
This section presents convergence results of Algorithm~\ref{algoStoScalTR} using ideas inspired by~\cite{audet2019stomads,billingsley1995ThirdEdition,chen2018stochastic,dzahini2020expected,dzahini2020constrained}. Section~\ref{Section4p1} presents preliminary results necessary for the proof of the main results. Section~\ref{Section4p2} proves that the sequence of random trust-region radii converges to zero almost surely. Section~\ref{Section4p3} demonstrates the existence of a subsequence of random iterates generated by the proposed method, which drives the norm of $\nabla f$ to zero almost surely. Proofs that are not presented in this section can be found in the Appendix.

\subsection{Preliminary results.}\label{Section4p1}
In the remainder of the manuscript, the following inspired by~\cite[Assumptions~4.1 and~4.3]{chen2018stochastic}, respectively, will be assumed.
\begin{assumption}\label{fGradfAssumption}
For given $\bm{\bm{x}_0}$, $\delta_{\max}$, and $\bm{Q_0}$, let $\mathscr{L}(\bm{\bm{x}_0};\bm{Q_0})\subset\mathbb{R}^n$ be the set containing all iterates of Algorithm~\ref{algoStoScalTR}. Define the region $\mathscr{L}^{\star}(\bm{\bm{x}_0};\bm{Q_0})$ considered by the algorithm realizations as $\mathscr{L}^{\star}(\bm{\bm{x}_0};\bm{Q_0})=\displaystyle{\underset{\accolade{k\, \geq\, 0,\, \bm{{x}_k}\in\mathscr{L}(\bm{\bm{x}_0};\bm{Q_0})}}{\bigcup}}\mathcal{B}(\bm{{x}_k},\delta_{\max};\Qk)$. Then $f(\bm{x})\geq f_{\min}$ for all $\bm{x}\in\mathbb{R}^n$ and for some constant $f_{\min}>-\infty$. Moreover, $f$ and its gradient are Lipschitz continuous on $\mathscr{L}^{\star}(\bm{\bm{x}_0};\bm{Q_0})$.
\end{assumption}

\begin{assumption}\label{boundedHessian}
There exists some $\kappa_h\geq 1$ such that for all $k\geq 0$, the Hessian $\bm{\hat{H}}_k$ of all realizations of $\Mhatk$ satisfies $\norme{\bm{\hat{H}}_k}\leq \kappa_h.$
\end{assumption}
The following subspace variant of~\cite[Lemma~4.5]{chen2018stochastic} shows that if $\delta_k$ is small enough compared with the size $\normii{\bm{\hat{g}_k}}$ of a $(\kef,\keg;\Qk)$-fully linear model gradient, then the trial step $\bm{s_k}$ provides a decrease in $f$ proportional to $\normii{\bm{\hat{g}_k}}$.
\begin{lemma}\label{LemmaA}
Suppose that the function $\mhatk$ is a $(\kef,\keg;\Qk)$-fully linear model for~$f$ in $\B{\bm{{x}_k}}{\delta_k;\Qk}$. Then the trial step $\bm{s_k}$ leads to an improvement in $f$ such that
\begin{equation}\label{LemmaAimprovement}
f(\bm{{x}_k}+\Qk\bm{s_k})-f(\bm{{x}_k})\leq -\frac{\kappa_{fcd}}{4}\normii{\bm{\hat{g}_k}}\delta_k,
\quad
\mbox{if }
\delta_k\leq \min\accolade{\frac{1}{\kappa_h},\frac{\kappa_{fcd}}{8\kappa_{ef}}}\normii{\bm{\hat{g}_k}}.
\end{equation}
\end{lemma}

The following subspace variant of~\cite[Lemma~4.6]{chen2018stochastic} shows that the guaranteed decrease in~$f$ provided by $\bm{s_k}$ is proportional to $\normii{\QkTranspose\nabla f(\bm{{x}_k})}$ if $\delta_k$ is small enough compared with $\normii{\QkTranspose\nabla f(\bm{{x}_k})}$.
\begin{lemma}\label{LemmaB2}
Let Assumption~\ref{boundedHessian} hold, and suppose the model $\mhatk$ is $(\kef,\keg;\Qk)$-fully linear. Define $C_1=\frac{\kappa_{fcd}}{4} \max\accolade{\frac{\kappa_h}{\kappa_h+\kappa_{eg}}, \frac{8\kappa_{ef}}{8\kappa_{ef}+\kappa_{fcd}\kappa_{eg}}}$. Then the trial step $\bm{s_k}$ leads to an improvement in $f$ such that

\begin{equation}\label{LemmaB2improvement}
f(\bm{{x}_k}+\Qk\bm{s_k})-f(\bm{{x}_k})\leq -C_1\normii{\QkTranspose\nabla f(\bm{{x}_k})}\delta_k,
\end{equation}
\begin{equation}\label{deltaLemmaB2}
\mbox{if}\quad
\delta_k\leq \min\accolade{\frac{1}{\kappa_h+\kappa_{eg}}, \frac{1}{\frac{8\kappa_{ef}}{\kappa_{fcd}}+\kappa_{eg}}}\normii{\QkTranspose\nabla f(\bm{{x}_k})}.
\end{equation}
\end{lemma}
The next subspace variant of~\cite[Lemma~4.7]{chen2018stochastic} shows that if the model and the estimates are sufficiently accurate and  $\delta_k$ is small enough compared with $\normii{\bm{\hat{g}_k}}$, then the iteration is successful.
\begin{lemma}\label{LemmaC}
Let Assumption~\ref{boundedHessian} hold. Assume that $\mhatk$ is $(\kef,\keg;\Qk)$-fully linear and the estimates $\accolade{f_k^0,f_k^s}$ are $\varepsilon_f$-accurate with $\varepsilon_f\leq\kappa_{ef}$. Then the $k$th iteration is successful if
\begin{equation}\label{deltaLemmaC}
\delta_k\leq\min\accolade{\frac{1}{\kappa_h}, \frac{1}{\eta_2}, \frac{\kappa_{fcd}(1-\eta_1)}{8\kappa_{ef}}}\normii{\bm{\hat{g}_k}}.
\end{equation}
\end{lemma}

\begin{proof}
The proof immediately follows from that of~\cite[Lemma~4.7]{chen2018stochastic} with minor modifications and is not presented here again.
\end{proof}
The next result, which  is a subspace variant of~\cite[Lemma~4.8]{chen2018stochastic}, guarantees an amount of decrease in $f$ on true successful iterations.
\begin{lemma}\label{LemmaD}
Suppose that Assumption~\ref{boundedHessian} holds and the estimates $\accolade{f_k^0,f_k^s}$ are $\varepsilon_f$-accurate with $\varepsilon_f<\frac{1}{4}\eta_1\eta_2\kappa_{fcd}\min\accolade{\frac{\eta_2}{\kappa_h}, 1}$. If the $k$th iteration is successful, then the improvement in $f$ is such that
\begin{equation}\label{LemmaDimprovement}
f(\bm{{x}_{k+1}})-f(\bm{{x}_k})\leq -C_2\delta_k^2,
\end{equation}
where $C_2=\frac{1}{2}\eta_1\eta_2\kappa_{fcd}\min\accolade{\frac{\eta_2}{\kappa_h}, 1} -2\varepsilon_f>0$.
\end{lemma}

\begin{proof}
Again, the proof immediately follows from that of~\cite[Lemma~4.8]{chen2018stochastic} with minor modifications and is not presented here.
\end{proof}
Inspired by a result from the proof of~\cite[Theorem~4.11]{chen2018stochastic}, the next lemma quantifies the maximum possible amount of increase in~$f$ on false successful iterations.
\begin{lemma}\label{LemmaE}
Under Assumptions~\ref{QmaxAssumption} and~\ref{fGradfAssumption}, assume that Algorithm~\ref{algoStoScalTR} erroneously accepts a step $\bm{s_k}$ leading to an increase in the objective function when $\norme{\QkTranspose\nabla f(\bm{{x}_k})}\geq \zeta\delta_k$ for some constant $\zeta>0$, and let $C_3(\zeta)=1+\frac{3L_g}{2\zeta}\Qmax^2$. Then the increase in $f$ is such that
\begin{equation}\label{LemmaEincrease}
f(\bm{{x}_k}+\Qk\bm{s_k})-f(\bm{{x}_k})\leq C_3(\zeta)\norme{\QkTranspose\nabla f(\bm{{x}_k})}\delta_k.
\end{equation}
\end{lemma}

Later, to prove the key result of Theorem~\ref{zerothOrderScal}, we favor techniques derived in~\cite{paquette2018stochastic}, unlike~\cite{blanchet2016convergence,chen2018stochastic}. These techniques, also used in~\cite{audet2019stomads,dzahini2020expected,dzahini2020constrained}, make use of event indicator functions. Here, by means of Lemma~\ref{indicatorMeasurable}, we introduce a measurability result that is related to one of these indicator functions and is crucial in the present analysis. But first, the following result is required.
\begin{lemma}\label{partialMeasurable}
Any partial derivative
$\frac{\partial f}{\partial x^i}: \left(\mathbb{R}^n, \mathscr{B}(\mathbb{R}^n)\right)\to \left(\R, \mathscr{B}(\R)\right),\ i\in\intbracket{1,n} $,
of a differentiable function $f$ is a Borel measurable function.
\end{lemma}

\begin{proof}
Let $(\bm{e_1},\bm{e_2},\dots,\bm{e_n})$ be the canonical basis of $\mathbb{R}^n$. For any $i\in\intbracket{1,n} $ and $\bm{x}=(x^1,x^2,\dots,x^n)$, define, for all $j\in\mathbb{N}$,
$D_{ij}(\bm{x})=\frac{f(\bm{x}+h_j\bm{e_i})-f(\bm{x})}{h_j}$,
where $\accolade{h_j}_{j\in\mathbb{N}}$ is a sequence of positive real numbers converging to~$0$; $f$ is differentiable and hence continuous, and thus is Borel measurable~\cite[Theorem~13.2]{billingsley1995ThirdEdition}. Therefore, since for all $i\in\intbracket{1,n}$, $\underset{j\to\infty}{\lim}D_{ij}(\bm{x})=\frac{\partial f}{\partial x^i}(\bm{x}),$
then $\frac{\partial f}{\partial x^i}$ is also Borel measurable~\cite[Theorem~13.3 and Theorem~13.4-(ii)]{billingsley1995ThirdEdition}.
\end{proof}

\begin{lemma}\label{indicatorMeasurable}
Let $\zeta\geq 0$ be a constant, and consider the event $\Gamma_k:=\accolade{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}_k\geq 0}$, where ${{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}_k:=\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}-\zeta{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k$. Then ${{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}_k$, $\mathds{1}_{\Gamma_k}$, and $\mathds{1}_{\bar{\Gamma}_k}$ are $\mathcal{F}_{k-1}^Q$-measurable.
\end{lemma}

\begin{proof}
By construction, all the entries $\left(\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}\right)_{ij},\ (i,j)\in\intbracket{1,n}\times\intbracket{1,p}$ of $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$ are $\mathcal{F}_{k-1}^Q$-measurable, while ${\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k$ and the components $\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k}^i, \ i\in\intbracket{1, n}$ of $\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k}$ are $\mathcal{F}_{k-1}$-measurable and hence $\mathcal{F}_{k-1}^Q$-measurable since $\mathcal{F}_{k-1}\subseteq\mathcal{F}_{k-1}^Q$. Thus, for any $i\in\intbracket{1, n}$, $\frac{\partial f}{\partial x^i}(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})$ is $\mathcal{F}_{k-1}^Q$-measurable as the composition of the Borel measurable function $\frac{\partial f}{\partial x^i}$ (from Lemma~\ref{partialMeasurable}) and the random map $\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k}:(\Omega,\mathcal{F}_{k-1}^Q)\to\left(\mathbb{R}^n, \mathscr{B}(\mathbb{R}^n)\right)$~\cite[Theorem~13.3 and Theorem~13.1-(ii)]{billingsley1995ThirdEdition}. Recall that  sums and products of measurable functions are also measurable~\cite[Theorem~13.3]{billingsley1995ThirdEdition},
and note that the random variable
\begin{equation}\label{ukRandomMeasurable}
\ubar{u}_k:=\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}^2=\sum_{j=1}^p\left(\sum_{i=1}^n\left(\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}\right)_{ij}\frac{\partial f}{\partial x^i}(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k}) \right)^2
\end{equation}
is consequently $\mathcal{F}_{k-1}^Q$-measurable. Thus ${{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}_k=\sqrt{{\ubar{u}_k}}-\zeta{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k$ is also $\mathcal{F}_{k-1}^Q$-measurable. The $\mathcal{F}_{k-1}^Q$-measurability of ${{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}_k$ ensures in particular that $\Gamma_k={{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}_k^{-1}([0, \infty))\in\mathcal{F}_{k-1}^Q$ since $[0, \infty)$ is a Borel set of $\R$, whence the {\it simple} real functions $\mathds{1}_{\Gamma_k}$ and $\mathds{1}_{\bar{\Gamma}_k}$~\cite[Equation~(13.3)]{billingsley1995ThirdEdition} are also $\mathcal{F}_{k-1}^Q$-measurable.
\end{proof}

\subsection{Zeroth-order convergence result.}\label{Section4p2}
To prove in Theorem~\ref{zerothOrderScal} that the sequence of random trust-region radii converges to zero almost surely, we assume the following.
\begin{assumption}\label{assumptionsOnEstimatesAndModels}
The sequences of estimates $\accolade{{\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^0, {\stackunder[0.6pt]{$f$}{\rule{.8ex}{.075ex}}}_k^s}$ and models $\accolade{\Mhatk}$ generated by Algorithm~\ref{algoStoScalTR} are, respectively, $\beta_f$-probabilistically $\varepsilon_f$-accurate for
$\varepsilon_f<\min\accolade{\kappa_{ef}, \frac{1}{4}\eta_1\eta_2\kappa_{fcd}\min\accolade{\frac{\eta_2}{\kappa_h}, 1}}$ and $\beta_m$-probabilistically $(\kef,\keg;\Qk)$-fully linear for some $\beta_f, \beta_m\in (0,1)$.
\end{assumption}

Recall Definitions~\ref{wellAlignedDeterministic} and~\ref{wellAlignedDefinition}. The result presented next shows that the sequence of trust-region radii converges to zero whether the sequence of subspace selection matrices is well aligned or not. The corresponding proof improves on the analyses presented in~\cite{blanchet2016convergence,chen2018stochastic}, for example, by reducing the number of subcases of {\bf Case~1}, discussed below, from four to three. This improvement is the outcome of a simple motivation to circumvent the derivation of a lower bound on the probability of the event that ``either the model is good and the estimates are bad, or the model is bad and the estimates are good,'' as was done in~\cite{blanchet2016convergence,chen2018stochastic}. Another remarkable difference is the fact that both subcases do not straightforwardly use the same $\sigma$-algebras for conditioning on the past, as is the case in~\cite{blanchet2016convergence,chen2018stochastic}, which results from the need to introduce a random variable $\mathds{1}_{G_k}$ related to the model gradient in {\bf Case~2}. We also note that unlike prior similar works, the remainder of the present analysis explicitly emphasizes the way all the $\sigma$-algebras $\mathcal{F}_{k-1}$, $\mathcal{F}_{k-1}^Q$, and $\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$ work together for the proofs of the proposed results.

\begin{theorem}\label{zerothOrderScal}
Let all assumptions that were made in Lemmas~\ref{LemmaA}-\ref{LemmaD} hold with the same constants $C_1, C_2$, and $C_3:=C_3(\zeta)$, for some fixed $\zeta>0$ satisfying
\begin{equation}\label{zetaChoice}
\zeta\geq \kappa_{eg}+\max\accolade{\eta_2,\kappa_h, \frac{8\kappa_{ef}}{\kappa_{fcd}(1-\eta_1)}}.
\end{equation}
Let $\nu\in (0,1)$ be chosen according to
\begin{equation}\label{nuChoiceScal}
\frac{\nu}{1-\nu} \geq \frac{8\gamma^2}{\min\accolade{C_2, \zeta C_1}}.
\end{equation}
Assume further that Assumption~\ref{assumptionsOnEstimatesAndModels} holds with $\beta_f, \beta_m\in (0,1)$ satisfying\footnote{In~\eqref{probChoices}, the result is intentionally presented with $(1-\beta_f)+(1-\beta_m)$ instead of $2-\beta_f-\beta_m$ in order to emphasize one of its differences compared with~\cite{blanchet2016convergence,chen2018stochastic} in which a term similar to $(1-\beta_f)(1-\beta_m)$ was derived.}
\begin{equation}\label{probChoices}
\frac{\beta_f\beta_m-\frac{1}{2}}{(1-\beta_f)+(1-\beta_m)}\geq \frac{C_3}{C_1}\quad \mbox{ and }\quad \frac{\beta_f}{1-\beta_f}\geq \frac{2\nu\zeta\gamma^2C_3}{(1-\nu)(\gamma^2-1)}+2\gamma^2.
\end{equation}
Define $\varrho=\frac{1}{2}\beta_f(1-\nu)(1-\frac{1}{\gamma^2})>0$. Then, the random function ${\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k:=\nu (f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})-f_{\min})+(1-\nu){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2$ satisfies
\begin{equation}\label{PhikIncrementBound}
\E{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k|\mathcal{F}_{k-1}^Q}\leq -\varrho{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2 \quad \mbox{for all}\ k\in\mathbb{N},
\end{equation}
which implies that
\begin{equation}\label{seriesResult}
\sum_{k=0}^{\infty}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2<\infty\quad \mbox{almost surely.}
\end{equation}
\end{theorem}

\begin{proof}
The proof is inspired by those of~\cite[Theorem~3]{blanchet2016convergence} and~\cite[Theorem~4.11]{chen2018stochastic}.
The overall goal is to prove~\eqref{PhikIncrementBound}. Indeed noticing that ${\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k>0$ and then taking expectations on both sides of the inequality in~\eqref{PhikIncrementBound}, we obtain~\eqref{seriesResult} (see, e.g.,~\cite[Theorem~3]{dzahini2020expected} for details). Let $\phi_k$ denote realizations of ${\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k$, and recall that on all successful iterations, $\bm{{x}_{k+1}}=\bm{{x}_k}+\Qk\bm{s_k}$ and ${\delta_{k+1}}=\min\accolade{\gamma\delta_k,\delta_{\max}}$,  implying
\begin{equation}\label{un}
\phi_{k+1}-\phi_k\leq \nu\left(f(\bm{{x}_{k+1}})-f(\bm{{x}_k})\right)+(1-\nu)(\gamma^2-1)\delta_k^2,
\end{equation}
while on unsuccessful iterations, $\bm{{x}_{k+1}}=\bm{{x}_k}$ and ${\delta_{k+1}}=\gamma^{-1}\delta_k$, in which case
\begin{equation}\label{bOne}
\phi_{k+1}-\phi_k\leq(1-\nu)\left(\frac{1}{\gamma^2}-1\right)\delta_k^2=:b_1<0.
\end{equation}

Recall the event $\Gamma_k=\accolade{\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}\geq\zeta{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k}$ of Lemma~\ref{indicatorMeasurable} with $\zeta$ satisfying ~\eqref{zetaChoice}. The proof considers two cases: $\mathds{1}_{\Gamma_k}=1$ and $\mathds{1}_{\bar{\Gamma}_k}=1$. Inspired by~\cite[proof of Theorem~4.11]{chen2018stochastic}, Case~1 aims to show that
\begin{equation}\label{huit}
\E{\mathds{1}_{\Gamma_k}\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k\right)|\mathcal{F}_{k-1}^Q}\leq-2\mathds{1}_{\Gamma_k}(1-\nu)(\gamma^2-1){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2
\leq -\frac{1}{2}\mathds{1}_{\Gamma_k}\beta_f(1-\nu)(1-\frac{1}{\gamma^2}){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2,
\end{equation}
where the last inequality follows from $1-\frac{1}{\gamma^2}<\gamma^2-1$ and $\beta_f\leq 1$. On the other hand, inspired by~\cite[proof of Theorem~3]{blanchet2016convergence}, Case~2 shows that
\begin{equation}\label{last}
\E{\mathds{1}_{\bar{\Gamma}_k}\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k\right)|\mathcal{F}_{k-1}^Q}\leq-\frac{1}{2}\mathds{1}_{\bar{\Gamma}_k}\beta_f(1-\nu)(1-\frac{1}{\gamma^2}){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2.
\end{equation}
Thus, combining~\eqref{huit} and~\eqref{last} leads to~\eqref{PhikIncrementBound}.

$ $\\
\textbf{Case~1: $\mathds{1}_{\Gamma_k}=1$}\\
$ $\\
Unlike the proof of~\cite[Theorem~4.11, Case~1]{chen2018stochastic} where four subcases were considered, only three are analyzed next.

\noindent \textit{(Subcase 1i)}
Good model ($\mathds{1}_{I_k^Q}=1$) and good estimates ($\mathds{1}_{J_k^Q}=1$). By~\eqref{zetaChoice},
\[\normii{\QkTranspose\nabla f(\bm{{x}_k})}\geq \max\accolade{\eta_2+\kappa_{eg}, \kappa_h+\kappa_{eg}, \frac{8\kappa_{ef}}{\kappa_{fcd}}+\kappa_{eg}}\delta_k,\]
implying in particular~\eqref{deltaLemmaB2} and hence a decrease in $f$ according to~\eqref{LemmaB2improvement}. Moreover, by $(\kef,\keg;\Qk)$-full linearity,
\[\normii{\bm{\hat{g}_k}}\geq \normii{\QkTranspose\nabla f(\bm{{x}_k})}-\kappa_{eg}\delta_k\geq \max\accolade{\eta_2,\kappa_h, \frac{8\kappa_{ef}}{\kappa_{fcd}(1-\eta_1)}}\delta_k, \]
thus implying~\eqref{deltaLemmaC}; and since $\varepsilon_f\leq \kappa_{ef}$ per assumptions of Lemma~\ref{LemmaC}, the iteration is successful. Consequently,~\eqref{un} together with~\eqref{LemmaB2improvement} yields
\begin{equation}\label{bTwo}
\begin{split}
	\phi_{k+1}-\phi_k &\leq -\nu C_1\normii{\QkTranspose\nabla f(\bm{{x}_k})}\delta_k+(1-\nu)(\gamma^2-1)\delta_k^2=:b_2\\
	&\leq \left(-\nu C_1\zeta+(1-\nu)(\gamma^2-1)\right)\delta_k^2<0,
\end{split}
\end{equation}
where the last inequality follows from~\eqref{nuChoiceScal}. Denoting by ${\ubar{b}}_2({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)$ the random variable with realizations $b_2<0$, with ${\ubar{u}}_k$ defined in~\eqref{ukRandomMeasurable}, we have from~\eqref{bTwo} that
\begin{equation}\label{bTwoRandom}
\mathds{1}_{\Gamma_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k\right)\leq \mathds{1}_{\Gamma_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}{\ubar{b}}_2({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k).
\end{equation}
\noindent \textit{(Subcase 1ii)}
Bad model ($\mathds{1}_{\bar{I}_k^Q}=1$) and good estimates ($\mathds{1}_{J_k^Q}=1$). In this case, regardless of the iteration type (i.e., successful or unsuccessful), the change in $\phi_k$ can always be bounded by using~\eqref{un} and~\eqref{LemmaEincrease} as follows:
\begin{equation}\label{bThree}
\phi_{k+1}-\phi_k\leq\nu C_3\norme{\QkTranspose\nabla f(\bm{{x}_k})}\delta_k+(1-\nu)(\gamma^2-1)\delta_k^2=:b_3.
\end{equation}
Denoting by ${\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)$ the random variable with realizations $b_3>0$, from~\eqref{bThree} we have that
\begin{equation}\label{bThreeRandom2}
\mathds{1}_{\Gamma_k}\mathds{1}_{\bar{I}_k^Q}\mathds{1}_{J_k^Q}\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k\right)\leq \mathds{1}_{\Gamma_k}\mathds{1}_{\bar{I}_k^Q}\mathds{1}_{J_k^Q}{\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)\leq \mathds{1}_{\Gamma_k}\mathds{1}_{\bar{I}_k^Q}{\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k).
\end{equation}
\noindent \textit{(Subcase 1iii)}
Bad estimates ($\mathds{1}_{\bar{J}_k^Q}=1$). From \textit{Subcase 1ii}, it always holds that
\begin{equation}\label{bThreeRandom3}
\mathds{1}_{\Gamma_k}\mathds{1}_{\bar{J}_k^Q}\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k\right)\leq \mathds{1}_{\Gamma_k}\mathds{1}_{\bar{J}_k^Q}{\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k).
\end{equation}

With the subcases thus complete, since the random variables $\mathds{1}_{\Gamma_k}, {\ubar{b}}_2({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)$, and ${\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)$ are $\mathcal{F}_{k-1}^Q$-measurable thanks to the proof of Lemma~\ref{indicatorMeasurable}, combining~\eqref{bTwoRandom},~\eqref{bThreeRandom2}, and~\eqref{bThreeRandom3} and taking expectations with respect to $\mathcal{F}_{k-1}^Q$, we have
\begin{equation}\label{bTwoThree}
\begin{split}
\E{\mathds{1}_{\Gamma_k}\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k\right)|\mathcal{F}_{k-1}^Q}&\leq \mathds{1}_{\Gamma_k}{\ubar{b}}_2({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)\E{\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}|\mathcal{F}_{k-1}^Q}\\
&+\mathds{1}_{\Gamma_k}{\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)\left(\E{\mathds{1}_{\bar{I}_k^Q}|\mathcal{F}_{k-1}^Q}+\E{\mathds{1}_{\bar{J}_k^Q}|\mathcal{F}_{k-1}^Q}\right)\\
&\leq \mathds{1}_{\Gamma_k}\left(\beta_f\beta_m{\ubar{b}}_2({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)+{\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)(1-\beta_m+1-\beta_f)\right)=:\mathds{1}_{\Gamma_k} {\ubar{b}}_{2,3}({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k),
\end{split}
\end{equation}
where the last inequality follows from~\eqref{ikQProb} and~\eqref{prodBetafBetam} and the fact that ${\ubar{b}}_2({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)<0$ while ${\ubar{b}}_3({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)>0$.  Note that
\begin{eqnarray*}
{\ubar{b}}_{2,3}({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)&=&\nu \normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\left[-\beta_f\beta_m C_1+(2-\beta_f-\beta_m)C_3\right]\\
& &+(1-\nu)(\gamma^2-1)(2-\beta_f-\beta_m+\beta_f\beta_m){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2\\
&\leq& \nu \normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\left[-\beta_f\beta_m C_1+(2-\beta_f-\beta_m)C_3\right] + 2(1-\nu)(\gamma^2-1){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2,
\end{eqnarray*}
where the last inequality follows from $2-\beta_f-\beta_m+\beta_f\beta_m=1+(1-\beta_f)(1-\beta_m)\leq 2$.

For $\beta_f$ and $\beta_m$ chosen according to the first condition in~\eqref{probChoices}, and for $\nu$ satisfying~\eqref{nuChoiceScal}, the following holds:
\[\beta_f\beta_m C_1 - (2-\beta_f-\beta_m)C_3\geq \frac{1}{2}C_1\geq 2\frac{2(1-\nu)(\gamma^2-1)}{\nu\zeta}.\]
Hence, since $\normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}\geq \zeta{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k$, then
\begin{eqnarray*}
{\ubar{b}}_{2,3}({\ubar{u}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k)&\leq& -\frac{1}{2}\left[\beta_f\beta_m C_1-(2-\beta_f-\beta_m)C_3\right]\nu \normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\\
&\leq& -\frac{1}{4}C_1\nu \normii{\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\leq -2(1-\nu)(\gamma^2-1){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2,
\end{eqnarray*}
which, together with~\eqref{bTwoThree}, leads to~\eqref{huit}.

$ $\\
\textbf{Case~2: $\mathds{1}_{\bar{\Gamma}_k}=1$}
$ $\\
Consider the event $G_k:=\accolade{\normii{\bm{\hat{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}}_k}}\geq \eta_2{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k}$, and note that since $\Mhatk:={\ubar{f}_k}+\bm{\hat{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}}_k}^{\top}\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}} +\frac{1}{2}\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}}^{\top}\hat{{\stackunder[1.0pt]{$H$}{\rule{1.3ex}{.075ex}}}}_k\bm{{\stackunder[1.0pt]{$s$}{\rule{.8ex}{.075ex}}}}$ is $\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$-measurable by construction, then in particular $\normii{\bm{\hat{{{\stackunder[1.0pt]{$g$}{\rule{1.0ex}{.075ex}}}}}_k}}$ and hence $\mathds{1}_{G_k}, \mathds{1}_{\bar{G}_k}$ are also $\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$-measurable. Observe that if $\mathds{1}_{\bar{G}_k}=1$, then the iteration is unsuccessful since $\normii{\bm{\hat{g}_k}}<\eta_2\delta_k$, leading to~\eqref{bOne}. Hence,
\begin{equation}\label{GkbarCase2}
\begin{split}
\E{\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{\bar{G}_k}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k)|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}&\leq -(1-\nu)\left(1-\frac{1}{\gamma^2}\right)\E{\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{\bar{G}_k}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}\\
&\leq -\frac{1}{2}\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{\bar{G}_k}\beta_f(1-\nu)\left(1-\frac{1}{\gamma^2}\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2,
\end{split}
\end{equation}
where the last inequality is due to the $\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$-measurability of $\mathds{1}_{\bar{\Gamma}_k}$ and ${\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k$, given that $\mathcal{F}_{k-1}\subset\mathcal{F}_{k-1}^Q\subset\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$.

Inspired by~\cite[Theorem~3]{blanchet2016convergence}, which improved on the proof of~\cite[Theorem~4.11, Case~2]{chen2018stochastic} in which four subcases were considered, only two are analyzed next. The overall goal is to prove that
\begin{equation}\label{GkCase2}
\E{\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k)|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}\leq -\frac{1}{2}\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}\beta_f(1-\nu)\left(1-\frac{1}{\gamma^2}\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2,
\end{equation}
which, combined with~\eqref{GkbarCase2}, leads to
\begin{equation*
\E{\mathds{1}_{\bar{\Gamma}_k}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k)|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}\leq -\frac{1}{2}\mathds{1}_{\bar{\Gamma}_k}\beta_f(1-\nu)\left(1-\frac{1}{\gamma^2}\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2
\end{equation*}
and consequently
\begin{equation*
\begin{split}
\E{\mathds{1}_{\bar{\Gamma}_k}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k)|\mathcal{F}_{k-1}^Q}&=\E{\E{\mathds{1}_{\bar{\Gamma}_k}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k)|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}|\mathcal{F}_{k-1}^Q}\\
&\leq \E{-\frac{1}{2}\mathds{1}_{\bar{\Gamma}_k}\beta_f(1-\nu)\left(1-\frac{1}{\gamma^2}\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2 \, |\mathcal{F}_{k-1}^Q}=-\frac{1}{2}\mathds{1}_{\bar{\Gamma}_k}\beta_f(1-\nu)\left(1-\frac{1}{\gamma^2}\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2,
\end{split}
\end{equation*}
 showing~\eqref{last}, where the last equality follows from the $\mathcal{F}_{k-1}^Q$-measurability of $\mathds{1}_{\bar{\Gamma}_k}$ and ${\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k$.  What remains to be shown is~\eqref{GkCase2} as achieved next, assuming that $\mathds{1}_{G_k}=1$.

\noindent \textit{(Subcase 2i)}
Good estimates ($\mathds{1}_{J_k^Q}=1$) and $\mathds{1}_{G_k}=1$. While $f$ is decreased on successful iterations because  of good estimates thanks to Lemma~\ref{LemmaD} and in this case~\eqref{LemmaDimprovement} holds,  $\delta_k$ is reduced on unsuccessful iterations, leading to~\eqref{bOne}. However, combining~\eqref{LemmaDimprovement} and~\eqref{un} yields
$\phi_{k+1}-\phi_k\leq \left[-\nu C_2+(1-\nu)(\gamma^2-1)\right]\delta_k^2\leq b_1,$
where the last inequality is due to~\eqref{nuChoiceScal}, which shows that~\eqref{bOne} holds in any case. Consequently,
\begin{equation}
\label{GammakBarGkJk}
\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}\mathds{1}_{J_k^Q}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k) \leq -\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}\mathds{1}_{J_k^Q}(1-\nu)\left(1-\frac{1}{\gamma^2}\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2.
\end{equation}

\noindent \textit{(Subcase 2ii)}
Bad estimates ($\mathds{1}_{\bar{J}_k^Q}=1$) and $\mathds{1}_{G_k}=1$. In this case, a successful step can lead to an increase in $f$
according to~\eqref{LemmaEincrease}, which combined with~\eqref{un} yield~\eqref{bThree} and hence
\begin{equation}\label{GammakBarGkJkBar}
\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}\mathds{1}_{\bar{J}_k^Q}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k) \leq -\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}\mathds{1}_{\bar{J}_k^Q}\left(\nu C_3\zeta+(1-\nu)(\gamma^2-1)\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2.
\end{equation}

With the subcases complete, recall that $\E{\mathds{1}_{J_k^Q}|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}\geq\beta_f$. Then  combining~\eqref{GammakBarGkJk} and~\eqref{GammakBarGkJkBar} and taking expectations with respect to~$\mathcal{F}^{\hat{m}\cdot Q}_{k-1}$ lead to
\begin{eqnarray*
\E{\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k)|\mathcal{F}^{\hat{m}\cdot Q}_{k-1}}&\leq& \mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}\left(-\beta_f(1-\nu)(1-1/\gamma^2)\right.\\
& &\left. +\ (1-\beta_f)\left(\nu C_3\zeta+(1-\nu)(\gamma^2-1)\right) \right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2\\
&\leq& -\frac{1}{2}\mathds{1}_{\bar{\Gamma}_k}\mathds{1}_{G_k}\beta_f(1-\nu)\left(1-\frac{1}{\gamma^2}\right){\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2,
\end{eqnarray*}
thus proving~\eqref{GkCase2}, where the last inequality follows from the second condition in~\eqref{probChoices}, and the proof is complete.
\end{proof}

\subsection{Liminf-type convergence.}\label{Section4p3}
The next result demonstrates the existence of a subsequence of random iterates generated by Algorithm~\ref{algoStoScalTR}, which drives $\nabla f$ to zero almost surely. While the corresponding proof is inspired by that of~\cite[Theorem~4.16]{chen2018stochastic}, we point out the additional difficulty introduced in the present work by the random matrices $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}$. Unlike~\cite[Section~4]{chen2018stochastic} where ``auxiliary lemmas'' similar to those of Section~\ref{Section4p1} are directly related to $\normii{\nabla f(\bm{{x}_k})}$ and the size $\normii{\bm{g_k}}$ of a full space model gradient, thus easing the proof of the liminf-type result, this is not the case here. Instead, recovering $\normii{\nabla f(\bm{{x}_k})}$ through $\normii{\QkTranspose \nabla f(\bm{{x}_k})}$ by means of the well-alignment assumption becomes crucial.
\begin{theorem}\label{liminfTheorem}
Let all the assumptions made in Theorem~\ref{zerothOrderScal} hold. Assume further that Assumption~\ref{wellAlignedAssumption} holds with $0\leq\beta_Q<1-\frac{1}{2\beta_f\beta_m}$. Then, the sequence $\accolade{\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k}}_{k\in\mathbb{N}}$ of random iterates generated by Algorithm~\ref{algoStoScalTR} satisfies
\[\underset{k\to\infty}{\liminf}\normii{\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}=0\quad\mbox{ almost surely.}\]
\end{theorem}

\begin{proof}
The result is proved by contradiction conditioned on the almost sure event $E_0:=\accolade{{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\to 0}$ (thanks to Theorem~\ref{zerothOrderScal}), inspired by the proof of~\cite[Theorem~4.16]{chen2018stochastic} and also using ideas from~\cite{audet2019stomads,dzahini2020constrained}. Assume that with nonzero probability there exists a random variable ${\ubar{\varepsilon}}'>0$ such that
$\normii{\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}\geq {\ubar{\varepsilon}}' \mbox{ for all } k\in\mathbb{N}.$

Let $\accolade{\bm{{x}_k}}, \accolade{\delta_k}$, and ${{\varepsilon}}'$ be realizations of $\accolade{\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k}}, \accolade{{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k}$, and ${\ubar{\varepsilon}}'$, respectively, for which $\normii{\nabla f(\bm{{x}_k})}\geq {{\varepsilon}}'$ for all $k$. Since $\delta_k\to 0$, there exists $k_0\in\mathbb{N}$ such that
\begin{equation}\label{bChoice}
\delta_k<b:=\alpha_Q\cdot\min\accolade{\frac{{{\varepsilon}}'}{2(\kappa_h+\kappa_{eg})}, \frac{{{\varepsilon}}'}{2(\eta_2+\kappa_{eg})}, \frac{{{\varepsilon}}'}{\frac{16\kappa_{ef}}{\kappa_{fcd}(1-\eta_1)}+2\kappa_{eg}}, \frac{\delta_{\max}}{\gamma}}
\end{equation}
for all $k\geq k_0$.
Let $\ubar{r}_k$ be the random variable with realizations $r_k=\log_{\gamma}\left(\frac{\delta_k}{b}\right)$. Then $r_k<0$ for all $k\geq k_0$. The main idea is to show that such realizations occur only with probability zero, leading to a contradiction.

To prove that $\accolade{\ubar{r}_k}$ is a submartingale, recall the events $I^Q_k$ and $J^Q_k$ satisfying~\eqref{otherProbIneqs} and $\mathcal{A}_k$ satisfying~\eqref{Akineq}.
Consider some iteration $k\geq k_0$ for which $\mathcal{A}_k$, $I^Q_k$, and $J^Q_k$ all occur, which happens with probability at least $\beta_f\beta_m(1-\beta_Q)$ conditioned on $\mathcal{F}_{k-1}$, as was shown in~\eqref{prodBetaIneq}. It follows from~\eqref{bChoice} that
\begin{equation}\label{ineq4}
\normii{\nabla f(\bm{{x}_k})}\geq \frac{1}{\alpha_Q}\max\accolade{\kappa_h+\kappa_{eg}, \eta_2+\kappa_{eg}, \frac{8\kappa_{ef}}{\kappa_{fcd}(1-\eta_1)}+\kappa_{eg}}\delta_k,
\end{equation}
which, together with the $(\kef,\keg;\Qk)$-full linearity and $\alpha_Q$-well alignment, yield
\begin{equation}\label{ineq5}
\normii{\bm{\hat{g}_k}}\geq\normii{\QkTranspose\nabla f(\bm{{x}_k})}-\kappa_{eg}\delta_k\geq \alpha_Q \normii{\nabla f(\bm{{x}_k})}-\kappa_{eg}\delta_k
\geq \max\accolade{\kappa_h, \eta_2, \frac{8\kappa_{ef}}{\kappa_{fcd}(1-\eta_1)}}\delta_k.
\end{equation}
Thus, the $k$th iteration is successful according to Lemma~\ref{LemmaC}. Hence, ${\delta_{k+1}}=\gamma\delta_k$, and hence $r_{k+1}=r_k+1$. For all other outcomes of $\mathcal{A}_k$, $I^Q_k$, and $J^Q_k$, which occur with a total probability of at most $1-\beta_f\beta_m(1-\beta_Q)$, it always holds that ${\delta_{k+1}}\geq\gamma^{-1}\delta_k$, which implies that $r_{k+1}\geq r_k-1$. Hence,
\begin{eqnarray*}
\E{\mathds{1}_{\mathcal{A}_k\cap I^Q_k\cap J^Q_k} \left(\ubar{r}_{k+1}-\ubar{r}_k\right)|\mathcal{F}_{k-1}}&\geq& \beta_f\beta_m(1-\beta_Q)\quad \quad\mbox{and}\\
\E{\mathds{1}_{\overline{\mathcal{A}_k\cap I^Q_k\cap J^Q_k}} \left(\ubar{r}_{k+1}-\ubar{r}_k\right)|\mathcal{F}_{k-1}}&\geq& \beta_f\beta_m(1-\beta_Q)-1,
\end{eqnarray*}
which shows that $\accolade{\ubar{r}_k}$ is a submartingale if $0\leq\beta_Q<1-\frac{1}{2\beta_f\beta_m}$.

Using the random walk defined in~\eqref{randomWalk}, one  can easily  show (following, e.g.,~\cite[Theorem~4]{audet2019stomads}, \cite[Theorem~4.16]{chen2018stochastic} and~\cite[Theorem~3.6 and Lemmas~4.7 and~4.12]{dzahini2020constrained}) that $\accolade{\ubar{w}_k}$ is also a submartingale with bounded increments and, as such, cannot converge to a finite value, and therefore $\prob{\underset{k\to\infty}{\limsup}\ \ubar{w}_k=\infty}=1$ (thanks to~\cite[Theorem~4.4]{chen2018stochastic}). However, since $r_k-r_{k_0}\geq w_k-w_{k_0}$ by construction,  the sequence of realizations~$r_k$ such that $r_k<0\ \forall k\ \geq k_0$ occurs with probability zero, leading to a contradiction, which achieves the proof.
\end{proof}

\section{Expected complexity analysis.}\label{Section5}
Section~\ref{Section5p1} introduces relevant assumptions, definitions, and theorems derived in the analysis of a general renewal-reward stochastic process and its associated stopping time introduced in~\cite{blanchet2016convergence} for the expected complexity analysis of a stochastic trust-region method. That renewal-reward process was also used in~\cite{dzahini2020expected,paquette2018stochastic} for the  analyses of stochastic direct-search and line-search methods. We show  in Section~\ref{ExpectedSubsection}  how the aforementioned assumptions are satisfied for Algorithm~\ref{algoStoScalTR} and we bound the expected number of iterations required to achieve  $\norme{\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}\leq \varepsilon$.

\subsection{A renewal-reward martingale process.}\label{Section5p1}

A formal definition from~\cite{blanchet2016convergence} of the stopping time related to a discrete time stochastic process is as follows.
\begin{definition}
A random variable~${{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}$ is called a stopping time with respect to a given discrete time stochastic process $\accolade{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k}_{k\in\mathbb{N}}$ if the event $\accolade{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}=k}$ belongs to the $\sigma$-algebra $\sigma\left({\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_1,\dots,{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k \right)$ generated by ${\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_1,\dots,{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k$, for each $k\in\mathbb{N}$.
\end{definition}

Consider a stochastic process $\accolade{\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\right)}_{k\in\mathbb{N}}$, where ${\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k, {\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\in [0,\infty)$, and introduce on the same probability space as $\accolade{\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\right)}_{k\in\mathbb{N}}$ a biased random walk process $\accolade{{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}}_{k\in\mathbb{N}}$ obeying the following dynamics:
\begin{equation}\label{dynamicsWk}
\prob{{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_{k+1}}=1|\mathcal{F}_k}=q\qquad \mbox{ and }\qquad \prob{{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_{k+1}}=-1|\mathcal{F}_k}=1-q,
\end{equation}
where $q\in (1/2, 1)$ and $\mathcal{F}_k=\sigma\left( \left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_0,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_0,{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_0}\right),\dots,\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k,{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}\right) \right)$, with ${{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_0}=1}$.

Let $\accolade{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}}_{\epsilon>0}$ be a family of stopping times with respect to $\accoladekinN{\mathcal{F}_k}$, parameterized by $\epsilon$. A bound on $\E{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}}$ is derived in~\cite{blanchet2016convergence,dzahini2020expected,paquette2018stochastic} under the following assumption.

\begin{assumption}\label{assumptionRenewal}
The following hold for the stochastic process $\accoladekinN{\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k,{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}\right)}$.
\begin{itemize}
\item[(i)] There exist $\lambda\in (0, \infty)$ and $\delta_{\max}=\delta_0 e^{\lambda j_{\max}}$ for some $j_{\max}\in\mathbb{Z}$ such that ${\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\leq \delta_{\max}$ for all $k$.
\item[(ii)] There exists $\delta_{\epsilon}=\delta_0 e^{\lambda j_{\epsilon}}$, for some $j_{\epsilon}\in\mathbb{Z}$ with $j_{\epsilon}\leq 0$ such that for all $k$
\begin{equation}\label{dynamicsRenewal}
\mathds{1}_{\{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon} > k\}}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_{k+1}\geq \mathds{1}_{\{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon} > k\}}\min\accolade{{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k e^{\lambda{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_{k+1}}}, \delta_{\epsilon}},
\end{equation}
where ${{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_{k+1}}$ satisfies~\eqref{dynamicsWk}.
\item[(iii)] There exists a nondecreasing function $h:[0,\infty)\to (0,\infty)$ and a constant $\varrho>0$ such that for all $k$
\[\E{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k|\mathcal{F}_k}\mathds{1}_{\{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon} > k\}}\leq-\varrho h\left({\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\right)\mathds{1}_{\{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon} > k\}}.\]
\end{itemize}
\end{assumption}

Noticing that the event $\accolade{{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\geq\delta_{\epsilon}}$ occurs sufficiently frequently on average, often $\E{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k}$ can be bounded by some negative fixed constant \cite{blanchet2016convergence}, thus allowing a bound on the expected stopping time $\E{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}}$, as stated next.

\begin{theorem}\label{theorComplexity1}
Under Assumption~\ref{assumptionRenewal},
$\E{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}}\leq \frac{q}{2q-1}\cdot\frac{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_0}{\varrho h\left(\delta_{\epsilon}\right)}+1.$
\end{theorem}

\subsection{Expected complexity result.}\label{ExpectedSubsection}

Consider the process $\accolade{\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k,{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}\right)}_{k\in\mathbb{N}}$, where~${\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k$ is defined in Theorem~\ref{zerothOrderScal}, ${\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k$ is the trust-region radius, and ${{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}$ is  defined by ${{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}=2\left(\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}-\frac{1}{2} \right)$.  Given $\epsilon\in (0,1)$, let ${{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}$ be the number of iterations required by Algorithm~\ref{algoStoScalTR} to first drive the norm of the gradient of~$f$ below $\epsilon$:
\[{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}:=\inf \accolade{k\in\mathbb{N}: \normii{\nabla f(\bm{{\stackunder[1.0pt]{$x$}{\rule{1.0ex}{.075ex}}}_k})}\leq\epsilon}.\]
Then ${{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}$ is a stopping time for the stochastic process generated by Algorithm~\ref{algoStoScalTR}, and hence for $\accolade{\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k,{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}\right)}_{k\in\mathbb{N}}$ \cite{blanchet2016convergence,dzahini2020expected,paquette2018stochastic}. More precisely, the occurrence of $\accolade{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}=k}$ can be determined by observing $\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_0,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_0,{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_0}\right)$, $\dots$ $({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k-1},{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_{k-1},{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_{k-1}})$, which means that ${{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}$ is a stopping time with respect to the filtration $\accoladekinN{\mathcal{F}_{k-1}}$. To bound $\E{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}}$ using Theorem~\ref{theorComplexity1}, we show next that Assumption~\ref{assumptionRenewal} holds for $\accolade{\left({\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k,{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k,{{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}\right)}_{k\in\mathbb{N}}$.

From~\eqref{PhikIncrementBound}, since $\mathcal{F}_{k-1}\subset\mathcal{F}_{k-1}^Q$, we observe that\footnote{Unlike the analysis in~\cite{blanchet2016convergence}, the first equality of~\eqref{PhikIncrSecond} is essential and due to the fact that in the present framework ${{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}$ is a stopping time with respect to $\sigma$-algebras $\mathcal{F}_{k-1}$ that are smaller than $\mathcal{F}_{k-1}^Q$ with respect to which $\accoladekinN{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k}$ was proved in Theorem~\ref{zerothOrderScal} Eq.~\eqref{PhikIncrementBound} to be a supermartingale.} for all $k\in\mathbb{N}$:
\begin{equation}\label{PhikIncrSecond}
\E{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k|\mathcal{F}_{k-1}}=\E{\E{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k|\mathcal{F}_{k-1}^Q}|\mathcal{F}_{k-1}}\leq\E{-\varrho{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2|\mathcal{F}_{k-1}}=-\varrho{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2.
\end{equation}
Note that~\eqref{PhikIncrSecond} holds for any realizations of the random variables $\E{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_{k+1}-{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_k|\mathcal{F}_{k-1}}$ and ${\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k^2$ and hence on the event $\accolade{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon} >k}$ in particular, which shows that Assumption~\ref{assumptionRenewal}-$(iii)$ holds with the constant $\varrho=\frac{1}{2}\beta_f(1-\nu)(1-\frac{1}{\gamma^2})$ of Theorem~\ref{zerothOrderScal} and $h(t)=t^2$. Assumption~\ref{assumptionRenewal}-$(i)$ automatically follows from the initialization strategy in Algorithm~\ref{algoStoScalTR} with $\lambda=\ln(\gamma)$. Before showing by means of Lemma~\ref{smallLemma}, inspired by~\cite[Lemma~7]{blanchet2016convergence}, that~\eqref{dynamicsRenewal} is satisfied, we first define the constant
\begin{equation}\label{deltaEpsilonDefinition}
\delta_{\epsilon}:=\frac{\epsilon}{\xi}, \quad\mbox{with } \xi\geq  \frac{1}{\alpha_Q}\max\accolade{\kappa_h+\kappa_{eg}, \eta_2+\kappa_{eg}, \frac{8\kappa_{ef}}{\kappa_{fcd}(1-\eta_1)}+\kappa_{eg}}=:\hat{\xi},
\end{equation}
inspired by the proof of Theorem~\ref{liminfTheorem}, especially~\eqref{ineq4}. Then, following ~\cite{blanchet2016convergence} exactly, one can assume without loss of generality that $\delta_{\epsilon}=\gamma^i\delta_0$ for some integer $i=:j_{\epsilon}\leq 0$, whence ${\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k=\gamma^{{\ubar{i}}_k}\delta_{\epsilon}$ for any $k$ and some integer~${\ubar{i}}_k$.

\begin{lemma}\label{smallLemma}
Let all the assumptions made in Theorem~\ref{liminfTheorem} hold. Then~\eqref{dynamicsRenewal} is satisfied for\\ ${{\stackunder[1.0pt]{$\tilde{w}$}{\rule{0.8ex}{.075ex}}}_k}=2\left(\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}-\frac{1}{2} \right)$, with $\lambda=\ln(\gamma)$ and some fixed $q\in (\tilde{\beta}, 1)$ with $\tilde{\beta}=\beta_f\beta_m(1-\beta_Q)$.
\end{lemma}

\begin{proof}
The result is proved by suitably adapting the proof of~\cite[Lemma~7]{blanchet2016convergence} as was done in~\cite{dzahini2020expected}. First, we notice that~\eqref{dynamicsRenewal} trivially holds when $\mathds{1}_{\{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon} > k\}}=0$. Next, we show that when $\mathds{1}_{\{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon} > k\}}=1$, then
\[{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_{k+1}\geq \min\accolade{\delta_{\epsilon}, \min\accolade{\delta_{\max},\gamma{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k}\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}+\gamma^{-1}{\stackunder[1.0pt]{$\delta$}{\rule{.8ex}{.075ex}}}_k\left(1-\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}\right)}. \]
Recall that $\delta_k=\gamma^{i_k}\delta_{\epsilon}$ for some integer $i_k$. Therefore, if $\delta_k>\delta_{\epsilon}$, then $\delta_k\geq \gamma\delta_{\epsilon}$, which implies that ${\delta_{k+1}}\geq\gamma^{-1}\delta_k\geq\delta_{\epsilon}$. Now assume that $\delta_k\leq \delta_{\epsilon}$. Since ${{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}>k$, we have $\normii{\nabla f(\bm{{x}_k})}>\epsilon=\xi\delta_{\epsilon}\geq\xi\delta_k$. Thus, it follows from the definition of~$\xi$ that~\eqref{ineq4} holds. If $\mathds{1}_{\mathcal{A}_k}=1$, $\mathds{1}_{I_k^Q}=1$, and $\mathds{1}_{J_k^Q}=1$, then~\eqref{ineq5} holds, and consequently the $k$th iteration is successful, as was  explained in the proof of Theorem~\ref{liminfTheorem}. Hence, $\bm{{x}_{k+1}}=\bm{{x}_k}+\Qk\bm{s_k}$ and ${\delta_{k+1}}=\min\accolade{\delta_{\max},\gamma\delta_k}$. If $\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}=0$, then it always holds that ${\delta_{k+1}}\geq \gamma^{-1}\delta_k$. The proof is completed by  observing that $\prob{\mathds{1}_{\mathcal{A}_k}\mathds{1}_{I_k^Q}\mathds{1}_{J_k^Q}=1}=q$ for some fixed $q\geq \tilde{\beta}$.
\end{proof}
The main complexity result is provided by the next  theorem.

\begin{theorem}
Let all the assumptions made in Theorem~\ref{liminfTheorem} hold. Then
\[\E{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}}\leq \frac{\tilde{\beta}}{2\tilde{\beta}-1}\cdot\frac{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_0\tilde{\xi}^2}{\varrho \epsilon^2}+1\]
for some $\tilde{\xi}\geq\hat{\xi}$, where $\varrho$ is the constant of Theorem~\ref{zerothOrderScal}, $\tilde{\beta}$ is the same constant of Lemma~\ref{smallLemma}, and $\hat{\xi}$ is defined in~\eqref{deltaEpsilonDefinition}.
\end{theorem}

\begin{proof}
Recall the choice of $q$ in Lemma~\ref{smallLemma}, and pick some $\tilde{\xi}\geq\hat{\xi}$ such that $\delta_{\epsilon}=\frac{\epsilon}{\tilde{\xi}}=\gamma^{i}\delta_0$ for some integer $i\leq 0$ (without loss of generality), as discussed above. The proof follows by employing Theorem~\ref{theorComplexity1} with $h(t)=t^2$, which yields
\[\E{{{\stackunder[1.0pt]{$\tau$}{\rule{0.6ex}{.075ex}}}_\epsilon}}\leq \frac{q}{2q-1}\cdot\frac{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_0\tilde{\xi}^2}{\varrho \epsilon^2}+1\leq \frac{\tilde{\beta}}{2\tilde{\beta}-1}\cdot\frac{{\stackunder[1.0pt]{$\phi$}{\rule{0.8ex}{.075ex}}}_0\tilde{\xi}^2}{\varrho \epsilon^2}+1.\]
\end{proof}

\input{numerical}

\section*{Conclusion.}
This work introduces STARS, the first DFO algorithm developed for stochastic objective functions that achieves scalability using random models constructed in low-dimensional random subspaces. The analysis of STARS extends an existing framework of model-based stochastic DFO (where randomness comes only from the stochasticity of the objective function), to settings where additional randomness stems from the mechanics of the algorithm. Making use of a supermartingale-based framework, we prove that the expected complexity of STARS, which uses subspace models, is similar to that of stochastic DFO algorithms in a smooth nonconvex setting.
Numerical experiments demonstrate the performance of STARS on large-scale problems using linear interpolation models in subspaces of various dimensions.

We note that for ease of exposition we have focused in
\eqref{probl1} on $f_{\ubar{\theta}}$ being an unbiased estimator of $f$, but this can be relaxed in all ways that the STORM framework can address. Similarly, for concreteness, the stated algorithm and numerical results have focused respectively on linear and quadratic models, but the analysis readily applies to more general probabilistically fully linear models in subspaces.
Apart from the use of particular random interpolation models, the present work also employs a classical Monte Carlo sampling strategy for computing estimates. Future works can improve this sampling strategy while also using more general random models.

\bibliographystyle{siamplain}

\section{Numerical results.}\label{Section6}
We now illustrate the performance of STARS for various random subspace and full space (i.e., STORM-like) forms. We consider stochastically noisy variants of $40$ deterministic unconstrained problems ranging in dimension from $n=98$ to $n=125$; see Table~\ref{tab:probs}. All objective functions are sums of squares, that is, $f(\bm{x})=\sum_{i=1}^m f_i(\bm{x})^2$,
and are corrupted with either {\it additive} or {\it multiplicative} stochastic noise. In the former case, the noisy $f_{\ubar{\theta}}$ is given by $f_{\ubar{\theta}}(\bm{x})=f(\bm{x})+{\ubar{\theta}}$;  in the latter, $f_{\ubar{\theta}}(\bm{x})=f(\bm{x})(1+{\ubar{\theta}})$. In both cases, the centered random variable ${\ubar{\theta}}$ with standard deviation~$\sigma=10^{-3}$ is either normally or uniformly distributed.
In order to take into account the variability due to the stochastic noise or random subspaces, $20$ replications (corresponding to random seeds common across the tests) were performed for each of the~$40$ problems; combined with the two distributions of ${\ubar{\theta}}$, we thus have a total of 1,600 {\it problem instances}.

For all tested variants, we
use $p$-dimensional linear interpolation models $\mhatk$, where $\nabla \mhatk = \bm{\hat{g}_k}$ is obtained from the forward finite-difference approximation from Corollary~\ref{CorQkFully}. We employ the forward finite-difference parameter $h=\min\accolade{\hopt,\delta_k}$, with $\hopt$ obtained following \cite{more2011edn} and provided in Table~\ref{tab:probs} and $\delta_k$ denoting the current trust-region radius.
We employ such linear interpolation models so as to focus on differences due to the size and form of $\QkTranspose$ and since the corresponding trust-region subproblems can be solved exactly:
$\bm{s_k}^{\star}=-\delta_k\frac{\bm{\hat{g}_k}}{\normii{\bm{\hat{g}_k}}}$ if $\bm{\hat{g}_k}\neq \bm{0}$ and $\bm{s_k}^{\star}=\bm{0}$ otherwise.
When computing estimates of unknown function values at each iteration by means of a Monte Carlo approach using $n_k=25$ noisy function evaluations for all~$k$, available samples from previous iterations are reused, following the strategy described in the last paragraph of~\cite[Section~2.3]{audet2019stomads}, which was also used in~\cite[Section~5.2]{dzahini2020constrained}. For $\QkTranspose$, motivated by Theorems~\ref{GaussianJLTmatrix} and~\ref{sHashingTheorem}, we tested both Gaussian and $(r=1)$-Hashing strategies.
We found these to perform comparably for the tested settings and hence report
the results for the Gaussian case (labeled {\sf G-STARS-}$p$) whereby the entries of $\bm{{\stackunder[0.6pt]{$Q$}{\rule{1.3ex}{.075ex}}}_k}^{\top}\in \mathbb{R}^{p\times n}$ are distributed as $\mathcal{N}(0,1/p)$.
We also test a STORM-like instance of STARS (labeled {\sf I-STARS-}$n$) where $\Qk=\bm{I_n}\in\mathbb{R}^{n\times n}$. All algorithmic variants used the parameters $\gamma=2$, $\eta_2= 90\eta_1=0.9$, $\delta_{\max}=5$, $c_1=1$ and~$\delta_0=1$.

All STARS variants are assessed by using data
profiles~\cite{MoWi2009}. For each of the 1,600 noisy problem instances, let~$\bm{x_N}$ be the point with the best~$f$ function value obtained by an algorithm after $N$ evaluations of~$f_{\ubar{\theta}}$, denote by~$f^{\star}$ the least such value found by all the algorithms, and let $\bm{x_0}$ be the starting point. A problem is considered successfully solved within a convergence tolerance~$\tau\in [0,1]$ after $N$ evaluations if
\begin{equation*
	f(\bm{x_N})\leq f^{\star} + \tau(f(\bm{x_0}) - f^{\star}).
\end{equation*}
The horizontal and vertical axes of the data profiles show, respectively, the number of noisy function evaluations divided by~$(n+1)$ and the proportion of problems solved. During the experiments a budget of 1,500$(n+1)$ noisy function evaluations is allocated to all the algorithms. For multiplicative noise, Figure~\ref{fig:dataprofsig3M} (left) shows that when the overall noise (i.e., $f(\bm{x}){\ubar{\theta}}$) is large (which is the case for many of these problems: for nearly half of these problem instances $f(\bm{x}_0)\ge 10^3$), the {\sf G-STARS} variants with small values of $p$ perform well relative to the full space {\sf I-STARS-}$n$ variant. On the other hand, for this sampling budget and additive noise level, Figure~\ref{fig:dataprofsig3M} (right) shows that {\sf I-STARS-}$n$ is competitive with {\sf G-STARS} with $p=2, 5$. Unsurprisingly, in both cases, {\sf G-STARS-}$n$ (which uses $p=\min\accolade{n,100}$) is outperformed by the {\sf I-STARS-}$n$, which is using a better conditioned, deterministic matrix throughout.
\begin{figure}
	\centering
	\includegraphics[width=0.49\linewidth]{Data_Tau2_Multiplicative} \hfill
	\includegraphics[width=0.49\linewidth]{Data_Tau2_Additive}
	\caption[]{\small{Data profiles for convergence tolerance $\tau=10^{-2}$ on 1,600  problem instances for multiplicative noise (left) and additive noise (right) with standard deviation $\sigma = 10^{-3}$.}}
	\label{fig:dataprofsig3M}
\end{figure}