Stochastic approximation

1. If ${\textstyle f(\theta )}$ is twice continuously differentiable, and strongly convex, and the minimizer of ${\textstyle f(\theta )}$ belongs to the interior of ${\textstyle \Theta }$, then the Robbins–Monro algorithm will achieve the asymptotically optimal convergence rate, with respect to the objective function, being ${\textstyle \operatorname {E} [f(\theta _{n})-f^{*}]=O(1/n)}$, where ${\textstyle f^{*}}$ is the minimal value of ${\textstyle f(\theta )}$ over ${\textstyle \theta \in \Theta }$.[5][6]
2. Conversely, in the general convex case, where we lack both the assumption of smoothness and strong convexity, Nemirovski and Yudin[7] have shown that the asymptotically optimal convergence rate, with respect to the objective function values, is ${\textstyle O(1/{\sqrt {n}})}$. They have also proven that this rate cannot be improved.

While the Robbins–Monro algorithm is theoretically able to achieve ${\textstyle O(1/n)}$ under the assumption of twice continuous differentiability and strong convexity, it can perform quite poorly upon implementation. This is primarily due to the fact that the algorithm is very sensitive to the choice of the step size sequence, and the supposed asymptotically optimal step size policy can be quite harmful in the beginning.[6][8]

Chung (1954)[9] and Fabian (1968)[10] showed that we would achieve optimal convergence rate ${\textstyle O(1/{\sqrt {n}})}$ with ${\textstyle a_{n}=\bigtriangledown ^{2}f(\theta ^{*})^{-1}/n}$ (or ${\textstyle a_{n}={\frac {1}{(nM'(\theta ^{*}))}}}$). Lai and Robbins[11][12] designed adaptive procedures to estimate ${\textstyle M'(\theta ^{*})}$ such that ${\textstyle \theta _{n}}$ has minimal asymptotic variance. However the application of such optimal methods requires much a priori information which is hard to obtain in most situations. To overcome this shortfall, Polyak (1991)[13] and Ruppert (1988)[14] independently developed a new optimal algorithm based on the idea of averaging the trajectories. Polyak and Juditsky[15] also presented a method of accelerating Robbins–Monro for linear and non-linear root-searching problems through the use of longer steps, and averaging of the iterates. The algorithm would have the following structure:

$${\displaystyle \theta _{n+1}-\theta _{n}=a_{n}(\alpha -N(\theta _{n})),\qquad {\bar {\theta }}_{n}={\frac {1}{n}}\sum _{i=0}^{n-1}\theta _{i}}$$

The convergence of ${\displaystyle {\bar {\theta }}_{n}}$ to the unique root ${\displaystyle \theta ^{*}}$ relies on the condition that the step sequence ${\displaystyle \{a_{n}\}}$ decreases sufficiently slowly. That is A1)

$${\displaystyle a_{n}\rightarrow 0,\qquad {\frac {a_{n}-a_{n+1}}{a_{n}}}=o(a_{n})}$$

Therefore, the sequence ${\textstyle a_{n}=n^{-\alpha }}$ with ${\textstyle 0<\alpha <1}$ satisfies this restriction, but ${\textstyle \alpha =1}$ does not, hence the longer steps. Under the assumptions outlined in the Robbins–Monro algorithm, the resulting modification will result in the same asymptotically optimal convergence rate ${\textstyle O(1/{\sqrt {n}})}$ yet with a more robust step size policy.[15] Prior to this, the idea of using longer steps and averaging the iterates had already been proposed by Nemirovski and Yudin[16] for the cases of solving the stochastic optimization problem with continuous convex objectives and for convex-concave saddle point problems. These algorithms were observed to attain the nonasymptotic rate ${\textstyle O(1/{\sqrt {n}})}$.

A more general result is given in Chapter 11 of Kushner and Yin[17] by defining interpolated time ${\textstyle t_{n}=\sum _{i=0}^{n-1}a_{i}}$, interpolated process ${\textstyle \theta ^{n}(\cdot )}$ and interpolated normalized process ${\textstyle U^{n}(\cdot )}$ as

$${\displaystyle \theta ^{n}(t)=\theta _{n+i},\quad U^{n}(t)=(\theta _{n+i}-\theta ^{*})/{\sqrt {a_{n+i}}}\quad {\mbox{for}}\quad t\in [t_{n+i}-t_{n},t_{n+i+1}-t_{n}),i\geq 0}$$

Let the iterate average be ${\displaystyle \Theta _{n}={\frac {a_{n}}{t}}\sum _{i=n}^{n+t/a_{n}-1}\theta _{i}}$ and the associate normalized error to be ${\displaystyle {\hat {U}}^{n}(t)={\frac {\sqrt {a_{n}}}{t}}\sum _{i=n}^{n+t/a_{n}-1}(\theta _{i}-\theta ^{*})}$.

With assumption A1) and the following A2)

A2) There is a Hurwitz matrix ${\textstyle A}$ and a symmetric and positive-definite matrix ${\textstyle \Sigma }$ such that ${\textstyle \{U^{n}(\cdot )\}}$ converges weakly to ${\textstyle U(\cdot )}$, where ${\textstyle U(\cdot )}$ is the statisolution to

