# How would you minimize the sum of squares if the predictive function is a black box?

I'm solving an optimization problem, using the mean squared error:

$$\arg\min_{\mathcal{M}} ||y - \hat{y}||$$

$$y$$ is the true value and $$\hat{y}$$ is obtained from some black box function. $$\mathcal{M}$$ is the set of weights/parameters that I am trying to optimize. Since the function that determines $$\hat{y}$$ is a black box, I obviously can't compute the gradient, so any gradient-based optimization methods like GD, Newton's, is out the window.

How would you approach solving this problem? What criteria do you use to determine which gradient-free optimization method is best?

Edit 1: In my particular case, $$\hat{y}$$ is determined from a scientific computing simulation code. It's essentially a measure of temperature at various locations in a material. The simulation code can be quite expensive to evaluate (maybe 5 minutes per iteration). Because I don't have a closed form solution for the black box function, I don't know if the objective function is convex or not. My suspicion is it's not because I think there are multiple parameters $$\mathcal{M}$$ that can result in the same objective function value.

$$\mathcal{M}$$ here is a set. The size of the set is about 150-180. Each variable $$\in \mathcal{M}$$ is continuous. $$\hat{y}, y$$ are vectors with about 10,000 values. (The black box simulation code outputs a 10,000-sized vector $$\hat{y}$$ for a given input $$\mathcal{M}$$)

• Can you give us some context? However you might want to look into Bayesian Optimisation, it works very well and doesnt require any gradient information Jun 15, 2020 at 5:20
• @jcken Sure. In my particular case, $\hat{y}$ is determined from a scientific computing simulation code. It's essentially a measure of temperature at various locations in a material. The simulation code can be quite expensive to evaluate (maybe 5 minutes per iteration). Because I don't have a closed form solution for the black box function, I don't know if the objective function is convex or not. My suspicion is it's not because I think there are multiple parameters $\mathcal{M}$ that can result in the same objective function value. Jun 15, 2020 at 6:02
• @jcken I've never used Bayesian optimization before. Are you familiar with particle swarm? It seems pretty similar? Jun 15, 2020 at 6:03
• Particle swarm looks similar to BO. However BO builds an approximate version of your blaxbox function and optimised that. It would also be useful to know how big your input space is, that can help you decide which method is best Jun 15, 2020 at 6:45
• @jcken Got it. Just edited the information into the OP under "#EDIT 1" Jun 15, 2020 at 7:48

Let me expand a bit over what has been discussed in the comments. The bottleneck of your problem is an expensive evaluation of an unknown black-box function $$f$$ and somewhat high dimensionality of a problem (if I understood correctly we're looking at a set of weights, $$\mathcal{M} = [0, 1]^{150}$$), which seems like an ideal problem to be solved by Bayesian Optimisation, which quantitatively represents the uncertainty of "unseen" regions in the search space, allowing for efficient selection of next evaluation candidates, making search very time-efficient.

However, Bayesian Optimisation might struggle with high-dimensional spaces, see A Tutorial on Bayesian Optimization, Peter I. Frazier or, as discussed in the abstract in High-dimensional Bayesian optimization using low-dimensional feature spaces, Riccardo Moriconi, Marc P. Deisenroth, K. S. Sesh Kumar:

Bayesian optimization (BO) is a powerful approach for seeking the global optimum of expensive black-box functions and has proven successful for fine tuning hyper-parameters of machine learning models. However, BO is practically limited to optimizing 10-20 parameters. To scale BO to high dimensions, we usually make structural assumptions on the decomposition of the objective and/or exploit the intrinsic lower dimensionality of the problem, e.g. by using linear projections.

Which indicates that most likely you will need to approach the problem from one of the two perspectives:

1. Heuristics
2. Proxy

Heuristic methods would include, as already mentioned Particle Swarm Optimisation, Genetic Algorithms, Simulated Annealing etc. but don't give you any guarantees about the optimum. However, if you're an expert (or have sufficient knowledge) in the field you're studying you might be better off by trying to define some assumptions about the shape/form of the underlying black box and then trying to fit a proxy function (see Surrogate model) that approximates your $$f$$, which would allow you to use e.g. gradient/hessian methods and find optimum quickly (with respect to the proxy function).

Lastly - again assuming you have enough expert knowledge - maybe you can reduce the dimensionality manually such that Bayesian Optimisation is feasible for your problem.

• Bayesian optimization doesn't give you any guarantees about the optimum either right? My understanding is that if the predictive function is a black box, hence you can't find the hessian and therefore can't get an idea of convexity, then it's impossible to know whether the solution to any optimization problem is actually the optimum? You can find that a solution is a local optimum if you run multiple trials that converged to different values, but it's not possible to know if the solution is a global optimum? Jun 15, 2020 at 20:40
• Entirely correct. BO requires a selection of an aquisition function that trades off exploration and exploitation so you're also not guaranteed to find a global optimum. I should have stressed that in my answer. Jun 15, 2020 at 20:47
• Does BO work if $\mathcal{M}$ here was inside some continuous function rather than discrete variables? In other words, I want to find an approximately function for $\mathcal{M}$? I've never solved this kind of problem before, but I was thinking of reposing my problem to be such that $\mathcal{M}$ is continuous. Jun 15, 2020 at 21:06

This is going to be a fairly general purpose solution to the problem, but I'm going to name drop some ideas.

Your computer model is essentially $$\mathbf{y} = f(\mathbf{x})$$ Where $$\mathbf{x}$$ has approximately a dimension of $$160$$ and $$\mathbf{y}$$ is of dimension $$10,000$$ (approx).

Your problem is quite high dimensional, I'm assuming your code is deterministic. The first think you should do is perform PCA on the $$\mathbf{y}$$ space to reduce it's dimension dramatically. There is a lot of info on PCA online, once you have performed the PCA call these new reduce dimension outputs $$\mathbf{z}$$ where $$dimension(\mathbf{z}) << 10,000$$. I suspect you could do some kind of dimension reduction of $$\mathbf{x}$$ too, but <$$200$$ dimensions might not be too difficult.

Now the simulation code is reasonably expensive, you're going to need some kind of surrogate model to make the computation feasible, for a general overview of surrogates see wikipedia or this recent open source book by Bobby Gramacy, he is one of the world's leading experts on surrogates. Since your problem is quite high dimensional you're probably going to want to build something like a Neural Network, a polynomial fit or perhaps a generalised additive model (GAM). A Gaussian process surrogate might not work very well here (although they're my go-to).

