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.

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.