$${\displaystyle dU=AU\,dt+\Sigma ^{1/2}\,dw}$$

where ${\textstyle w(\cdot )}$ is a standard Wiener process.

satisfied, and define ${\textstyle {\bar {V}}=(A^{-1})'\Sigma (A')^{-1}}$. Then for each ${\textstyle t}$,

$${\displaystyle {\hat {U}}^{n}(t){\stackrel {\mathcal {D}}{\longrightarrow }}{\mathcal {N}}(0,V_{t}),\quad {\text{where}}\quad V_{t}={\bar {V}}/t+O(1/t^{2}).}$$

The success of the averaging idea is because of the time scale separation of the original sequence ${\textstyle \{\theta _{n}\}}$ and the averaged sequence ${\textstyle \{\Theta _{n}\}}$, with the time scale of the former one being faster.

Suppose we want to solve the following stochastic optimization problem

$${\displaystyle g(\theta ^{*})=\min _{\theta \in \Theta }\operatorname {E} [Q(\theta ,X)],}$$

where ${\textstyle g(\theta )=\operatorname {E} [Q(\theta ,X)]}$ is differentiable and convex, then this problem is equivalent to find the root ${\displaystyle \theta ^{*}}$ of ${\displaystyle \nabla g(\theta )=0}$. Here ${\displaystyle Q(\theta ,X)}$ can be interpreted as some "observed" cost as a function of the chosen ${\displaystyle \theta }$ and random effects ${\displaystyle X}$. In practice, it might be hard to get an analytical form of ${\displaystyle \nabla g(\theta )}$, Robbins–Monro method manages to generate a sequence ${\displaystyle (\theta _{n})_{n\geq 0}}$ to approximate ${\displaystyle \theta ^{*}}$ if one can generate ${\displaystyle (X_{n})_{n\geq 0}}$ , in which the conditional expectation of ${\displaystyle X_{n}}$ given ${\displaystyle \theta _{n}}$ is exactly ${\displaystyle \nabla g(\theta _{n})}$, i.e. ${\displaystyle X_{n}}$ is simulated from a conditional distribution defined by

$${\displaystyle \operatorname {E} [H(\theta ,X)|\theta =\theta _{n}]=\nabla g(\theta _{n}).}$$

Here ${\displaystyle H(\theta ,X)}$ is an unbiased estimator of ${\displaystyle \nabla g(\theta )}$. If ${\displaystyle X}$ depends on ${\displaystyle \theta }$, there is in general no natural way of generating a random outcome ${\displaystyle H(\theta ,X)}$ that is an unbiased estimator of the gradient. In some special cases when either IPA or likelihood ratio methods are applicable, then one is able to obtain an unbiased gradient estimator ${\displaystyle H(\theta ,X)}$. If ${\displaystyle X}$ is viewed as some "fundamental" underlying random process that is generated independently of ${\displaystyle \theta }$, and under some regularization conditions for derivative-integral interchange operations so that ${\displaystyle \operatorname {E} {\Big [}{\frac {\partial }{\partial \theta }}Q(\theta ,X){\Big ]}=\nabla g(\theta )}$, then ${\displaystyle H(\theta ,X)={\frac {\partial }{\partial \theta }}Q(\theta ,X)}$ gives the fundamental gradient unbiased estimate. However, for some applications we have to use finite-difference methods in which ${\displaystyle H(\theta ,X)}$ has a conditional expectation close to ${\displaystyle \nabla g(\theta )}$ but not exactly equal to it.