To build your surrogate (this might be a Gaussian process, a polynomial, neural network) by running the model at lots of different inputs (you will need to choose these carefully, e.g. by a Maximin Latin Hypercube design). We will now run the computer model lots of times and obtain data $$(\mathbf{x}_i,\mathbf{y}_i)$$; reduce the dimension of the $$\mathbf{y}_i$$ using the exact same algorithm as you did for $$\mathbf{y}$$. Our aim is to predict $$\mathbf{z}$$ using some kind of surrogate, we have data $$(\mathbf{x}_i, \mathbf{z}_i)$$ train your surrogate on this data. Denote predictions from the surrogate to be $$\hat{\mathbf{z}}(\mathbf{x})$$

We then want to minimise $$\Omega(\mathbf{x}) = ||\mathbf{z}_i - \hat{\mathbf{z}}(\mathbf{x})||$$ where $$|| \cdot ||$$ is some metric in the $$\mathbf{z}$$ space, e.g. euclidean distance.

I guess we are now at the point of answering your question: how to acutally minimise this thing.

In the past I've used the Nelder-Mead method with good success. There is an R implementation of Nelder-Mead and it's probably available in whatever programming language you're using. The optimisation will give you $$\hat{\bf{x}}_z =\text{argmin}_{\mathbf{x} \in \mathcal{M}} || \mathbf{z}_i - \hat{\mathbf{z}}(\mathbf{x}) ||$$ This will not be the ''true'' minimum $$\hat{\bf{x}} =\text{argmin}_{\mathbf{x} \in \mathcal{M}} || \mathbf{y}_i - \mathbf{y}(\mathbf{x}) ||$$ but we frequently have to make sacrifices in these high-dim settings.

As with any complex optimisation, run the optimisation a few times from different starting points to assess convergence. Finally, check that your optimal value $$\hat{\mathbf{x}}_z$$ is appropriate by computing $$\mathbf{y}(\hat{\mathbf{x}}_z)$$ against $$\mathbf{y}$$; the ''true'' values.

• What is the purpose of performing PCA on the output, $\hat{y}$? Isn't PCA typically performed on the input space, $\mathcal{M}$? Jun 15, 2020 at 20:21
• Basically performing PCA on the $\mathbf{y}$ space means I have to predict fewer things. I.e. I only need to predict the elements of $\mathbf{z}$, if the PCA is useful here this might only be a handful of numbers, rather than predicting the $10000$ elements of $\mathbf{y}$ Jun 15, 2020 at 20:26
• Is this to intended prevent overfitting, or do reduce computation time of the black box function? The latter actually won't be reduced because the black box function's evaluation time is independent of the size of the output vector, $\hat{y}$. It could be size 1 or size 10,000, and the computation time would be essentially the same. Jun 15, 2020 at 20:28
• It is just much simpler to predict a low dimension object than a high dim object. You're going to have to train a very complex surrogate to predict the $\mathbf{y}$ whereas predicting $\mathbf{z}$ will be relatively simple. If the dimension of y is 10K then I need to predict 10K things, if the dim of z in 100 I only need to predict 100 things. Training the surrogate will also be greatly simplified Jun 15, 2020 at 20:34
• Ah I think I see what you're getting at now. Just to confirm, typically PCA is used on the input variable, $\boldsymbol{x}$, in your notation, right? Jun 15, 2020 at 22:23