# Gradient descent-like optimization on a convex landscape with noisy sampling

This is a rewrite of the original positing (below), and is crossposted to (https://cstheory.stackexchange.com/questions/21496/gradient-descent-like-optimization-on-a-convex-landscape-with-noisy-sampling):

I have a strictly convex function $f(x,y)$, i.e. something like $f(\theta,\phi) \approx t_1 Sin[\theta] + t_2 Sin[\phi]$ (where $(\theta,\phi) \in [0,\pi]$ and $(t_1,t_2) \in \mathbb{R}^+$), with a global minimum at some $h_{min}=(x_{min},y_{min})$.

We now play a kind of guessing game where you provide me some $h^*$ on a plane that is guaranteed to be within a distance $R$ of $h_{min}$, and I provide you with a coordinate $c_i = (x_i,y_i)$ that corresponds to a "guess" for the value of $h_{min}$. You then provide me the value $j_i = f(x_i,y_i) + q_i$ where $q_i \in \mathbb{R}$ corresponds to Gaussian noise drawn from a distribution with parameters $(\mu_n,\sigma_n)$, e.g. we could have $(\mu_n,\sigma_n) = (0,1)$. If it simplifies things, I suppose we could instead draw $q_i$ from some uniform distribution s.t. $q_i \in [A,B]$.

If you only allow me to perform $k$ guesses: what is my optimal strategy for guessing a set of coordinates $(c_1, c_2, ..., c_k)$ s.t. $c_k$ minimizes $\delta = (h_{min} - c_k)$? On average, how well can I do? Provided this optimal strategy, can we write down something like a probability distribution for $\delta$ as a function of $k$?

Original posting:

I have a bivariate Gaussian distribution $f(x,y)$ with mean $\mu$, covariance matrix $M$, and square roots of the eigenvalues of the covariance matrix $(\sigma_1,\sigma_2)$. We now play a kind of guessing game where you provide me some $\mu^*$ on a plane that is guaranteed to be within a distance $R$ of $\mu$, and I provide you with a coordinate $c_i = (x_i,y_i)$ that corresponds to a "guess" for the value of $\mu$. You then provide me the value $j_i = f(x_i,y_i) + q_i$ where $q_i \in \mathbb{R}$ corresponds to Gaussian noise drawn from a distribution with parameters $(\mu_n,\sigma_n)$, e.g. we could have $(\mu_n,\sigma_n) = (0,1)$. If it simplifies things, I suppose we could instead draw $q_i$ from some uniform distribution s.t. $q_i \in [A,B]$.

If you only allow me to perform $k$ guesses: what is my optimal strategy for guessing a set of coordinates $(c_1, c_2, ..., c_k)$ s.t. $c_k$ minimizes $\delta = (\mu - c_k)$? On average, how well can I do? Provided this optimal strategy, can we write down something like a probability distribution for $\delta$ as a function of $k$?

Update -- I understand that the bivariate Gaussian distribution $f(x,y)$ is not a convex function, so perhaps I shouldn't be using the "convex-optimization" tag, though this feels most relevant to this problem. That said, I would be happy to recast $f(x,y)$ as a convex or strictly convex function and ask the same question above, with the hope that convex optimization methods can be relevant and helpful.

Update 2 -- To simplify things, let's specifically call $f(x,y)$ a strictly convex function, i.e. something like $f(\theta,\phi) \approx c_1 Sin[\theta] + c_2 Sin[\phi]$ where $(\theta,\phi) \in [0,\pi]$ and $(c_1,c_2) \in \mathbb{R}^+$. The question now is, what is the best gradient descent-like strategy provided a strictly convex landscape with the aforementioned sampling noise, and how well can we do with this algorithm after some number of steps (guesses) $k$?

• You may want to try stats.stackexchange.com – Bjørn Kjos-Hanssen Mar 11 '14 at 22:51
• @BjørnKjos-Hanssen Do you there or something like cstheory.stackexchange.com? Is this really a statistics question? – O1155 Mar 11 '14 at 23:04
• @BjørnKjos-Hanssen Well, we can give it a shot. – O1155 Mar 11 '14 at 23:35
• search for: stochastic programming, noisy gradients, noisy function evaluations, etc. – Suvrit Mar 12 '14 at 1:04
• @Suvrit I've been going through a lot of the machine learning literature, but I've been having trouble finding instances were we're in the limit of only a few samplings (guesses) $k$. – O1155 Mar 12 '14 at 2:33

You may try starting a local search from $h^*$ bounding the region by constraint t the distance from $h^*$ to be less than $R$.

For the noise, there is quite a large number of papers on that. You could start from the "implicit filtering" method

http://www4.ncsu.edu/~ctk/iffco.html

and then search accordingly.

The optimization problem might look like: $$\min f(h) + q \quad s.t. ||h-h^*||_2\leq R$$

The guess provided by $c$ is more tricky to use. If it is trustable, i.e. it does not underestimate the minimum value of $f$, you can include an additional constraint as $f(h) \leq c$ which s still convex.

Note that this approach does not use any information on the involved distributions. The implicit filtering" only assume the noise to be small compared to the function you want to optimize.

General Discussion: Sampling

I'm inclined to take a "sampling" approach to the problem, which requires that $k$ is not especially small. Now, I am not going to present a proof that this is necessarily the "optimal" strategy here. My main motivation for presenting this method is that I know it is an effective technique for tolerating a stochastic objective function, and there are some tricks we can use to estimate the distribution of the resulting error.

If $k$ is large, we can treat $-f$ as though it were an unnormalized probability density and use sampling methods (e.g. MCMC) to draw $k$ samples from it. The samples will generally lie near the max of the function (since, in the context of a density, this corresponds to the region of highest probability) and so we can estimate the mode of the distribution from the samples via kernel density estimation, or alternatively take the expectation of our samples to estimate the coordinates of the minimum of $f$.

An alternative but related approach would be to use simulated annealing with cooling time equal to $k$, but again: this really only works if we're permitted a large enough $k$. I'll constrain my attention for the rest of the post to discussing the MCMC approach.

These strategies are targeted towards precision and generally require a large number of samples to produce good results, but they tolerate the "noisiness" of the objective function well. They are not especially "efficient" in small $k$.

Optimizaiton via Monte Carlo Integration

(a lot of general results about MCMC methods to follow)

Let's assume for the sake of argument that the surface we are optimizing over is structured such that when treated as a probability density, the mode is collocated with the mean (sort of a big assumption but just roll with it). Therefore, we can equivalently solve for this expectation to find the minimimum of the objective function. The method I proposed earlier of drawing samples and taking their mean then amounts to estimation of this quantity (which we are assuming is a suitable proxy for the max of $-f$) via monte carlo integration. This works because the law of large numbers says the sample mean will converge almost surely to the true mean as the number of samples approaches infinity.

So how good can we do? Let $c_t$ denote the $t^th$ candidate sampled from our markov chain. Then our estimator is:

$m_c = \frac{1}{k} \sum \limits_{t=1}^{k} c_t$

Then the error of our estimator is:

$error = m_c - h_{min} = \frac{1}{k} \sum \limits_{t=1}^{k} ( c_t - h_{min} ) = \delta(k)$

Which demonstrates that our error is improved by either increasing $k$ or "burning" the early samples (since samples drawn outside of the stationary distribution will disproportionately influence the error calculation).

Going a bit further

We can make a slightly more useful statement about the distribution of the error here by invoking the central limit theorem and essentially "translating it" into the contexts of a markov chain. For brevity I'm going to provide you with several results without their proofs because it's really late (I may come back to this post later and build out the proofs if there's demand for it) and will just assert that the following theorem holds given that we are already in the stationary distribution at $T=0$:

$\lim \limits_{T \rightarrow \infty} \frac{1}{\sqrt{T}} \sum \limits_{t=0}^{T} (c_t - h_{min} ) = N(0,\sigma_{TV}^2)$

where $\sigma_{TV}^2$ is the "time averaged variance", which is a function of the natural variance of the objective function, $\sigma_f^2$, and the autocorrelation of the sampling chain. Let $\rho(t)$ denote the $t$-lag autocorrelation of the markov chain. Then:

$\sigma_{TV}^2 = \frac{\sigma_f^2}{T} \sum \limits_{t=1}^{T} \sum \limits_{t'=1}^{T} \rho(|t-t'|)$

Now, the autocorrelation of the chain collapses on the scale of the mixing of our sampler, so now we can see that we can further improve our estimate by choosing a good proposal distribution and thereby giving us a low relaxation time. The equation for the time-averaged variance can be further reduced to:

$\sigma_{TV}^2 = \sigma_f^2 \alpha \tau$

Where $\tau$ is the relaxation time of the chain and $\alpha$ is some constant. Note: if the samples were drawn using an i.i.d. process (instead of a markov process, which introduces autocorrelation in our samples) then the time-averaged variance would reduce to the variance of the objective function, and the "MCMC translated CLT" just reduces to the CLT.

With these tools on hand, we can rewrite the distribution of the error as:

$m_c - h_{min} \sim N(0, \frac{\sigma_f^2 \alpha \tau}{k})$

So we have the result that the error using this method is distributed as a function of $k$ and the relaxation time of our chain.

Making this all practical

Obviously, this is of limited utility since we don't know what $\alpha$ is, but you can use the following method to estimate $\sigma_{TV}^2$:

let $\Delta t$ be an approximation for the chain's relaxation time $\tau$. Let

$w_i = \frac{1}{\Delta t} \sum \limits_{t=i \Delta t}^{(i+1) \Delta t} c_t$.

Then

$\frac{\sigma_{TV}^2}{k} \cong \frac{1}{j-1} \sum \limits_{i=0}^{j-1} (w_i - m_c)^2$

Where a recommended value for $j$ is generally around 20, or at least small enough to force $\Delta t >> \tau$. This method is known as "batch means variance estimation." You can read more about it here:

James M. Flegal, Galin L. Jones. "Batch means and spectral variance estimators in Markov chain Monte Carlo." Annals of Statistics 2010, Vol. 38, No. 2, 1034-1070