\begin{document}

 
\title{Markovian Foundations for Quasi-Stochastic Approximation
\
with Applications to Extremum Seeking Control}
 
	
 
 \begin{abstract}
 

This paper concerns  quasi-stochastic approximation (QSA) to solve root finding problems commonly found in applications to optimization and reinforcement learning.     The general constant gain algorithm may be expressed as the time-inhomogeneous ODE $ \ddt \Uptheta_t = \alpha f_t (\Uptheta_t)$,  with state process $\Uptheta$ evolving on $\field{R}^n$.    Theory is based on an almost periodic vector field, so that in particular the time average of $f_t(\theta)$ defines the time-homogeneous mean vector field ${\widebar{f}}\colon\field{R}^d\to\field{R}^d$ with ${\widebar{f}}(\theta^*) = {\sf Z}ero$.    Under smoothness assumptions on the functions involved, the following exact representation is obtained:
\[
	\ddt \Uptheta_t    =   \alpha [{\widebar{f}} (\Uptheta_t)   -\alpha \widebar{\Upupsilon}_t +    \alpha^2 {\cal W}_t^0  +   \alpha \frac{d}{dt}   {\cal W}_t^1 +\frac{d^2}{dt^2}   {\cal W}_t^2]

\]
along with formulae for the smooth signals $\{ \widebar{\Upupsilon}_t, {\cal W}_t^i :  i=0, 1, 2\}$.    This representation is based on the application of techniques from Markov processes,  for which Poisson's equation plays a central role.    This new representation, combined with new conditions for ultimate boundedness, has many applications for furthering the theory of QSA and its applications, including the following implications that are developed in this paper:  
\begin{alphanum}
\item 
A proof that the estimation error $\| \Uptheta_t - \theta^* \|$ is of order $O(\alpha)$, but can be reduced to $O(\alpha^2)$  using a second order linear filter.   

\item
In application to extremum seeking control  
	[an approach to gradient free optimization], 
it is found that the results do not apply because the standard algorithms are not Lipschitz continuous.   A new approach is presented to ensure that the required Lipschitz bounds hold, and from this we obtain stability, transient bounds,   asymptotic bias of order  $O(\alpha^2)$,
 and asymptotic variance of order  $O(\alpha^4)$.    

\item It is in general possible to obtain better than 
$O(\alpha)$  bounds on error in traditional stochastic approximation when there is Markovian noise.

\end{alphanum}

\noindent
AMS-MSC: 68T05 (Primary), 93E35, 49L20 (Secondary)
\end{abstract}

\tableofcontents

\section{Introduction}
\label{s:intro}

The basic problem of interest in this paper is root finding:  find or approximate a solution $\theta^\text{\rm s.t.\,}arp\in\field{R}^d$  to ${\widebar{f}}(\theta^\text{\rm s.t.\,}arp) = 0$,  where  ${\widebar{f}}\colon\field{R}^d \to\field{R}^d$ may not be available in closed form, but noisy measurements of ${\widebar{f}}(\theta)$ are available for any $\theta$.      There is a large library of solution techniques.   

\subsection{Root finding with noisy observations}

The field of stochastic approximation (SA) was born in the work of Robbins and Monro \cite{robmon51a}. The goal is to solve the root finding problem in which the function
 ${\widebar{f}}\colon\field{R}^d\to\field{R}^d$ is expressed as the expectation  
	\begin{equation}
		{\widebar{f}}(\theta) \mathbin{:=} {\sf E}[f(\theta,\qsaprobe)]
		\label{e:barf}
\end{equation}
with $\qsaprobe$ a random vector taking values in $\field{R}^m$,   and $f\colon\field{R}^d\times\field{R}^m\to\field{R}^d$.   

The basic SA algorithm is expressed as the $d$-dimensional recursion,
\begin{equation}
		\theta_{n+1} = \theta_{n}+ \alpha_{n+1} f(\theta_n,\qsaprobe_{n+1}) \, , \qquad n\ge 0,
 \label{e:SA_recur}
\end{equation}
in which $\{\alpha_n\}$ is the non-negative  step-size sequence, and $\qsaprobe_{n+1} \buildrel{\rm d}\over\longrightarrow \qsaprobe$ as $n \to \infty$.  
Writing $\widetilde Xi_n =  f(\theta_n,\qsaprobe_{n+1})  -  {\widebar{f}}(\theta_n)$, the SA recursion can be expressed in the suggestive form,
\begin{equation}
\theta_{n+1} = \theta_{n}+ \alpha_{n+1} \{ {\widebar{f}}(\theta_n) +\widetilde Xi_n  \}\, , \qquad n\ge 0  
\label{e:SA_recurDelta}
\end{equation}
This is interpreted as a noisy Euler approximation of the ODE
\begin{equation}
	\ddt \upvartheta_t  = {\widebar{f}}  ( \upvartheta_t  )
	\label{e:barf_ODE_gen}
\end{equation}
Convergence   theory for SA is based on the recognition that the identity ${\widebar{f}}(\theta^\text{\rm s.t.\,}arp) =0$ means that $\theta^\text{\rm s.t.\,}arp$ is the stationary point for this ODE.   Convergence of the SA recursion is established under the assumption that  	\eqref{e:barf_ODE_gen} is globally asymptotically stable,  along with some regularity conditions on the noise,  and the standard assumptions ensuring success of noiseless Euler approximations \cite{bor20a}.   Note that in early work it was assumed 
that the sequence ${\bf q}saprobe \mathbin{:=} \{ \qsaprobe_n  : n\ge 1\}$ is i.i.d.\  (independent and identically distributed).   It is made clear in  \cite{bor20a} that convergence does not require such strong assumptions;  it is only when we turn to rates of convergence that we require strong assumptions and far more elaborate analysis.  

While there have been exciting advances in SA theory over recent decades,    these advances  do not allow us to escape an unfortunate truth:    in the majority of applications of SA, the mean square error decays no faster than  $1/n$.   \textit{There appears to be no way to avoid the Central Limit Theorem}.
	
While this curse of variance is true in broad generality,   there is a remedy made possible by the flexibility we have in many  application domains, such as optimization and reinforcement learning:  in these applications it is  \textit{we who design the noise} for purposes of ``exploration''.    We obtain much more efficient algorithms by abandoning randomness in the algorithm.

 
	It is also convenient to perform analysis in continuous time, so that the recursion \eqref{e:SA_recur} is replaced with an ordinary differential equation.   
\wham{Quasi-Stochastic Approximation} 
	is a deterministic analog of stochastic approximation.   The  \textit{QSA ODE} is defined by the ordinary differential equation
\begin{equation}
		\ddt\Uptheta_t = a_t f(\Uptheta_t,\qsaprobe_t) 
\label{e:QSA_ODE_gen}
\end{equation}
where $f\colon\field{R}^d\times \field{R}^m \to \field{R}^d$ is continuous.  The $m$-dimensional deterministic continuous time process ${\bf q}saprobe$ is called the \textit{probing signal}.  

Subject to mild assumptions,   it was shown in \cite{laumey22} that the algorithm can be designed so that the  squared error decays at rate $1/t^{4\rho}$ for the vanishing gain algorithm, using
\begin{equation}
a_t = \alpha/(1+t/t_e)^\rho  
\label{e:arho}
\end{equation}
with  $\rho \in (\half, 1)$,   and $\alpha>0$,  $t_e\ge 1$ arbitrary constants.  

The present article focuses on the case   $\rho=0$,  meaning the gain is constant:  $a_t \equiv \alpha$. 
	
	
	
	
 
It is assumed throughout the paper that the  probing signal is a nonlinear function of sinusoids:  
$\qsaprobe_t  = G_0(\qsaprobe _t^0  )$,  with $G_0\colon\field{R}^K \to\field{R}^m$ continuous,  and
\begin{equation}
	\qsaprobe _t^0  =  [  \cos (2\pi [\,  \omega_1 t  +  \phi_1 ] ) ;\,  \dots  ;\,    \cos (2\pi [\,  \omega_K t  +  \phi_K ] )  ]
\label{e:qsaprobe0}
\end{equation}
Motivation for a nonlinearity may be to create rich probing signals from simple ones.   For example,  with $K=1$ we can obtain 
		 $\qsaprobe_t  =    [  \cos (2\pi \omega_1 t ) ; \, \cos (4\pi \omega_1 t ) ;\, \dots  ; \,  \cos (2m\pi \omega_1 t   )  ]$ for any $m\ge 1$,  with  $G_0\colon\field{R} \to\field{R}^m$ polynomial.    
On choosing $G_0$ linear we obtain, for vectors $\{v^i\}\subset \field{R}^m$,
\begin{equation}
		\qsaprobe _t  = \sum_{i=1}^K  v^i   \cos (2\pi [  \omega_i t  +  \phi_i  ] )
\label{e:qSGD_probe}
\end{equation}
 
 
For analysis it is crucial   to abandon $\qsaprobe _t^0  $ in favor of the   $K$-dimensional clock-process   ${\bf P}hi$ with entries $\Phi_t^i =   \exp(2\pi j[\omega_i t + \phi_i])$.
On defining 
 $ G(z) = G_0( (z+1/z)/2  )$, where $1/z \mathbin{:=} (1/z_1,\dots,1/z_K)$ for $z\in \{\field{C}\setminus \{0\}\}^K$, we obtain 
\begin{equation}
\qsaprobe_t  = G( \Phi_t) 
\label{e:qsaprobeG}
\end{equation}
The function $G$ will be assumed continuous on the restricted domain $ \{\field{C}\setminus \{0\}\}^K$.  	

The clock-process   ${\bf P}hi$ is Markovian.    It can also be expressed as the solution to the linear system,
\begin{equation}
		\ddt \Phi_t = W \Phi_t, \quad \text{with } W =2\pi j \hbox{\rm diag\thinspace}(\omega_i) 
		\label{e:LinearMarkovQSA}
\end{equation}
It evolves on a compact subset of Euclidean space denoted $\Upomega \subset\field{C}^K$,  and it is evident that the uniform distribution on $\Upomega$ is the unique invariant measure for ${\bf P}hi$;   this is denoted $ \uppi$. 

The \textit{mean vector field} associated with \eqref{e:QSA_ODE_gen} is defined as the expectation,
\begin{equation}
		{\widebar{f}}(\theta) =   {\sf E}[ f(\theta, G(\Phi))  ]\,,\qquad \theta\in\field{R}^d\,,  \qquad \textit{in which  $ \Phi\sim \uppi$. 
}
\label{e:barfQSA}
\end{equation}

\notes{
}

\subsection{Perturbative mean flow}

The central conclusion on which most of the concepts in this paper are based upon is the following representation for the QSA ODE, which holds under smoothness conditions on $f$ and $G_0$.

\wham{Perturbative Mean Flow:}  The  solution to the QSA ODE admits the exact description
 \begin{equation} 
		\begin{aligned}
			\ddt \Uptheta_t   & =  \alpha[ {\widebar{f}} (\Uptheta_t)  
						  -  
						\alpha \widebar{\Upupsilon}_t
						  +   {\cal W}_t]
			\\
			& \qquad    {\cal W}_t =   \alpha^2 {\cal W}_t^0   +  \alpha \frac{d}{dt} {\cal W}_t^1 +    \frac{d^2}{dt^2} {\cal W}_t^2
		\end{aligned} 
\label{e:BigGlobalODE}
\end{equation}	
The details: 
\begin{romannum}
\item
The deterministic processes  $\{{\cal W}_t^i         : i=0,1,2\}$   have explicit representations,  
	 given in \eqref{e:AllTheNoise0}--\eqref{e:AllTheNoise2},
	 as smooth functions of the larger state process ${\bf P}si =  (\bm{\Uptheta},{\bf P}hi)$.
		
		\item The function $\widebar{\Upupsilon}_t = \widebar{\Upupsilon}(\Uptheta_t)  $ appears when there is multiplicative noise in the QSA ODE.    It can contribute significantly to the estimation error $\|\Uptheta_t - \theta^*\|$, resulting in large bias and variance. Fortunately, it can be eliminated with careful design.
		 
\end{romannum}
The implications of the perturbative mean flow (or P-mean flow) representation to algorithm design are focuses of this paper.   

However,  there is one catch:  while the P-mean flow representation holds in broad generality, we cannot use it to establish stability (in the sense of ultimate boundedness) since to-date we have not found global Lipschitz bounds on the functions  
$\{ \widebar{\Upupsilon},  {\cal W}^i        : i=0,1,2\}$.      

\wham{Algorithm Performance:} The pair ${\bf P}si =  (\bm{\Uptheta},{\bf P}hi)$ is a Feller Markov process.   Existence of an invariant measure $\upvarpi \sim {\bf P}si  $ is guaranteed for QSA whenever the sample path $\bm{\Uptheta}$ is bounded from at least one initial condition. It is not known if $\upvarpi$ is unique, so bias and variance are defined in terms of sample path averages
\begin{equation}
	\begin{aligned}
		\beta_\Uptheta  &= \mathop{\rm lim{\,}sup}_{T\to\infty} {\sf B}igl\| \frac{1}{T}\int_0^T [\Uptheta_t - \theta^*] \, dt  {\sf B}igr\| \,,
		\\
		\sigma^2_\Uptheta  &=   {\sf B}igl(  \mathop{\rm lim{\,}sup}_{T\to\infty}    \frac{1}{T}\int_0^T \| \Uptheta_t -\theta^\text{\rm s.t.\,}arp \|^2 \, dt  {\sf B}igr) 
		-   \|\beta_\Uptheta\|^2  \,.
	\end{aligned} 
	\label{e:biasdef}
\end{equation}

The standard $L_p$ norms are also considered in their sample path forms: 
\begin{equation}
	\begin{aligned}
		\| \Uptheta -\theta^*\|_{L_1} = \mathop{\rm lim{\,}sup}_{T\to\infty}    \frac{1}{T}\int_0^T \| \Uptheta_t -\theta^*\| \, dt    \, ,
		\\
		\| \Uptheta -\theta^*\|_{L_2} = \hbox{\rlap{$\sqcap$}$\sqcup$}rt{\mathop{\rm lim{\,}sup}_{T\to\infty}    \frac{1}{T}\int_0^T \| \Uptheta_t -\theta^*\|^2 \, dt }   \, .
	\end{aligned}
	\label{e:norms}
\end{equation}
The $L_1$ norm is also referred to as the absolute average error (AAD). These quantities are related via 
\begin{equation}
	\| \Uptheta -\theta^*\|_{L_1}   \le  \| \Uptheta -\theta^*\|_{L_2}  = \hbox{\rlap{$\sqcap$}$\sqcup$}rt{   \sigma^2_\theta  +   \beta_\Uptheta^2  } \, .
	\label{e:cov_from_bias}
\end{equation}
We say that the  AAD is  $O(\eta_\alpha)$ if there is a finite constant $B_0$ such that
 $\| \Uptheta -\theta^*\|_{L_1} \le B_0 \eta_\alpha$   for all $\alpha>0$ in a neighborhood of the origin; $\eta_\alpha$ will be polynomial in $\alpha$.     

We find that the AAD of QSA may be order $O(\alpha)$ even when $\widebar{\Upupsilon}\equiv 0$.    \Cref{fig:QSArep} illustrates how linear filtering techniques can reduce this AAD to $O(\alpha^2)$.  The concept is simple:  if $\ddt u_t$ is passed through a low pass filter with bandwidth $\gamma>0$,   then the output $y_t$  satisfies $|y_t|   \le \gamma \|{\bf m}u\|_\infty + o(1)$.
In general we require a second order filter to attenuate the term   $\frac{d^2}{dt^2} {\cal W}_t^2$
 appearing in \eqref{e:BigGlobalODE}.
 
 
 We are also interested in the  \textit{target bias},  defined as the sample average
\begin{equation}
	b_{\widebar{f}}two \mathbin{:=}       \lim_{T\to\infty} \frac{1}{T} \int_0^T   {\widebar{f}} (\Uptheta_t)   \, dt
	\label{e:meanTargetBias}
\end{equation}

	

	
	\begin{figure}
\caption{Perturbative mean flow and its implication to filter design. }
\label{fig:QSArep}
\end{figure}
	
	

\subsection{Stability theory for QSA}
\label{s:QSAstable}

We opt for the usual control-theoretic definition for the full state process:
We say that ${\bf P}si  = (\bm{\Uptheta},{\bf P}hi)$ is \textit{ultimately bounded} if there is a fixed constant $B<\infty$ satisfying the following:   for each initial condition $\Uptheta_0 =\theta$,  $\Phi_0 = z$,  there is a finite time $t_0 = t_0(\theta,z)$ such that
\begin{equation}
	\|  \Uptheta_t\|   \le B\,, \qquad t\ge t_0
	\label{e:UB}
\end{equation}
Two criteria to establish this property are available based entirely on the mean vector field ${\widebar{f}}$:    
\begin{romannum}
\item[1.]   \textbf{Lipschitz Lyapunov Function:}      For a $C^1$ function $V\colon\field{R}^d\to\field{R}_+$ and a constant $\delta_0>0$,
\begin{equation} 
	\frac{d}{dt} V(\upvartheta_t)  \leq  -  \delta_0    V(\upvartheta_t)    \,, \qquad   \textit{when $\|\upvartheta_t\|\ge \delta_0^{-1}$}
	\label{e:ddt_bound_LyapfunTmp}
\end{equation}
This is equivalently expressed $\nabla  V(\theta) \cdot {\widebar{f}}(\theta)  \le  -  \delta_0    V(\theta)$ for all  $\|\theta\|\ge \delta_0^{-1}$.   This implies a similar bound for the QSA ODE under very general conditions.   Most crucial is that $f$,  ${\widebar{f}}$,   and $V$ each satisfy a global Lipschitz bound.   

\item[2.]   \textbf{ODE@$\infty$:}     This is the vector field obtained by scaling 
\begin{equation}
	{\widebar{f}}_\infty (\theta) 
	\mathbin{:=}
	\lim_{r \to \infty }   {\widebar{f}}^r (\theta)
	\quad 
	\textit{with}
	\quad {\widebar{f}}^r (\theta)= r^{-1} {\widebar{f}} ( r \theta)  \,,\qquad \theta\in\field{R}^d 
	\label{e:barfinfty_QSA}
\end{equation} 
If this exists, and if the ODE with this vector field is locally asymptotically stable, then we establish a relaxation of	\eqref{e:ddt_bound_LyapfunTmp} that also implies that ${\bf P}si$ is ultimately bounded. 
\end{romannum}

The first criterion might look odd to those accustomed to quadratic Lyapunov functions.     Suppose that $V_1$ is quadratic, of the form $V_1 (\theta) = \theta^{\rm tr\, }anspose P \theta$ with $P>0$  $(d\times d)$, and solves 
$\nabla  V_1(\theta) \cdot {\widebar{f}}(\theta)  \le  -  \delta_1    V_1(\theta)$ for some $\delta_1>0$ and all $\|\theta\| \ge \delta_1^{-1}$.   It follows from the chain rule that 	\eqref{e:ddt_bound_LyapfunTmp} holds using $V = \hbox{\rlap{$\sqcap$}$\sqcup$}rt{1+ V_1}$, and $\delta_0 \in (\delta_1/2,1)$.
 

The second approach is based on the very notion of ultimate boundedness:  if this condition fails to hold, then by definition the state blows up.   Hence we look at the evolution on a compressed spatial scale.   
An example is shown in \Cref{fig:ODEinfty},  with details postponed to \Cref{s:ESC}.

 
 \begin{wrapfigure}[16]{r}{0.32\textwidth}

 	\centering
	

	
	\includegraphics[width=0.95\hsize]{./Figures/ODE-Inf1e10newVanishingJustOne.pdf}
	\caption{Solutions to \textit{extremum seeking control} for the Rastrigin objective, starting from very large initial conditions}
	\label{fig:ODEinfty}
\end{wrapfigure} 

\textit{The ODE@$\infty$ is also a useful technique to establish ultimate boundedness of the mean flow}:   consider a large initial condition,  with magnitude $r = \| \upvartheta_0\|$.    The scaled state process $\upvartheta^r_t  = r^{-1} \upvartheta_t$ solves the scaled ODE,
\begin{equation}
\ddt \upvartheta^r_t    =  r^{-1} {\widebar{f}}(\upvartheta_t)     = {\widebar{f}}^r(\upvartheta_t^r)   \,,\qquad   \|  \upvartheta^r_0 \|=1
\label{e:ODEinf-r}
\end{equation}
This simple conclusion leads to a proof of ultimate boundedness for the mean flow, following two observations:

\wham{(i)}
 The convergence in	\eqref{e:barfinfty_QSA}
is uniform on compact subsets of $\field{R}^d$  (an application of Lipschitz continuity).

\wham{(ii)}
 if ${\widebar{f}}_\infty$ is locally asymptotically stable, then it must be globally exponentially asymptotically stable, with the origin the unique stationary point, $ {\widebar{f}}_\infty({\sf Z}ero) ={\sf Z}ero$.
 

 Based on these conclusions we obtain semi-exponential stability, in the following sense:  for some $\delta_e>0$ and initial conditions satisfying $\| \upvartheta_0\| \ge \delta_e^{-1}$,  solutions to 	\eqref{e:barf_ODE_gen}  satisfy
 \begin{equation}
V(\upvartheta_t)  \le V(\upvartheta_0) e^{-\delta_e t} \,,\qquad 0\le t\le  T(\upvartheta_0) 
\label{e:ExpStableAlmost}
\end{equation}
where $T$ is the first entrance time: $
T(\upvartheta_0) \mathbin{:=}  \min \{  t\ge 0 : \| \upvartheta_t\|  \le \delta_e^{-1}  \} $.

The bound \eqref{e:ExpStableAlmost} is extended to the QSA ODE in  \Cref{t:ODEstate_is_bdd},  
under either of these two stability criteria.  

 
 

\subsection{Bridge building}      
\label{s:ESCetc}

As motivation, to provide examples, 
and to provide a short literature survey,
we present a brief survey of gradient free optimization (GFO) and how it is related to QSA and
extremum seeking control (ESC).   
The definition of ESC is postponed to \Cref{s:ESC}, based on interpretation of its standard block diagram description shown in \Cref{fig:ESC}.

		
	\begin{figure}
\caption{A Typical Architecture for Extremum Seeking Control for Gradient-Free Optimization}
\label{fig:ESC}
\end{figure}
	
	
\subsubsection{SPSA}

\begin{subequations}

Within the area of GFO two algorithms of Spall are most easily described and justified.   Each are SA algorithms of the form \eqref{e:SA_recur}, differing only in the definition of $f$:
\begin{align}

\text{\sf 1SPSA:} \quad 
f(\theta,\qsaprobe) &
=
		- \frac{1}{\varepsilon} 
		\qsaprobe \Upgamma(\theta + \varepsilon \qsaprobe)
\label{e:1qSPSA}
		 \\

\text{\sf 2SPSA:} \quad 
f(\theta,\qsaprobe) &=
		-
		\frac{1}{2\varepsilon}   \qsaprobe\bigl[\Upgamma(\theta + \varepsilon \qsaprobe) -\Upgamma(\theta - \varepsilon \qsaprobe) \bigr]   
\label{e:2qSPSA}
\end{align}
It is assumed that the probing sequence $\{\qsaprobe_n\}$  is i.i.d., and in much of the theory it is assumed that the entries of $\qsaprobe_n$ have support in $\{-1,1\}$.   
In either case, the mean vector field is defined via \eqref{e:barf}:
\[
 {\widebar{f}}(\theta) \mathbin{:=}   {\sf E}[  f(\theta,\qsaprobe_{n+1})  ] \,,  \qquad \theta\in\field{R}^d  
\]
This is independent of $n$ under the assumption that the probing sequence is i.i.d..   

\label{e:SPSA_Gen}
\end{subequations}

\begin{subequations}

Motivation for either is usually based on a first-order Taylor series approximation.  The following exact representations follow from   the Fundamental Theorem of Calculus:
\begin{align}

\text{\sf 1SPSA:} \quad 
f(\theta,\qsaprobe) &=
			 		-  \qsaprobe\qsaprobe^{\rm tr\, }anspose  
						 \frac{1}{\varepsilon}  \int_{0}^\varepsilon   \nabla \Upgamma\,( \theta + t \qsaprobe)  \, dt   
						  -  \frac{1}{\varepsilon}  \qsaprobe \Upgamma( \theta)  
\label{e:1qSPSA_FTC}
		 \\

\text{\sf 2SPSA:} \quad 
f(\theta,\qsaprobe) &
		=  -   \qsaprobe\qsaprobe^{\rm tr\, }anspose  \frac{1}{2\varepsilon}  \int_{-\varepsilon}^\varepsilon   \nabla \Upgamma\,( \theta + t \qsaprobe)  \, dt  
\label{e:2qSPSA_FTC}
\end{align}
These representations in terms of average gradients provide the first step in the proof of the following:

 
\label{e:SPSA_Gen_FTC}
\end{subequations}

 

  \begin{proposition}
\label[proposition]{t:SPSAbias}
Suppose that the $d$-dimensional sequence ${\bf q}saprobe$ is i.i.d.\ with zero mean,  and   that $ \qsaprobe_{n+1} \buildrel{\rm dist}\over = - \qsaprobe_{n+1}$  (equality in distribution).   Then the mean vector fields for 1SPSA and 2SPSA coincide.

Suppose in addition that the objective $\Upgamma$ is continuously differentiable ($C^1$)  and strongly convex in a neighborhood of the optimizer $\theta^\text{\tiny\sf opt}$.
Then,   there exists $\varepsilon^0>0$ such that the function ${\widebar{f}}$ has a  root $\theta^\text{\rm s.t.\,}arp$ for $0<\varepsilon\le \varepsilon^0$, and for either 1SPSA or 2SPSA,  
\begin{equation}
\| \theta^\text{\rm s.t.\,}arp - \theta^\text{\tiny\sf opt} \| = O(\varepsilon^2)  \qquad \textit{and}\qquad 
\Upgamma( \theta^\text{\rm s.t.\,}arp ) -  \Upgamma(\theta^\text{\tiny\sf opt} )  = O(\varepsilon^4)  
\label{e:SPSAbias}
\end{equation}
\qed
\end{proposition}

SPSA is a perfect example of a successful algorithm that can be improved using deterministic exploration.

\subsubsection{ESC and QSA}   
\label{s:ESC_QSA}

Spall's SPSA recursions admit natural analogs as QSA ODEs, with vector fields
\begin{subequations} 
\begin{align}

\text{\sf 1qSGD:} \quad 
f(\Uptheta_t,\qsaprobe_t)  &=
- \frac{1}{\varepsilon} 
		\qsaprobe_t \Upgamma(\Uptheta_t + \varepsilon \qsaprobe_t)
\label{e:1qSGD}
		 \\

\text{\sf 2qSGD:} \quad 
f(\Uptheta_t,\qsaprobe_t) &=
-
		\frac{1}{2\varepsilon}   \qsaprobe_t\bigl[\Upgamma(\Uptheta_t + \varepsilon \qsaprobe_t) -\Upgamma(\Uptheta_t - \varepsilon \qsaprobe_t) \bigr]   
\label{e:2qSQD}
\end{align}
in which {\sf qSGD} refers to \textit{quasi-stochastic gradient descent}.   In examples that follow we typically take the sinusoidal probing signal \eqref{e:qSGD_probe}.   

\label{s:qSGD}
\end{subequations} 

In \Cref{s:ESC} we show that 1qSGD is a simple instance of extremum seeking control.

And we also find a challenge:  in both domains,  ESC and 1qSGD,  convergence theory is limited because the vector field $f$ is not Lipschitz continuous in its first variable.   \textit{Lipschitz continuity is a  critical component of theory for both SA and QSA}.   
A remedy is introduced in this paper, through the introduction of a state-depending probing gain.   Details are found in    \Cref{s:SDprobe}.

\subsubsection{Averaging}

