Classic Neural Network Layout for regression: is the score function derivable? I need to solve an optimization problem based on a function f(X). This function is not known, but it can be estimated from a training set. So first I train a model, then I get the score function f(X), and then I can use this function to solve my optimization problem.

So far, in my problem, I used a simple polynomial regression, and hence the score function is a polynomial. With this score function, the optimization problem works very well since it is infinitely derivable.

Now I want to improve f using a more complex model, i.e. a classic neural network for regression. I wonder: what kind of properties does the score function have in neural networks? Is it at least continuous? Is it derivable?
It depends on what activation function you use.  With all standard activation functions, neural networks are continuous.  With ReLU activations (the most typical choice of activation function), they are piecewise differentiable.  The partial derivatives (the gradient) can be computed efficiently (using so-called "backpropagation"), and standard neural network frameworks will automate this process for you, as computing the derivative is part of how neural networks are trained.  Consequently, it is possible to use gradient descent methods to find $x$ that (approximately) minimizes $f(x)$, when $f(x)$ is instantiated as a neural network.  $f(x)$ is not convex so you are not guaranteed to find the global minimum but for many problems this may be adequate; in others, you can adjust gradient descent (e.g., using random restarts) to mitigate those problems.