Stochastic approximation explained

Stochastic approximation methods are a family of iterative methods typically used for root-finding problems or for optimization problems. The recursive update rules of stochastic approximation methods can be used, among other things, for solving linear systems when the collected data is corrupted by noise, or for approximating extreme values of functions which cannot be computed directly, but only estimated via noisy observations.

In a nutshell, stochastic approximation algorithms deal with a function of the form $f(\backslash theta)\; =\; \backslash operatorname\; E\_\; [F(\backslash theta,\backslash xi)]$which is the expected value of a function depending on a random variable $\backslash xi$. The goal is to recover properties of such a function $f$ without evaluating it directly. Instead, stochastic approximation algorithms use random samples of $F(\backslash theta,\backslash xi)$ to efficiently approximate properties of $f$ such as zeros or extrema.

Recently, stochastic approximations have found extensive applications in the fields of statistics and machine learning, especially in settings with big data. These applications range from stochastic optimization methods and algorithms, to online forms of the EM algorithm, reinforcement learning via temporal differences, and deep learning, and others.^[1] Stochastic approximation algorithms have also been used in the social sciences to describe collective dynamics: fictitious play in learning theory and consensus algorithms can be studied using their theory.^[2]

The earliest, and prototypical, algorithms of this kind are the Robbins–Monro and Kiefer–Wolfowitz algorithms introduced respectively in 1951 and 1952.

Robbins–Monro algorithm

The Robbins–Monro algorithm, introduced in 1951 by Herbert Robbins and Sutton Monro,^[3] presented a methodology for solving a root finding problem, where the function is represented as an expected value. Assume that we have a function $M(\backslash theta)$, and a constant $\backslash alpha$, such that the equation $M(\backslash theta)\; =\; \backslash alpha$ has a unique root at $\backslash theta^*$. It is assumed that while we cannot directly observe the function $M(\backslash theta)$, we can instead obtain measurements of the random variable $N(\backslash theta)$ where $\backslash operatorname\; E[N(\backslash theta)]\; =\; M(\backslash theta)$. The structure of the algorithm is to then generate iterates of the form:

\theta _{n+1} = \theta _{n} - a _{n(N(\theta} _{n)} - \alpha)

Here,

a _{1,} a _{2,} ... is a sequence of positive step sizes. Robbins and Monro proved^{, Theorem 2} that \theta _{n} converges in L ^{2} (and hence also in probability) to \theta ^{*} , and Blum^[4] later proved the convergence is actually with probability one, provided that:

• $N(\backslash theta)$ is uniformly bounded,
• $M(\backslash theta)$ is nondecreasing,
• $M\text{'}(\backslash theta^*)$ exists and is positive, and
• The sequence $a\_n$ satisfies the following requirements:

\sum ^{inf ty} _{n =0} a _{n} = inf ty and \sum ^{inf ty} _{n =0} a ^{2} _{n} < inf ty

A particular sequence of steps which satisfy these conditions, and was suggested by Robbins–Monro, have the form: $a\_n=a/n$, for $a\; >\; 0$. Other series are possible but in order to average out the noise in $N(\backslash theta)$, the above condition must be met.

Complexity results

1. If $f(\backslash theta)$ is twice continuously differentiable, and strongly convex, and the minimizer of $f(\backslash theta)$ belongs to the interior of $\backslash Theta$, then the Robbins–Monro algorithm will achieve the asymptotically optimal convergence rate, with respect to the objective function, being $\backslash operatorname\; E[f(\backslash theta\_n)\; -\; f^*]\; =\; O(1/n)$, where $f^*$ is the minimal value of $f(\backslash theta)$ over $\backslash theta\; \backslash in\; \backslash 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 $O(1/\backslash sqrt)$. They have also proven that this rate cannot be improved.

Subsequent developments and Polyak–Ruppert averaging

While the Robbins–Monro algorithm is theoretically able to achieve $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.^[8]

Chung (1954)^[9] and Fabian (1968)^[10] showed that we would achieve optimal convergence rate $O(1/\backslash sqrt)$ with $a\_n=\backslash bigtriangledown^2f(\backslash theta^*)^/n$ (or $a\_n=\backslash frac$). Lai and Robbins^{[11] [12]} designed adaptive procedures to estimate $M\text{'}(\backslash theta^*)$ such that $\backslash 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:$$\backslash theta\_\; -\; \backslash theta\_n\; =\; a\_n(\backslash alpha\; -\; N(\backslash theta\_n)),\; \backslash qquad\; \backslash bar\_n\; =\; \backslash frac\; \backslash sum^\_\; \backslash theta\_i$$The convergence of

\bar{\theta} _{n} to the unique root \theta ^{*} relies on the condition that the step sequence \{a _{n\}} decreases sufficiently slowly. That is

A1) $$a\_n\; \backslash rightarrow\; 0,\; \backslash qquad\; \backslash frac\; =\; o(a\_n)$$

Therefore, the sequence $a\_n\; =\; n^$ with satisfies this restriction, but $\backslash 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 $O(1/\backslash sqrt)$ yet with a more robust step size policy. 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 $O(1/\backslash sqrt)$.

A more general result is given in Chapter 11 of Kushner and Yin^[17] by defining interpolated time $t\_n=\backslash sum\_^a\_i$, interpolated process $\backslash theta^n(\backslash cdot)$ and interpolated normalized process $U^n(\backslash cdot)$ as

where $g(\backslash theta)\; =\; \backslash operatorname[Q(\backslash theta,X)]$ is differentiable and convex, then this problem is equivalent to find the root

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

$$\backslash operatorname[H(\backslash theta,X)|\backslash theta\; =\; \backslash theta\_n]\; =\; \backslash nabla\; g(\backslash theta\_n).$$

Here

H(\theta, X) is an unbiased estimator of \nabla g(\theta) . If X depends on \theta , there is in general no natural way of generating a random outcome 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 H(\theta, X) . If X is viewed as some "fundamental" underlying random process that is generated independently of \theta , and under some regularization conditions for derivative-integral interchange operations so that \operatorname{E} [ (\partial)/(\partial\theta) Q(\theta,X) ] = \nabla g(\theta) , then H(\theta, X) = (\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 H(\theta, X) has a conditional expectation close to \nabla g(\theta) but not exactly equal to it.

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

$$\backslash theta\_\; =\; \backslash theta\_n\; -\; \backslash varepsilon\_n\; H(\backslash theta\_n,X\_).$$

Convergence of the algorithm

The following result gives sufficient conditions on

\theta _{n} for the algorithm to converge:^[18]

C1)

\varepsilon _{n} \geq 0 , \forall n\geq 0 .

C2)

\sum ^{inf ty} _{n =0} \varepsilon _{n} = inf ty

C3)

\sum ^{inf ty} _{n =0} \varepsilon ^{2} _{n} < inf ty

C4)

|X _{n|} \leq B, forafixedbound B.

C5)

g(\theta) isstrictlyconvex,i.e.

Then

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

Here are some intuitive explanations about these conditions. Suppose

H(\theta _{n,} X _{n+1} ) is a uniformly bounded random variables. If C2) is not satisfied, i.e. \sum ^{inf ty} _{n =0} \varepsilon _{n} < inf ty , then$$\backslash theta\_n\; -\; \backslash theta\_0\; =\; -\backslash sum\_^\; \backslash varepsilon\_i\; H(\backslash theta\_i,\; X\_)$$is a bounded sequence, so the iteration cannot converge to \theta ^{*} if the initial guess \theta _{0} is too far away from \theta ^{*} . As for C3) note that if \theta _{n} converges to \theta ^{*} then

$$\backslash theta\_\; -\; \backslash theta\_n\; =\; -\backslash varepsilon\_n\; H(\backslash theta\_n,\; X\_)\; \backslash rightarrow\; 0,\; \backslash text\; n\backslash rightarrow\; \backslash infty.$$ so we must have

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

Example (where the stochastic gradient method is appropriate)

Suppose

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

Kiefer–Wolfowitz algorithm

The Kiefer–Wolfowitz algorithm was introduced in 1952 by Jacob Wolfowitz and Jack Kiefer,^[19] and was motivated by the publication of the Robbins–Monro algorithm. However, the algorithm was presented as a method which would stochastically estimate the maximum of a function.

Let

M(x) be a function which has a maximum at the point \theta . It is assumed that M(x) is unknown; however, certain observations N(x) , where \operatorname E[N(x)] = M(x) , can be made at any point x . The structure of the algorithm follows a gradient-like method, with the iterates being generated asx _{n+1} = x _{n} + a _{n ⋅ \left( (N(x _{n} + c _{n)} - N(x _{n} - c _{n)})/(2 c _{n})} \right) where N(x _{n+c} _{n)} and N(x _{n-c} _{n)} are independent. At every step, the gradient of M(x) is approximated akin to a central difference method with h =2 c _{n} . So the sequence \{c _{n\}} specifies the sequence of finite difference widths used for the gradient approximation, while the sequence \{a _{n\}} specifies a sequence of positive step sizes taken along that direction.

