# Min-max stochastic function of few variables

I'm facing an optimization problem I'll have to run (with a variety of different hyperparameters) and I'd like some tips on how best to optimize it. I have a function $f$ which I can only evaluate in a noisy fashion, that is, there is a random variance to it. I have some 7 variables in group "A", and another 4 variables in group "B". I would like to maximize the expectation of $f$ with regards to the A variables, while simultaneously minimizing with regards to "B". Each evaluation of $f$ takes a few seconds up to a minute, so I'm happy to run a relatively complicated algorithm using data from all my previous evaluations, in order to converge in less calls to $f$. The only constraint is that the variables are all positive. Since they're also operating at very different scales (some are ~1, some are ~0.001, some are ~1000), so I'm making the problem a bit easier by taking log of all the parameters first. That is, I'm actually evaluating $f(\exp(\bf x))$, and tweaking my parameter vector $\bf x$. I expect that shifting my destination optimum will be within $\pm 3$ in each coordinate of $\bf x$, then. (So, within an order of magnitude of my original parameters.) This also means that the problem is essentially completely unconstrained now.

Here is one idea I had: Start by evaluating $f$ at one initial point $\bf x_0$, and then evaluate $f$ at other values for $x_0$, with each coordinate $\pm 1$. Do a quadratic fit to these points (I already have a library set up to do that), and then inspect the quadratic term. If I had only group A variables, then I would check that the quadratic term is positive definite, and then run a test at the estimated global optimum of that quadratic. As I take more and more test points close to the optimum, the quadratic fit will begin to pay more attention to fitting correctly around there (since I have more samples), and so this should converge to the correct value. The downside is that it might end up weighted "incorrectly" by values far from the optimum for quite a while; and I'm not sure how I would handle when it's not positive definite (just run a random test point near an old expected optimum, maybe??); and I'm not sure how I would handle adding in the B-group variables would handle this.

One possibility would be writing it as the sum of a quadratic A in and a quadratic in B, and then moving to the maximum of one term and minimum of the other. But this neglects the A/B variable cross terms, which I expect to matter here. Again, ideally, as more data accumulate around the true optimum, it should converge (i.e. all the B data is collected near the A optimum, and vice versa), but that convergence could be very slow. And I'm still not sure how to approach the positive definiteness.

I could find some information on stochastic optimization algorithms, but they all seem to be much more designed for the case with very tight memory constraints, so that they avoid using earlier evaluations too much. For instance, this is the problem with the Python noisyopt algorithms. Since I'm min-maxing on different variables, also, I basically need to optimize with regards to A, and then holding A constant, optimize the other with regards to B, which isn't something that the most packages seem to like very much, so I think I need to roll my own algorithm here (and am quite happy to -- I just want these optimizations to be fast!)

The final code where I run many such optimization problems, I expect to execute over the course of a month or so as part of a research project, so I'm really happy writing my own algorithm if I think it will be more stable and require fewer invocations.

• I think this question would benefit from a reference to Generalized Nash Equilibrium Problems. hal.archives-ouvertes.fr/hal-00813531/document Commented Sep 26, 2017 at 20:13
• Upon further reading, this actually is a standard (as opposed to a generalized) nash equilibium problem, because the feasible set for each player's variables doesn't depend on the strategy chosen by the other player. Commented Sep 26, 2017 at 20:26
• You're right that this is a standard Nash equilibrium problem. I might edit the title/tags to reflect this. However, most references on that topic seem to assume perfect observation of the expectation, or -- even worse! -- gradient methods. Are you perhaps aware of robust methods for tackling Nash eq problems? It seems popular to rephrase them as Variational Inequality problems but I know almost nothing about those. Commented Sep 26, 2017 at 20:51
• There exists an example in which the Nelder Mead method has been applied to minimize the difference between the best response function and the current response. The Nelder Mead method is favourable for not requiring gradients or calculation of many data points. However it does require, in this particular application, the computation of the best response function which can be costly. Still you might be able to replace it with some other (simpler) measure that goes to zero as you approach the Nash equilibrium (e.g. a polynomial fit minimum). ieeexplore.ieee.org/abstract/document/4224254 Commented Sep 27, 2017 at 14:21
• This sounds like a variation on V-optimal DOE with online experiment repair. Commented Sep 27, 2017 at 14:40

I suggest trying Simultaneous Perturbation Stochastic Approximation . The link has a very clear beginner tutorial and a link to Matlab code; I'll alter the basic algorithm slightly to deal with the min / max issue and present some R code on a toy saddlepoint problem.

SPSA is designed for problems which involve minimizing or maximizing an unconstrained, continuous-in-expectation, stochastic function. It can be pretty easily extended to problems with a small number of constraints, and there is a root-finding variant as well. It has been extended to discrete problems as well. I've used it on problems with hundreds of parameters with success.

The basic idea is that, at every iteration $k$, we form an estimate of the gradient in a randomly-chosen direction by sampling the function $f$ at two points $\theta_k \pm c_k\Delta$ around the current parameter value $\theta_k$. $\Delta$ is typically chosen by generating a vector of $\pm 1$ values; $c_k$ is a scaling parameter that decreases slowly as the number of iterations increases. Our gradient estimate is formed in the obvious way:

