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

Question

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}$)

Answer 1

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.

Answer 2

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.