How to find the gradient when a black box I/O function is involved in evaluation of the loss? I am trying to learn a neural network $NN_\pi$ to minimize the loss function
$$  L_{\pi} = || Y_{true} - F(X_{true},    NN_{\pi}(X_{true}) ) ||^2 $$
where $F$ is a black box (I/O) function (we only have access to the output given an input to $F$, but we have no way to compute the gradient of $F$), and $X_{true},Y_{true}$ are training data.
The issue is we cannot backpropagate as we cannot find gradient of $L$ (as $F$ is an oracle function, we cannot differentiate it).
Can it be possible to learn $NN_{\pi}$ without knowing $F$?  Or, can we find the gradient with any other methods?
In my application, $X \in \mathbb{R}^2$, $NN_\pi(X) \in \mathbb{R}$, and $F(\cdots),Y \in \mathbb{R}^2$, but I am interested in general solutions.  None of these variables will be in a high dimensional space.
 A: It's not possible to learn $NN_{pi}$ without (at least approx.) knowing $\mathbf F$. You can approximate the derivative of $\mathbf 
 F$ wrt the output of the neural network, say $o$, numerically. Then, apply back-propagation steps to the rest as usual. The numerical derivative can be calculated in various ways. One simple formula is as follows:
$$\mathbf F'(o)\approx\frac{\mathbf F(o+h)-\mathbf F(o-h)}{2h}$$
you should choose $h$ small.
A: Here are two methods you can use.
Numerical derivative
As @gunes explains, you can estimate the gradient of $F$ using numerical differentiation methods, such as finite differences.  In particular, the partial derivative $\frac{\partial F}{\partial x_i}(x)$ can be approximated as
$$\frac{\partial F}{\partial x_i}(x) \approx \frac{F(x + \delta e_i) - F(x - \delta e_i)}{2\delta},$$
where $\delta>0$ is sufficiently small and $e_i=(0,\dots,0,1,0,\dots,0)$ is a vector that is zeros in all coordinates except in the $i$th dimension.  If $F$ is a function of $n$ variables, doing this $n$ times will give you an approximation of the gradient of $F$ at $x$.  This then lets you compute the gradient of the loss function $L$.
Now you can use gradient descent to train the neural network $NN$.  You will use gradient descent to minimize the loss function.  This requires knowledge of the gradient of $F$, which you can compute using the method above.
How well this will work is likely to depend on how "locally smooth" $F$ is.
Approximation with a neural network
A different option is to train a neural network $\tilde{F}$ to approximate the function $F$.  Then, $\tilde{F}$ will be differentiable, so you can minimize the loss
$$\tilde{L}_{\pi} = || Y_{true} - \tilde{F}(X_{true},    NN_{\pi}(X_{true}) ) ||^2 $$
directly using gradient descent.
How do you train $\tilde{F}$?  Well, you pick a random input $x$ and some plausible $\pi$, compute $y=F(x, NN_{\pi}(x))$, and then add $(w,y)$ to the training set for $\tilde{F}$, where $w=(x, NN_{\pi}(x))$.  Repeat many times until you have a sufficiently large training set for $\tilde{F}$, then use standard methods for training neural nets to train $\tilde{F}$ on this training set.
A tricky bit is that the there is a circular dependency here: the training set for $\tilde{F}$ depends on $\pi$, but $\pi$ is obtained by training $NN$ on a training set and you have to have $\tilde{F}$ to do that.  So, it might be better to use optimization to jointly optimize both $\pi$ and $\tilde{F}$.  In particular, we might construct the following loss function:
$$L^* = \| Y - \tilde{F}_{\rho}(X, NN_{\pi}(X) )\|^2 + \lambda \cdot \| \tilde{F}_{\rho}(X, NN_{\pi}(X)) - F(X, NN_{\pi}(X))) \|^2,$$
where here $(X,Y)$ are chosen randomly from the training distribution (i.e., the training set).  Then we can use gradient descent to simultaneously find parameters $\pi,\rho$
that minimize $\mathbb{E}[L^*]$.  Note that, for a fixed $X,Y$, you can compute the gradient of $L^*$ with respect to $\rho,\pi$ without difficulty.
Here $\lambda>0$ is some hyperparameter that you can set using cross-validation.
How well this will work likely depends on how well-behaved $F$ is and how easy it is to approximate, given some input-output samples for it.
A: It is always possible to reframe model-fitting problems as reinforcement learning problems. The actions are the choice of model parameters and the reward is the negative of the error (here loss) function.
A possible implementation would be to use baseline comparison policy-gradient (aka the REINFORCE algorithm) . This requires that the likelihood of an action (here choice of neural-net parameters) is a differential function of its parameters.
A possible method

*

*Sample the weights/parameters to be used in the neural-net from a convenient distribution. An independent normal distribution for each parameter would be a convenient choice for this action distribution.$w_i\sim\left(\mu_i,\sigma^2\right)$.

*Use those parameters in the neural-net with the inputs of an example drawn from the training set and the weights drawn from the action distribution i.e. calculate $F\left(X,NN\left(X;\mathbf{w}\right)\right)$

*Calculate the loss function $L$ for the sampled example.

*Perform updates for all the parameters $\Delta \mu_i=\alpha\left(b-L \right)\left(w_i-\mu_i \right)$

*And repeat.

There are three hyper-parameters to choose with this method. 1) The variance of the parameters tried ($\sigma^2$): It is possible to make this a learnable parameter as well. But this can make choosing the learning rate trickier. It's probably safest to go small. Also, it is important to scale the inputs so that they have roughly the same range. 2) The learning rate ($\alpha$): important in all stochastic gradient algorithms. To ensure convergence, a small value is usually best but doesn't have to be too small if the variance is also small. 3) The baseline that the observed loss is compared to ($b$): It is known that the mean of $L$ isn't optimal but some fixed value near to the expected mean is the safest choice.
Note while this is possible, it may not be very efficient. This would be a problem if the neural-net is very complex and takes a long time to calculate for each new set of parameters.