$$g_k = (f(\theta + c_k\Delta)-f(\theta - c_k\Delta)/2c_k\Delta$$

where the division is done elementwise so $g$ is a vector of the appropriate length. Once we have a gradient, we step in the appropriate direction (depending on whether we are minimizing or maximizing) to get the next iteration's estimate of $\theta$

$$\theta_{k+1} = \theta_k - \alpha_kg_k$$

where $a_k$ is a stepsize that decreases as $k$ increases. Obviously, I've assumed we are facing a minimization problem in the line above. Guidelines for choosing the sequences $a_k, c_k$ abound, but it's really pretty easy and in my experience not much experimentation is needed.

In your case, since you are minimizing with respect to some parameters and maximizing with respect to others, you'd simply switch the sign of the derivative for the parameters you're maximizing on.

On to an example! We have the function $f(x) = x_1^2 - x_2^2$ with some noise (distributed $N(0,1)$) added to the return value. We will try to minimize with respect to the first parameter and maximize with respect to the second parameter. Using the excellent package plot3D in R, we plot the de-noised function in the vicinity of the saddlepoint:

The saddlepoint, which we are trying to find, is at $x = (0,0)$. We will start at $x=(2,-2)$, for which $f(x)=0$, just as it does at the saddlepoint. The R code, written for clarity instead of efficiency, is:

f <- function(x) {
x[1]*x[1] - x[2]*x[2] + rnorm(1)
}

# Parameters for SPSA
alpha <- 0.606  # Often used as a default value
gamma <- 0.16   # Often used as a default value
A <- 10         # Set at 5-10% of rough guess at max iteration count
a0 <- 1
c0 <- 1

x_history <- matrix(0,51,2)

# Starting values (optimum is (0,0))
# Note f at starting values = f at optimum = 0
x_history[1, ] <- c(2,-2)

# SPSA loop a fixed number of times, as this is a demo
for (k in 1:(nrow(x_history)-1)) {
# Step size parameters
ck <- c0 / (k^gamma)
ak <- a0 /(A+k)^alpha

# Random step direction
delta <- ck * (2*rbinom(2, 1, 0.5) - 1)

# Function evaluation + derivative calculation
f_plus <- f(x_history[k,] + delta)
f_minus <- f(x_history[k,] - delta)
deriv <- (f_plus - f_minus) / (2*delta)

# We are trying to maximize on the second parameter,
# so switch the sign of its "derivative"
deriv[2] <- -deriv[2]

x_history[k+1,] <- x_history[k,] - ak*deriv
}


The trace of the parameter values at each iteration looks like:

> matplot(x_history, xlab='Iteration count', ylab='Parameter values')


Or, superimposed on the de-noised function:

> f_values <- function(x,y) x*x - y*y
> M <- mesh(seq(-2.2, 2.2, by=0.1), seq(-2.2, 2.2, by=0.1))
> surf3D(x=M$x, y=M$y, z=f_values(M$x, M$y), theta=30)
> scatter3D(x=x_history[,1], y=x_history[,2], z=f_values(x_history[,1], x_history[,2]),


As we can see, the algorithm is, give or take some randomness, moving towards $x = (0,0)$ along the line $f(x) = 0$. In this iteration of the answer, it actually jumped almost all the way there on iteration 4.

One common approach to forming the final estimate of the optimum parameter values is to average over some number of iterations. This isn't useful if alpha = 1 for technical reasons, but is more useful as alpha -> 0.606 (the lowest it should be.) For our sample run, we get:

> colMeans(x_history[4:51,])
[1] 0.1106594 0.0896351
> f_values(0.1106594, 0.0896351)
[1] 0.004211052


which isn't too bad, given that we have done no parameter tuning. Note that because the values at each step are going to be autocorrelated, the sample standard deviation will be biased, although you could remove a lot of the bias by applying a simple ARMA model to the series over the iterations where it seems to have stabilized and using the resulting estimates.

Of course, your mileage will vary. It may be that you require several hundred iterations to achieve convergence, which given your function evaluation times might take several hours of computation. Parameter tuning may also take some considerable time. There are second-derivative versions that converge more quickly when you get in the neighborhood of the optimum, but require somewhat more coding and parameter tuning. There is a structured way of selecting the $\Delta$ at each iteration that may speed convergence a little. But, generally speaking, I haven't found any off-the-shelf techniques for optimizing stochastic functions that are as good across as wide a range of problems as SPSA and its variants.

• I appreciate this a good amount! But is there any nice way to extend it naturally to use more of the history? e.g. if it has sampled one point (or very close) 10 times already, to use all of that in conjunction with new samplings? Awarding the bounty anyway because expiry, but I'm really hoping there's a nice way to extend this, again... also, if you have tips/references on how to pick c_n / a_n that would be great. Commented Oct 2, 2017 at 8:28
• Thanks! The references are in the link. Typically, though, there's some trial and error about the initial values $c_0, a_0$. Insofar as using more of the history is concerned, if you use iterate averaging (which Spall doesn't like, but others do) you form your final estimate by averaging over the last $M$ iterations where $M$ is the range where you appear to be just bouncing around - as in the example, where I average over the last 60% of the observations. There are also extensions that provide updates of the gradient / Hessian based on more history. I'll extend the answer later today. Commented Oct 2, 2017 at 14:39
• I've updated the answer with some information about default values and iterate averaging. Since I altered one of the parameters in line with the default values, the results graphs have changed too. Commented Oct 3, 2017 at 1:03