We then define a recursion analogously to Newton's Method in the deterministic algorithm:

$${\displaystyle \theta _{n+1}=\theta _{n}-\varepsilon _{n}H(\theta _{n},X_{n+1}).}$$

The following result gives sufficient conditions on ${\displaystyle \theta _{n}}$ for the algorithm to converge:[18]

C1) ${\displaystyle \varepsilon _{n}\geq 0,\forall \;n\geq 0.}$

C2) ${\displaystyle \sum _{n=0}^{\infty }\varepsilon _{n}=\infty }$

C3) ${\displaystyle \sum _{n=0}^{\infty }\varepsilon _{n}^{2}<\infty }$

C4) ${\displaystyle |X_{n}|\leq B,{\text{ for a fixed bound }}B.}$

C5) ${\displaystyle g(\theta ){\text{ is strictly convex, i.e.}}}$

$${\displaystyle \inf _{\delta \leq |\theta -\theta ^{*}|\leq 1/\delta }\langle \theta -\theta ^{*},\nabla g(\theta )\rangle >0,{\text{ for every }}0<\delta <1.}$$

Then ${\displaystyle \theta _{n}}$ converges to ${\displaystyle \theta ^{*}}$ almost surely.

Here are some intuitive explanations about these conditions. Suppose ${\displaystyle H(\theta _{n},X_{n+1})}$ is a uniformly bounded random variables. If C2) is not satisfied, i.e. ${\displaystyle \sum _{n=0}^{\infty }\varepsilon _{n}<\infty }$ , then

$${\displaystyle \theta _{n}-\theta _{0}=-\sum _{i=0}^{n-1}\varepsilon _{i}H(\theta _{i},X_{i+1})}$$

is a bounded sequence, so the iteration cannot converge to ${\displaystyle \theta ^{*}}$ if the initial guess ${\displaystyle \theta _{0}}$ is too far away from ${\displaystyle \theta ^{*}}$. As for C3) note that if ${\displaystyle \theta _{n}}$ converges to ${\displaystyle \theta ^{*}}$ then

$${\displaystyle \theta _{n+1}-\theta _{n}=-\varepsilon _{n}H(\theta _{n},X_{n+1})\rightarrow 0,{\text{ as }}n\rightarrow \infty .}$$

so we must have ${\displaystyle \varepsilon _{n}\downarrow 0}$ ，and the condition C3) ensures it. A natural choice would be ${\displaystyle \varepsilon _{n}=1/n}$. Condition C5) is a fairly stringent condition on the shape of ${\displaystyle g(\theta )}$; it gives the search direction of the algorithm.

Suppose ${\displaystyle Q(\theta ,X)=f(\theta )+\theta ^{T}X}$, where ${\displaystyle f}$ is differentiable and ${\displaystyle X\in \mathbb {R} ^{p}}$ is a random variable independent of ${\displaystyle \theta }$. Then ${\displaystyle g(\theta )=\operatorname {E} [Q(\theta ,X)]=f(\theta )+\theta ^{T}\operatorname {E} X}$ depends on the mean of ${\displaystyle X}$, and the stochastic gradient method would be appropriate in this problem. We can choose ${\displaystyle H(\theta ,X)={\frac {\partial }{\partial \theta }}Q(\theta ,X)={\frac {\partial }{\partial \theta }}f(\theta )+X.}$

1. The Kiefer Wolfowitz algorithm requires that for each gradient computation, at least ${\displaystyle d+1}$ different parameter values must be simulated for every iteration of the algorithm, where ${\displaystyle d}$ is the dimension of the search space. This means that when ${\displaystyle d}$ is large, the Kiefer–Wolfowitz algorithm will require substantial computational effort per iteration, leading to slow convergence.
   1. To address this problem, Spall proposed the use of simultaneous perturbations to estimate the gradient. This method would require only two simulations per iteration, regardless of the dimension ${\displaystyle d}$.[20]
2. In the conditions required for convergence, the ability to specify a predetermined compact set that fulfills strong convexity (or concavity) and contains the unique solution can be difficult to find. With respect to real world applications, if the domain is quite large, these assumptions can be fairly restrictive and highly unrealistic.