As Paul states, without more information, it is hard to give advice without assumptions.
With 10-20 variables and expensive function evaluations, the tendency is to recommend derivative-free optimization algorithms. I am going to disagree strongly with Paul's advice: you generally need a machine-precision gradient unless you're using some sort of special method (for instance, stochastic gradient descent in machine learning will exploit the form of the objective to come up with reasonable gradient estimates).
Each quasi-Newton step is going to be of the form:
\begin{align}
\tilde{H}(x_{k})d_{k} = -\nabla{f}(x_{k}),
\end{align}
where $\tilde{H}$ is some approximation of the Hessian matrix, $d_{k}$ is the search direction, $x_{k}$ is the value of the decision variables at the current iterate, $f$ is your objective function, and $\nabla{f}$ is the gradient of your objective, and the decision variables are updated like $x_{k+1} = x_{k} + \alpha_{k}d_{k}$, where $\alpha_{k}$ is a step size determined in some fashion (like a line search). You can get away with approximating the Hessian in certain ways and your iterations will converge, although if you use something like finite difference approximations of the Hessian via exact gradients, you might suffer from issues due to ill-conditioning. Typically, the Hessian is approximated using the gradient (for instance, BFGS type methods with rank-1 updates to the Hessian).
Approximating the Hessian and the gradient both via finite differences is a bad idea for a number of reasons:
- you're going to have error in the gradient, so the quasi-Newton method you're applying is akin to finding the root of a noisy function
- if function evaluations are expensive and you're trying to evaluate a gradient with respect to $N$ variables, it's going to cost you $N$ function evaluations per iteration
- if you have error in the gradient, you're going to have more error in your Hessian, which is a big problem in terms of the conditioning of the linear system
- ...and it's going to cost you $N^{2}$ function evaluations per iteration
So to get one bad iteration of quasi-Newton, you're doing something like up to 420 function evaluations at 30 minutes per evaluation, which means you're either going to be waiting a while for each iteration, or you're going to need a big cluster just for the function evaluations. The actual linear solves are going to be of 20 by 20 matrices (at most!), so there's no reason to parallelize those. If you can get gradient information by, for instance, solving an adjoint problem, then it might be more worthwhile, in which case, it might be worth looking at a book like Nocedal & Wright.
If you're going to do a lot of function evaluations in parallel, you should look instead at surrogate modeling approaches or at generating set search methods before considering quasi-Newton approaches. The classic review articles are the one by Rios and Sahinidis on derivative-free methods, which was published in 2012 and provides a really good, broad comparison; the benchmarking article by More and Wild from 2009; the 2009 textbook "Introduction to Derivative-Free Optimization" by Conn, Scheinberg, and Vicente; and the review article on generating set search methods by Kolda, Lewis, and Torczon from 2003.
As linked above, the DAKOTA software package will implement some of those methods, and so will NLOPT, which implements DIRECT, and a few of Powell's surrogate modeling methods. You might also take a look at MCS; it's written in MATLAB, but maybe you can port the MATLAB implementation to the language of your choice. DAKOTA's essentially a collection of scripts you can use to run your expensive simulation and collect data for optimization algorithms, and NLOPT has interfaces in a large number of languages, so choice of programming language shouldn't be a huge issue in using either software package; DAKOTA does take a while to learn, though, and has a massive amount of documentation to sift through.