Kiefer and Wolfowitz proved that, if

M(x) satisfied certain regularity conditions, then x _{n} will converge to \theta in probability as n\to inf ty , and later Blum in 1954 showed x _{n} converges to \theta almost surely, provided that:\operatorname{Var}(N(x))\le S< inf ty for all x .

• The function

M(x) has a unique point of maximum (minimum) and is strong concave (convex)

• • The algorithm was first presented with the requirement that the function

M( ⋅ ) maintains strong global convexity (concavity) over the entire feasible space. Given this condition is too restrictive to impose over the entire domain, Kiefer and Wolfowitz proposed that it is sufficient to impose the condition over a compact set C _{0} \subset R ^{d} which is known to include the optimal solution.

• The function

M(x) satisfies the regularity conditions as follows:

• • There exists

\beta> 0 and B> 0 such that $$|x\text{'}-\backslash theta|+|x$$-\theta|<\beta \quad \Longrightarrow \quad |M(x')-M(x)||

• • There exists

\rho> 0 and R> 0 such that $$|x\text{'}-x$$|<\rho \quad \Longrightarrow \quad |M(x')-M(x)|

• • For every

\delta> 0 , there exists some \pi(\delta)> 0 such that $$|z-\backslash theta|>\backslash delta\; \backslash quad\; \backslash Longrightarrow\; \backslash quad\; \backslash inf\_\backslash frac$$>\pi(\delta)

• The selected sequences

\{a _{n\}} and \{c _{n\}} must be infinite sequences of positive numbers such thatc _{n} → 0 as n\to inf ty \sum ^{inf ty} _{n =0} a _{n} =inf ty \sum ^{inf ty} _{n =0} a _{nc} _{n} < inf ty \sum ^{inf ty} _{n =0} a ^{-2} _{n} < inf ty

A suitable choice of sequences, as recommended by Kiefer and Wolfowitz, would be

a _{n} = 1 /n and c _{n} = n ^{-1/3} .

Subsequent developments and important issues

1. The Kiefer Wolfowitz algorithm requires that for each gradient computation, at least

d +1 different parameter values must be simulated for every iteration of the algorithm, where d is the dimension of the search space. This means that when d is large, the Kiefer–Wolfowitz algorithm will require substantial computational effort per iteration, leading to slow convergence.

1. 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

d .^[20]

1. 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.

Further developments

An extensive theoretical literature has grown up around these algorithms, concerning conditions for convergence, rates of convergence, multivariate and other generalizations, proper choice of step size, possible noise models, and so on.^{[21] [22]} These methods are also applied in control theory, in which case the unknown function which we wish to optimize or find the zero of may vary in time. In this case, the step size

a _{n} should not converge to zero but should be chosen so as to track the function.^{[21] , 2nd ed., chapter 3}

C. Johan Masreliez and R. Douglas Martin were the first to apply stochastic approximation to robust estimation.^[23]

The main tool for analyzing stochastic approximations algorithms (including the Robbins–Monro and the Kiefer–Wolfowitz algorithms) is a theorem by Aryeh Dvoretzky published in 1956.^[24]

See also

• Stochastic gradient descent
• Stochastic variance reduction

Notes and References