The theory of averaging for variance reduction is well-established in the SA literature,  with application to optimization  \cite{frapol20,bacmou11,moujunwaibarjor20} and to  GFO  \cite{spa03}.   The idea is very simple:  given the estimates $\{\theta_n : 0\le n  \le N\}$ from any SA recursion,  obtain a new estimate via  
\begin{equation}
\theta^{\text{\tiny\sf PR}}_N \mathbin{:=} \frac{1}{N-N_0}  \sum_{n=N_0+1}^N   \theta_n
\label{e:PR}
\end{equation}
where $1\le N_0 <N$ is chosen to discard large transients.  
This apparently simple ``hack'' provides enormous benefit:  under mild conditions, the variance of this estimate of $\theta^\text{\rm s.t.\,}arp$ is minimal in a strong sense  \cite{bor20a,kusyin97}.   The technique goes by the name of Polyak-Ruppert averaging in honor of the creators of the technique and original analysis  \cite{rup88,pol90,poljud92} (note Polyak's  contributions in \cite{pol90} prior to his collaboration with Juditsky).
The article \cite{dipren97} is an early application of these ideas to accelerate GFO.

This is similar to the low pass filter shown in \eqref{fig:ESC}, generating  $\UpthetaF_t$ from $\Uptheta_t$.
We will find better options for fixed gain QSA   in \Cref{s:filterDesign}.

 \subsubsection{Multiplicative noise}
 
One objective of this paper is to bring to light the  challenges posed by multiplicative noise.     Many optimistic results in prior work consider additive noise models for which the nuisance term $\Upupsilon$ is not present:
\begin{romannum}
\item
Estimation error bounds of order $O(\alpha^2)$  are obtained in \cite{arikrs02} for  ESC when the objective is quadratic.    
 
\item
There are recent optimistic results in the stochastic approximation literature for  fixed gain algorithms \cite{moujunwaibarjor20,durmounausamscawai21}.   Through design it is possible to obtain zero asymptotic bias and optimize variance (via averaging).
These results are obtained for models with additive white noise.   
 \end{romannum}
 For QSA we obtain  both $O(\alpha^2)$ bounds for AAD and target bias through careful design of filters and the probing signal.  
 Theory shows that we  cannot expect    $O(\alpha^2)$ AAD for general models, unless a filter is  designed following the specifications in \Cref{t:Couple_theta}.
 
In stochastic settings with Markovian noise it is not at all clear how to obtain such strong conclusions---averaging cannot remove the inherent $O(\alpha)$ AAD that is a consequence of statistical memory.  Details are found in  \Cref{s:SAbias}.
 
 

  
  

  
 
\wham{Overview}
The remainder of the body of the paper is organized into five additional sections:     
 The main technical results of the paper are contained in \Cref{s:QSAtheory}, including a proof of the P-mean flow representation \eqref{e:BigGlobalODE}.    
 Implications to extremum seeking control are summarized in \Cref{s:ESC}, 
with examples in \Cref{s:Experiments},  
  and conclusions and directions for future research are summarized in
\Cref{s:conc}. 
A literature survey is contained in \Cref{s:lit}, and proofs of some technical results are contained in the Appendix.  
 

\section{QSA Theory}
\label{s:QSAtheory}

The theory described in this section is based on the fixed-gain QSA ODE with $a_t\equiv \alpha>0$:    
\begin{equation}
	\ddt\Uptheta_t = \alpha f(\Uptheta_t,\qsaprobe_t)
\label{e:QSAgen_alpha}
\end{equation}
Convergence of $\{\Uptheta_t\}$ to $\theta^*$  cannot be expected;  instead with careful design of the algorithm we obtain bounds on asymptotic bias of order $O(\alpha^2)$ and variance of order $O(\alpha^4)$.   

	

\subsection{Markovian framework for analysis}
\label{s:markov}

	In the theory of SA the  stochastic sequence $\{\qsaprobe_{n+1}\}$ appearing in \eqref{e:SA_recur} is not always assumed to be i.i.d..
Sharp results on asymptotic variance are available when it can be expressed as a function of a well behaved Markov chain.  
It may come as a big surprise to learn that the techniques extend to QSA analysis, and the application of techniques rooted in Markov chain theory are the only tools available to obtain \eqref{e:BigGlobalODE} and its progeny. 

	

	For the fixed gain QSA ODE it is clear that the pair ${\bf P}si =  (\bm{\Uptheta},{\bf P}hi)$ is itself the state process for a time-homogeneous dynamical system.      It can also be regarded as a Feller Markov process on the closed state space $\Uppi =  \field{R}^d \times \Upomega$. 
	Provided  ${\bf P}si $  is ultimately bounded,   it admits at least one invariant measure \cite[Thm.~12.1.2]{MT}, which we denote $\upvarpi$.
	This is valuable conceptually, and in terms of notation.   We can unambiguously write ${\widebar{g}} \mathbin{:=} {\sf E}[ g(\Uptheta,\qsaprobe)] $ for a function $g\colon\field{R}^d\times \field{R}^K\to\field{R}$,  where the expectation is with respect to $\Psi\sim\upvarpi$.

The following simple result will aid in interpretation of other expectations.   
For any function $h\colon\Uppi\to\field{R}$ that is continuously differentiable, define the continuous function $\mathcal{D} h$ via
\begin{equation}
 \mathcal{D} h\, (\theta,z)  \mathbin{:=}  \partial_\theta h\,  (\theta,z )   \cdot f (\theta, G(z))    
	+   \partial_z   h \, (\theta, z) \cdot W  z    \,,\qquad (\theta, z) \in\Uppi  
\label{e:diffGen}
\end{equation}
In Markov terminology, the  functional $\mathcal{D}$ is known as the \textit{differential generator} for ${\bf P}si$.

\begin{proposition}
\label[proposition]{t:QSAgenerator}
Let $h\colon\Uppi\to\field{R}$ be a $C^1$ function, and denote $g=\mathcal{D} h$.  The following then hold: 
\begin{romannum}
\item   
$ g(\Psi_t)  =\ddt h(\Psi_t) $.

\item
Suppose that for some initial condition $\Psi_0$, the resulting trajectory $\{ \Psi_t : t\ge 0\}$ is uniformly bounded.   Then,
\[
\lim_{T\to\infty}   \frac{1}{T}   \int_0^T  g(\Psi_t)  \, dt   = {\sf Z}ero
\]

\item
Suppose that an invariant measure  $\upvarpi$ exists with compact support.
Then ${\sf E}[g(\Psi)]=0$  when $(\Uptheta,\Phi) = \Psi\sim \upvarpi$.   
In particular,  on choosing $h(\theta,z) \equiv \theta$,
\begin{equation}
 			{\sf E}[f(\Uptheta, G(\Phi))]= {\sf Z}ero
 \label{e:fMeanZero}
\end{equation}
\qed
\end{romannum}
\end{proposition}

\begin{proof}
Part (i) follows from the chain rule,  and (ii) from the fundamental theorem of calculus:
\[
  h(\Psi_T)  =  h(\Psi_0)     +   \int_0^T  g(\Psi_t)  \, dt    \,,\qquad T>0
\]
Part (iii) also follows from the fundamental theorem of calculus, using the expression above,  but with   a random initial condition.
Choose  $\Psi_0 =  (\Uptheta_0, \Phi_0)$ to be a $\Uppi$-valued random variable with distribution $ \upvarpi$.   We then have $\Psi_t\sim \upvarpi $ for all $t$.   Taking expectations of each side  completes the proof:
 \[
 {\sf E}[  h(\Psi_0)  ]   =  {\sf E}[  h(\Psi_0)  ]     + T  {\sf E}[  g(\Psi_0) ]
\] 
\end{proof}
	

Our main motivation for a Markovian framework comes from the application of solutions to \textit{Poisson's equation}:  
	for functions   $g,\, {\hat g}\colon\Upomega\to\field{R}$, this is expressed
\begin{equation}
		{\hat g}(\Phi_0) = 
		{\hat g}(\Phi_T) +   \int_0^T   \tilde g(\Phi_t)\, dt\,,\qquad  T\ge 0\,, 
		\label{e:PoissonDefn}
\end{equation}
	where $\tilde g(z) = g(z) -{\widebar{g}}$ for $z\in\Upomega$,  with ${\widebar{g}}$ the steady-state mean.  We say that ${\hat g}$ is the \textit{solution} to Poisson's equation, and that $g$ is the \textit{forcing function}.    If ${\hat g}$ is continuously differentiable, then from \eqref{e:LinearMarkovQSA} we obtain $\nabla^{\rm tr\, }anspose {\hat g}(z) W z = - \tilde g(z)$ for $z\in\Upomega$.
	
	Justification for smooth solutions to Poisson's equation 
	is
	 obtained in \cite{laumey22} when  the forcing function   is analytic,  and subject to assumptions on the frequencies defining ${\bf P}hi$. 	
	 A review is contained in the Appendix, and a summary of required assumptions in Assumption (A0) below.
	\begin{romannum}
		\item[]\textit{Please note:}
		\item
		The solution ${\hat g}$ is not unique.   We always normalize so that ${\sf E}[{\hat g}(\Phi)] =0$.
		
\item
In the applications of Poisson's equation used in this paper, it is frequently true that the function $g$ depends on $\Phi_t$ only through $\qsaprobe_t$.
\textit{This is this is not in general true for ${\hat g}$}. 
		
		\item 
		Finally, on notation:  we often  write ${\hat g}_t$ instead of  $ {\hat g}(\Uptheta_t,\Phi_t)$  (similar notation for other functions of ${\bf P}si$).
		
	\end{romannum}
	
	
	
	
	We will not require solutions for the joint process ${\bf P}si$, but   require a slight abuse of notation:  for a vector-valued function on the joint state space  $g\colon\Uppi\to \field{R}$,  we denote by ${\hat g} (\theta,\,\cdot\,)$   the solution to Poisson's equation with $\theta\in\field{R}^d$ fixed:
	\begin{equation}
		{\hat g}(\theta,  \Phi_0) = 
		{\hat g}(\theta, \Phi_T) +   \int_0^T   \tilde g(\theta,  \Phi_t)\, dt\,,\qquad  T\ge 0\, .
		\label{e:PoissonDefnAbuse}
	\end{equation}
	This is applied with $g=f$  [the QSA ODE vector field \eqref{e:QSA_ODE_gen}];
	the solution ${\hat f}$ plays a part in the    definition of $\widebar{\Upupsilon}$.
	
		
Further implications require a few assumptions and notation:
	
\wham{Assumptions}

The following assumptions are in force throughout much of the paper.

The first assumption sets restrictions on frequencies:
we begin with a  ``basis'' of irrationally related frequencies specified as follows:
For an integer $K_0$,  frequencies $\{\omega_1^0\,,\dots\,,  \omega_{K_0}^0\}$ are chosen of the following form:

\begin{equation}
\begin{aligned}
&\omega_i^0   = \log(a_i/b_i) > 0\,, \ \  1\le i\le K_0\,, \ \ \textit{where   $\{a_i,b_i\}$ are positive integers.   }
   \\
            &  \textit{$\{\omega_i^0 : 1\le i\le K_0\}$ are linearly independent over the field of rational numbers.}
\end{aligned} 
\label{e:logFreq}
\end{equation}

 
\begin{subequations}

\wham{(A0a)}   
$\qsaprobe_t  = G_0(\qsaprobe _t^0  )$ for all $t$,   with $\qsaprobe _t^0$ defined in \eqref{e:qsaprobe0}.  The function 
$G_0\colon\field{R}^K \to\field{R}^m$ is assumed to be analytic, with the coefficients in the Taylor series expansion for $G_0(\qsaprobe _t^0  )$ absolutely summable.     

\wham{(A0b)}   
The frequencies $\{\omega_1\,,\dots\,,  \omega_K\}$   satisfy one of the two conditions:

\begin{flushleft} 
\textbf{\textit{Periodic}:}  
For a fixed but arbitrary frequency $\omega_1>0$ and an  increasing sequence of positive integers $\{n_i\}$,
\begin{equation}
\omega_i = n_i \omega_1\,,\quad 2\le i\le K
\label{e:freqPeriodic}
\end{equation}

\textbf{\textit{Transcendental Almost Periodic}:}      
The frequencies $\{\omega_i : 1\le i\le K\}$ are positive,  distinct,  and in the positive integer span of the $\{ \omega_i^0 : 1\le i\le K_0 \}$:
\begin{equation}
\omega_i = \sum_j n_{i,j} \omega_j^0
\quad
\textit{
 for integers  $\{n_{i,j} \} \subset \field{Z}_+$.}
\label{e:logFreqSpan}
\end{equation}
\end{flushleft}
\label{e:TwoFreqChoices}
\end{subequations}

\begin{subequations}

\wham{(A1)}    
The functions ${\widebar{f}}$ and $f$ are
	 Lipschitz continuous:  for   a constant $L_f  <\infty$,
\begin{equation}
\begin{aligned} 
		\|{\widebar{f}}(\theta') - {\widebar{f}}(\theta)\| &\le L_f \|\theta' - \theta\|, 
		\\
		\|f(\theta',\qsaprobe) - f(\theta,\qsaprobe)\| 
		+
		\|f(\theta,\qsaprobe') - f(\theta,\qsaprobe)\| 
		&\le L_f \bigl[ \|\theta' - \theta\|  +   \|\qsaprobe' - \qsaprobe\|  \bigr]
\end{aligned} 
\label{e:fLip}
\end{equation}
for all $ \theta' \,, \theta\in\field{R}^d$,
$\qsaprobe\,, \qsaprobe' \in\field{R}^m $.

\wham{(A2)}   
The vector fields $f$ and ${\widebar{f}}$   are each twice continuously  differentiable,  with derivatives denoted 
\begin{align}
A(\theta, z) & = \partial_\theta f \, (\theta, z)  \,,
\qquad
 {\bar{A}}(\theta) = \partial_\theta {\widebar{f}}\, (\theta)
		\label{e:barfDer}
\end{align}

\wham{(A3)} 
There exists a  solution ${\hat f}$ to Poisson's equation that is continuously differentiable on $\Uppi$,  and
normalized so that $ {\sf E}[ {\hat f}(\theta,  \Phi) ] = {\sf Z}ero$ for each $\theta$,  with $ \Phi\sim \uppi$.   Its Jacobian with respect to $\theta$ is denoted
\begin{equation}
\widehat A(\theta ,z)  \mathbin{:=}  \partial_\theta {\hat f}\, (\theta, z)
\label{e:der_haA}
\end{equation}
Moreover,  there are continuously differentiable solutions to Poisson's equation for each of the two forcing functions ${\hat f}$ and $\Upupsilon$,  with
	\begin{equation}
		\Upupsilon(\theta,z) =   - \widehat A   (\theta,z )    f (\theta, G(z))    
\label{e:Upupsilon}
\end{equation}
The solutions are denoted, respectively,    ${\hat h}af$  and ${\widehat U}pupsilon$:
\begin{align}
{\hat h}af(\theta,  \Phi_{t_0} )  &= 
			\int_{t_0} ^{t_1}  {\hat f}( \theta ,\Phi_t)    \, dt    + {\hat h}af( \theta,   \Phi_{t_1} )  
\nonumber
   \\
{\widehat U}pupsilon(\theta,  \Phi_{t_0} )  &= 
			\int_{t_0} ^{t_1}  [ \Upupsilon( \theta ,\Phi_t)-  \widebar{\Upupsilon}(\theta ) ]     \, dt    + {\widehat U}pupsilon( \theta,   \Phi_{t_1} )  			 \,,   && 0\le t_0\le t_1
			\label{e:haEverythingElse}
\\
\textit{with}\quad 
\widebar{\Upupsilon}(\theta ) & =     {\sf E}[   \Upupsilon(\theta ,\Phi )  ]
				   = -
				\int_\Upomega   \widehat A(\theta ,z )      f(\theta,G(z))   \,  \uppi(dz) \,,  && \theta\in\field{R}^d 
\nonumber
\end{align}
They are also normalized: $    {\sf E}[ {\hat h}af(\theta,  \Phi) ] =    {\sf E}[ {\widehat U}pupsilon(\theta,  \Phi) ] =   {\sf Z}ero$ for each $\theta$. 
\qed
	\end{subequations}	

\smallbreak

Assumption~(A0b) is justified in \Cref{s:BakerPoisson}, but deserves some explanation here.
First, it would appear that \eqref{e:freqPeriodic} is a special case of \eqref{e:logFreqSpan} in which $K_0 =1$.  Note however that no restriction is placed on $\omega_1$ in \eqref{e:freqPeriodic}, except that it is strictly positive.   

The assumption is imposed for two reasons:  
\begin{romannum} 
\item   It is a sufficient condition ensuring that $\widebar{\Upupsilon}\equiv {\sf Z}ero$,   assuming that Assumption~(A3) holds so that  $\Upupsilon$ is well defined.    See the example following \Cref{t:Couple_theta} to see that constraints on the frequencies are indeed needed for this. 

\item  If   $f(\theta, G(z))$ is an analytic function of $(\theta,z)$   on an appropriate domain, then  Assumption~(A3) holds.
\end{romannum}
Claim (i) is an application of \Cref{t:orthCorBigPre};
this result doesn't require the full (A0b)---it is sufficient that the frequencies are distinct.   
Claim~(ii) follows from \Cref{t:hah}.

As for the remaining assumptions, we require  (A1)--(A3) and further assumptions to obtain sharp bounds on bias and variance. 
In particular, the functions introduced in (A3) determine all of the terms in  \eqref{e:BigGlobalODE}, with  $\widebar{\Upupsilon}$ is defined in \eqref{e:haEverythingElse}.

The assumptions on $G_0$ in Assumption~(A0a) are imposed so that the probing signal is almost periodic \cite{amepro13}.  
This combined with the Lipschitz conditions in  (A1) imply a uniform version of the Law of Large Numbers:
\begin{proposition}
\label[proposition]{t:uniLLN}
Under (A0a) and (A1) the uniform Law of Large Numbers holds:
\begin{equation}
\lim_{T\to\infty }   \sup_{\theta, z_0} \frac{1}{ 1 +  \| \theta\|  }  \left |
\frac{1}{T }\int_0^T  \bigl[f(\theta,\qsaprobe_t) - {\widebar{f}}(\theta) \bigr]  \, dt
				\right|   =0
\label{e:LLNf}
\end{equation}
where the supremum is over  $\theta_0\in\field{R}^d$,
			  and   initial conditions $\Phi_0 = z_0\in\Upomega$.
\qed
\end{proposition}

 

 
    

\subsection{Three steps towards the perturbative mean flow }
 
 The notation in \eqref{e:Upupsilon}  is complex, which brings us to the following compact alternative:
	For $h\colon \Uppi \to\field{R}^d$ that is continuously differentiable in its first variable we let $D^f h$ denote its directional derivative in the direction $f$,  where $f$ is the vector field for the QSA ODE:	
\begin{equation}
[ D^f h ] (\theta,z )  = \partial_\theta h\,  (\theta,z )   \cdot f (\theta, G(z))  \,, \qquad (\theta,z )  \in \Uppi
		\label{e:Dfg}
\end{equation} 
This defines a component of the differential generator \eqref{e:diffGen};
 recall from  \Cref{t:QSAgenerator}  that  $g=\mathcal{D} h$ implies that  $ g(\Psi_t)  =\ddt h(\Psi_t) $.

 

The following companion to  \Cref{t:QSAgenerator} will be used in the following:
	for any smooth solution to Poisson's equation with forcing function $h\colon\Uppi\to\field{R}$,
	\begin{align}
		\ddt {\hat h}(\Uptheta_t,z )  &=  \alpha [D^f {\hat h} ] (\Uptheta_t,z )  \,,   \ z\in\Upomega  
		\label{e:chainRule_f}
		\\
		\ddt {\hat h}(\Uptheta_t, \Phi_t )  &=  \alpha [D^f {\hat h} ] (\Uptheta_t, \Phi_t )     
		-  [ h(\Uptheta_t, \Phi_t )  -  {\overline {h}}(\Uptheta_t) ]   
		\label{e:chainRule_f2}
	\end{align}	
Using $h={\hat f}$,  the definition \eqref{e:Upupsilon} gives
\begin{equation}
\Upupsilon_t  =          - [D^f {\hat f} ] ( \Uptheta_t, \Phi_t)     =  - \frac{1}{\alpha}  \{ \ddt {\hat f}_t +  [ f(\Uptheta_t, \qsaprobe_t )  -  {\widebar{f}}(\Uptheta_t) ]\}    
\label{e:UpupsilonDer}
\end{equation}

	
We now proceed through the three steps, with this representation in view:
	\begin{equation}
		\ddt\Uptheta_t = \alpha \bigl[  {\widebar{f}} (\Uptheta_t )     +  \widetilde Xi_t  ] \,,
		\label{e:QSA_ODE_fixed}
	\end{equation} 
	where $  \widetilde Xi_t   =  f(\Uptheta_t,\qsaprobe_t) -  {\widebar{f}}(\Uptheta_t) $ is called the  \textit{apparent noise}.  Understanding \eqref{e:BigGlobalODE} is equivalent to determining the functions $\{ {\cal W}^i \}$ in the representation
	\begin{equation}
		\begin{aligned} 
			\widetilde Xi_t  &  =
			-  \alpha \widebar{\Upupsilon}_t 
				+    \alpha^2 {\cal W}_t^0 + \alpha \ddt {\cal W}_t^1  +   \tfrac{d^2}{dt^2} {\cal W}_t^2  
		\end{aligned} 
		\label{e:EffNoiseDecomp}
	\end{equation}
	
	\wham{Step 1:}     Apply \eqref{e:chainRule_f2}  with $g={\hat f}$:
	\[
	\ddt  {\hat f}(\Uptheta_t,  \Phi_t)   =   \partial_\theta {\hat f}(\Uptheta_t,  \Phi_t)  \ddt \Uptheta_t    - 
	[ f(\Uptheta_t,  \qsaprobe_t)  -{\widebar{f}}(\Uptheta_t)] 
	\]
	This gives the first transformation of the apparent noise:
	\[
	\widetilde Xi_t  =    \underbrace{-\ddt  {\hat f}(\Uptheta_t,  \Phi_t)  }_{\text{High pass}}
	+    
	\underbrace{\alpha   \partial_\theta {\hat f}(\Uptheta_t,  \Phi_t) f(\Uptheta_t,  \qsaprobe_t)   \phantom{\ddt }}_{\text{Attenuation}}
	\]
Recalling \eqref{e:Upupsilon}  gives  in shorthand notation,
\begin{equation}
		\widetilde Xi_t  =   -\ddt  {\hat f}_t   -\alpha \Upupsilon_t
		\label{e:EndOfStep1}
\end{equation}
 

	\wham{Step 2:}    Repeat previous argument with ${\hat f}$:
	\[
	\begin{aligned}
		\ddt {\hat h}af  (\Uptheta_t,\Phi_t)    &=    \alpha [D^f {\hat h}af ] (\Uptheta_t, \Phi_t )      -  {\hat f}_t     
		\\
		&\Rightarrow  {\hat f}_t   =   \alpha [D^f {\hat h}af ] (\Uptheta_t, \Phi_t )   
		-  \ddt {\hat h}af_t  				
		\\					
		&\Rightarrow 
		\ddt {\hat f}_t   =   \alpha \ddt  [D^f {\hat h}af ] (\Uptheta_t, \Phi_t ) 
		-  \tfrac{d^2}{dt^2} {\hat h}af_t  
	\end{aligned} 
	\] 
	
\wham{Step 3:}    Repeat   with $\Upupsilon$:
	\[
	\begin{aligned}
		\ddt {\widehat U}pupsilon (\Uptheta_t,\Phi_t)    &=    \alpha [D^f {\widehat U}pupsilon ] (\Uptheta_t, \Phi_t )      - [  \Upupsilon_t     -   \widebar{\Upupsilon}(\Uptheta_t)]
		\\
		&\Rightarrow  \Upupsilon_t   = \widebar{\Upupsilon}_t +  \alpha [D^f {\widehat U}pupsilon ] (\Uptheta_t, \Phi_t ) 
		-\ddt {\widehat U}pupsilon_t
	\end{aligned} 
	\] 
	Steps 2 and 3 combined with  \eqref{e:EndOfStep1} lead to the P-mean flow  representation:
	
 

\begin{theorem}
\label[theorem]{t:PMF}
Suppose that continuously differentiable solutions to Poisson's equation exist, for each of the three
forcing functions $f$,   ${\hat f}$,  and $\Upupsilon$.   Then, 
\begin{romannum}
\begin{subequations}
\item   The pre-P mean flow representation holds:
\begin{align} 
		\ddt Y_t  & =   \alpha\bigl[   {\widebar{f}}( Y_t    )   
			- \alpha  \bigl(B_t  {\hat f}_t  +  \Upupsilon_t  \bigr)  \bigr]
		\,,  \qquad B_t    =   \int_0^1 {\bar{A}} ( Y_t   -  r \alpha {\hat f}_t)    \, dr  
 \label{e:prematureMeanflow}
\\
 \Uptheta_t & =	Y_t    - \alpha {\hat f}_t
\end{align}
\end{subequations}
\begin{subequations}
\item
The P mean flow  representation \eqref{e:BigGlobalODE}  holds with   
		\begin{align}
			{\cal W}_t^0  = {\cal W}^0  (\Uptheta_t,\Phi_t)    & \mathbin{:=}  -  [D^f {\widehat U}pupsilon ] (\Uptheta_t,\Phi_t)    
			\label{e:AllTheNoise0}   
			\\
			{\cal W}_t^1  = {\cal W}^1  (\Uptheta_t,\Phi_t)    & \mathbin{:=} - [D^f {\hat h}af ] (\Uptheta_t,\Phi_t)   +    {\widehat U}pupsilon  (\Uptheta_t,\Phi_t)       
			\label{e:AllTheNoise1}
			\\
			{\cal W}_t^2  = {\cal W}^2  (\Uptheta_t,\Phi_t)    & \mathbin{:=}   {\hat h}af(\Uptheta_t,\Phi_t)  
			\label{e:AllTheNoise2}
\end{align}  
\qed
 	\label{e:AllTheNoise}
\end{subequations}
\end{romannum}
\end{theorem}
 
 

\Cref{t:hah} provides conditions ensuring existence of smooth solutions to Poisson's equation.
Analysis in \Cref{s:BakerPoisson} also leads to the proof that $\widebar{\Upupsilon}(\theta) = {\sf Z}ero$ when the frequencies satisfy (A0).  
Below is a summary of conclusions:  

	\begin{subequations}		
	\begin{lemma}
		\label{t:BiasVariance}
		Suppose that $\theta^*\in\field{R}^d$ is the unique solution to ${\widebar{f}}(\theta)={\sf Z}ero$.  
		Suppose moreover that Assumption (A2) holds, and denote ${\tilde \theta} = \theta-\theta^\text{\rm s.t.\,}arp$.
		Then,
		\begin{romannum}
		\item  There is a function ${\cal E}_A\colon\field{R}^d\to\field{R}^d$ satisfying
		\begin{equation}
			{\widebar{f}}(\theta) = A^\text{\rm s.t.\,}arp {\tilde \theta}  +   {\cal E}_A(\theta)  \,,\qquad \theta\in\field{R}^d  \, .
			\label{e:MVTtarget}
		\end{equation}
		The error term is Lipschitz continuous, and admits the quadratic bound $ {\cal E}_A(\theta) \le L_A  \|{\tilde \theta}\|^2$.
		
		\item  If $A^\text{\rm s.t.\,}arp$ is invertible, 
		\begin{equation}
			\Uptheta_t  -\theta^\text{\rm s.t.\,}arp  =     [A^\text{\rm s.t.\,}arp ]^{-1}  \bigl[    {\widebar{f}}(\Uptheta_t)  -  {\cal E}_A(\Uptheta_t)  \bigr]  \, .
			\label{e:biasTargetBiasVariancePre}
		\end{equation} 
		And, provided the target bias and variance are finite,   
		\begin{equation}
			\beta_\Uptheta  \le  \|  [A^\text{\rm s.t.\,}arp ]^{-1} \|_F \bigl[  \beta_{\widebar{f}}two     + L_A \,   \sigma_\Uptheta^2 \bigr]  \, .    
			\label{e:biasTargetBiasVariance}
		\end{equation}
		where   the subscript $F$ indicates the Frobenius norm.
	\end{romannum}
\qed
	\end{lemma}
	
\end{subequations}	

  
\begin{proposition}
\label[proposition]{t:orthCorBigPre}   
Consider the QSA ODE \eqref{e:QSAgen_alpha} subject to  (A0a) and (A1), but with arbitrary choices of frequencies $\{\omega_i\}$.  
Then, subject to the assumptions of \Cref{t:PMF},
 \begin{romannum}
\item  
For any invariant measure $ \upvarpi$ exists,  with compact support,
\[
 b_{\widebar{f}}two = {\sf E}  \bigl[ {\widebar{f}}(\Uptheta)  \bigr]
		 =   \alpha  {\sf E}  \bigl[  \Upupsilon(\Uppsi)  \bigr]   - \alpha^2 {\sf E}  \bigl[  {\cal W}^0    (\Uppsi)  \bigr]  
\]
where $\Uppsi = (\Uptheta, \Upphi)  $ denotes a $\Uppi$-valued random vector with distribution $ \upvarpi$. This implies $b_{\widebar{f}}two = O(\alpha)$.

		 
		 

\item
If in addition (A0) holds, then     $\widebar{\Upupsilon}(\theta) = 0 $  for all $\theta$, and
\[
 b_{\widebar{f}}two = {\sf E}  \bigl[ {\widebar{f}}(\Uptheta)  \bigr]
		 =   - \alpha^2 {\sf E}  \bigl[  {\cal W}^0    (\Uppsi)  \bigr]  
\]
with   ${\cal W}^0   $ defined below  \eqref{e:AllTheNoise0}, implying $b_{\widebar{f}}two = O(\alpha^2)$ in this case.
\qed
\end{romannum}
\end{proposition}
 
 \begin{proof}
 Part (i) follows  from \eqref{e:EndOfStep1}.
For (ii), the function classes $S$ and ${\widehat S}$ are \textit{orthogonal}: 
for
$g = h\circ G\in S$ and $\ell\in{\widehat S}$ we must have
\begin{equation}
	\int h(G(z)) \ell(z)  \, \uppi(dz)   =0 \,.
	\label{e:S-haSorth}
\end{equation}
In view of \eqref{e:Upupsilon},  the $i$th entry of $\Upupsilon(\theta)$ may be expressed
\[
\Upupsilon_i= \sum_{j=1}^d  \widehat A _{i,j }     g_j     \, ,
\]
with $g_j(\theta,z)  = f_j(\theta, G(z))$, so that $g_j\in S$.
For each $\theta\in\field{R}^d$ we have $\widehat A _{i,j }(\theta,\,\cdot\,) \in{\widehat S}$, so the result follows from \eqref{e:S-haSorth}.

The result that the target bias is of order $O(\alpha^2)$ follows from (ii) in \Cref{t:orthCorBigPre}.
\end{proof}
 The conclusion $\widebar{\Upupsilon}(\theta) = 0 $ in part~(ii) is most surprising;  it is based on geometry illustrated in \Cref{fig:FishOrthogonality}.
   \begin{figure}
\caption{Orthogonality of two function classes: 1.~Analytic functions of the probing signal, and 2.~Corresponding solutions to Poisson's equation.}
\label{fig:FishOrthogonality}
\end{figure}  

\subsection{Implications to filter design}
\label{s:filterDesign}

Filter design proceeds under the assumption that the state process is ultimately bounded, in the sense of 	\eqref{e:UB}, but we require something slightly stronger for the family of QSA ODEs with fixed gain: 

\wham{\boldmath{$\alpha^0$}-Ultimate Boundedness:}   With  $\alpha^0>0$ constant, the family of QSA ODE models is  
\textit{$\alpha^0$-ultimately bounded} if there is a fixed constant $B$ such that  for each $\alpha \in (0, \alpha^0]$,
 initial condition $\Uptheta_0 =\theta$,  and $\Phi_0 = z$,   there is a finite time $t_0 = t_0(\theta,z,\alpha)$ such that \eqref{e:UB}  holds for \eqref{e:QSAgen_alpha}.
 It is also assumed that $t_0$ is continuous on its domain.
 
 

 

 
Criteria for $\alpha^0$-ultimate boundedness are surveyed  in \Cref{s:someStability}---see  \Cref{t:ODEstate_is_bdd}.

Assumption~(A4) will be put in place throughout the remainder of the paper:

 
\wham{(A4)} 
The family of QSA ODE models is  
$\alpha^0$-ultimately bounded,  and the mean flow satisfies the following two conditions:
\begin{romannum}
\item   	The ODE $\ddt \upvartheta_t  = {\widebar{f}}  ( \upvartheta_t  )$ is globally asymptotically stable with unique equilibrium $\theta^\text{\rm s.t.\,}arp$.

\item  The matrix   $A^\text{\rm s.t.\,}arp = {\bar{A}} (\theta^\text{\rm s.t.\,}arp)$ is Hurwitz. 
\end{romannum}

 
 	\begin{subequations}

Subject to $\alpha^0$-ultimate boundedness we obtain the following consequence of the   P-mean flow representation: there  are finite constants $B_\ell \le B_\ell^+$, independent of $\alpha$, such that  for each $ \theta =\Uptheta_0\in\field{R}^d$  and  $\alpha \in (0, \alpha^0]$, the following first entrance time is finite:
\begin{equation}
T_\ell(\theta) 
=  \min \{  t\ge 0 :   \max_{z_0} \| \Uptheta_t -\theta^* \|  \le  B_\ell  \}  
\label{e:FirstEntranceLocal}
\end{equation}
where the maximum is over $\Phi_0 = z_0$.  
Moreover,  for $t\ge T_\ell$ the trajectory is constrained to the larger region:   For each $\Psi_0 = (\theta,\Phi_0)$,  
\begin{equation}
 \| \Uptheta_t -\theta^* \|  \le  B_\ell^+  \,,\qquad t\ge T_\ell(\theta) 
\label{e:FinalEntranceLocal}
\end{equation}
The time $ T_\ell(\theta) $   depends on $\alpha$,  but the constants $B_\ell $,  $ B_\ell^+$ do not. 

\end{subequations}

\Cref{t:ODEstate_is_bdd} contains a criterion for    $\alpha^0$-ultimate boundedness in a strong form. 
From \eqref{e:littleV4bddODEstate} it follows that we can construct $\alpha^0$, together with a   function  ${\bar{T}}^0_\ell$,  
satisfying  
\[
T^\alpha_\ell(\theta)  \le {\bar{T}}^\alpha_\ell(\theta)   \mathbin{:=}   \frac{\alpha^0}{\alpha} {\bar{T}}^0_\ell  (\theta)
\,,
\qquad
\textit{for all $\theta$ and all $0<\alpha\le \alpha^0$.}
\]

The matrix $A^*  =  {\bar{A}}(\theta^*)$ is invertible under (A4).
In this case we define a new function ${\bar{Y}}\colon\field{R}^d\to\field{R}^d$ via
\begin{equation}
 {\bar{Y}}  =[A^*]^{-1} \widebar{\Upupsilon}
\label{e:barY}
\end{equation}
along with the fixed vectors   $\widebar{\Upupsilon}^* = \widebar{\Upupsilon}(\theta^*)$  and ${\bar{Y}}^* = [A^*]^{-1} \widebar{\Upupsilon}^*$. 
This notation is used for approximations obtained in  \cite{laumey22,CSRL} for vanishing gain algorithms.

 
\begin{proposition}
\label[proposition]{t:PMFlin}
If  Assumptions~(A1)--(A4) hold, then  for each initial condition $\Psi_0 = (\theta, z_0)\in\Uppi$,
\begin{equation}
			\ddt \Uptheta_t    =  \alpha A^*  \bigl[ \Uptheta_t -\theta^\text{\rm s.t.\,}arp    -     \alpha   {\bar{Y}}^\text{\rm s.t.\,}arp 
			 +    {\cal W}_t   \bigr] +O(\alpha^3)     \,,\qquad \textit{for all   $ t\ge T_\ell(\theta) $}
\label{e:BigGlobalODElin}
\end{equation}
\end{proposition}

This linearization is useful only if  the matrix   $A^\text{\rm s.t.\,}arp = {\bar{A}} (\theta^\text{\rm s.t.\,}arp)$ is Hurwitz, which is why this assumption appears in (A4).     Subject to this assumption,  the approximation \eqref{e:BigGlobalODElin} strongly suggests that we apply linear filtering techniques to reduce volatility.  

 

For reasons that will become clear in the proof of \Cref{t:Couple_theta} that follows,   a filter that obtains AAD of order $O(\alpha^2) $  can be chosen to be of second order, with bandwidth approximately equal to $\alpha$.   
It must also have unity DC gain, which brings us to the following special form for the filtered estimates:
\begin{equation}
 \tfrac{d^2}{dt^2}  \UpthetaF_t 
+
2\gamma\zeta\tfrac{d}{dt}  \UpthetaF_t
+\gamma^2\UpthetaF_t  = \gamma^2 \Uptheta_t
\label{e:GVSOS}
\end{equation}
It is the unity gain that prevents us from eliminating the error inherited from $\widebar{\Upupsilon}$,  illustrated using ${\bar{Y}}^\text{\rm s.t.\,}arp$ in
\Cref{fig:QSArep}.   

If ${\bar{Y}}^*=0$, then the coefficients $\gamma$ and $\zeta$ can be designed to reduce AAD dramatically.
Theory predicts that a first order filter is not effective in general.   

In the following companion to \Cref{t:BigQSA} we consider the family of QSA ODEs with fixed gain  $\alpha$,
and let the natural frequency $\gamma$ scale with $\alpha$.

 
\begin{subequations}
\begin{theorem}
\label[theorem]{t:Couple_theta}
Consider the solution to the QSA ODE with fixed gain.  
For each $\alpha>0$,   the filtered estimates are obtained using \eqref{e:GVSOS},  in which
the value of $\zeta \in (0,1)$ is fixed,  and $\gamma =\eta \alpha$ with $\eta>0$ also fixed.    

 

Subject to  Assumptions~(A1)--(A4), the following conclusions hold for $0<\alpha\le \alpha^0$:
\begin{align}
	\|\Uptheta_t - \theta^*\|  & \leq O(\alpha)  + o(1)  
\label{e:fixedGain1}
\\	
\|\UpthetaF_t - \theta^*\|  & \leq O(\alpha \gamma + \alpha^2)  + O(\alpha \|\widebar{\Upupsilon}^*\|) + o(1)   
\label{e:fixedGain1PR}
\end{align}
with $\widebar{\Upupsilon}^* =\widebar{\Upupsilon}(\theta^*)$. 
Consequently,   the AAD for   $\UpthetaF$  is  bounded by $O(\alpha^2)$ provided    $\widebar{\Upupsilon}^* = 0$.
\end{theorem}
 \label{e:BigFilter}
\end{subequations} 

 The full proof is postponed to \Cref{s:AppStability}.  
The proof of \eqref{e:fixedGain1} follows from two conclusions from the given assumptions:  that the ODE \eqref{e:barf_ODE_gen} is exponentially asymptotically stable,  by \Cref{t:qsaprobeExpAS},   which implies the existence of a smooth Lyapunov function for the mean flow.
One solution is
\begin{equation}
	V(\theta) = \int_0^{T} e^{\delta_0 t} \| \upvartheta_t -\theta^* \|^2\,dt  \,,\qquad   \upvartheta_0 =\theta\in\field{R}^d
	\label{e:notBMV}
\end{equation}
with $\delta_0>0$ chosen so that the integrand is vanishing,  and $T>0$ sufficiently large so that  $e^{\delta_0 T} \| \upvartheta_T -\theta^* \|  \le \half  \| \theta -\theta^* \| $ for any $\theta$.  
This Lyapunov function is then applied to  the representation  \eqref{e:prematureMeanflow}.       The proof of \eqref{e:fixedGain1PR} is then simpler since we can apply \eqref{e:fixedGain1} to justify linearization around $\theta^*$,  applying an entirely linear systems theory analysis.  

The approximations \eqref{e:BigFilter}  imply bounds on the absolute deviation of parameter estimates,  and hence the AAD.
Bounds on bias and variance also  follow as corollaries to \Cref{t:Couple_theta}.  

\begin{corollary}
	\label{t:asympt_bias} 
	Under the assumptions of \Cref{t:Couple_theta},
	
	\begin{subequations}
		\begin{romannum}
		\item
		The asymptotic bias and variance \eqref{e:biasdef} admit the bounds,  
		\begin{align}

			\beta_{\Uptheta}  &= O(\alpha)\,,
			\quad 
			\sigma^2_{\Uptheta} = O(\alpha^2)  \,,
			\label{e:part1corr}
			\\
			\beta_{\UpthetaF} &= O(\alpha) \, ,
			\quad \sigma^2_{\UpthetaF} = O(\alpha^2)  \, ,
			\label{e:part1corrF}

		\end{align}	
		\item
		If in addition (A0) holds, then  
		\begin{align}
			\beta_{\Uptheta} &= O(\alpha^2)  \,,
			\quad  \sigma^2_{\Uptheta} = O(\alpha^2)  \, , 
			\label{e:part2corr}
			\\
			\beta_{\UpthetaF} &= O(\alpha^2)\,,
			\quad  \sigma^2_{\UpthetaF} = O(\alpha^4) \, .
			\label{e:part2corrF}
		\end{align}
	\end{romannum}
	\end{subequations}
\qed
\end{corollary}

Assumption~(A0) has the largest impact on bias and variance.    Equation~\eqref{e:part2corr}
tells us that the variance is of order $O(\alpha^2)$ subject to this restriction on frequencies,
which is remarkable when compared with standard results from SA theory 
(see  \eqref{e:BM} and discussion that follows).    
Filtering brings the variance down to   $O(\alpha^4)$;   a restatement of the second bound in \eqref{e:part2corrF}.

  
  

  \begin{subequations}
	
\paragraph{Example }
Two linear QSA ODEs will be used to illustrate the impact of $\widebar{\Upupsilon}^*$ on AAD:  
	\begin{align}
		\ddt \Uptheta_t &= \alpha[A_t\Uptheta_t + 2\sin(\omega t)+1]  
		\label{e:lin_barY0}
		\\
		\ddt \Uptheta_t &= \alpha[   A_t\Uptheta_t + 2\cos(\omega t)+1]  
		\label{e:lin_barYnot}
\end{align}  
where  $A_t = -(1+\sin(\omega t))$ with $\omega>0$. 
Assumption (A0b) holds for \eqref{e:lin_barY0} with $K=1$  in \eqref{e:qsaprobe0}, using
 $\qsaprobe^0_t = \cos(\omega t)$. 
 In the second case \eqref{e:lin_barYnot}, Assumption (A0b) is \textit{violated}:  any realization of $\qsaprobe^0_t $
  will require $K=2$ with $\omega_1=\omega_2$.   

	
We   obtain a formula for  $\widebar{\Upupsilon}^* $ in each case, based on the following observations:  
	\begin{romannum}
		\item They share the common mean vector field ${\widebar{f}}(\theta) =-\theta+1$, so that $\theta^\text{\rm s.t.\,}arp =1$.
		
		\item  They also share the same linearizations:  $A(\theta,\qsaprobe) = \partial_\theta f(\theta,\qsaprobe) =  -(1+\sin(\omega t)) $  and $A^\text{\rm s.t.\,}arp=-1$.   
		
Consequently, ${\bar{Y}}(\theta) = [A^\text{\rm s.t.\,}arp]^{-1} \widebar{\Upupsilon}(\theta) = - \widebar{\Upupsilon}(\theta) $,  and
for any $\theta$ and $t$,
\[ 
 \widehat A( \theta, \Phi_t )   =  - \omega^{-1} \cos(\omega t)   
\]
Observe that  $ \int \widehat A( \theta, z )  \,\uppi(dz) = {\sf Z}ero$ 
and $\ddt  \widehat A( \theta, \Phi_t ) = -[  A (\theta,\qsaprobe_t) - {\bar{A}}(\theta)]$, as required.
 
		
		

		   
				
		
		\item  A crucial difference is the steady-state apparent noise:  $f(\theta^\text{\rm s.t.\,}arp, \qsaprobe_t) =  \sin(\omega t)$ for  \eqref{e:lin_barY0},  and 
		$f(\theta^\text{\rm s.t.\,}arp, \qsaprobe_t) = - \sin(\omega t) + 2\cos( \omega t)$ for  \eqref{e:lin_barYnot}. 
	\end{romannum}
From the \eqref{e:haEverythingElse}  we find that $\widebar{\Upupsilon}(\theta) = \widebar{\Upupsilon}^\text{\rm s.t.\,}arp$ (independent of $\theta$):  
\[
	\begin{aligned}
-{\bar{Y}}(\theta)  = 
		\widebar{\Upupsilon}(\theta)  & \mathbin{:=} 
		- \lim_{T\to\infty}\frac{1}{T} \int_0^T  \widehat A( \theta, \Phi_t ) \cdot  f(\theta, \qsaprobe_t)  \, dt
		\\
		&
		= 
		\begin{cases}
			\displaystyle
			\ \ 

			  \   0    &   \text{Model    \eqref{e:lin_barY0}}
			\\
			\displaystyle
			\ \ 
			-\lim_{T\to\infty}\frac{1}{T} \int_0^T \{ -\omega^{-1}  \cos(\omega t) \} \cdot \{ 2\cos(\omega t) \} \, dt  =  \frac{1}{\omega}  \quad  &   \text{Model    \eqref{e:lin_barYnot}}
		\end{cases}
\end{aligned} 
\] 
In each case,  we have used the fact that  $ \widehat A( \theta, \Phi_t )$ and $ \widehat A( \theta, \Phi_t )    \cdot \sin(\omega t)$ each have time average equal to zero.
	
\label{e:TwoLinExamples}
\end{subequations}

\wham{Simulation Setup} 
We are interested in comparing AAD for $\bm{\Uptheta}$ with and without filtering.
The special case $\omega =  0.1$ was considered, 
many  gains were tested in the range $\alpha \in [10^{-3}, 0.9]$,
and in each experiment the solution to the QSA ODE was approximated
using an Euler scheme with sampling time of $0.1$~sec. 
 
 
 Data $\{\Uptheta_t\}$ for each choice of $\alpha$ was then used to obtain the   two filtered estimates:
\[
{\Uptheta}^{1 \text{\tiny\sf F}}_t = \int^t_0 h^1_{t - \tau} \Uptheta_\tau\,d\tau, \qquad {\Uptheta}^{2 \text{\tiny\sf F}}_t = \int^t_0 h^2_{t - \tau} \Uptheta_\tau\,d\tau
\]
where $h^1_t$ and $h^2_t$ have Laplace transforms:
\begin{equation}
	H_1(s) = \frac{\gamma}{s+\gamma}, \qquad  H_2(s) = \frac{\gamma^2}{s^2+2\zeta \gamma s + \gamma^2 }
	\,,
	\qquad  \textit{ with $\gamma = \alpha$, $\zeta =0.8$. }
	\label{e:filter_numerics}
\end{equation}
The steady-state $L_1$ norm of the error was approximated using the Monte-Carlo estimate,
\[
\| \Uptheta_{\alpha}\|_{L_1} \mathbin{:=}  \frac{1}{T - T_0} \int^T_{T_0}  |\Uptheta_\tau -\theta^* | \, d\tau
\]
taking  $T_0 = 0.8 T$, so that only the final $20\%$ of the run was included.  This was repeated to obtain $ \|{\Uptheta}^{1 \text{\tiny\sf F}}_\alpha \|_{L_1}$ and $\|{\Uptheta}^{2 \text{\tiny\sf F}}_\alpha \|_{L_1}$ for each value of $\alpha$ considered. 

\begin{figure}
\caption{AAD as a function of $\alpha$ for different filtering techniques (left).   Impact of ${\bar{Y}
\label{fig:alpha2bias}
\end{figure} 

\wham{Results}
The final estimates for each case are shown plotted against $\alpha$ on a logarithmic scale in \Cref{fig:alpha2bias}~(a).
Also shown for comparison with expected bounds are plots of the linear functions $r_1(\alpha) = k_1 \alpha$, $r_2(\alpha) = k_2 \alpha^2$. The constants $k_1, k_2$ were chosen to aid comparison. The results agree with \Cref{t:Couple_theta}. In particular, we see in \Cref{fig:alpha2bias} that the term $\alpha{\bar{Y}}^*$ dominates the estimation error when ${\bar{Y}}^* \neq 0$. In this case, filtering has no improvement on reducing AAD of order $O(\alpha)$. When ${\bar{Y}}^*=0$, \Cref{fig:alpha2bias}~(a) shows that  filtering reduces AAD to $O(\alpha^2)$ for  $\alpha<0.1$,   and \textit{this is observed for both first and second order low pass filters}.  An explanation for the similar performance is contained in \Cref{s:FiltFoot}.

The evolution of the process $\{{\Uptheta}^{2 \text{\tiny\sf F}}_t\}$ for the specific case of $\alpha =0.01$ is shown in \Cref{fig:alpha2bias}~(b).    The volatility of $\bm{\Uptheta}$ is massive in each case, and reduced dramatically with filtering.

These results illustrate the importance of carefully designing the probing signal: the inclusion of a phase shift in one of the entries of ${\bf q}saprobe$ might appear harmless. In fact,  this small change results in significant estimation error---10\%\ in this example.  

\subsection{Some stability theory}
\label{s:someStability}

Here we establish ultimate boundedness of the QSA ODE under either of the two stability criteria surveyed in the introduction.   The Lyapunov bound is recalled here in slightly modified form, and given a new name:
\begin{romannum}  
\item[\textbf{(V4)}] 
There exists a $C^1$ function $V:\field{R}^d\rightarrow[1,\infty)$ and  a constant $\delta_0>0$ such that  the following bound holds for \eqref{e:barf_ODE_gen} for any time $\tau\ge 0$ for which $ \|\upvartheta_\tau\| \ge   \delta_0^{-1} $:
	\begin{equation}

\frac{d}{d\tau} V(\upvartheta_{\tau})  \leq  -  \delta_0  V(\upvartheta_{\tau}) 
	\label{e:ddt_bound_Lyapfun}
\end{equation}
Moreover, the  function is Lipschitz continuous with linear growth:  
there exists a constant $L_V  <\infty$ such that  
\begin{equation}
\begin{aligned}
|  V(\theta') - V(\theta)|  &  \le L_V  \|\theta' - \theta\|\,, &&   \textit{ for all $\theta$,  $\theta'$. }
   \\
 V(\theta)  &\ge  \|\theta \|\,,  && \textit{when $\|\theta \|  \ge \delta_0^{-1} $}
 \end{aligned} 
\label{e:VLip}
\end{equation}	  
\end{romannum}
The use of ``(V4)'' is because of its close connection with a similar ``drift condition'' in the theory of Markov chains and processes \cite{MT}.  

It is convenient to consider a relaxation:
\begin{romannum}
\item[\textbf{(V4')}] 
There is a function $V\colon\field{R}^d\to\field{R}_+$ satisfying the bound \eqref{e:VLip},

 together with constants $T<\infty$ and $\delta_1>0$ such that for each initial condition $\upvartheta_0$ and  $\tau\ge 0$,
 \begin{equation}
	V(\upvartheta_{\tau+T}) - V(\upvartheta_\tau) 
	\leq   - \delta_1 \|\upvartheta_\tau \| , \quad \text{for all } \tau \geq  0\,, \ \|\upvartheta_\tau \|> \delta_1^{-1}
\label{e:V4T}
\end{equation}
\end{romannum}
This unifies the theory since it is obviously implied by (V4),  and we will see that it is  implied by asymptotic
stability of the ODE@$\infty$.

\begin{theorem}
\label[theorem]{t:ODEstate_is_bdd}
Suppose that  (A1) and the Lyapunov bound   (V4') hold for the fixed gain QSA ODE \eqref{e:QSAgen_alpha}.
Then,   there is $\alpha^0>0$ and   positive constants $b$ and $\delta$ such that the following bounds hold for any $\alpha\in (0,\alpha^0]$ and any initial condition $\Uptheta_0,\Phi_0$:
\begin{equation}
\| \Uptheta_t \|   \le  b \| \Uptheta_0  \| \exp( -\alpha \delta t ) \,,
	\quad
	 \textit{for $ t \le T_1  $,  where $T_1 = \min\{ t : \| \Uptheta_t \|    \le \delta^{-1} \}$.}
\label{e:littleV4bddODEstate}
\end{equation} 
Consequently,  the family of QSA ODE models is  $\alpha^0$-ultimately bounded.
\end{theorem}

The full proof of the theorem is postponed to the Appendix, based on the following two steps:
\begin{romannum}
\item     Bounds on the difference between solutions to
QSA  ODE \eqref{e:QSAgen_alpha}
and
 the solutions to the ODE, 
\begin{equation}
	\ddt \upvartheta_t  =   \alpha {\widebar{f}}  ( \upvartheta_t  )
	\label{e:barf_ODE_gen_alpha}
\end{equation}
  
\item  A proof that (V4') implies  a version of \eqref{e:littleV4bddODEstate}  for the mean flow.
\end{romannum}
A similar result is established in \cite{laumey22} for the QSA ODE with vanishing gain.

  
  

  
The implications of (V4') for the mean flow are summarized here:  

\begin{proposition}
\label[proposition]{t:qsaprobeExpAS}
Suppose that (V4') holds, along with (A1) and (A4).   Then the ODE \eqref{e:barf_ODE_gen} is exponentially asymptotically stable:  for positive constants $b$ and $\delta$,  and any initial  condition $\upvartheta_0$,
\[
\| \upvartheta_t -  \theta^\text{\rm s.t.\,}arp \|   \le  b \| \upvartheta_0-  \theta^\text{\rm s.t.\,}arp \|    e^{-\delta t}  \,,\qquad t\ge 0
\]
\end{proposition}

\begin{proof}  
The proof is obtained by interpreting 	\Cref{fig:qsaprobeExpAS}.   The time $T_1$ is  the first $t\ge 0$ such that $\upvartheta_t\in B_1$,  and $T_2$   the first $t\ge T_1$ such that $\upvartheta_t\in B_2$.   The set $B_1$ is defined in the figure, and $B_2$ is a region of exponential asymptotic stability. 

\begin{figure}
\caption{Roadmap for the proof of exponential asymptotic stability for $	\ddt \upvartheta_t  = {\widebar{f}
\label{fig:qsaprobeExpAS}
\end{figure} 

\begin{subequations}
\Cref{t:littleV4bdd} tells us that,
under (V4'),  the constant $K$ defining
$B_1$ can be chosen sufficiently large so that there  are positive constants
 $b_1$ and $\delta_1$ such that the following bound holds for $\upvartheta_0\in B_1^c$:    
\begin{equation}
\| \upvartheta_{T_1} -  \theta^\text{\rm s.t.\,}arp \|   \le  b_1 \| \upvartheta_0-  \theta^\text{\rm s.t.\,}arp \|    e^{-\delta_1 T_1}
\label{e:EAS1}
\end{equation}
The set $B_1$ is bounded,  so that $T_2 - T_1$ is uniformly bounded over all initial conditions for $\upvartheta_0$---this is where we apply globally asymptotically stability from (A4).   Consequently, the constant $b_1$ can be increased if necessary to obtain
\begin{equation}
\| \upvartheta_{T_2} -  \theta^\text{\rm s.t.\,}arp \|   \le  b_1 \|  \upvartheta_{T_1} -  \theta^\text{\rm s.t.\,}arp \|     
\label{e:EAS2}
\end{equation}
   Finally we specify $B_2$:  it is chosen so that it is absorbing, that is $\upvartheta_t\in B_2$ for all $t\ge T_2$,   and for some positive constants
 $b_2$ and $\delta_2$,
\begin{equation}
\| \upvartheta_t -  \theta^\text{\rm s.t.\,}arp \|   \le  b_2 \| \upvartheta_{T_2} -  \theta^\text{\rm s.t.\,}arp \|    e^{-\delta_2 (t- T_2) } \,,\qquad t\ge T_2.
\label{e:EAS3}
\end{equation}
The proof is completed on combining the three bounds \eqref{e:EASproof}, with $\delta = \min(\delta_1,\delta_2)$.
\label{e:EASproof}
\end{subequations}
\end{proof}

 
 

\subsubsection{QSA solidarity}

	
Here we establish bounds between the solutions of the QSA ODE and the mean flow \eqref{e:barf_ODE_gen_alpha}.    However, it  is convenient to eliminate the $\alpha$ in  \eqref{e:barf_ODE_gen_alpha} through a change of time-scale:  
If $\bm{\upvartheta}$ is a solution to  \eqref{e:barf_ODE_gen_alpha},   then $\{ \upvartheta_{t/\alpha} : t\ge 0\}$ is a solution to \eqref{e:barf_ODE_gen}.    

Similarly, on denoting $\Uptheta_t^\alpha = \Uptheta_{t/\alpha} $,    the scaling is removed from \eqref{e:QSAgen_alpha}:
\begin{equation}
	\ddt\Uptheta^\alpha_t =   f(\Uptheta^\alpha_t, \qsaprobe^\alpha_t   )
\label{e:QSAgen_alpha_time-scaled}
\end{equation}
with $\qsaprobe^\alpha_t \mathbin{:=} \qsaprobe_{t/\alpha} $.
Hence time is speeded up for the probing signal when $\alpha>0 $ is small.   
 
\begin{proposition}
\label[proposition]{t:Couple_theta_local}
Consider the QSA ODE  subject to (A1).
There is $\alpha^0>0$, $B_0<\infty $,  and $\varepsilon_0 \colon\field{R}_+\to\field{R}_+$
 such that the following uniform bounds hold for respective solutions to \eqref{e:QSAgen_alpha_time-scaled} and \eqref{e:barf_ODE_gen}, with common initial condition  $\Uptheta_0 =\upvartheta_0=\theta_0\in\field{R}^d$.
 For each
$\Phi_0\in \Upomega$, 
$T>0$,
 and $\alpha\in (0,\alpha^0]$: 
\begin{equation}
\|\Uptheta_T^\alpha - \upvartheta_T\| \leq      B_0 \exp(2 L_f T)    \bigl[  1     +      \|\theta_0\|  \bigr]  \varepsilon_0( \alpha)
\label{e:Couple_theta_local}
\end{equation}
with  $  \varepsilon_0( \alpha  ) \to 0 $ as $\alpha {\hbox{\tiny down}}arrow 0$.  
\end{proposition}
 
 
 
	
	\begin{figure}
\caption{Solidarity between the time-scaled QSA ODE and the mean flow. }
\label{fig:SolidarityFixedGain}
\end{figure}
	
 
The proof of the proposition is contained in the Appendix. 
The main idea is illustrated in \Cref{fig:SolidarityFixedGain}:  the time axis is partitioned into intervals of width $\Delta>0$ on which $\Uptheta_t^\alpha$ is approximately constant.     With 
		$\tau_n = n\Delta$,   we let   $\{ \upvartheta^n_t :    t \ge   \tau_{n} \}$ denote the solution to \eqref{e:barf_ODE_gen}  with $\upvartheta^n_{\tau_n} = \Uptheta^\alpha_{\tau_n}$.    The proof follows in three steps:  
		1.\  bounds on the error $\|   \upvartheta^n_t - \Uptheta^\alpha_t \|$ for $t\in [\tau_n, \tau_{n+1}]$;
		2.\     bounds on  $\|   \upvartheta^n_t - \upvartheta_t \|$ for all $t$ and $n$,  and
		3.\  combine 1 and 2 (with minor additional approximations) to complete the proof.

 
\subsubsection{ODE@\texorpdfstring{$\infty$}{infinity}}

We now explain how to establish (V4') based on the ODE@$\infty$ introduced in \Cref{s:QSAstable}, with vector field $
	{\widebar{f}}_\infty $ defined in 	\eqref{e:barfinfty_QSA}.  
The limit exists in many applications, and is often much simpler than ${\widebar{f}}$.  In particular, ${\widebar{f}}_\infty(0) =0$ (the origin is a stationary point),  and the function is radially linear:  ${\widebar{f}}_\infty (c\theta)  = c{\widebar{f}}_\infty (\theta) $ for any $\theta\in\field{R}^d$ and $c\ge 0$.  Based on these properties it is known that the origin is locally stable in the sense of Lyapunov if and only if it is globally exponentially asymptotically stable \cite{bormey00a}.

The following result follows from arguments leading to the proof of   \cite[Prop.~4.22]{CSRL}:
\begin{theorem}
\label[theorem]{t:V4TfromBM}
	Suppose that the limit \eqref{e:barfinfty_QSA} exists, and suppose that the origin is  asymptotically stable for the  ODE@$\infty$.   Then, there is a Lipschitz continuous solution to \eqref{e:V4T}.  \qed
\end{theorem}

The proof of the theorem is contained in the Appendix, where it is shown that the following function is a solution to \eqref{e:V4T},
 similar to \eqref{e:notBMV}:
\begin{equation}
	V(\theta) = \int_0^{T}  \| \upvartheta_t^\infty \|\,dt  \,,\qquad   \upvartheta_0^\infty =\theta\in\field{R}^d
\label{e:BMV}
\end{equation}
  Lipschitz continuity of both ${\widebar{f}}_\infty$ and $V$ is inherited from ${\widebar{f}}$.
 

To illustrate application of the theorem, consider the gradient flow for which ${\widebar{f}} =   -  \nabla \Upgamma$:
\begin{equation}
\ddt \upvartheta =     -  \nabla \Upgamma\, (\upvartheta)
\label{e:GradFlow}
\end{equation}
with $\Upgamma\colon\field{R}^d\to\field{R}$,   with minimum denoted $\Upgamma^\text{\tiny\sf opt}$.
 The existence of ${\widebar{f}}_\infty$ implies that the following limit also holds: 
 \begin{equation}
\Upgamma_\infty (\theta) \mathbin{:=} \lim_{r \to \infty }  r^{-2}  \Upgamma(r\theta)    \,, \qquad \theta\in\field{R}^d 
\label{e:ObjInfinity}
\end{equation} 
This follows from the representation
\[
\Upgamma(\theta) -\Upgamma( 0) = \int_0^1 \ddt \Upgamma(t\theta) \, dt  = \int_0^1 \theta^{\rm tr\, }anspose\nabla \Upgamma\, (t\theta) \, dt   
\]
which implies the representation ${\widebar{f}}_\infty(\theta) = - \nabla \Upgamma_\infty(\theta)$ and
\begin{equation}
\Upgamma_\infty (\theta)   =  - \int_0^1 \theta^{\rm tr\, }anspose  \nabla \Upgamma_\infty(t\theta) \, dt   \,, \qquad \theta\in\field{R}^d 
\label{e:ObjInfinity2}
\end{equation}

 

The gradient flow is ultimately bounded under mild assumptions:

\begin{proposition}
	\label[proposition]{t:gradFlowODE}
	Suppose that $\Upgamma_\infty$ is  $C^1$,   and $\Upgamma_\infty(\theta)>0$   for $\theta\neq 0$.     Then the ODE@$\infty$ for the gradient flow is exponentially asymptotically stable:  for some $B<\infty$, $\delta>0$, and any initial $\upvartheta^\infty_0\in\field{R}^d$,
	\begin{equation}
		\|\upvartheta^\infty_t \|^2 
		\le  B   \|\upvartheta^\infty_0 \|^2  e^{ -  \delta    t} \,,\qquad \textit{for $t\ge 0$.}
		\label{e:expASgradFlowODE}
\end{equation}
Consequently,  the gradient flow itself is ultimately bounded.
\end{proposition}

\begin{proof}
Both $\Upgamma_\infty$  and $\| \nabla \Upgamma_\infty \|^2$ are bounded above and below by quadratic functions---this follows from continuity and the radial homogeneity properties
\[
\Upgamma_\infty( r \theta)  = r^2\Upgamma_\infty(  \theta) \,,\qquad   \nabla \Upgamma_\infty\, ( r \theta)  = r\nabla\Upgamma_\infty\, (  \theta)\,,\qquad \theta\in\field{R}^d\,, \  r>0
\]
Hence the following constants exist as finite positive numbers,
\[
\kappa_+ = \max \frac{\Upgamma_\infty(\theta)} {\| \theta\|^2} 
 \,,\quad 
\kappa_- = \min \frac{\Upgamma_\infty(\theta)} {\| \theta\|^2}  \,, 
\qquad
\chi_+ = \max \frac{ \| \nabla\Upgamma_\infty\, (  \theta) \|^2} {\| \theta\|^2}  \,,\quad 
\chi_- = \min \frac{ \| \nabla\Upgamma_\infty\, (  \theta) \|^2} {\| \theta\|^2}  \,,  
\]
where the min and max are over non-zero $\theta\in\field{R}^d$.

We have $\ddt \Upgamma_\infty( \upvartheta^\infty) = -  \| \nabla \Upgamma_\infty\, ( \upvartheta^\infty) \|^2$,  and from the above definitions we conclude that  $\ddt \Upgamma_\infty( \upvartheta^\infty) = - \delta  \Upgamma_\infty( \upvartheta^\infty)$ with $\delta = \chi_- /\kappa_+$.
Hence $ \Upgamma_\infty( \upvartheta^\infty_t) \le e^{-\delta t}  \Upgamma_\infty( \upvartheta^\infty_0)$,  which completes the proof using 
$B = \kappa_+/{\kappa_-}$.  
\end{proof}
 
 
 
 

\subsection{Extensions}
\label{s:extensions}

Three extensions follow.  

\begin{romannum}
\item 
The first is to almost periodic linear systems to illustrate the connections, and because one conclusion is useful in establishing finer bounds for general QSA ODEs.   
\item
Then, we review the main conclusions from \cite{laumey22} for vanishing gain QSA.

\item
Finally, we explain how bias bounds extend to stochastic approximation.    
\end{romannum}

\subsubsection{Almost periodic linear systems}
\label{s:QPLS}

A \textit{quasi periodic linear system} is a special case of QSA with constant gain:
\begin{equation}
\ddt \Uptheta_t =  \alpha A_t \Uptheta_t\,,\qquad A_t\mathbin{:=} A(\qsaprobe_t) \,,
\label{e:linQSA}
\end{equation}
The matrix-valued  function of time $D_T = \int_0^T  A_t\,  dt$,  $T\ge 0$, emerges in the theory.   In the notation of \Cref{t:QPL} this becomes $D_T = T{\bar{A}}^\text{\rm s.t.\,}arp   + \widehat A_0  -   \widehat A_T$.    If $A_T$ and $D_T$ commute for each $T$, one obtains a simple expression for the state transition matrix \cite{luk82,subwasred21}. 

In this linear case,  a bound on the Lyapunov exponent is equivalent to a bound on the following for each non-zero initial condition $\theta$:  
\begin{equation}
\Lambda_{\Uptheta}(\theta) = \lim_{t\to\infty}\tfrac{1}{t} \log (\| \Uptheta_t \|)  \,,\qquad  \Uptheta_0 =\theta\in\field{R}^d
\label{e:LyapExpODEstate}
\end{equation}
 
\begin{proposition}[Application to Quasi Periodic Linear Systems]
\label[proposition]{t:QPL}
Consider the linear QSA ODE \eqref{e:linQSA}, subject to the following subset of Assumptions (A1)--(A4):   
\begin{romannum}
\item 
  $A^\text{\rm s.t.\,}arp \mathbin{:=} {\sf E}[A(\qsaprobe)]$ is Hurwitz.    

\item  ${\hat f}(\theta, z)  = \widehat A(z) \theta $, which solves
\[
\widehat A(  \Phi_{t_0} ) = 
			\int_{t_0} ^{t_1}   [ A( \qsaprobe_t)  - A(\theta) ] \, dt    + \widehat A(   \Phi_{t_1} )   \,,\qquad 0\le t_0\le t_1
\]
\end{romannum}
 For each $\alpha>0$, let $\Lambda_{\Uptheta}^\alpha$ denote the Lyapunov exponent associated with \eqref{e:linQSA}.
 Then  $\Lambda_{\Uptheta}^\alpha <0$ for sufficiently small $\alpha>0$, and
\begin{equation}
\lim_{\alpha{\hbox{\tiny down}}arrow 0}  \frac{1}{\alpha} \Lambda_{\Uptheta}^\alpha  =   \text{\rm Re}(\lambda_1) <0
\label{e:linQSAasymptotics}
\end{equation}
where $\lambda_1$ is an eigenvalue of  $A^\text{\rm s.t.\,}arp$ with maximal real part.
\end{proposition}

\begin{proof}
We establish only the   upper bound,
\begin{equation}
\mathop{\rm lim{\,}sup}_{\alpha{\hbox{\tiny down}}arrow 0}  \frac{1}{\alpha} \Lambda_{\Uptheta}^\alpha  \le      \text{Re}(\lambda_1) 
\label{e:linQSAasymptoticsUpper}
\end{equation}
leaving the technical details for the lower bound to the reader. 

The following calculation is identical to the first step in the proof of  \Cref{t:PMF}, 
using the shorthand $\widehat A_t = \widehat A(\Phi_t)$:
\[
\begin{aligned}
\ddt {\hat f}(\Uptheta_t,  \Phi_t)    \mathbin{:=}  \ddt \bigl( \widehat A_t \Uptheta_t  \bigr)
            & = \bigl(  \ddt \widehat A_t  \bigr)  \Uptheta_t  +  \widehat A_t  \bigl(  \ddt \Uptheta_t \bigr)  
     \\
     & =    - [ A_t  - A^\text{\rm s.t.\,}arp ] \Uptheta_t    +  \alpha \widehat A_t  A_t \Uptheta_t   
\end{aligned} 
\] 
Note that    $\Upupsilon_t \mathbin{:=} - \widehat A_t  A_t \Uptheta_t  $ in the notation of
 \Cref{t:PMF}.

  
The resulting identity $A_t \Uptheta_t  =   A^\text{\rm s.t.\,}arp +   \alpha \widehat A_t  A_t  -  \ddt \bigl( \widehat A_t \Uptheta_t  \bigr)$  is substituted into  \eqref{e:linQSA}  to obtain
\[
\ddt \Uptheta_t =  \alpha   \bigl[  A^\text{\rm s.t.\,}arp +   \alpha \widehat A_t  A_t  \bigr] \Uptheta_t   - \alpha \ddt \bigl( \widehat A_t \Uptheta_t  \bigr)
\]
This motivates the introduction of a new state variable $Y_t =  [I + \alpha  \widehat A_t] \Uptheta_t $ to obtain the state space model
\[
\ddt Y_t =  \alpha   \bigl[  A^\text{\rm s.t.\,}arp +   \alpha \widehat A_t  A_t  \bigr]     \bigl[  I + \alpha  \widehat A_t  \bigr]   ^{-1}  Y_t
\]
The inverse exists for all sufficiently small $\alpha>0$ since $\{  \widehat A_t \}$ is bounded,   from which we obtain 
\begin{equation}
\ddt Y_t =  \alpha   \bigl[  A^\text{\rm s.t.\,}arp +   \alpha {\cal E}^\alpha_t  \bigr]      Y_t
\label{e:linQSArep}
\end{equation}
in which  $\{{\cal E}^\alpha_t\}$ is uniformly bounded in $t$, and $\alpha$ in a neighborhood of the origin.  
For each $t$,
\[
\lim_{\alpha{\hbox{\tiny down}}arrow 0} \tfrac{1}{\alpha} 
{\cal E}^\alpha_t   = \widehat A_t  A_t   -    A^\text{\rm s.t.\,}arp\widehat A_t   
\] 

The remainder of the proof is standard.  First, for any positive definite matrix $P$, on taking $V(\theta) = \theta^{\rm tr\, }anspose P \theta $
we have 
\[
2\Lambda_{\Uptheta}(\theta) 
 = \lim_{t\to\infty}\tfrac{1}{t} \log (V(\Uptheta_t) )   
 = \lim_{t\to\infty}\tfrac{1}{t} \log (V(Y_t) )  
\]
This is simply a consequence of the fact that  all norms are equivalent on $\field{R}^d$,  along with uniform bounds on the norm of $[I + \alpha  \widehat A_t] $ and its inverse.

The matrix $\delta I +  A^\text{\rm s.t.\,}arp $ is Hurwitz for any $0\le \delta <  | \text{Re}(\lambda_1)|$, so that there is a solution $P>0$ to the Lyapunov equation  ${A^\text{\rm s.t.\,}arp}^{\rm tr\, }anspose P + P A^\text{\rm s.t.\,}arp   + 2 \delta P =  - I $.   Applying \eqref{e:linQSArep} gives
\[
\ddt V(Y_t)  =   -\alpha \bigl[  2\delta V(Y_t)  +  \| Y_t\|^2 \bigr]    + \alpha^2 Y_t^{\rm tr\, }anspose \bigl[  {{\cal E}^\alpha_t  }^{\rm tr\, }anspose P + P {\cal E}^\alpha_t  \bigr]  Y_t.
\]
For $\alpha>0$ sufficiently small we obtain the upper bound,
\[
\begin{aligned}
2\Lambda_{\Uptheta}(\theta)    =  \lim_{t\to\infty}\tfrac{1}{t} \log (V(Y_t) )    & \le -2\alpha \delta +  O(\alpha^2)
   \\
   \textit{giving}\quad       \mathop{\rm lim{\,}sup}_{\alpha{\hbox{\tiny down}}arrow 0}  \frac{1}{\alpha} \Lambda_{\Uptheta}^\alpha & \le  -   \delta
\end{aligned} 
\]
This establishes \eqref{e:linQSAasymptoticsUpper} since  $\delta <  | \text{Re}(\lambda_1)|$ is arbitrary.
\end{proof}

The representation \eqref{e:linQSArep}  suggests that we might go further, using repeated applications of Poisson's equation as in our treatment of general QSA ODEs.   On writing  
\[
\ddt Y_t =  \alpha   A_t^1     Y_t  \,,\qquad    A_t^1  =   A^\text{\rm s.t.\,}arp +   \alpha {\cal E}^\alpha_t  
\]
successive applications of this approach would lead to a family of almost periodic linear systems,
\[
\ddt Y_t^i =  \alpha   A_t^i     Y_t^i    \,,\qquad    i\ge 2
\]
The hope is that we would obtain $ A_t^i  \to A^*$ as $i\to \infty  $ for each $t$, leading to the
 time-homogeneous linear system, combined with time varying coordinate transformation $\mathcal{T}_t$:  
\[
\ddt  Y_t^\infty = A^* Y_t^\infty\,,\qquad    \Uptheta_t = \mathcal{T}_t Y_t^\infty  
\]

\subsubsection{Vanishing gain}

The main result  of \cite{laumey22},  concerns rates of convergence for the QSA ODE,  and averaging techniques to improve these rates.  The gain is taken of the form \eqref{e:arho} with  $\rho \in (\half, 1)$ and $t_e=1$.    
In this setting, the use of  Polyak-Ruppert averaging is sometimes effective.   At time $T$, the averaged parameter estimate is defined as 
\begin{equation}
\UpthetaPR_T \mathbin{:=} \frac{1}{T-T_0}  \int_{T_0}^T   \Uptheta_t\, dt
\label{e:PRQSA}
\end{equation}
in which the starting time scales linearly with $T$:  for fixed $\kappa>1$ we chose   $T_0 $ to solve  $1/(T-T_0) = \kappa/T$.

 \begin{subequations}
\begin{theorem}
\label[theorem]{t:BigQSA}
Suppose that  Assumptions~(A1)--(A4) hold,   and that the solutions to the QSA ODE \eqref{e:QSA_ODE_gen} have bounded sample paths.  Then the following conclusions hold:
\begin{align}
\Uptheta_t  & = \theta^\text{\rm s.t.\,}arp  -  a_t   {\hat f}^*_t  + O(a_t \|{\bar{Y}}^*\|)  + o(a_t)    &&  t\ge t_0
\label{e:BigCouplingQSAcor}
\\	
		 \UpthetaPR_T  &= \theta^*+   a_T  [c(\kappa,\rho)  +o(1)]   {\bar{Y}}^* +  O(T^{-2\rho}) \,,\qquad &&T_0\ge t_0
		\label{e:PR_2rhobound} 
\end{align}
where $c(\kappa,\rho) =c(\kappa,\rho) = \kappa[1-(1-1/\kappa)^{1-\rho}]/(1-\rho)$,
and  ${\bar{Y}}^* = [A^\text{\rm s.t.\,}arp]^{-1} \widebar{\Upupsilon}^* $.  
 
Consequently,    $\UpthetaPR_T$ converges to $\theta^\text{\rm s.t.\,}arp$ with rate bounded by $O(T^{-2\rho})$ if and only if ${\bar{Y}}^* = 0$.
\end{theorem}

 \end{subequations}

This is a remarkable conclusion, since we can choose $\rho$ to obtain a rate of convergence arbitrarily close to $O(T^{-2})$.   

The key step in the proof is to obtain a version of the perturbative mean flow, this time for the scaled error,
 \begin{equation}
Z_t = \frac{1}{a_t}  \bigl(  \Uptheta_t -  \theta^\text{\rm s.t.\,}arp \bigr) \, ,   \qquad t\ge t_0
\label{e:scaled_error}
\end{equation}
The analysis is entirely  local, so that a linearization around $\theta^\text{\rm s.t.\,}arp$ is justified:
\[
\ddt Z_t   =    a_t       [o(1) I +  A^\text{\rm s.t.\,}arp  ]     Z_t   +      \widetilde Xi_t   
\]
in which the definition of the apparent noise $\{  \widetilde Xi_t    \}$  is unchanged.    Applying 	\eqref{e:EndOfStep1} then gives,
\begin{equation} 
\ddt Z_t   =    a_t       [o(1) I +  A^\text{\rm s.t.\,}arp  ]     Z_t    -\ddt  {\hat f}_t   -a_t \Upupsilon_t
\label{e:scaledErrorRep_cor_rho}
\end{equation} 
The change of variables $Y_t = Z_t   +      {\hat f}_t $ then gives an approximation similar to  \eqref{e:prematureMeanflow}:
\[
\ddt Y_t   =    a_t    \bigl(     [o(1) I +  A^\text{\rm s.t.\,}arp  ]     Y_t     \bigr)
	  -a_t  \bigl(         [o(1) I +  A^\text{\rm s.t.\,}arp  ]      {\hat f}_t   +\Upupsilon_t \bigr)
\]
This leads to two consequences:   a rate of convergence of order $O(a_t)$ for $\Uptheta_t$,  and a bound of order $O(a_t^2)$  for $\UpthetaPR_t$.  However, \textit{a rate of order $O(a_t^2)$ is only possible if $\widebar{\Upupsilon}^*={\sf Z}ero$.}
 
 
\subsubsection{Implications to stochastic approximation}
\label{s:SAbias}

To obtain extensions to the stochastic setting recall  the definition of $\Upupsilon$ in \eqref{e:Upupsilon},  
and the representations for its steady-state mean value:
\begin{equation}
\begin{aligned}
		\Upupsilon(\theta,z)   &=   -  [\partial_\theta {\hat f}   (\theta,z ) ]\cdot   f (\theta, G(z))       
   \\
\alpha {\sf E}[ \Upupsilon(\Psi)  ]   &            
			=   - {\sf E} [ \widetilde Xi]  \,,\qquad\qquad \widetilde Xi = f(\Psi )  -  {\widebar{f}}(\Uptheta)   
\end{aligned} 
\label{e:UpupsilonDer2}
\end{equation}
The second identity is obtained on combining \eqref{e:EndOfStep1} with \Cref{t:QSAgenerator}  (with expectations in steady-state).

 

Something similar appears in a stochastic setting,   and presents a similar challenge for 
	variance
 reduction via averaging.  Here we establish an analog of \Cref{t:orthCorBigPre}~(i).

Consider the general SA recursion \eqref{e:SA_recur} with constant gain,
\[
 \theta_{n+1} = \theta_{n}+ \alpha  \bigl[  {\widebar{f}}(\theta_n)  +  \widetilde Xi_n \bigr]
		\,,\qquad    \widetilde Xi_n   =  f(\theta_n,\qsaprobe_{n+1})     - {\widebar{f}}(\theta_n)
\]
It is assumed throughout that $f$ satisfies the uniform Lipschitz condition imposed in (A1).

Assume a finite state space Markovian realization for the probing sequence:   there is an irreducible Markov chain on a finite state space ${\sf X}$ and function $G\colon{\sf X} \to \field{R}^{d\times d}$ that specifies $\qsaprobe_n = G(\Phi_n)$.   
Then $\Psi_n \mathbin{:=}(\theta_n,\Phi_n)$ is a Markov chain, and it satisfies the Feller property under our standing assumption that $f$ is continuous in $\theta$.  

Assume that the sequence is bounded in $L_2$ for at least one initial condition:  $\sup_n {\sf E}[\| \theta_n\|^2 ]<\infty$.   The Feller property combined with tightness implies the existence of  at least one invariant probability measure $\upvarpi$ on  ${\cal B}(\field{R}^d\times {\sf X})$.  Moreover,  under the invariant distribution,  the random vector $\theta_n$ has a finite $p$th moment for $p<2$.
Existence of $\upvarpi$  is justified in Section 12.1.1 of  \cite{MT}, and the uniform second moment implies uniform integrability  of $\{ \| \theta_n\|^p : n\ge 1\}$ from one initial condition, implying $ {\sf E}[\| \theta_n\|^p] <\infty$ 
when $(\theta_n, \Phi_n)\sim \upvarpi$.

Suppose that $\Psi_0 \sim \upvarpi$ so that $\{ \Psi_n : n\ge 0\}$ is a stationary process.  Taking expectations of both sides of the recursion gives:
\[
{\sf E}[ \theta_n ] ={\sf E}[ \theta_{n+1} ] 
			={\sf E}[  \theta_{n} ]
			+ \alpha {\sf E} \bigl[  {\widebar{f}}(\theta_n) \bigr] 
			+ \alpha {\sf E} \bigl[ \widetilde Xi_n \bigr] 
\]
From this we obtain the identity  
\[
{\sf E} \bigl[  {\widebar{f}}(\theta_n) \bigr] 
			= -{\sf E} \bigl[ \widetilde Xi_n \bigr] 
\]
When $\{ \widetilde Xi_n \}$ is a martingale difference sequence, such as when ${\bf P}hi$ is i.i.d., the right hand side is zero,  which is  encouraging.   In the Markov case this conclusion is no longer valid.

We can see this by bringing in Poisson's equation for the Markov chain:

\[
{\sf E}[ {\hat f}(\theta ,  \Phi_{n+1})   \mid {\cal F}_n]  =  {\hat f}(\theta ,  \Phi_{n})  -  
\bigl[ f(\theta,\qsaprobe_{n})     - {\widebar{f}}(\theta)\bigr]  
\]
where ${\cal F}_n = \sigma(\Psi_0\,,\dots\,, \Psi_n)$.   This implies a useful representation for the apparent noise, 
\begin{align}
\widetilde Xi_{n}   &=  f(\theta_n,\qsaprobe_{n+1})     - {\widebar{f}}(\theta_n)
=
  \bigl[ {\hat f}(\theta_n ,  \Phi_{n+1})   -  {\hat f}(\theta_n ,  \Phi_{n+2})  \bigr]  + {\cal W}_{n+2} 
 \label{e:SAdecomp}
\end{align}
where $ {\cal W}_{n+1}  =  {\hat f}(\theta ,  \Phi_{n+2})  - {\sf E}[ {\hat f}(\theta ,  \Phi_{n+2})   \mid {\cal F}_{n+1}]  $ is a martingale difference sequence.  

It is assumed that the solution to Poisson's equation is normalized to have zero mean:  ${\sf Z}ero = \sum_z {\hat f}(\theta , z) \, \uppi(z)$ for each $\theta$,  with $\uppi$ the unique invariant probability measure for ${\bf P}hi$.
Under the Lipschitz assumption on $f$ we have $\| {\hat f}(\theta , z) \|  \le L_{{\hat f}} [1 + \| \theta\| ]$ for some constant $ L_{{\hat f}} $ and all $(\theta, z)$   (apply any of the standard representations of Poisson's equation in standard texts, such as \cite{MT}).

\begin{proposition}
\label[proposition]{t:ApparentNoiseSA}
The  tracking bias may be expressed
\[
{\sf E}[{\widebar{f}}(\theta_n)]  =     \alpha {\sf E}[\Upupsilon_n ]  \,,
\qquad
\textit{with}
\quad
 \Upupsilon_n      \mathbin{:=}   -  \frac{1}{\alpha} \bigl[ {\hat f}(\theta_{n+1} ,  \Phi_{n+2})   -  {\hat f}(\theta_n ,  \Phi_{n+2}) \bigr] \,,
\]
where the expectations are taken in steady-state, so independent of $n$.
\end{proposition}
This is entirely consistent with the conclusion in the continuous time case,     in which the derivative appearing in \eqref{e:UpupsilonDer2} is replaced by a difference. 
 

 There are deeper connections explained in the Appendix---see in  particular the formula \eqref{e:newUpupsilon}.   As in the deterministic setting, we expect $\Upupsilon_n$ to be bounded over $0 < \alpha\le\alpha^0$ for suitably small $\alpha^0>0$,  but it may be large if the Markov chain has significant memory.
 

\begin{proof}
Applying \eqref{e:SAdecomp},  the   partial sums of the apparent noise may be expressed 
\[
\sum_{n=1}^{N} \widetilde Xi_{n}  = 
{\cal M}_{N+2}     
+
\bigl[ {\hat f}(\theta_0 ,  \Phi_2)  -  {\hat f}(\theta_{N} ,  \Phi_{N+2}) \bigr]
+
\sum_{n=0}^{N-1} \bigl[ {\hat f}(\theta_{n+1} ,  \Phi_{n+2})   -  {\hat f}(\theta_n ,  \Phi_{n+2}) \bigr]
\]
where $\{  {\cal M}_{N+2}  =   \sum_{n=1}^{N}{\cal W}_{n+2}    \}$ is a zero-mean martingale.
Consider now the stationary regime with  $\Psi_0 \sim \upvarpi$.  
Taking expectations of both sides, dividing each side by $N$, and letting $N\to\infty$ gives the desired result.  \end{proof}

 

\section{Implications to Extremum Seeking Control}
\label{s:ESC}

\subsection{What is extremum seeking control?}

	While much of ESC theory concerns tracking the minimizer of a time-varying objective  (that is,  $\Upgamma$ depends on both the parameter $\theta$ and time $t$),   it is simplest to first explain the ideas in the context of global optimization of the static objective $\Upgamma \colon\field{R}^d\to\field{R}$.

We begin with an explanation of the appearance of $- M  {\widetilde{\nabla}\!}_t\Upgamma $ in Fig~\ref{fig:ESC}.     The low pass filter with output  $\Uptheta_t$ is designed so that the derivative $\ddt \Uptheta_t$ is small in magnitude, to justify a quasi-static analysis.   
	An example is
	\begin{equation}
		\ddt \Uptheta_t =   - \sigma [ \Uptheta_t  -\theta^\text{\tiny{\sf ctr}} ]  +   \alpha   U_t  \,,   \qquad U_t = -  M {\widetilde{\nabla}\!}_t\Upgamma 
		\label{e:ESC_LP}
\end{equation}
with parameters satisfying $0< \sigma < \alpha$. 

Motivation for the high pass filter   is provided using two extremes.  In each case it is assumed that the probing gain is constant:   $\upepsilon_t \equiv \varepsilon$, independent of time.  
	
	\wham{Pure differentiation:}   In this case, the figure is interpreted as follows
	\begin{equation}
		M {\widetilde{\nabla}\!}_t\Upgamma  =  \bigl(\ddt \qsaprobe_t  \bigr) \bigl(   \tfrac{1}{\varepsilon}  \ddt  \Upgamma(\Uptheta_t + \varepsilon \qsaprobe_t)  {\sf B}igr)
		\label{e:ESC1}
	\end{equation}
	Adopting the notation from the figure $  \check{\qsaprobe}_t = \ddt \qsaprobe_t  $,  we obtain by the chain rule 
	\begin{equation}
		M {\widetilde{\nabla}\!}_t\Upgamma  =   \check{\qsaprobe}_t {\check{\qsaprobe}_t}^{\rm tr\, }anspose \nabla \Upgamma\, (\Uptheta_t + \varepsilon \qsaprobe_t)   + {\cal W}_t
		\label{e:ESC1explained}
	\end{equation}
	where ${\cal W}_t =    \check{\qsaprobe}_t     \nabla^{\rm tr\, }anspose \Upgamma\, (\Uptheta_t + \varepsilon \qsaprobe_t) \ddt  \Uptheta_t  $ is small by design of the low pass filter---consider \eqref{e:ESC_LP}, with $\alpha>0$  small.
	
	This justifies the diagram with $M_t= \check{\qsaprobe}_t {\check{\qsaprobe}_t}^{\rm tr\, }anspose $  time varying, but it is the expectation of $M_t$ that is most important in analysis.
	
	
	The pair of equations (\ref{e:ESC_LP}, \ref{e:ESC1}) is an instance of the QSA ODE \eqref{e:QSA_ODE_gen}, with $2d$-dimensional probing signal $( \qsaprobe_t,\check{\qsaprobe}_t)$.
	
	
\wham{All pass:} This is the special case in which the high pass filter is removed entirely, giving 
	\begin{equation}
		M {\widetilde{\nabla}\!}_t\Upgamma  =    \qsaprobe_t  \tfrac{1}{\varepsilon} \Upgamma\, (\Uptheta_t + \varepsilon \qsaprobe_t)    
		\label{e:ESC-0}
\end{equation}
This takes us back to the beginning:  based on \eqref{e:1qSGD} we conclude that 1qSGD is the ESC algorithm in which $\sigma=0$ in the low pass filter, and the high pass filter is unity gain.

In \Cref{s:ESCbig} we explain how the more general ESC ODE can be cast as QSA.    
	
	

\subsection{State dependent probing}
\label{s:SDprobe}

\begin{subequations}

It is time to explain the time varying probing gain, denoted $\upepsilon_t $ in  \Cref{fig:ESC}.    This will be state dependent for two important reasons:
\begin{romannum}
\item   The vector fields for ESC and the special case 1qSGD are \textit{not Lipschitz continuous} unless the objective is Lipschitz.    For 1qSGD,  the vector field includes  $    {\cal Y}n_t   = \frac{1}{\varepsilon} \Upgamma(\Uptheta_t + \varepsilon \qsaprobe_t)$,  and we show below that every ESC algorithm admits a state space representation driven by the same signal.

\item  If the observed cost $\Upgamma({\cal Z}_t)$ is large, then it makes sense to increase the exploration gain to move more quickly to a more desirable region of the parameter space.   
\end{romannum}
 Two choices for $\upepsilon_t \equiv \upepsilon(\Uptheta_t)$ are proposed here:
 \begin{align}
 \upepsilon(\theta) &= \varepsilon \hbox{\rlap{$\sqcap$}$\sqcup$}rt{1 + \Upgamma(\theta) - \Upgamma^-} 
 \label{e:upepsy1}
\\
 \upepsilon(\theta) &= \varepsilon\hbox{\rlap{$\sqcap$}$\sqcup$}rt{1+ \| \theta - \theta^\text{\tiny{\sf ctr}}\|^2 / \upsigma_p^2 }
\label{e:upepsy2}
\end{align}
where in  \eqref{e:upepsy1} the constant  $ \Upgamma^- $  chosen so that $ \Upgamma(\theta) \ge \Upgamma^-$ for all $\theta$.
In the second option,   $\theta^\text{\tiny{\sf ctr}}$ is interpreted as an a-priori estimate of $\theta^\text{\tiny\sf opt}$,  and $\upsigma_p$ plays the role of standard deviation around this prior.   
 \label{e:upepsy}
 \end{subequations}

 
The first is most intuitive, since it directly addresses (ii):  the exploration gain $\upepsilon_t$  is large when $ \Upgamma(\Uptheta_t) $ is far from its optimal value.   However, it does not lead to an online algorithm since $\Upgamma(\Uptheta_t) $ is not observed.   In a discrete-time implementation we would use the online version:
\[
 \upepsilon_{t_n}  = \varepsilon \hbox{\rlap{$\sqcap$}$\sqcup$}rt{1 + \Upgamma({\cal Y}_{t_{n-1}} ) - \Upgamma^-}
\]

In cases  \eqref{e:upepsy1} or  \eqref{e:upepsy2} we introduce the new notation,
\begin{equation}
    {\cal Y}n( \theta,\qsaprobe)    = \frac{1}{\upepsilon} \Upgamma(\theta + \upepsilon \qsaprobe)\,, 
        \label{e:normalizedObservations}
\end{equation}
with the usual shorthand  $   {\cal Y}n_t  = {\cal Y}n(\Uptheta_t,\qsaprobe_t) $.  
The use of \eqref{e:upepsy2} leads to a Lipschitz vector field for ESC:
 \begin{proposition}
 \label[proposition]{t:Lip_qSGD} 
 The function $   {\cal Y}n$ defined in \eqref{e:normalizedObservations}
is uniformly Lipschitz continuous in $\theta$ under either of the following assumptions on $\upepsilon = \upepsilon(\theta)$,  and with   objective $\Upgamma$ whose gradient is uniformly Lipschitz continuous:  
 \begin{romannum}
\item 
 $\upepsilon$ is defined by \eqref{e:upepsy2}.
 
\item 
$\upepsilon$ is defined by \eqref{e:upepsy1},   
and
there is $\delta>0$ such that  $\|\nabla \Upgamma (\theta)\|  \ge \delta \| \theta\| $ whenever $\|\theta\| \ge \delta^{-1}$.     
 \end{romannum}
Moreover, under either (i) or (ii), the following approximation holds,
\begin{equation}
	{\cal Y}n( \theta,\qsaprobe)    =    \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta)   +   \qsaprobe^{\rm tr\, }anspose \nabla\, \Upgamma(\theta)   +  O(\upepsilon)     \, ,
	\label{e:clYnTS}
\end{equation}
where the error term $ O(\upepsilon)$ is bounded by a fixed constant times $\upepsilon(\theta)$.  
 \end{proposition}
 
 \begin{proof}
 The proof follows from a variation on \eqref{e:1qSPSA_FTC} and the definition \eqref{e:normalizedObservations}:
\[
    {\cal Y}n( \theta,\qsaprobe)    = 
\qsaprobe^{\rm tr\, }anspose   \int_{0}^1   \nabla \Upgamma\,( \theta + t  \upepsilon(\theta)\qsaprobe)  \, dt   
						 +				 
    \frac{1}{\upepsilon(\theta)} \Upgamma(\theta )  
 \]
 The integrand in the first term is Lipschitz in $\theta$,  since the composition of Lipschitz functions is Lipschitz.  The second term is Lipschitz since its gradient is bounded under the given assumptions.
\end{proof}

\subsection{ESC and QSA}
\label{s:ESCbig}

We demonstrate here that the general ESC ODE can be interpreted as QSA,    driven by the observation process   $   {\cal Y}n_t  = {\cal Y}n(\Uptheta_t,\qsaprobe_t) $.  We find that a more informative interpretation of ESC is described as a \textit{two time-scale} variant of the QSA ODE.
  
  
For simplicity we opt for the first order low pass filter,
\begin{equation}
\ddt \Uptheta_t  = -\sigma \Uptheta_t  -   \alpha  {\widetilde{\nabla}\!}_t\Upgamma\,,\qquad{\widetilde{\nabla}\!}_t\Upgamma  =  \check{\qsaprobe}_t \check{\mathcal{Y}}^{\text{\footnotesize\sf n}}_t
\label{e:LPss}
\end{equation}
 The high pass filter is a general linear filter of order $q\ge 1$,   with state space realization defined by matrices  $(F,G,H,J)$  of compatible dimension.  For a scalar input $u_t$,  with the output of the high pass filter denoted $y_t$,  and the $q$-dimensional state  at time $t$ is denoted $Z_t$, we then have by definition,  
\begin{equation}
\begin{aligned}
\ddt Z_t   &=  F Z_t + G u_t
   \\
     y_t        & =   H^{\rm tr\, }anspose  Z_t + J u_t
\end{aligned} 
\label{e:HPss}
\end{equation}
In the architecture \Cref{fig:ESC} the high pass filter is used for  $d+1$ different choices of input and zero initial conditions are assumed: 
when $u_t={\cal Y}n_t$, the output is $y_t = \check{\mathcal{Y}}^{\text{\footnotesize\sf n}}_t$,   and the $i$th component   $\check{\qsaprobe}_t^i$  of the filtered probing signal is the output with $u_t = \qsaprobe_t^i$.

\begin{proposition}
\label[proposition]{t:ESC=QSA}
The ESC ODE obtained using the low pass filter \eqref{e:LPss} and high pass filter \eqref{e:HPss} is itself a state space model with $(d+q)$-dimensional state $X_t = [\Uptheta_t^{\rm tr\, }anspose,  Z_t^{\rm tr\, }anspose]^{\rm tr\, }anspose$.  The nonlinear dynamics can be expressed 
	\begin{equation}
	\begin{aligned}
		\ddt X_t  &=  
		\alpha  \begin{bmatrix}    
			- \frac{\sigma}{\alpha}  I   &  -     \check{\qsaprobe}_t \text{\rm H} ^{\rm tr\, }anspose
			\\
			0   & \frac{1}{\alpha}   \text{\rm F}   
		\end{bmatrix}  X_t   
		+
		\alpha \begin{bmatrix}    
			-  \text{\rm J}  \check{\qsaprobe}_t   
			\\
			\frac{1}{\alpha}\text{\rm G} 
		\end{bmatrix}   {\cal Y}n_t  \, ,
	\end{aligned} 
	\label{e:ESC=QSA}
\end{equation}
It is a controlled nonlinear state space model with input $(\qsaprobe, \check{\qsaprobe})$ for any choice of $\upepsilon$ in  \eqref{e:upepsy}. 
\end{proposition}

\Cref{t:PMF} can be freely applied to the state space representation \eqref{e:ESC=QSA} because the theorem makes no assumptions on the magnitude of $\alpha$, or even stability of the QSA ODE.

Recall that $\alpha$ is a fixed constant, so the fact that $f$ depends on this gain is irrelevant in the definition for the QSA vector field,
\begin{equation}
	f(x,\check{\qsaprobe}, \qsaprobe) =    
	\begin{bmatrix}    
		- \frac{\sigma}{\alpha}  I   &  -     \check{\qsaprobe} \text{\rm H} ^{\rm tr\, }anspose
		\\
		0   & \frac{1}{\alpha}   \text{\rm F}   
	\end{bmatrix}  x
	+
	\begin{bmatrix}    
		-  \text{\rm J}  \check{\qsaprobe}   
		\\
		\frac{1}{\alpha}\text{\rm G} 
	\end{bmatrix}   {\cal Y}n (\theta,\qsaprobe) \, .
	\label{e:ESCvectorField}
\end{equation}
Three solutions to Poisson's equation are required to write down the P-mean flow:
\whamrm{(i)}  The solution 
$\widehat{\mathcal{Y}}^{\text{\tiny\sf n}}$ with forcing function ${\cal Y}n$.
\whamrm{(ii)} $ \widehat{\cqsaprobe}$   with forcing function $ \check{\qsaprobe}$ (similar to ${\widehat G}$ in \eqref{e:ESCapprox_haf}).  
\whamrm{(iii)}    $ \widehat{Q}$   with forcing  function $Q(\theta,\Phi) = -  \text{\rm J}  \check{\qsaprobe}{\cal Y}n (\theta,\qsaprobe)$.

\begin{subequations}
	
	
	\begin{theorem}
		\label{t:PMFESC}
		If $\Upgamma$ is analytic, the  P-mean flow representation holds,
		\begin{equation} 
			\ddt X_t    =  \alpha[ {\widebar{f}} (X_t)    - \alpha  \widebar{\Upupsilon}_t  +   {\cal W}_t] \,,
			\label{e:ESC=QSA-PMF}
		\end{equation}
		in which
		\begin{equation}
			{\widebar{f}}(x) =    
			\begin{bmatrix}    
				- \frac{\sigma}{\alpha}  I   & 0
				\\
				0   & \frac{1}{\alpha}   \text{\rm F}   
			\end{bmatrix}  x
			+
			\begin{bmatrix}    
				-  \text{\rm J}  {\sf E}[  \check{\qsaprobe} (\Phi)  {\cal Y}n (\theta,\qsaprobe]
				\\
				\frac{1}{\alpha}\text{\rm G} {\sf E}[  {\cal Y}n (\theta,\qsaprobe)  ]
			\end{bmatrix}    \, ,
			\label{e:SuperBarf}
		\end{equation}
		with  expectations in steady-state.
		The functions   ${\hat f}$    and $\widebar{\Upupsilon}$ admit the representations,  
		\begin{align}
			{\hat f}(x,z) & =    
			\begin{bmatrix}    
				0	   &  -     \widehat{\cqsaprobe}(z )\text{\rm H} ^{\rm tr\, }anspose
				\\
				0   & 0   
			\end{bmatrix}  x
			+
			\begin{bmatrix}    
				-  \text{\rm J}  \widehat{Q}(\theta,z)
				\\
				\frac{1}{\alpha}\text{\rm G} \widehat{\mathcal{Y}}^{\text{\tiny\sf n}} (\theta,z)
			\end{bmatrix}   \,  ,
			\label{e:BigESCfish}
			\\
			\Upupsilon(x,z) & = \frac{1}{\alpha} 
			\begin{bmatrix} \widehat{\cqsaprobe}(z) \text{\rm H}^{\rm tr\, }anspose \{ \text{\rm F} \zeta + \text{\rm G}{\cal Y}n(\theta,\qsaprobe(z) )\}
				\\
				0
			\end{bmatrix} \, .
			\label{e:UpOopsESC}
		\end{align}
	\end{theorem}
\end{subequations}

\begin{proof}
 The expression  \eqref{e:BigESCfish} follows directly from \eqref{e:ESCvectorField}.  There is simplification because terms not involving  $   \qsaprobe$ or $\check{\qsaprobe}$ vanish.  The formula \eqref{e:UpOopsESC} then follows from the definition   $\Upupsilon(x,z)  = -\partial_x{\hat f}(x,z) f(x,z)   $.

\end{proof}

The interpretation of the P-mean field representation is entirely different here because $\Upupsilon$ is no longer a nuisance term, but a critical part of the dynamics.
Applications to design remains a topic for future research.

\wham{Interpretation as two time-scale QSA:}
The ODE \eqref{e:ESC=QSA} falls into this category,  where   $(Z_t, \Phi_t)$ represents the fast state variables. The basic principle is that we can approximate 
the solution to the $(d+q)$-dimensional state space model through the following quasi-static analysis.    For a given time $t$,   let $\{  {\bar{Z}}_r  : r\ge t\}$ denote the solution to the state space model defining ${\bf m}Z$ with $\Uptheta_r \equiv \theta$ for all $-\infty < r<\infty$:
\[
 {\bar{Z}}_r (\theta)  = \int_{-\infty}^r  e^{F (r-\tau)} G Y_\tau\, d\tau\,, \quad 
 					\textit{with} \ \  Y_\tau =  \frac{1}{\upepsilon(\theta) }   {\cal Y}(\theta + \upepsilon(\theta) \qsaprobe_\tau )   
 \]
The next step is to substitute the solution to obtain the approximate dynamics,
\begin{equation}
\ddt \Uptheta_t  \approx       
			-\sigma \Uptheta_t     - \alpha     \bigl[ \check{\qsaprobe}_t H^{\rm tr\, }anspose {\bar{Z}}_t(\Uptheta_t)  			   +   J \check{\qsaprobe}_t    {\cal Y}n_t \bigr]
\label{e:QSAtwo-time}
\end{equation}
A more useful approximation is obtained on applying \Cref{t:Lip_qSGD}:
\[
\begin{aligned}
	{\cal Y}n_t
	& \approx  \frac{1}{\upepsilon(\theta)}   \Upgamma(\Uptheta_t)   +   \qsaprobe_t^{\rm tr\, }anspose \nabla \Upgamma\, (\Uptheta_t) \, ,
	\\
	\text{\rm H} ^{\rm tr\, }anspose {\bar{Z}}(\theta,  \Phi_t)  
	&\approx  \text{\rm H} ^{\rm tr\, }anspose  \int_{-\infty}^t  e^{\text{\rm F}  (t-\tau)} \text{\rm G} \bigl\{  
	\frac{1}{\upepsilon(\theta)}   \Upgamma(\theta)   +   \qsaprobe_\tau^{\rm tr\, }anspose \nabla\, \Upgamma(\theta)  \bigr\}\, d\tau   
	\\
	& =h_0  \frac{1}{\upepsilon(\theta)}  \Upgamma(\theta)   +    \check{\qsaprobe}_t^{\rm tr\, }anspose \nabla\, \Upgamma(\theta) \,  ,
\end{aligned} 
\]
where $h_0 = -\text{\rm H} ^{\rm tr\, }anspose \text{\rm F}^{-1} \text{\rm G}$ is the DC gain of the high pass filter.

Substitution in \eqref{e:QSAtwo-time} gives
\[
\begin{aligned} 
	\ddt \Uptheta_t  \approx       
	-\sigma \Uptheta_t   &  - \alpha M_t \nabla  \Upgamma(\Uptheta_t) 
	- \alpha \check{\qsaprobe}_t    \frac{1}{\upepsilon_t} (h_0+\text{\rm J}) \Upgamma(\Uptheta_t)     \,,
	\\
	&\textit{with}\quad M_t = \check{\qsaprobe}_t   [\check{\qsaprobe}_t  + \text{\rm J} \qsaprobe_t ]^{\rm tr\, }anspose  \, .
\end{aligned} 
\]

Once tight error bounds in these approximations are established,  design and analysis proceeds based on the $d$-dimensional QSA approximation.   
In particular, without approximation this QSA ODE is stable provided the high pass filter is \textit{passive}, such as a lead compensator.   Passivity combined with the positivity of $\Sigma_{\text{\tiny$\qsaprobe$}}$ implies that $M+M^{\rm tr\, }anspose>0$ with
$M = {\sf E}_\uppi  [ M_t]$.

Theory for two time-scale QSA is not yet available to justify this approximation.   
The extension of stochastic theory is a topic of current research.

\subsection{Perturbative mean flow for qSGD}

\Cref{t:bigQSA+ESC} below is a summary of the conclusions of  \Cref{t:PMF} and  \Cref{t:Couple_theta} for qSGD ODEs.

We obtain some additional local stability structure, in terms of Lyapunov exponents, defined in terms of the sensitivity process
\begin{equation}
{\cal S}_t \mathbin{:=} \frac{\partial }{ \partial\Uptheta_0} \Uptheta_t
\label{e:SensitivityODEstate}
\end{equation}
It is a
 solution to the time varying linear system
 \begin{equation}
\ddt {\cal S}_t  =   \alpha A(\Uptheta_t, \qsaprobe_t) {\cal S}_t\,,\qquad {\cal S}_0 = I  \ (d\times d)
\label{e:SensitivityODE}
\end{equation}
 with $ A( \theta ,  \qsaprobe)$ defined in \eqref{e:barfDer}.   The
Lyapunov exponent is defined as  
\begin{equation}
\Lambda_{\Uptheta} = \lim_{t\to\infty}\tfrac{1}{t} \log (\| {\cal S}_t \|)
\label{e:LyapExp}
\end{equation}
If $\Lambda_{\Uptheta} $ is negative, then the trajectories of the QSA ODE  from distinct initial conditions couple  in a topological sense.  It will be seen in \Cref{t:bigQSA+ESC}  that this holds true for 2qSGD subject to conditions on the objective function and the gains.  

  In \eqref{e:LTIESClin} the  time $t_0$ is chosen so that $\| \Uptheta_t - \theta^*\| \le B \alpha $ for some $B$, and all $t\ge t_0$.   

\begin{subequations}

 
\begin{theorem}[QSA Theory for ESC]
\label[theorem]{t:bigQSA+ESC}
Consider an objective function $\Upgamma$ satisfying the following:  
\begin{romannum}
\item    $\Upgamma$ is analytic with Lipschitz gradient satisfying  $\|  \nabla\Upgamma (\theta)\| \ge \delta\|\theta\| $  for some $\delta>0$ and all $\theta$ satisfying $\|\theta\|\ge \delta^{-1}$.   
   
\item  The objective has a  unique minimizer $\theta^\text{\tiny\sf opt}$,  and it is the only solution to $\nabla\Upgamma (\theta) ={\sf Z}ero$.

\item  $P = \nabla^2\Upgamma\, (\theta^\text{\tiny\sf opt})$ positive definite.

\end{romannum} 

Consider 1qSGD or 2qSGD with constant gain $\alpha$.  Suppose that the probing signal is chosen of the form 		\eqref{e:qSGD_probe} with $K=d$,    frequencies satisfying \eqref{e:logFreq}, and
   $\Sigma_{\text{\tiny$\qsaprobe$}}>0$.

Then, for either algorithm, there are positive constants  $\varepsilon^0$ and $\alpha^0$ such that the following approximations are valid for $0<\varepsilon\le \varepsilon^0$,  $0<\alpha\le \alpha^0$:     

\begin{align}
   \ddt \Uptheta_t   &=    \alpha \bigl[ -  \Sigma_{\text{\tiny$\qsaprobe$}} \nabla\Upgamma\, (\Uptheta_t)  + {\cal W}_ t  + O(  \varepsilon^2)  \bigr]  \,,  &&  t\ge 0
\label{e:LTIESC}
\\
   &=    \alpha \bigl[ - \Sigma_{\text{\tiny$\qsaprobe$}} P  [\Uptheta_t  -\theta^\text{\tiny\sf opt}]  + {\cal W}_ t  + O(\alpha^2 +   \varepsilon^2) +o(1)  \bigr]  \,, && t\ge t_0
\label{e:LTIESClin}
\end{align} 
Moreover, the
 following hold for $0<\varepsilon\le \varepsilon^0$,  $0<\alpha\le \alpha^0$:
\begin{romannum}
\item   $\|  \Uptheta_t -\theta^\text{\tiny\sf opt} \| \le O( \alpha+\varepsilon^2) + o(1)$  

\item   $\|  \UpthetaF_t -\theta^\text{\tiny\sf opt} \| \le O( \alpha^2+\varepsilon^2) + o(1)$,  with filtered estimate obtained using the criteria of \Cref{t:Couple_theta}.
\end{romannum}
In addition,  the Lyapunov exponent \eqref{e:LyapExp} is negative in the special case of   2qSGD.
\end{theorem}
 
\end{subequations}
\begin{proof}
The analytic assumption is imposed so that (A3) holds, thanks to  \Cref{t:hah}. 

\Cref{t:Lip_implies_V4}  tells us that  the QSA ODE is $\alpha^0$   ultimately bounded  for some $\alpha^0>0$.   
The estimation error bounds in (i) and (ii) are then immediate from   \Cref{t:Couple_theta}.

Justification of \eqref{e:LTIESC} follows from \Cref{t:PMF} and the approximation $  
{\widebar{f}}( \theta)   = - \Sigma_{\text{\tiny$\qsaprobe$}}\nabla\, \Upgamma(\theta)   +  O(\varepsilon^2)   $ obtained for 1qSGD in \Cref{t:ESCapproximation}    and for 2qSGD in \Cref{t:ESCapproximation2}.    
 The term $\widebar{\Upupsilon}_t$ appearing in \eqref{e:BigGlobalODE} is zero on application of \Cref{t:orthCorBigPre}.

The approximation \eqref{e:LTIESClin} then follows from  \eqref{e:LTIESC}  and the Taylor series approximation
\[
\nabla\Upgamma\, (\Uptheta_t)     =   \nabla\Upgamma\, ( \theta^\text{\tiny\sf opt})   + P  [\Uptheta_t  -\theta^\text{\tiny\sf opt}]   +   O (\| \Uptheta_t  -\theta^\text{\tiny\sf opt} \|^2)
\]  
The first term above is zero by assumption, $   \nabla\Upgamma\, ( \theta^\text{\tiny\sf opt})   ={\sf Z}ero$,   and $ \| \Uptheta_t  -\theta^\text{\tiny\sf opt} \|  \le \| \Uptheta_t  -\theta^\text{\rm s.t.\,}arp \|  +  \| \theta^\text{\rm s.t.\,}arp  -\theta^\text{\tiny\sf opt} \|  = O(\alpha + \varepsilon^2)$ for   $t\ge t_0$,  giving $ \| \Uptheta_t  -\theta^\text{\tiny\sf opt} \|^2  \le  O(\alpha^2 + \varepsilon^4)$  (recall that $\| \Uptheta_t - \theta^*\| \le B \alpha $ for $t\ge t_0$).
These approximations imply \eqref{e:LTIESClin}.

 
  
The proof that $\Lambda_{\Uptheta} <0$ for 2qSGD  follows from identifying the linearization matrix used in the LTI approximation \eqref{e:BigGlobalODElin}.   The following is established in  \Cref{t:ESCapproximation2}:  
\[
A( \theta,\qsaprobe)    =   A^0( \qsaprobe) +  O(\varepsilon)  \,,\qquad \textit{with }   A^0( \qsaprobe)  =
 - \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla^2\, \Upgamma(\theta)   
\]
in which the error $O(\varepsilon)$ depends smoothly on $\theta$.       
Hence  the dynamics  \eqref{e:SensitivityODE} for the 
 sensitivity process can be expressed
\begin{equation}
\ddt {\cal S}_t  = -  \alpha  [ (\alpha+ \varepsilon^2) M_t +  \qsaprobe_t  \qsaprobe^{\rm tr\, }anspose_t P  ] {\cal S}_t\, ,
						 \qquad t\ge t_0\,,
\label{e:SensitivityODE2aSGD}
\end{equation}
where  the norm of the matrix process $\{ M_t : t\ge t_0 \}$ is bounded by a constant that is
 independent of $\varepsilon$ or $\alpha$,  though $t_0$ is dependent on these parameters.
 
 The remainder of the proof is then identical to the proof of  \Cref{t:QPL}.     

\end{proof}

Approximations in \Cref{s:ESCtheory} suggest that a negative Lyapunov exponent is not likely for any variant of 1qSGD.    We obtain in this case,
\[
A( \theta,\qsaprobe)   =  
 -  \partial_\theta{\sf B}igl( \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta) {\sf B}igr) \qsaprobe
 - \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla^2\, \Upgamma(\theta)   +  O(\varepsilon)  
\]
If $\upepsilon =\varepsilon$, independent of $\theta$,  this simplifies to 
\[
A( \theta,\qsaprobe)   =  
 -   {\sf B}igl( \frac{1}{ \varepsilon}    \qsaprobe   \nabla^{\rm tr\, }anspose  \Upgamma\, (\theta) 
+ \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla^2\, \Upgamma(\theta)  {\sf B}igr)  +  O(\varepsilon)  
\]
It is likely that the Lyapunov exponent is negative for small $\varepsilon>0$ and  $\alpha\ll \varepsilon$, but this is not an interesting setting.   

 
 
 
 
\section{Examples}
\label{s:Experiments}

In each of the following experiments\footnote{Publicly available code for experiments were obtained under a GNU General Public License v2.0.}, the QSA ODEs were implemented using an Euler approximation with time discretization of $1~$sec. The probing signal respected \eqref{e:qSGD_probe} with $\{\omega_i\}$ irrationally related. 

\subsection{Optimization of Rastrigin's Objective}
\label{s:Opt_Rast_NIPS}
We illustrate the fast rates of convergence of QSA with vanishing gain by surveying results from an experiment in \cite{laumey22}. 
\paragraph{Simulation Setup}

The $1$qSGD algorithm \eqref{e:1qSGD} and its stochastic counterpart \eqref{e:1qSPSA} were used for the purpose of optimizing the Rastrigin Objective in $\field{R}^2$
\begin{equation}

\Upgamma(\theta) = 20 +\sum _{i=1}^2\left[\theta_{i}^{2}-10\cos(2\pi \theta_{i})\right]  
\label{e:Rastrigin}
\end{equation}  
The experiments were conducted before discovery of the Lipschitz variant of $1$qSGD so $\upepsilon(\theta) = \varepsilon=0.25$ and the sample paths of $\bm{\Uptheta}$ were projected onto $[-B,B]^2$ with $B=5.12$, which is the evaluation domain commonly used for this objective \cite{simulationlib}. The QSA ODE was implemented through an Euler scheme with time sampling of $1~$sec with $a_t = \min\{0.5,(t+1)^{-0.85}\}$. $M=200$ independent experiments were run with $\Uptheta^m_0$ uniformly sampled from $[-B,B]^2$ for $\{1 \leq m \leq M\}$. For each experiment, the probing signal was selected of the form
\begin{equation}
	\qsaprobe^m_t = 2 [\sin(t/4+ \phi^m), \sin(t/e^2+ \phi^m)]^{\rm tr\, }anspose
	\label{e:exp_qsaprobe}
\end{equation}
where for each $m$, $\phi^m$ was uniformly sampled from $[-\pi/2, \pi/2]$. 
The noise for stochastic algorithm was a scaled and shifted Bernoulli with $p=1/2$ and had support in the set $\{-\hbox{\rlap{$\sqcap$}$\sqcup$}rt{2}, \hbox{\rlap{$\sqcap$}$\sqcup$}rt{2}\}$. The values $\pm \hbox{\rlap{$\sqcap$}$\sqcup$}rt{2}$ were chosen so that the covariance matrices in both algorithms were equal.

In order to test theory for rates of convergence with acceleration through PR averaging, the empirical covariance was computed for $1$qSGD across all initial conditions,
\begin{equation}
	\begin{aligned}
	{\bar{S}}igma_T & = \frac{ 1}{M} \sum_{i=1}^{M} \Uptheta_T^i  {\Uptheta_T^i }^{\rm tr\, }anspose   -  \widebar{\Uptheta}_T \widebar{\Uptheta}_T^{\rm tr\, }anspose\,,\qquad \widebar{\Uptheta}_T= \frac{ 1}{M} \sum_{i=1}^{M} \Uptheta_T^i \\ 
		\hat{\sigma}_T & = \hbox{\rlap{$\sqcap$}$\sqcup$}rt{{\rm tr\, }({\bar{S}}igma_T)}
	\end{aligned}
	\label{e:empCov}
\end{equation}
In each experiment the trajectories $\{\Uptheta_t\}$ were accelerated through PR averaging with $\kappa =5$. If the scaled empirical covariance $a_T^{-2} \hat{\sigma}_T$ is bounded in $T$, we have $\|\UpthetaPR_T -\theta^* \| = O(T^{-2\rho})$.

\begin{figure}
\caption{Rastrigin objective (top left), scaled empirical covariance (bottom left), histograms of estimation error for $1$SPSA with PR averaging (top middle and top right), histograms of estimation error for $1$qSGD with PR averaging (bottom middle and bottom right).}
\label{fig:Rastrigin_results}
\end{figure} 

\paragraph{Results}
\Cref{fig:Rastrigin_results} shows plots for the scaled empirical covariance for $1$qSGD as well as histograms for the estimation error $\UpthetaPR_T - \theta^*$ for both algorithms. We see in \Cref{fig:Rastrigin_results} that the convergence rates of $O(T^{-2\rho})$ are achieved by both filtering techniques. 

The variance of the estimation error for the deterministic algorithm is much smaller than for its stochastic counterpart: \Cref{fig:Rastrigin_results} shows that the reduction is roughly two orders of magnitude. 

Outliers were identified with Matlab's \textit{isoutlier} function and removed from the histograms in \Cref{fig:Rastrigin_results}. Around $40\%$ of the estimates were considered outliers for the stochastic algorithm, while none was observed for its deterministic counterpart. 

Recall that bias is inherent in both 1qSGD and 1SPSA when using a non-vanishing probing gain.
It would appear from \Cref{fig:Rastrigin_results} that this bias is larger for 1qSGD  [consider the histogram for estimates of $\theta_1^\text{\rm s.t.\,}arp=0$, which shows that final estimates typically exceed this value].    In fact, there is no theory  to predict if 1qSGD is better or worse than 1SPSA in terms of bias.    The bias is imperceptible in the stochastic case due to the higher variance combined with the removal of outliers.  

\subsection{Vanishing vs Fixed Gain Algorithms -  Rastrigin's Objective}
\label{s:vanishing_vs_fixed_gain}
Experiments were performed to test   the performance of the Lipschitz version of $1$qSGD for both constant and vanishing gain algorithms.   Two take-aways from the numerical results surveyed below:
\begin{romannum}
	\item 
	The value of a second order filter can be substantial when using a fixed gain algorithm.
	Recall that this conclusion was not at all clear for either of the linear examples 	(\ref{e:lin_barY0}, \ref{e:lin_barYnot}),  even though in the first case the example satisfies the assumptions of \Cref{t:Couple_theta}.
	
	
	\item  
	However, in this example the flexibility of a vanishing gain algorithm is evident: 
	with high gain during the start of the run  (when $a_t\approx \alpha$) there is a great deal of exploration,  especially when $\Uptheta_t$ is far from the optimizer  due to the use of $\upepsilon(\theta)$.
	The vanishing gain combined with the use of distinct frequencies  then justifies the use of PR-averaging, which results in extremely low AAD in this example.
\end{romannum}

\paragraph{Simulation Setup}
For this objective, $M=5$ individual experiments were carried out in a similar fashion to the experiments summarized in \Cref{s:Opt_Rast_NIPS}. The Lipschitz variant of $1$qSGD was implemented using each of the initial conditions so projection of the sample paths $\{\Uptheta^m_t\} $ was not necessary. Each experiment used $\upepsilon(\theta) = \varepsilon \hbox{\rlap{$\sqcap$}$\sqcup$}rt{1 + \|\theta \|^2}$ with $\varepsilon = 0.6$ and the same probing signal $\qsaprobe^m_t$ as in \eqref{e:exp_qsaprobe}. The following three choices of gain were tested for each initial condition:
\[
\text{(a) } a_t = 0.1 (t+1)^{-0.65}  
\quad \text{(b) } a_t = \alpha_b = 3 \times 10^{-3} 
\quad \text{(c) } a_t = \alpha_s = 7 \times 10^{-4}
\]
For the vanishing gain case (a) the process $\{\Uptheta^m_t\}$ was accelerated through PR averaging with $\kappa =5$, while for the fixed-gain case two filters of the form \eqref{e:filter_numerics} where used with $\zeta =0.8$ and $\gamma=\eta  \alpha$ for various values of $\eta \geq 1$. These experiments were repeated for $\Uptheta^m_0$ uniformly sampled from $10^{10}[-B,B]^2$ with $B =5.12$.

\paragraph{Results}
The top row of \Cref{fig:Rastrigin_comp} shows the evolution of $\{\Uptheta^m_t\}$ for each $m$ with $\eta=5$. The second row shows the evolution of $\{\Upgamma(\UpthetaPR_T),\Upgamma({\Uptheta}^{1 \text{\tiny\sf F}}_T),\Upgamma({\Uptheta}^{2 \text{\tiny\sf F}}_T),\Upgamma(\Uptheta_T)\}$ for the single path yielding the best performance for each gain choice across all 5 initial conditions. 

We see in \Cref{fig:Rastrigin_comp} a clear advantage of vanishing gain algorithms: the algorithm explores much more in the beginning of the run, eventually deviating from $\theta^\text{\tiny\sf opt}$ by the bias inherent to the $1$qSGD algorithm.

For the runs that used $\alpha_b$ (case (b)), we have a good amount of exploration, but the steady state behavior is poor. Case (c) using the smaller value of $\alpha$ yielded better results.

\Cref{fig:Rastrigin_comp} shows the benefit of variance reduction from second order filter as opposed to a first order filter, based on runs that used $\alpha_s$. As opposed to the results in \Cref{fig:alpha2bias}, this example shows that we can't always obtain the best AAD with a first order filter. For these experiments, the filtered final estimation of $\theta^*$ yields comparable objective values to what is seen in the vanishing gain experiment. However, we do not always observe desirable exploration as in the vanishing gain case: we see several sample paths hovering around local minima of the objective rather than $\theta^\text{\tiny\sf opt}$.

Results for the trajectories with initial conditions of order $10^{10}$ were presented in \Cref{fig:ODEinfty}.

\begin{figure}
\caption{Comparison of the performance between vanishing and constant gain $1$qSGD for the Rastrigin Objective.}
\label{fig:Rastrigin_comp}
\end{figure} 

\subsection{The Walking Camel} 
\label{s:Camel_Walk}

The next numerical results are based on a tracking problem using the $1$qSGD algorithm. Although the results in this subsection will show that the algorithm successfully tracks the signal $\{\theta^\text{\tiny\sf opt}_t\}$, we are not advocating the use of $1$qSGD for tracking applications. The algorithm was employed to illustrate general algorithm design principles;
better performance is likely in algorithms with a carefully designed low pass filter such as \eqref{e:ESC_LP} with $\sigma>0$. Several experiments were conducted to address the following questions: 
\begin{romannum}
	\item 
	How does the frequency content of the probing signal affect performance?
	
	\item
	How does the rate of change of the signal $\{\theta^\text{\tiny\sf opt}_t\}$ impact tracking performance?
	
	\item  
	How does filtering reduce estimation error and variance in tracking?
	
\end{romannum}

\begin{figure}
\caption{Tracking for the travelling camel with ESC-0}
\label{fig:moving_camel_lotus}
\end{figure} 

\paragraph{Simulation Setup}

The \textit{Three-Hump Camel} objective is the sixth-order polynomial function on $\field{R}^2$:
\begin{equation}
  \Upgamma(\theta) =2\theta_1^{2}     +\theta_1\theta_2+\theta_2^{2}
			+\tfrac{1}{6} \theta_1^6 -1.05\theta_1^{4}
\label{e:3Hump}
\end{equation}
A plot of the negative of this objective is shown on the left hand side of \Cref{fig:moving_camel_lotus}.   

Observe that the infimum of $\Upgamma$ over $\field{R}^2$ is not finite. 
The domain of this function is usually taken to be the bounded region $[-5,5]^2$   \cite{simulationlib}. 

We considered a time varying version, of the form
\[
\Upgamma_t = \Upgamma(\theta - \theta^\text{\tiny\sf opt}_t)
\]
with the following choices of $\{\theta^\text{\tiny\sf opt}_t\}$:
\begin{equation}
	\begin{aligned}
		&\text{(a) }  \theta^\text{\tiny\sf opt}_t = \begin{bmatrix}
		 m \cos(\omega^\circ t) - h\cos(m\omega^\circ t/b^\circ) \\                                                                                      
		 m \sin(\omega^\circ t) -h\sin(m\omega^\circ t/b^\circ)
		\end{bmatrix}
	 \quad
		\text{where }
		 \omega^\circ = 2 \times 10^{-3}~\text{rad/sec} , b^\circ = 3/5, m =8/5, h = 0.6	\\
		&\text{(b) }\theta^\text{\tiny\sf opt}_t =  
		\begin{bmatrix}
			1\\1
		\end{bmatrix} 
	g_t, 
	\quad \text{where } 
	g_t=	\begin{cases} 
			\Delta(t/T_0),& t\leq T/2 \\
			\Pi(t/T_0) ,& t>T/2 
		\end{cases} \quad
		\text{and }  
		T_0=  \frac{c}{b} 
	~\text{sec with } b \in \{ 1, 3 , 5   \}, c = 5 \times 10^{-3} 
\end{aligned}
\label{e:theta_opt}
\end{equation}   
The choice (a) is an example of an epitrochoid path.  In (b) the notation $\Delta,\Pi$ is used to denote the unit triangle and square waves, respectively. The objective $\Upgamma_t$ was normalized so that $\theta^\text{\tiny\sf opt} = {\sf Z}ero$.

Experiments were conducted for a single initial condition $\Uptheta_0$, chosen to be 
		a
 second (non-optimal) local minima of the objective, and the following probing signals:
\begin{equation}
\qsaprobe^a_t = 2 \begin{bmatrix}\sin(t/9)\\\sin(t/e^3) \end{bmatrix}, \quad \qsaprobe^b_t = 2 \begin{bmatrix} \sin(t/4) \\ \sin(t/e^2) \end{bmatrix},\quad \qsaprobe^c_t = 2 \begin{bmatrix} \sin(\hbox{\rlap{$\sqcap$}$\sqcup$}rt{20}t/10) \\ \sin(\pi t/10)  \end{bmatrix}
\label{e:probingcamel}
\end{equation}
The assumptions of \Cref{t:Lip_qSGD}  are violated in this example, since the gradient of this objective is not Lipschitz continuous.     For this reason, the state-dependent gain  $\upepsilon$ was abandoned in favor of a constant gain in these experiments.
Choices for the parameters $\upepsilon(\theta) \equiv \varepsilon$ and $\alpha$ were problem specific: (a) $\varepsilon=0.2$ and $\alpha = 6 \times10^{-3}$, (b) $\varepsilon=0.5$ and $\alpha = 5\times10^{-3}$. 

The process $\{\Uptheta_t\}$ was projected onto the set $[-5,5]^2$.
The filters in \eqref{e:filter_numerics} were used for variance reduction,
 with $\zeta =0.8$ and $\gamma= \eta \alpha$ for several $\eta \geq 1$.

\begin{figure}
\caption{Tracking of $\{\theta^\text{\tiny\sf opt}
\label{fig:moving_camel_qsaprobe}
\end{figure} 

\bigbreak

\wham{Results}  

\wham{Case (a)}
Phase plots of the trajectories $\{\Uptheta_t\},\{{\Uptheta}^{1 \text{\tiny\sf F}}_t \} $ are shown in \Cref{fig:moving_camel_lotus}. These plots display results with the probing signal $\qsaprobe^b_t$ and $\eta =5$. It is possible to see in \Cref{fig:moving_camel_lotus} that the output of $1$qSGD tracks $\{ \theta^\text{\tiny\sf opt}_t\}$ after a transient period. Since the gain is non-vanishing, $\bm{\Uptheta}$ oscillates around the trajectory $\{\theta^\text{\tiny\sf opt}_t\}$. As expected by \Cref{t:Couple_theta}, the filtered process has much less variability.

\wham{Case (b)}

\Cref{fig:moving_camel_qsaprobe} shows the evolution of $\{\Uptheta_t\}$ using each probing signal in \eqref{e:probingcamel}. Here, the signal $\{\theta^\text{\tiny\sf opt}_t\}$ has period $T_0$ with $b=3$ \eqref{e:theta_opt}. The plots in \Cref{fig:moving_camel_qsaprobe} illustrate how the frequency content of the probing signal impacts tracking performance: higher frequencies yield less variability over most of the run.

\begin{figure}
\caption{Estimation error for different $b$ \eqref{e:theta_opt}
\label{fig:moving_camel_norm}
\end{figure} 

 

The evolution of the estimation error $\{\|\Uptheta_t - \theta^\text{\tiny\sf opt}_t \|\}$ is shown in \Cref{fig:moving_camel_norm} for each value of $b$ listed in \eqref{e:theta_opt}. Only $\qsaprobe^c_t$ in \eqref{e:probingcamel} was used as the probing signal for these experiments. We see in \Cref{fig:moving_camel_norm} that tracking performance deteriorates as the rate of change of $\{\theta^\text{\tiny\sf opt}_t\}$ increases. When $t>T/2$, the signal $\{\theta^\text{\tiny\sf opt}_t\}$ is nearly static for each period and hence we expect the rate of change of $\{\theta^\text{\tiny\sf opt}_t\}$ to have little effect on the AAD.
 This is confirmed by the results in shown in \Cref{fig:moving_camel_norm}: the steady state estimation errors are roughly the same for all values of $b$.

\Cref{fig:moving_camel_eta} shows the evolution $\{\Upgamma(\Uptheta_t), \Upgamma(  {\Uptheta}^{1 \text{\tiny\sf F}}_t  ), \Upgamma( {\Uptheta}^{2 \text{\tiny\sf F}}_t) \}$ as functions of time for $\eta \in \{5,15\}$. Again, the selected probing signal was $\qsaprobe^c_t$ and the signal $\{\theta^\text{\tiny\sf opt}_t\}$ had period $T_0$ with $b=3$ for both choices of $\eta$. We see in \Cref{fig:moving_camel_eta} that the filtered estimates $\{{\Uptheta}^{1 \text{\tiny\sf F}}_t\}$ yield lower objective values than $\{{\Uptheta}^{2 \text{\tiny\sf F}}_t\}$ for $t \leq  T/2$. When $t>T/2$, the estimates $\{{\Uptheta}^{1 \text{\tiny\sf F}}_t\}$ only yield lower objective values than $\{{\Uptheta}^{2 \text{\tiny\sf F}}_t\}$ for very small portions of the run. That is, a second order filter is preferable when the rate of change of $\{\theta^\text{\tiny\sf opt}_t\}$ is small. 

Recall that the parameter $\eta$ scales the bandwidth of the filter. We see in \Cref{fig:moving_camel_eta} an improvement in tracking for both filters when $\eta$ is increased as a consequence of the increased bandwidth.

\begin{figure}
\caption{Evolution of $\Upgamma$ for different filtering techniques.}
\label{fig:moving_camel_eta}
\end{figure} 

 

\section{Conclusions}
\label{s:conc}
 This paper introduces a new exact representation 
 of the QSA ODE, opening several doors for analysis.     
Major outcomes
of the perturbative mean flow representation include the clear path to obtaining transient bounds based on those of the mean flow, and guidelines on how to design filters to obtain bias 
	and AAD
of order $O(\alpha^2)$. There is much more to be unveiled:
 
 

\wham{\small$\triangle$} 
The use of filtering for acceleration of algorithms is not at all new.  
It will be exciting to investigate the implications of the acceleration techniques pioneered by  Polyak and Nesterov for nonlinear optimization,
particularly in their modern form  (see \cite[\S~2.2]{les22} and the references therein). 

The integration of these two disciplines may provide insight on how to design the high pass filters shown in  \Cref{fig:ESC}, or suggest entirely new architectures.   

\wham{\small$\triangle$} 
The introduction of normalization in the observations in the general form \eqref{e:normalizedObservations} was crucial to obtain global stability of ESC ODEs.   There are many improvements to consider.   First,   on considering the Taylor series approximation \eqref{e:clYnTS}, performance is most likely improved via a second normalization:
\[
{\cal Y}n_t =   \frac{1}{\upepsilon_t} {\sf B}igl[ \Upgamma(\Uptheta_t + \upepsilon_t \qsaprobe_t)  -   \Upgamma^\text{\rm s.t.\,}ar_t {\sf B}igr] \,, 
\]
in which $\{ \Upgamma^\text{\rm s.t.\,}ar_t \}$ are estimates of the minimum of the objective.   These might be obtained by passing $\{ {\cal Y}_r \mathbin{:=} \Upgamma(\Uptheta_r + \upepsilon_r \qsaprobe_r) \}$ through a low pass filter.  

\wham{\small$\triangle$} 
Far better performance might be obtained through an observation process inspired by  2SPSA.  Consider first a potential improvement of 2SPSA:   A state dependent exploration gain is introduced, so that
\eqref{e:2qSPSA} becomes
\[
\theta_{n+1} = \theta_n    -  \alpha_{n+1} \frac{1}{2\varepsilon}   \qsaprobe_{n+1} \bigl[\Upgamma(\theta_{n}  + \upepsilon_n \qsaprobe_{n+1} ) -\Upgamma(\theta_{n}  -  \upepsilon_n \qsaprobe_{n+1} ) \bigr]  \,
\] 
with $\upepsilon_n = \upepsilon( \theta_{n})$.  The division by $2\varepsilon$ (independent of state) remains, since 2SPSA in its original form satisfies the required Lipschitz conditions for SA provided $\nabla \Upgamma$ is Lipschitz continuous.

There are surely many ways to obtain an online version based on QSA.   One approach is through sampling:  denote $T_n = nT$ for  
a  given sampling interval $T>0$,  and
take ${\cal Y}n_t$  constant on each interval $[T_n,T_{n+1})  $, 
designed to mimic 2SPSA.   One option is the simple average,
\[
\begin{aligned}
	{\cal Y}n_{T_{n+1}}  \mathbin{:=} \frac{2}{T}  \frac{1}{2\varepsilon} 
	&\int_{T_n}^{T_n+T/2}   \qsaprobe_{t}  \bigl[ \Upgamma( \theta_n   +  \upepsilon_t \qsaprobe_t ) - \Upgamma(\theta_n - \upepsilon_t \qsaprobe_t ) ]  \, dt
\end{aligned} 
\]
with $\theta_n=  \Uptheta_{T_n} $.
This can be computed in real-time, based on two sets of observations:
\[
\begin{aligned}
	\Upgamma( \theta_n   +  \upepsilon_t \qsaprobe_t   )  \qquad  &    T_n \le t\le  T_n + T/2
	\\
	\Upgamma( \theta_n   -  \upepsilon_n \qsaprobe_{t-T/2}  ) \qquad   &     T_n + T/2  \le t\le T_{n+1}
\end{aligned} 
\]

\wham{\small$\triangle$}
The implications to reinforcement learning deserve much greater attention.  The applications of QSA and ESC in  \cite{mehmey09a,kildrs06,CSRL} are only a beginning.   

\wham{\small$\triangle$} 
It may be straightforward to extend the P-mean flow representation  \eqref{e:BigGlobalODE} to tracking problems.   This may start with the general topic of   time inhomogeneous QSA, of the form 
\[
\ddt\Uptheta_t = \alpha f(\Uptheta_t,\qsaprobe_t;   t)\,.
\]
Analysis would require consideration of solutions to Poisson's equation, such as ${\hat f}(\theta,\,\cdot\,; t)$ for each $\theta\in\field{R}^d$ and $t\in \field{R}$.  
The representation will be more complex than 	\eqref{e:BigGlobalODE}, but will likely lead to sharper bounds than are presently available.  

 

 

\section{Literature and Resources}
\label{s:lit}

Background on stability via Lyapunov criteria is standard  \cite{bor20a,kusyin97}.
The 
 ODE@$\infty$ approach was introduced in 
\cite{bormey00a}, and refined in a series of papers since  \cite{bha11,rambha17,chedevborkonmey21}.  

 
 
 
\wham{Quasi-stochastic approximation:} 

The first appearance of the QSA ODE is the almost  periodic linear system discussed in
\Cref{s:QPLS}---see \cite{luk82,subwasred21} for the century-old history.   

The first work on QSA as an alternative to SA appears in \cite{lappagsab90,larpag12} with applications to finance,  and \cite{bhafumarwan03} contains initial results on convergence rates that are better than that obtained for their stochastic counterparts, but nothing like the bounds obtained here.   The setting for this prior work is in discrete time.

The QSA ODE appeared in applications to reinforcement learning in \cite{mehmey09a},  which motivated the convergence rate bounds in  \cite{shimey11} for multi-dimensional linear models, of which 	\eqref{e:TwoLinExamples} are special cases (except that analysis is restricted to vanishing gain algorithms).  

The work of \cite{shimey11} was extended to nonlinear QSA in \cite{chedevbermey20b,chedevbermey21},  and initial results on PR averaging and the challenges with multiplicative noise is one topic of \cite[\S~4.9]{CSRL}.   Until recently it was believed that the best convergence rate possible is $O(1/T)$ for the vanishing gain algorithms considered.  
Recently, this was shown to be fallacy:  convergence rates arbitrarily close to $O(T^{-2})$ are obtained in  \cite{laumey22} via PR-averaging, but again this requires that a version of $\Upupsilon$ be zero. 

\wham{Extremum seeking control:} ESC is perhaps among the oldest gradient-free optimization methods. It was born in the 1920s when a extremum seeking control solution was developed to maintain maximum transfer of power between a transmission line and a tram car \cite{leb22,EShistory2010,liukrs12}. Successful stories of the application of ESC to various problems have been shared over the 20th century -- e.g. \cite{rastrigin1960extremum,ras73,drali51,mee67,oba67}, but theory was always behind practice: the first Lyapunov stability analysis for ESC algorithms appeared in the 70s for a very special case \cite{luxlee71}.    

General local stability and bounds on estimation error for ESC with scalar-valued probing signals were established 30 years later in \cite{krswan00} and later extended for multi-dimensional probing in \cite{arikrs02}, but the authors do not specify a domain of attraction for stability. Stronger stability results was than established in \cite{tannesmar06}, where the authors show that for an arbitrarily large set of initial conditions, it is possible to find control parameters so that the extermum seeking method results in a solution arbitrarily close the extrema to be found. This is known as semi-global pratical stability in the ESC literature. These analytical contributions sparked further research \cite{docpergua11,guadoc15,durstaebejoh13,pokeis22}.

It is pointed out in \cite{mil16} that there is a potential  ``curse of dimensionality'' when using sinusoidal probing signals,  which echos the well known curse in Quasi Monte Carlo \cite{asmgly07}.

\wham{Convergence rates for stochastic approximation:}

Consider first the case of vanishing step-size:
with step-size $\alpha_n = (1+n)^{-\rho}$,  $\half <\rho <1$, we obtain the following limits for the scaled  asymptotic mean-square error (AMSE):
	\begin{align}
		\lim_{n\to\infty} \frac{1}{\alpha_n} {\sf E}[ \| \theta_n -\theta^\text{\rm s.t.\,}arp\|^2]  & = {\rm tr\, }ace (\Sigma_\theta(\rho)) 
		\label{e:PRrate1}
		\\
		\lim_{n\to\infty}  n {\sf E}[ \| \theta^{\text{\tiny\sf PR}}_n -\theta^\text{\rm s.t.\,}arp\|^2]  & = {\rm tr\, }ace (\Sigma_\theta^*) 
		\label{e:PRrate2}
\end{align} 
where $\{\theta^{\text{\tiny\sf PR}}_n\}$ are obtained from $\{\theta_n\}$ via Polyak-Ruppert averaging \eqref{e:PR}.   The covariance matrix in the second limit has the explicit form $\Sigma_\theta =  G \Sigma_\Delta G^{\rm tr\, }anspose$ with $G =  {A^*}^{-1}$ and $\Sigma_\theta$ the covariance of $f(\theta^\text{\rm s.t.\,}arp,\qsaprobe_{n+1}) $   (see \eqref{e:SA_recur}).  The covariance matrix $\Sigma_\theta^*$ is  minimal in a strong sense  \cite{pasgho13,chedevborkonmey21,bor20a,kusyin97}.    

 
There is an equally long history of analysis for algorithms with constant gains.  For simplicity assume an SA algorithm with ${\bf q}saprobe$ i.i.d., so that  the stochastic process ${\bf t}heta$ is a Markov chain.   Stability of the ODE@$\infty$ implies a strong form of geometric ergodicity under mild assumptions on the algorithm
\cite{bormey00a}, which leads to bounds on the variance 
\begin{equation}
	{\sf E}[ \| \theta_n -\theta^\text{\rm s.t.\,}arp\|^2]    =
	{\sf E}[ \| \theta_\infty -\theta^\text{\rm s.t.\,}arp\|^2]    + O(\varrho^n)
	\label{e:BM}
\end{equation}
where the steady-state mean $ {\sf E}[ \| \theta_\infty -\theta^\text{\rm s.t.\,}arp\|^2] $ is of order $O(\alpha)$ with step-size $\alpha>0$,  and $\varrho<1$ since ${\bf t}heta$ is geometrically ergodic.   Averaging can reduce variance significantly \cite{moujunwaibarjor20,durmounausamscawai21}.

\wham{Zeroth order optimization:}

We let  $\theta^\text{\rm s.t.\,}arp$ denote the limit of the algorithm (assuming it exists), and  $\theta^\text{\tiny\sf opt}$ the global minimizer of $\Upgamma$. 

We begin with background on methods based on stochastic approximation,  which was termed gradient free optimization (GFO) in the introduction.
\begin{subequations}
	
\wham{1. GFO }    
Much of the theory is based on vanishing gain algorithms, and the probing gain is also taken to be vanishing in order to establish that the algorithm is \textit{asymptotically unbiased}. 

	
	

	
	
\end{subequations}

Polyak was also a major contributor to the theory of convergence rates for GFO algorithms that are asymptotically unbiased,  with \cite{poltsy90} a highly cited starting point.  If the objective function is $p$-fold differentiable at $\theta^\text{\rm s.t.\,}arp$ then the  best possible rate for the AMSE is $O(n^{-\beta} )$ with $\beta = \frac{p-1}{ 2p} $. 
This work was motivated in part by the upper bounds and algorithms of Fabian \cite{fab68b}.
See \cite{dipren97,dip03,pasgho13} for more recent history, and algorithms based on averaging that achieve this bound;  typically using   averaging combined with vanishing step-sizes.

 

\wham{2. ESC} Theory has focused on constant-gain algorithms because many authors are interesting in tracking rather than solving exactly a static optimization problem. 

It is shown in \cite{krswan00}, that it is possible to implement ESC to find a control law that drives a system near to an equilibrium state $x^\text{\tiny\sf opt}$, minimizing the steady-state value of its output at $-y^\text{\tiny\sf opt}$. It is assumed that $x^\text{\tiny\sf opt} = \Upgamma(\theta^\text{\tiny\sf opt})$ where $\Upgamma: \field{R} \to \field{R}^d$ is a smooth function. For a scalar-valued perturbation signal of the form $\qsaprobe_t = \alpha \sin (\omega t)$, the estimation error $\| y^* - y^\text{\tiny\sf opt} \|$ is of order $O(\alpha^2)$, assuming stability of the algorithm. These results are later extended to the case where $\theta,\qsaprobe \in \field{R}^d$ in \cite{arikrs02} satisfying: $\qsaprobe^j_t = \alpha \sin(\omega_j t)$, for $\{1 \leq j \leq d\}$. In this multivariate setting, the function $\Upgamma : \field{R}^d \to \field{R}^d$ was assumed to be quadratic and time-varying. It is shown that the estimation error for tracking is not only a function of $\alpha$ but also the frequencies of $\qsaprobe$, since $\| y^*- y^\text{\tiny\sf opt} \|$ was of order $O(1/\omega^2_1 +d^2\alpha^2)$ where $\omega_1$ is the smallest frequency among $\{\omega_j\}$.

 
 
 
 

 
 

 

 
 
  

 
 \appendix
 
 
 

\section{Appendix}

 \subsection{ODE QSA Solidarity}
 
 
The first step in the proof 
of \Cref{t:Couple_theta_local}
is to  re-interpret the LLN \eqref{e:LLNf}:
\begin{lemma}
\label[lemma]{t:LLNfast}
If (A1) holds then the following uniform Law of Large Numbers holds:  there is a function $\varepsilon_f \colon\field{R}\to\field{R}_+$ such that for any $\theta\in\field{R}^d$ and $\Phi_0\in\Upomega$,
\begin{equation}  
\left |
\frac{1}{T }\int_0^T  \bigl[f(\theta,\qsaprobe^\alpha_{t}) - {\widebar{f}}(\theta) \bigr]  \, dt
				\right|     \le    \bigl(1 +  \| \theta\| \bigr)  \varepsilon_f(  T /\alpha   ) 
\label{e:LLNfast}
\end{equation}
with  $  \varepsilon_f( r  ) \to 0 $ as $r\to\infty$.  
\end{lemma}

Armed with this lemma, we now compare solutions to \eqref{e:QSAgen_alpha_time-scaled} with solutions to the mean flow \eqref{e:barf_ODE_gen}.   If both are initialized at the same value $\upvartheta_0 =\Uptheta_0 =\theta_0$, we obtain for any $T>0$,
\[
\Uptheta^\alpha_T = \theta_0+  \int_0^T   f(\Uptheta_r^\alpha ,\qsaprobe^\alpha_r) \, dr   \,,
		\qquad
\upvartheta_T= \theta_0+  \int_0^T {\widebar{f}}( \upvartheta_r) \, dr
\]
 The challenge is that the bound \eqref{e:LLNfast} holds with $\theta$ fixed, so that it is not immediately applicable to comparison of the two ODEs.
 
 
 
To proceed we divide the positive real line into adjacent intervals $[\tau_n, \tau_{n+1})$,   and let $\{ \upvartheta^n_t :    t \ge   \tau_{n} \}$ denote the solution to \eqref{e:barf_ODE_gen}  with $\upvartheta^n_{\tau_n} = \Uptheta^\alpha_{\tau_n}$.   
We take   $\tau_n = \Delta  n$ for each $n$, with $\Delta \in (0,1]$ to be chosen as a function of $\alpha\le 1$.  The following corollary to   \Cref{t:LLNfast} will provide constraints on this choice. 

\begin{lemma}
\label[lemma]{t:shortSolidarity}
The following uniform bounds hold under (A1):  for a fixed constant $B_f$, and each $n\ge 0$,
\[
\|\upvartheta^{n+1}_{\tau_{n+1}}  -\upvartheta^n_{\tau_{n+1}}    \|   
=
\|  \Uptheta^\alpha_{\tau_{n+1}}   -  \upvartheta^n_{\tau_{n+1}}  \|  \le     
				\Delta B_f^n(\Delta,\alpha)  \mathbin{:=}
			B_f \Delta^2    +  \Delta [1 + \|  \upvartheta^n_{\tau_{n}}  \|  ]  \varepsilon_f(\Delta/\alpha)
\]
\end{lemma}

The proof of the lemma relies on \Cref{t:LLNfast} and an application of the Bellman-Gr{\"o}nwall inequality:
\begin{lemma}
\label[lemma]{t:BG}
Let $\bm{\upvartheta}^i$, $i=1,2$, denote two solutions to \eqref{e:barf_ODE_gen}, with different initial conditions.   Then, subject to (A1),
\[
\|  \upvartheta^1_t  -  \upvartheta^2_t \|  \le \|  \upvartheta^1_0  -  \upvartheta^2_0 \|  e^{L_f t}\,,\qquad t\ge 0.
\]
\end{lemma}

\begin{proof}[Proof of \Cref{t:shortSolidarity}]
The equality $
\|\upvartheta^{n+1}_{\tau_{n+1}}  -\upvartheta^n_{\tau_{n+1}}    \|    = \|  \Uptheta^\alpha_{\tau_{n+1}}   -  \upvartheta^n_{\tau_{n+1}}  \|  $ follows from the definitions.   
To bound the latter, we write
\[
\begin{aligned}
\Uptheta^\alpha_{\tau_{n+1}}   -  \upvartheta^n_{\tau_{n+1}}  
	&=  \int_{\tau_n}^{\tau_{n+1}}   \bigl[  f(\Uptheta^\alpha_r, \qsaprobe^\alpha_r)  -   {\widebar{f}}(\upvartheta^n_r)\bigr]\, dr
			&&    \textbf{Bounds:}
   \\
	&=  \int_{\tau_n}^{\tau_{n+1}}  
		 \bigl[  f(\Uptheta^\alpha_{\tau_n}, \qsaprobe^\alpha_r)  -   {\widebar{f}}(\upvartheta^n_{\tau_n})\bigr]\, dr
		 &&     O\bigl( \bigl( 1+ \|\upvartheta^n_{\tau_n}\|  \bigr)\varepsilon_f(  \Delta /\alpha   )   \bigr)
	\\
	&\quad
		 +   \int_{\tau_n}^{\tau_{n+1}}   \bigl[ f(\Uptheta^\alpha_r, \qsaprobe^\alpha_r)   -   f(\Uptheta^\alpha_{\tau_n}  , \qsaprobe^\alpha_r) \bigr]\, dr	
		 \qquad
		 	 &&     O\bigl(  \Delta^2\bigr)
		 \\
	&\quad
		 +   \int_{\tau_n}^{\tau_{n+1}}   \bigl[      {\widebar{f}}(\upvartheta^n_{\tau_n})   -   {\widebar{f}}(\upvartheta^n_r)\bigr]\, dr
		 		  &&     O\bigl(  \Delta^2\bigr)
\end{aligned} 
\]
That is, the error is the sum of three terms, in which the first is bounded via \Cref{t:LLNfast}, and the second two are bounded by $O(\Delta^2)$ for $\Delta\le 1$ under the Lipschitz assumptions on $f$ and ${\widebar{f}}$. 
\end{proof}

The next step is to obtain uniform bounds on the error ${\cal E}_n =   \|\upvartheta^{n}_{\tau_n}  - \upvartheta_{\tau_n}    \|   $ for each $n$,   with $\upvartheta_{\tau_n} $ the solution at time $\tau_n$ with initial condition $\upvartheta_0 = \Uptheta_0 =\theta_0$, so in particular ${\cal E}_0=0$.   
\begin{lemma}
\label[lemma]{t:shortSolidarity2}
The following uniform bounds hold under (A1),  provided  $\varepsilon_f(  \Delta /\alpha   )\le 1$:      for a fixed constant $B_f^1$, and each $n\ge 0$,
\[
    \|\upvartheta^{n}_{\tau_n}  - \upvartheta_{\tau_n}    \|    \le B_f^1\exp(2 L_f \tau_n)   \bigl[  1     +      \|\theta_0\|  \bigr]  \cdot  \bigl[ \Delta +  \varepsilon_f(  \Delta /\alpha   )   \bigr]
\]
\end{lemma}

\begin{proof}
The proof begins with the triangle inequality to obtain a recursive bound:
\[
\begin{aligned}
{\cal E}_{n+1} \mathbin{:=}   \|\upvartheta^{n+1}_{\tau_{n+1}}  - \upvartheta_{\tau_{n+1}}    \|  
  &   \le      \|\upvartheta^{n}_{\tau_{n+1}}  - \upvartheta_{\tau_{n+1}}    \|  
+ \|\upvartheta^{n+1}_{\tau_{n+1}}  - \upvartheta^n_{\tau_{n+1}}    \|
\\
  &   \le      \|\upvartheta^{n}_{\tau_{n+1}}  - \upvartheta_{\tau_{n+1}}    \|    
  +
\Delta B_f^n(\Delta,\alpha) 
\end{aligned} 
\]
The recursion is obtained on combining this bound with \Cref{t:BG}:
\[
{\cal E}_{n+1}  \le     e^{L_f \Delta} \|\upvartheta^{n}_{\tau_n}  - \upvartheta_{\tau_n}    \|    +\Delta B_f^n(\Delta,\alpha) 
 =     e^{L_f \Delta} {\cal E}_n    +\Delta B_f^n(\Delta,\alpha) 
\]
where $B_f^n$ appears in \Cref{t:shortSolidarity}.   
Applying the triangle inequality,
\[
B_f^n(\Delta,\alpha)  = 
			B_f \Delta   +    [1 + \|  \upvartheta^n_{\tau_{n}}  \|  ]  \varepsilon_f(\Delta/\alpha)
			\le
			\bigl(B_f \Delta   +  \varepsilon_f(\Delta/\alpha)  \bigr)   
			+   \bigl(   \|  \upvartheta_{\tau_{n}}  \|   +  {\cal E}_n\bigr) \varepsilon_f(\Delta/\alpha)
\]
An application of \Cref{t:BG} once more gives $ \|  \upvartheta_{\tau_{n}}  \| \le  e^{L_f \tau_n}   \| \theta_0\|$.

Combining these bounds gives
\[
\begin{aligned}
{\cal E}_{n+1}  &\le     \beta_\Delta {\cal E}_n    +  \Delta   \bigl[   {\cal U}^0   +    {\cal U}^1  e^{L_f \tau_n}  \bigr]
   \\
 \  \textit{with} \ \              &\beta_\Delta\mathbin{:=}   e^{L_f \Delta} +  \varepsilon_f(\Delta/\alpha)
            \\
            &{\cal U}^0  \mathbin{:=} \bigl(B_f \Delta   +  \varepsilon_f(\Delta/\alpha)  \bigr)    
            \,, \quad {\cal U}^1 \mathbin{:=}   \varepsilon_f(\Delta/\alpha)    \| \theta_0\|
\end{aligned} 
\]
The final bound is obtained by iteration, and recalling ${\cal E}_0=0$.     

It is simplest to obtain a comparison with the ordinary differential equation,
\[
\ddt X_t  =  \gamma_\Delta X_t   +   {\cal U}^0   +  {\cal U}^1  \exp(L_f t)  \,,\quad X_0 = {\cal E}_0 = 0
\] 
with $\gamma_\Delta = (\beta_\Delta  -1)/\Delta$.
\end{proof}

 \begin{proof}[Proof of \Cref{t:Couple_theta_local}]
 Combining \Cref{t:shortSolidarity} and
 \Cref{t:shortSolidarity2}
 we obtain the bound \Cref{t:Couple_theta_local} for the sampling times $\tau_n$, provided we choose $\Delta = \Delta(\alpha)$ so that 
 \[
 \varepsilon_0(\alpha) \mathbin{:=}  \Delta(\alpha) +  \varepsilon_f(  \Delta(\alpha) /\alpha   )      \to 0\,,\quad \textit{as $\alpha{\hbox{\tiny down}}arrow 0$}
 \]
For this we choose $\Delta =\alpha^{\rho}$ with $\rho\in (0,1)$.    

 Lipschitz continuity of $f$ and ${\widebar{f}}$ implies that the bound \eqref{e:Couple_theta_local}
 is preserved for all time,  with an increase in $B_0$ by a value of only $O(\Delta)$.  
\end{proof}

 

\subsection{Baker's Theorem and Poisson's Equation}
\label{s:BakerPoisson}

 
Here we explain how to obtain well behaved solutions to Poisson's equation for a function $h\colon\field{R}^d\to\field{R}$ that is
is analytic.    We denote the mixed partials by
\[
	h^{(\beta)} (z) 

	=  \frac{\partial^{\beta_1}}{ \partial z_1^{\beta_1}} \cdots  \frac{\partial^{\beta_d}}{ \partial z_d^{\beta_d}}  h\,  (z)\,,\quad  \beta \in \field{Z}_+^d\, .
\] 
For $\beta \in \field{Z}_+^d$ and a vector $v\in\field{R}^d$ we denote $v^\beta = v_1^{\beta_1}\cdots  v_d^{\beta_d}$,   and $|\beta|$ the $\ell_1$ norm:  $|\beta| =\sum |\beta_i|$. 

 A  candidate representation for ${\hat h}$ can be found on obtaining the solution to Poisson's equation for each of the primitive functions appearing in a Taylor series expansion.    Justification is then a challenge without assumptions on the frequencies.   

For any $ \beta\in \field{Z}^K$ denote $\omega_\beta = \sum \beta_i \omega_i$ and $\phi_\beta = \sum \beta_i \phi_i$.  
We will see that the terms in the Taylor series expansion for ${\hat h}$ resemble those for $h$, but with terms divided by $\omega_\beta$ (subject to the constraint that this is non-zero).    Consequently, a lower bound on  $\omega_\beta$ as a function of $|\beta|$ is required---such bounds are available under Assumption~(A0).

 

The proof of the following in the almost periodic setting is based on a version of Baker's Thereom  \cite[Thm. 1.8]{bug18}:   there is a constant $C>0$ such that for any $\beta$, provided $\omega_{\beta} \neq 0$,
\[
|\omega_{\beta} |   \ge k_\beta^{-C} \,,  \quad k_\beta = \max\{3, \beta_1,\dots,\beta_d\}
\]
The periodic case is much simpler since in this case  $|\omega_\beta| \ge \omega_1$ whenever $\omega_\beta$ is non-zero.  

\begin{theorem}
\label[theorem]{t:hah}
Suppose that  $h \colon \field{R}^K \to \field{R}$ is analytic.
  Then there are coefficients $\{ c_\beta  :  \beta\in \field{Z}^K \}$ such that the following representations hold for $h$ and its mean ${\overline {h}}$:
\begin{subequations}
\begin{align}
 h(\qsaprobe_t) & =  \sum_{ \beta\in \field{Z}^K} c_\beta  \cos(2\pi[\omega_\beta t + \phi_\beta])   \,,
\qquad
{\overline {h}}  =  \sum_{ \beta:  \omega_\beta=0} c_\beta  \cos(2\pi  \phi_\beta)
 \label{e:tilh}
\end{align}
If in addition  the frequencies $\{ \omega_i \}$ satisfy either of the conditions \eqref{e:TwoFreqChoices}, then
there is a solution to Poisson's equation ${\hat h}$ that has zero mean, 
and is analytic on the domain  $\{\field{C}\setminus \{0\}\}^K$.    
When restricted to $\Upomega$ it takes the form,
\begin{equation}
{\hat h}(\Phi_t) =  -   \frac{1}{2 \pi  }   \sum_{ \beta:  \omega_\beta \neq  0}  \frac{1}{ \omega_\beta}  c_\beta   \sin(2\pi[\omega_\beta t + \phi_\beta]) 
\label{e:hahTheoremEquations}
\end{equation}

\end{subequations}
\end{theorem}

\begin{proof}
The almost periodic case \eqref{e:logFreqSpan} is taken from \cite[Appendix~B]{laumey22}.  
 Based on  (\ref{e:tilh}, \ref{e:hahTheoremEquations}), the extensions of $\tilde{h} = h-{\overline {h}}$ and ${\hat h}$ to $z\in \{\field{C}\setminus \{0\}\}^K$ are easily identified:
 \begin{align*}
 \tilde{h}(G(z)) & =   \sum_{ \beta:  \omega_\beta \neq  0}   c_\beta   \half [z^\beta  +  z^{-\beta} ]   
 \\
 {\hat h}(z) & =   -  j \frac{1}{2 \pi  }   \sum_{ \beta:  \omega_\beta \neq  0}   \frac{1}{ \omega_\beta} c_\beta   \half  [z^\beta  -  z^{-\beta} ]   
\end{align*}
 \end{proof}

 

 
The theorem implies the following useful corollary, which is taken from   \cite[Thm.~2.4]{laumey22}.   The geometry is illustrated in \Cref{fig:FishOrthogonality}.

\begin{corollary}
\label[corollary]{t:orthCor}
Suppose that $g$ and $h$ are two real-valued analytic functions on $\field{R}^K$.  
If  the frequencies $\{ \omega_i \}$ satisfy either form of \eqref{e:TwoFreqChoices}, 
then  ${\sf E}[ g(\qsaprobe)  {\hat h}(\Phi)] = 0$,   where ${\hat h}$ is given in \eqref{e:hahTheoremEquations}.
\end{corollary}

The corollary is vital for ensuring that $\widebar{\Upupsilon}\equiv {\sf Z}ero$ when frequencies are chosen with care.   

This subsection is concluded with the remaining step in the proof of  \Cref{t:orthCorBigPre}.
\begin{proposition}
\label[proposition]{t:orthCorBig} 
Consider the QSA ODE \eqref{e:QSAgen_alpha} subject to (A0a) and (A1), but with arbitrary choices of frequencies $\{\omega_i\}$.   
Suppose that an invariant probability measure $ \upvarpi$ exists, and let $\Uppsi = (\Uptheta, \Upphi)  $ denote a $\Uppi$-valued random vector with distribution $ \upvarpi$.   
Then,  for any analytic function $g \colon \Uppi\to \field{R}$ with zero-mean solution to Poisson's equation ${\hat g}$,
\[
\alpha  {\sf E} \bigl[  [D^f {\hat g} ] (\Uppsi )     \bigr]  =  {\sf E} \bigl[  g( \Uppsi )  -  {\widebar{g}}(\Uptheta)        \bigr] 
\]
In particular,  if $f$ is analytic and there is a solution ${\hat f}$ to Poisson's equation,  then the   tracking bias may be expressed
\begin{equation}
{\sf E} \bigl[ {\widebar{f}}(\Uptheta)        \bigr]  =
\alpha {\sf E} \bigl[  \Upupsilon (\Uppsi )     \bigr]  
\label{e:newUpupsilon}
\end{equation}
\end{proposition}
 

 \begin{proof} 
 For (i) we  apply \Cref{t:QSAgenerator}.   
 To establish \eqref{e:newUpupsilon},   recall from below \eqref{e:Dfg} the 
formula   $\Upupsilon  = - D^f g  $,   using $g={\hat f}$.   Hence from (i),
 \[
 \alpha {\sf E} \bigl[  \Upupsilon (\Uppsi )     \bigr]  =  - {\sf E} \bigl[  f( \Uppsi )  -  {\widebar{f}}(\Uptheta)        \bigr] 
\]
The proof is completed on establishing $ {\sf E} \bigl[  f( \Uppsi )         \bigr] ={\sf Z}ero$,  which follows from \Cref{t:QSAgenerator}. 
    \end{proof}
   
 

 
 

\subsection{Filtering footnotes}
\label{s:FiltFoot}

\wham{Attenuation bounds}

We begin with justification of the bounds \eqref{e:clWattenuation}.
Consider a stable filter with transfer function denoted $H$ and impulse response $h$.

Let $g$ be any differentiable function with Laplace transform $G$, $g^{(m)}$ be its $m$th derivative, and denote
\[
{\hat g}^{(m)}_t = \int_0^t  h_{t-\tau} g^{(m)}_\tau\, d\tau\,,
\]
omitting the $m$ in the special case $m=0$.  
Let $\| g\|_1$, $\| g\|_\infty$ denote the respective $L_1$ and $L_\infty$-norms of the function $g\colon\field{R}_+\to \field{C}$.

Bounds on ${\hat g}^{(i)}$ are obtained in the following,  which are uniform for vanishing $\alpha$ for the choice of filter \eqref{e:GVSOS}.

\begin{proposition}
	\label[proposition]{t:FilterDesign}
	Consider the stable second order transfer function 
	\[
	H(s) =   H^1(s) H^2(s) \,,\quad 
	H^i(s) = \frac{\gamma_i }{ s+\gamma_i}\, ,  \qquad i=1,2,       
	\] 
	with  $\text{Re}(\gamma_i)<0$ for each $i$.  Then,
	\begin{equation}
		\begin{aligned}
			| {\hat g}^{(1)}_t|  & \le | \gamma_1|  \bigl( \| h^2 \|_1  +\| h_2 \|_1 \bigr)   \|g\|_\infty  + o(1)  
			\\
			| {\hat g}^{(2)}_t|   &\le | \gamma_1 \gamma_2|  \bigl( 1  + \| h^1 \|_1  +\| h^2 \|_1  +\| h_2 \|_1 \bigr)   \|g\|_\infty  + o(1)  
		\end{aligned} 
		\label{e:HattenDerivatives}
	\end{equation}
	For the special case \eqref{e:GVSOS} with $0<\zeta\le 1$ fixed,  
	\begin{align}
		\| h^2 \|_1 = \| h^1 \|_1 & \le  \frac{1}{\zeta}   \,,
		\qquad
		\| h\|_1   \le  \frac{1}{\zeta^2}
		\label{e:GVSOSgains}
	\end{align} 
\end{proposition}

\begin{proof}
	For any $m$, the Laplace transform of $ {\hat g}^{(m)}$ is    
	\[
	{\cal L}( {\hat g}^{(m)}) = 
	s^m H(s) G(s)  - H(s) \text{poly}_m(s)
	\]
	In particular, $ \text{poly}_1(s) =  g_0$ and  $ \text{poly}_2(s) = s g_0 + g'_0$.
	
	
	Consider a stable first order filter with transfer function $H^1(s) = \beta / (s+\beta)$, where $\text{Re}(\beta)>0$.
	If we pass the derivative of the function $g$ through this filter, then the Laplace transform of the output is  $s H^1(s) G(s)$, with   
	\[
	s H^1(s) =  \beta  [1-   H^1(s)  ]

	\]
	From this we conclude that we obtain attenuation by $\beta$:   
	$
	{\hat g}^{(1)}_t   =  \beta  (g_t - {\hat g}_t )  + o(1)  
	$,
	where the $o(1)$ term is a constant times $e^{-\gamma t}$.  Letting $\|g\|_\infty = \sup_t |g_t|$  gives
	\[
	| {\hat g}^{(1)}_t|   \le   |\beta|   (1+ \| h^1\|_1 ) \|g\|_\infty  + o(1)  
	\]  
	We have $\| h^1\|_1 = 1$ if $\beta$ is real.
	
	This easily leads to the desired bounds for $H = H^1 H^2$ using $\beta=\gamma_i$.    
	We have,
	\[
	\begin{aligned}
		s H(s) &=  {\gamma_1 } \bigl[  1-  H^1(s)   \bigr]    H^2(s)    =  {\gamma_1 } \bigl[    H^2(s)-  H(s)   \bigr]  
		\\
		s^2 H(s) &=  {\gamma_1\gamma_2} \bigl[  1-  H^1(s)   \bigr] \bigl[  1-   H^2(s)   \bigr] 
		={\gamma_1\gamma} \bigl[  1-  H^1(s)    -   H^2(s)  + H(s)  \bigr]  
	\end{aligned} 
	\] 
	which results in the bounds \eqref{e:HattenDerivatives}.
	
	
	We turn next to the bounds in \eqref{e:GVSOSgains} for the filter \eqref{e:GVSOS}. 
	The $L_1$ norm is the induced operator norm, with $h$ viewed as a linear operator on $L_\infty$. This implies it is sub-multiplicative.  In the notation above,
	\[
	\begin{aligned} 
		\| h \|_1 &=  \| h^1 *h^2 \|_1  \le \| h^1 \|_1 \| h^2 \|_1
		\\
		\textit{and}\quad
		\| h^1 \|_1 & \le {|\gamma_1|}  \int_0^\infty \exp(- \text{Re}\, (\gamma_1) t )\, dt  
		=   {|\gamma_1|}  \frac{1}{\text{Re}\, (\gamma_1)}  
	\end{aligned} 
	\]
	For these complex poles, $|\gamma_1| =\gamma$   and $ \text{Re}\, (\gamma_1)  =\zeta \gamma$, 
	which gives the desired bounds:
	\[
	\| h^2 \|_1=\| h^1 \|_1  \le  \frac{1}{\zeta}
	\qquad
	\| h \|_1   \le \| h^1 \|_1 \| h^2 \|_1  \le   \gamma^2 \frac{1}{\zeta^2 \gamma^2}  = \frac{1}{\zeta^2}
	\]
\end{proof}

\wham{Implications to design}

\Cref{t:Couple_theta}  is immediate from \Cref{t:PMFlin} and the elementary bounds \eqref{e:clWattenuation}.
However, these bounds are conservative, as we now show through a frequency domain analysis.

In view of \Cref{t:PMFlin},  filter design can be performed based on the linear system
\[
\ddt x_t =   \alpha A^*   x_t   		 +  \alpha {\cal W}_t  
\]
in which $x_t$ is an approximation of  $\Uptheta_t -\theta^\text{\rm s.t.\,}arp$ up to order $O(\alpha^2)$  (ignoring transients).  
The definition of ${\cal W}_t$ remains the same---a function of $(\Uptheta_t,\Phi_t)$.

Let $X(s)$, $W(s)$ denote the respective Laplace transforms of the state and input in this linear system,  and $W^i(s)$ the transforms of the components of ${\cal W}_t$ shown in \Cref{t:PMF}.
Taking Laplace transforms of each side gives,
\[
X(s)  =   \alpha  [Is -  \alpha A^* ] ^{-1}  W(s)   =   \alpha  [Is -  \alpha A^* ]^{-1}  \bigl(  \alpha^2  W^0(s)   + \alpha s  W^1(s)   + s^2  W^2(s)   \bigr) 
\]
Using a superscript ``{\sf F}'' for the filtered signals, we obtain
\[
X^{\text{\tiny\sf F}}(s)  =   \alpha  [Is -  \alpha A^* ] ^{-1}  \WF{}(s)   =   \alpha  [Is -  \alpha A^* ]^{-1}  \bigl(  \alpha^2   \WF{0}(s)   + \alpha s   \WF{1}(s)   + s^2   \WF{2}(s)   \bigr) 
\]
The  filter $H$ is designed so that $ s^2  \WF{2} (s) $ and $   s  \WF{1} (s)  $  are attenuated,  as explained in \Cref{t:FilterDesign}.      Hence a proof of \Cref{t:Couple_theta} can be conducted entirely in the frequency domain,  based on the fact that the induced operator norm of $ \alpha [Is - \alpha A^*]^{-1} $, viewed as a mapping on $L_\infty$, is uniformly bounded over $0<\alpha\le 1$.

The matrix-valued transfer function $[Is - \alpha A^*]^{-1}  $ may attenuate some of these signals if $\alpha$ is small, and the signals are band-limited.     An example is the experiment surveyed right after \Cref{t:Couple_theta} for which results are presented in \Cref{fig:alpha2bias}:  a second order filter is not needed in this special case, because $\alpha A^*= - \alpha$,  and the spectrum of ${\cal W}_t$ is discrete.  
 
\subsection{Stability}
\label{s:AppStability}

\begin{proposition}
\label[lemma]{t:CriteriaSolidarity}  
Either of the criteria (V4) or asymptotic stability of the ODE@$\infty$ imply (V4').  
\end{proposition}

\begin{proof}
The implication (V4) $\Rightarrow$ (V4') follows from integration:   for each $T$,  subject to \eqref{e:ddt_bound_Lyapfun},
\[
V(\upvartheta_T) - 
V(\upvartheta_0) \le    -  \delta_0 \int_0^T  V(\upvartheta_{\tau}) \, d\tau
\]
This combined with the bounds in \eqref{e:VLip} imply that (V4') also holds,  with $T>0$ sufficiently large.

 

We next establish (V4') subject to the assumption that the ODE@$\infty$ is asymptotically stable.   For this we  establish that the function $V$ defined in \eqref{e:BMV} is a Lipschitz continuous solution to \eqref{e:V4T}. 

 We have already remarked that Lipschitz continuity follows from (QSA2), which requires  Lipschitz continuity of ${\widebar{f}}$.

For the remainder of the proof, it is enough to establish that this function solves (V4') for the ODE@$\infty$, by  applying the following bound from \cite{bormey00a}:
There is a function $\varepsilon_0\colon\field{R}_+\to\field{R}_+$ for which 
 $ \varepsilon_0(r)\to 0$ as $r\to\infty$, and 
the following bound holds:
with $\theta = \upvartheta^\infty_0 = \upvartheta_0$,  and 
\[
\| \upvartheta^\infty_T  - \upvartheta_T  \|    =  \varepsilon_0(\| \theta\| )   \| \theta\|  
\]
This is a simple consequence of the representation \eqref{e:ODEinf-r}.

To obtain the inequality in \eqref{e:BMV} for $\bm{\upvartheta}^\infty$ requires that we take $T$ sufficiently large:
Choose $T>0$ so that  $\|\upvartheta_t^\infty \|  \le \half  \|\upvartheta_t^\infty \|$ for $t\ge T$  (existence follows from Lemma 4.23 of \cite{CSRL}).    
We have for any  initial condition $\upvartheta_0$, 
\[
	V(\upvartheta^\infty_T)   =
	\int_0^{T}     \| \upvartheta_{t+T} ^\infty \| \,dt
		\le \half 	V(\upvartheta^\infty_0) 
\]
It follows that this function $V$ solves (V4') for the ODE@$\infty$.    
\end{proof}

 \begin{lemma}
\label[lemma]{t:littleV4bdd}
Suppose that  (A1) and the Lyapunov bound   (V4') hold.   Then, the following semi-exponential stability bound holds for the mean flow:
for positive constants $\delta_1$, $b$ and $c$,  
\begin{equation}
\| \tilde\odestate_t \|   \le  b \| \tilde\odestate_0 \|   e^{-\delta_1 t}  \,,
	\quad
	 \textit{for $ t \le \tau_c  $,  where $\tau_c= \min\{ t : \| \tilde\odestate_t \|    \le c \}$  
	 and $\tilde\odestate_t=\upvartheta_t -\theta^\text{\rm s.t.\,}arp$.  }
\label{e:littleV4bdd}
\end{equation}
\end{lemma}

\begin{proof}

We can choose $\tau =0 $ in \eqref{e:V4T} without loss of generality since the system is autonomous. For $\delta_1,\beta>0$, \eqref{e:V4T} gives
\[
V(\upvartheta_{T}) \leq  V(\upvartheta_{0})  - \delta_1 \|\upvartheta_0\| \leq e^{ - \beta} V(\upvartheta_0) \quad \text{if } \|\upvartheta_\tau\| > c
\] 
The condition above can be extended to general $t$ by first defining $t \mathbin{:=} nT $. 
Then,
\[
V(\upvartheta_{t}) \leq e^{ - \beta_1 nT} V(\upvartheta_0), \quad \beta_1 = \beta/T  \quad \text{if } t = nT \leq \tau_c
\]
This establishes \eqref{e:littleV4bdd} at these discrete time points.   It extends to arbitrary $t\le \tau_c$ via Lipschitz continuity of $V$.   
\end{proof}

\begin{proof}[Proof of \Cref{t:Couple_theta}]
 
The first step in the proof is to recognize that $\alpha^0$-ultimate boundedness allows us to restrict analysis to $T\ge T_\ell(\theta)$, so that $\Uptheta_t$ is restricted to the region $\{ \|\Uptheta_t \| \le B \}$  [independent of $\theta=\upvartheta_0$ and $\alpha\in (0,\alpha^0)$].   In the proof we assume without loss of generality that $\|\Uptheta_0 \| \le B$.

Recall that under exponential asymptotic stability there is a Lyapunov function  $V\colon\field{R}^d \to\field{R}_+$ with Lipschitz gradient satisfying, for positive constants $\delta,\delta_0$,
 \[
 	\nabla V(x) \cdot {\widebar{f}}(x) \le -\delta_0 \|x - \theta^*\|^2  \le  -\delta  V(x) \,, \qquad x\in\field{R}^d \,.
 \]
 One solution is found in \eqref{e:notBMV}.
Based on this we compute the derivative of $V(Y_t)$ with respect to time, based on \eqref{e:prematureMeanflow}:
\[
\ddt V(Y_t)  =  \nabla V\, (Y_t) \cdot \ddt Y_t  =   \nabla V\, (Y_t) \cdot \bigl[ 
  \alpha\bigl[   {\widebar{f}}( Y_t    )   
			- \alpha  \bigl(B_t  {\hat f}_t  +  \Upupsilon_t  \bigr)  \bigr] \bigr] \] 
From the definitions this gives
\[
\ddt V(Y_t)   \le
			- \alpha  \delta  V( Y_t ) + O(\alpha^2) \| Y_t\| \bigl( \|  {\hat f}_t  \| + \|  {\widehat U}pupsilon_t  \| \bigr)
			 \le- \alpha  \delta  V( Y_t )   + O(\alpha^2) \bigl( 1 + V(Y_t)\bigr)
\] 
An application of the Bellman-Gr{\"o}nwall implies that both $\| Y_t  - \theta^*\|$  and $\| \Uptheta_t  - \theta^*\|$ are $O(\alpha)$ for $t\ge (\delta\alpha)^{-1}$, which establishes \eqref{e:fixedGain1}.

   
The proof of \eqref{e:fixedGain1PR} begins with  \eqref{e:fixedGain1}, which justifies restricting analysis to $t$ for which $\| \Uptheta_t  - \theta^*\|\le O(\alpha)$, so that a linearization around $A^*$ is justified.   The remainder of the proof  is immediate from \Cref{t:PMFlin}  and the following bounds obtained in \Cref{t:FilterDesign}:  Under the assumptions of the theorem,
\begin{equation}
| {\cal W}_t^{i \text{\tiny\sf F}} |   = O(\alpha^i) \| {\cal W}^i\|_\infty + o(1)   =   O(\alpha^3) + o(1)\,,\qquad t\ge 0\,,  \ i=1,2
\label{e:clWattenuation}
\end{equation}
 where ${\cal W}^{i \text{\tiny\sf F}}$ is the output of the filter when the input is $  {\cal W}^i$,
 and $\| {\cal W}^i\|_\infty $ denotes the $L_\infty $ norm.   We cannot obtain this bound for  ${\cal W}_t^{2 \text{\tiny\sf F}} $ using a first order filter unless we have information regarding the bandwidth of ${\cal W}_t^2$ [which may depend on $\alpha$].

\end{proof}
 

\subsection{Implications to ESC}
\label{s:ESCtheory}

 
\subsubsection{Stability}

Recall the following bound for a function $\Upgamma(\theta):\field{R}^d \to \field{R}$ with Lipschitz gradient:
 \begin{equation}
	\Upgamma(\theta) \leq \Upgamma(\theta^\prime) + [\theta - \theta^\prime]^{\rm tr\, }anspose \nabla \Upgamma(\theta^\prime) + \tfrac{1}{2}L \|\theta - \theta^\prime \|^2,\quad \theta^\prime, \theta \in \field{R}^d
	\label{e:Lip_Grad}
\end{equation}
where $L$ is the Lipschitz constant associated with $\nabla \Upgamma$.
This bound and a bit more allows us to construct a solution to (V4) for the qSGD ODEs:

  \begin{proposition}
\label[proposition]{t:Lip_implies_V4}
Suppose that $\Upgamma$ has a Lipschitz gradient, so that \eqref{e:Lip_Grad} holds.  Suppose moreover that   $\|  \nabla\Upgamma (\theta)\| \ge \delta\|\theta\| $  for some $\delta>0$ and all $\theta$ satisfying $\|\theta\|\ge \delta^{-1}$.   Suppose that the probing signal is chosen of the form \eqref{e:qSGD_probe} with $K=d$ and $\Sigma_{\text{\tiny$\qsaprobe$}}>0$. 
 Then,  $V(\theta) = \hbox{\rlap{$\sqcap$}$\sqcup$}rt{\Upgamma(\theta) - \Upgamma^\text{\tiny\sf opt} }$ satisfies assumption (V4) for both  1qSGD and 2qSGD. 
  \end{proposition}

\subsubsection{Vector field approximations}

To obtain bounds that are uniform in the scaling gain $\varepsilon>0$ it is useful to apply a Taylor series representation very different from what was used in \Cref{t:hah}.

 \begin{lemma}
\label[lemma]{t:ObjObsTaylor}
If the objective function $\Upgamma$ is analytic, then it and the normalized observation function ${\cal Y}n$ admit the representations,
\begin{equation}
\begin{aligned} 
\Upgamma(\theta + \upepsilon \qsaprobe) & =   \Upgamma(\theta) +  \sum_{| \beta|\ge 1}   c_\beta  	\Upgamma^{(\beta)} (\theta)  \upepsilon^{|\beta|}  \qsaprobe^{\beta}
\\
    {\cal Y}n( \theta,\qsaprobe)   & =  \frac{1}{\upepsilon}   \Upgamma(\theta) +  \qsaprobe^{\rm tr\, }anspose \nabla\, \Upgamma(\theta)  
    				+  \sum_{| \beta|\ge 2}   c_\beta  	\Upgamma^{(\beta)} (\theta)  \upepsilon^{|\beta|-1}  \qsaprobe^{\beta}
\end{aligned} 
\label{e:ObjObsTaylor}
\end{equation}
where in each case the sum   is restricted to $\beta\in\field{Z}_+^d$.
\end{lemma}

 

For general QSA theory there is no loss of generality in taking $G$ analytic in the definition  \eqref{e:qsaprobeG}.     
In applications to ESC it is convenient to take $G$ linear, since in this case we obtain better control of the terms in  \eqref{e:ObjObsTaylor},  and hence better approximations for ${\hat f}$ and its derivatives.   Before summarizing these results we require some notation.

Let ${\widehat G}$ denote the solution to Poisson's equation with forcing function $\qsaprobe_t = G(\Phi_t)$,   and ${\widehat{\Sigma}}$ the solution to Poisson's equation with forcing function $\qsaprobe_t \qsaprobe_t ^{\rm tr\, }anspose = G(\Phi_t)G(\Phi_t)^{\rm tr\, }anspose$.  
In the special case  \eqref{e:qSGD_probe}  we write    $\qsaprobe_t = V \qsaprobe_t^0$ using
\begin{equation}
	\qsaprobe _t^0  =  [  \cos (2\pi [\,  \omega_1 t  +  \phi_1 ] ) ;\,  \dots  ;\,    \cos (2\pi [\,  \omega_K t  +  \phi_K ] )  ]
\label{e:qSGD_probe0}
\end{equation}
 with $V$ the $m\times K$ matrix with columns equal to the $v^i$ appearing in  \eqref{e:qSGD_probe}.

\begin{subequations}

\begin{lemma}
\label[lemma]{t:FirstAndSecondOrderPoisson}
Zero mean solutions to Poisson's equation
 with respective forcing functions  $G(z)$,  $G(z)G(z)^{\rm tr\, }anspose$, 
may be expressed in the time domain using 
$\widehat qsaprobe_t = {\widehat G}  (\Phi_t) $ and   ${\widehat{\Sigma}}_t ={\widehat{\Sigma}}  (\Phi_t) $  with
\begin{align}
\!
\!
\!
\widehat qsaprobe_t  &=   V    \widehat qsaprobe^0_t   \,,
\qquad 
&&
  \widehat qsaprobe^{0i}_t  =  
   \frac{1}{2\pi \omega_i}  \sin (2\pi [\,  \omega_i t  +  \phi_i  ] ) \,,\qquad 1\le i\le d
\label{e:hapsaprobe}
\\
{\widehat{\Sigma}}_t  &=  V {\widehat{\Sigma}}^0_t  V^{\rm tr\, }anspose  \,,
&&
{\widehat{\Sigma}}^{0i,j}_t  =   \begin{cases}
				\phantom{+}\frac{1}{4\pi} \frac{1}{\omega_i +\omega_j}    \sin (2\pi [\,  (\omega_i + \omega_j)t  +  \phi_i +\phi_j   ] )  &  
\\  
				+\frac{1}{4\pi} \frac{1}{\omega_i -\omega_j}    \sin (2\pi [\,  (\omega_i - \omega_j)t  +  \phi_i -\phi_j    ] )  &   i\neq j
				\\[.5em]
				\phantom{+}
				\frac{1}{8\pi} \frac{1}{\omega_i }    \sin (4\pi [\,  \omega_i t  +  \phi_i   ] )    & i=j
				\end{cases}
\label{e:haSigmaqsaprobe}
\end{align}
\end{lemma}
\end{subequations}

\begin{proof}
The representation \eqref{e:hapsaprobe} is immediate from  $\qsaprobe_t = V \qsaprobe_t^0$ and the definition \eqref{e:PoissonDefn}.

The trigonometric identity $ \cos(2x  )  =  2  \cos( x)  ^2 -1$ and $\Sigma_{\text{\tiny$\qsaprobe$}} =  \half VV^{\rm tr\, }anspose $ gives   
\[
\begin{aligned} 
\qsaprobe_t \qsaprobe_t ^{\rm tr\, }anspose  	  -   \Sigma_{\text{\tiny$\qsaprobe$}}
	& = V   \hbox{\rm diag\thinspace}  [  \cos^2 (2\pi [\,  \omega_1 t  +  \phi_1 ] ) \,, \dots \,,   \cos^2 (2\pi [\,  \omega_K t  +  \phi_K ] )  ]   V^{\rm tr\, }anspose 	  -   \Sigma_{\text{\tiny$\qsaprobe$}}
	\\
	& = \half V  \hbox{\rm diag\thinspace} [  \cos (2\pi [\,  2\omega_1 t  +  2\phi_1 ] ) \,, \dots \,,   \cos (2\pi [\,  2\omega_K t  +  2\phi_K ] )  ]   V^{\rm tr\, }anspose 
\end{aligned} 
\]
The representation \eqref{e:haSigmaqsaprobe} follows from integration of each side, and applying the definition \eqref{e:PoissonDefn}.
\end{proof}

\Cref{t:ObjObsTaylor} is presented with $\upepsilon $ an independent variable.  It remains valid and useful when we substitute $\upepsilon = \upepsilon(\theta)$:

\begin{subequations} 
\begin{proposition}
\label[proposition]{t:ESCapproximation}   
If $\Upgamma$ is analytic,
	and  the frequencies $\{ \omega_i \}$ satisfy either version of \eqref{e:TwoFreqChoices},
 then we have the following approximations for 1qSGD when $\upepsilon$ is   a function of $\theta$:
\begin{align}
f( \theta,\qsaprobe)   & = - \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta) \qsaprobe
			- \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla\, \Upgamma(\theta)   +  O(\varepsilon)   
&&   
{\widebar{f}}( \theta)   = - \Sigma_{\text{\tiny$\qsaprobe$}}\nabla\, \Upgamma(\theta)   +  O(\varepsilon^2)   
\label{e:ESCapprox_f}
   \\
A( \theta,\qsaprobe)   & =  
 -  \partial_\theta{\sf B}igl( \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta) {\sf B}igr) \qsaprobe
 - \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla^2\, \Upgamma(\theta)   +  O(\varepsilon)  
&&
{\bar{A}}( \theta)    =   - \Sigma_{\text{\tiny$\qsaprobe$}} \nabla^2\, \Upgamma(\theta)   +  O(\varepsilon^2)   
  \label{e:ESCapprox_A} 
  \\
 {\hat f}( \theta, \Phi)   & = - \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta) {\widehat G}(\Phi)
			- {\widehat{\Sigma}}(\Phi) \nabla\, \Upgamma(\theta)   +  O(\varepsilon)   
&&   
\label{e:ESCapprox_haf} 
\end{align} 
where in each case the ratio $O(\varepsilon^p) /\varepsilon^p $ is uniformly bounded over $\varepsilon>0$, when $\theta$ is restricted to any compact subset of $\field{R}^d$.    
\end{proposition}
\end{subequations} 

\begin{proof}
The expressions for $f$ and $A$ are immediate from 
\Cref{t:ObjObsTaylor},  and the expression for ${\hat f}$ then follows from \Cref{t:hah} combined with  \Cref{t:FirstAndSecondOrderPoisson}.
The error is $O(\varepsilon^2)$ for the means ${\widebar{f}}$ and ${\bar{A}}$ because ${\sf E}[ \qsaprobe^\beta] = 0$ when $|\beta| =3$.
\end{proof}

 

Expressions for the other terms appearing in \Cref{t:PMF} immediately follow.  We have
 \[
  {\hat h}af( \theta,\Phi)    = - \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta) {\hat h}aqsaprobe(\Phi)
			- {\hat h}aSigma(\Phi) \nabla\, \Upgamma(\theta)   +  O(\varepsilon)   
\]
in which	the ``double hats'' are defined in analogy with ${\hat h}af$, and
\[
\Upupsilon ( \theta,\Phi) =  - [D^f {\hat f}]( \theta,\Phi)   = -  
 {\sf B}igl[  \partial_\theta{\sf B}igl( \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta) {\sf B}igr)   {\widehat G}(\Phi)
 		+
			{\widehat{\Sigma}}(\Phi) \nabla^2\, \Upgamma(\theta)  {\sf B}igr] {\sf B}igl[ \frac{1}{\upepsilon(\theta)}   \Upgamma(\theta) \qsaprobe
			+
			 \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla\, \Upgamma(\theta){\sf B}igr]
											 +  O(\varepsilon)  
\]
While \Cref{t:orthCor} tells us that $\widebar{\Upupsilon}    ( \theta)  \mathbin{:=}{\sf E}[\Upupsilon ( \theta,\Phi)  ]= 0$ for every $\theta$ under the assumptions of    \Cref{t:bigQSA+ESC},   the expression above  combined with
 \eqref{e:AllTheNoise0}   
and
 \eqref{e:AllTheNoise1}   
 suggest high volatility.  

\subsubsection{Conclusions for 2qSGD}  

The conclusions above simplify greatly in this case:  
Letting $f_1$ denote the QSA vector field for 1qSGD,  and $f_2$   the QSA vector field for 2qSGD,  we have
\[
f_2(\theta,\qsaprobe)  = \half[ f_1(\theta,\qsaprobe) + f_1(\theta, -\qsaprobe) ]
\]
The following corollary to  \Cref{t:ESCapproximation} is immediate:  

\begin{subequations} 
\begin{proposition}
\label[proposition]{t:ESCapproximation2}   
If $\Upgamma$ is analytic and  the frequencies $\{ \omega_i \}$ satisfy either version of \eqref{e:TwoFreqChoices}, 
we then have the following approximations for 2qSGD:  
\begin{align}
f( \theta,\qsaprobe)   & =  
			- \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla\, \Upgamma(\theta)   +  O(\varepsilon)   
&&   
{\widebar{f}}( \theta)   = - \Sigma_{\text{\tiny$\qsaprobe$}}\nabla\, \Upgamma(\theta)   +  O(\varepsilon^2)   
\label{e:ESCapprox_f2}
   \\
A( \theta,\qsaprobe)   & =   
 - \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla^2\, \Upgamma(\theta)   +  O(\varepsilon)  
&&
{\bar{A}}( \theta)    =   - \Sigma_{\text{\tiny$\qsaprobe$}} \nabla^2\, \Upgamma(\theta)   +  O(\varepsilon^2)   
  \label{e:ESCapprox_A2} 
  \\
 {\hat f}( \theta, \Phi)   & = 
			- {\widehat{\Sigma}}(\Phi) \nabla\, \Upgamma(\theta)   +  O(\varepsilon)   
&&     {\hat h}af( \theta,\Phi)    =  
			- {\hat h}aSigma(\Phi) \nabla\, \Upgamma(\theta)   +  O(\varepsilon)   
\\
\Upupsilon ( \theta,\Phi) &=   -   
			{\widehat{\Sigma}}(\Phi) \nabla^2\, \Upgamma(\theta) \qsaprobe \qsaprobe^{\rm tr\, }anspose \nabla\, \Upgamma(\theta) 
											 +  O(\varepsilon)   &&
\label{e:ESCapprox_haf2} 
\end{align} 
where in each case the ratio $O(\varepsilon^p) /\varepsilon^p $ is uniformly bounded over $\varepsilon>0$, when $\theta$ is restricted to any compact subset of $\field{R}^d$.    
\end{proposition}
\end{subequations} 

\end{document}