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.

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.