1. Toulis . Panos . Edoardo . Airoldi. Scalable estimation strategies based on stochastic approximations: classical results and new insights . Statistics and Computing . 25 . 4 . 2015 . 781–795. 10.1007/s11222-015-9560-y. 26139959 . 4484776 .
2. Web site: Le Ny. Jerome. Introduction to Stochastic Approximation Algorithms. Polytechnique Montreal. Teaching Notes. 16 November 2016.
3. Robbins . H. . Herbert Robbins. Monro . S. . 10.1214/aoms/1177729586 . A Stochastic Approximation Method . The Annals of Mathematical Statistics . 22 . 3 . 400 . 1951 . free .
4. Blum. Julius R.. 1954-06-01. Approximation Methods which Converge with Probability one. The Annals of Mathematical Statistics. EN. 25. 2. 382–386. 10.1214/aoms/1177728794. 0003-4851. free.
5. Sacks . J. . Asymptotic Distribution of Stochastic Approximation Procedures . 10.1214/aoms/1177706619 . The Annals of Mathematical Statistics . 29 . 2 . 373–405 . 1958 . 2237335. free .
6. Nemirovski . A. . Arkadi Nemirovski. Juditsky . A. . Lan . G. . Shapiro . A. . Robust Stochastic Approximation Approach to Stochastic Programming . 10.1137/070704277 . SIAM Journal on Optimization . 19 . 4 . 1574 . 2009 .
7. Problem Complexity and Method Efficiency in Optimization, A. Nemirovski and D. Yudin, Wiley -Intersci. Ser. Discrete Math 15 John Wiley New York (1983) .
8. https://books.google.com/books?id=f66OIvvkKnAC&q=%22Robbins-Monro%22 Introduction to Stochastic Search and Optimization: Estimation, Simulation and Control
9. Chung. K. L.. 1954-09-01. On a Stochastic Approximation Method. The Annals of Mathematical Statistics. EN. 25. 3. 463–483. 10.1214/aoms/1177728716. 0003-4851. free.
10. Fabian. Vaclav. 1968-08-01. On Asymptotic Normality in Stochastic Approximation. The Annals of Mathematical Statistics. EN. 39. 4. 1327–1332. 10.1214/aoms/1177698258. 0003-4851. free.
11. Lai. T. L.. Robbins. Herbert. 1979-11-01. Adaptive Design and Stochastic Approximation. The Annals of Statistics. EN. 7. 6. 1196–1221. 10.1214/aos/1176344840. 0090-5364. free.
12. Lai. Tze Leung. Robbins. Herbert. 1981-09-01. Consistency and asymptotic efficiency of slope estimates in stochastic approximation schemes. Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete. en. 56. 3. 329–360. 10.1007/BF00536178. 122109044. 0044-3719. free.
13. Polyak. B T. 1991. New stochastic approximation type procedures. (In Russian.). Automation and Remote Control. 7. 7.
14. Ruppert. David. Efficient estimators from a slowly converging robbins-monro process. Technical Report 781. Cornell University School of Operations Research and Industrial Engineering. 1988.
15. Polyak . B. T. . Juditsky . A. B. . 10.1137/0330046 . Acceleration of Stochastic Approximation by Averaging . SIAM Journal on Control and Optimization . 30 . 4 . 838 . 1992 .
16. On Cezari's convergence of the steepest descent method for approximating saddle points of convex-concave functions, A. Nemirovski and D. Yudin, Dokl. Akad. Nauk SSR 2939, (1978 (Russian)), Soviet Math. Dokl. 19 (1978 (English)).
17. Book: Stochastic Approximation and Recursive Algorithms and Harold Kushner Springer. www.springer.com. 2016-05-16. 9780387008943. Kushner. Harold. George Yin. G.. 2003-07-17.
18. Book: Numerical Methods for stochastic Processes. Bouleau. N.. Lepingle. D.. John Wiley. 1994. 9780471546412. New York.
19. Kiefer . J. . Wolfowitz . J. . 10.1214/aoms/1177729392 . Stochastic Estimation of the Maximum of a Regression Function . The Annals of Mathematical Statistics . 23 . 3 . 462 . 1952 . free .
20. Spall . J. C. . Adaptive stochastic approximation by the simultaneous perturbation method . 10.1109/TAC.2000.880982 . IEEE Transactions on Automatic Control . 45 . 10 . 1839–1853 . 2000 .
21. Book: Kushner . H. J. . Harold J. Kushner. Yin . G. G. . 10.1007/978-1-4899-2696-8 . Stochastic Approximation Algorithms and Applications . 1997 . 978-1-4899-2698-2 .
22. Stochastic Approximation and Recursive Estimation, Mikhail Borisovich Nevel'son and Rafail Zalmanovich Has'minskiĭ, translated by Israel Program for Scientific Translations and B. Silver, Providence, RI: American Mathematical Society, 1973, 1976. .
23. Martin . R. . Masreliez . C. . 10.1109/TIT.1975.1055386 . Robust estimation via stochastic approximation . IEEE Transactions on Information Theory . 21 . 3 . 263 . 1975 .
24. Dvoretzky . Aryeh . Aryeh Dvoretzky . Neyman . Jerzy . Jerzy Neyman . On stochastic approximation . https://projecteuclid.org/euclid.bsmsp/1200501645 . 84911 . 39–55 . University of California Press . Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1954–1955, vol. I . 1956.

This article is licensed under the GNU Free Documentation License. It uses material from the Wikipedia article "Stochastic approximation".

Except where otherwise indicated, Everything.Explained.Today is © Copyright 2009-2024, A B Cryer, All Rights Reserved. Cookie policy.