optimization (using python) how to find the most suitable solver?

Question

I have a fitting routine set up. It works, but pretty slow. I was wondering if there is a better method to use. I checked my (forward) code against some literature data and at least I do have no bugs and mistakes in my problem.

I fit the steric mass action model to some some isotherms like the ones shown below.

Additionally, i assume some arbitrary numbers of unknown components and fit just the measured concentrations of my known components (it is a multicomponent system). Additionally, the measured compound itself consists of some charge variants with different adsorption characteristics. Therefore, I am playing around right know with the number of components.

The model for a three component system takes following form. (the first component (the salt concentration) has to be calculated numerically)

\begin{align} qs &=& \frac{\Lambda \, c_1}{c_1 + (\sigma_2 + \nu_2)* \{ k_2 \, ( \frac{qs}{c_1} )\}^{\nu_2 - 1} + (\sigma_3 + \nu_3)* \{ k_3 \, ( \frac{qs}{c_1} )\}^{\nu_3 - 1} } \\ \end{align}

\begin{align} q2 &=& \frac{\Lambda \, {k \, ( \frac{qs}{c_1} )^{\nu_2 - 1} } }{c_1 + (\sigma_2 + \nu_2)* \{ k_2 \, ( \frac{qs}{c_1} )\}^{\nu_2 - 1} + (\sigma_3 + \nu_3)* \{ k_3 \, ( \frac{qs}{c_1} )\}^{\nu_3 - 1} } \\ \end{align}

\begin{align} q3 &=& \frac{\Lambda \, {k \, ( \frac{qs}{c_1} )^{\nu_3 - 1} } }{c_1 + (\sigma_2 + \nu_2)* \{ k_2 \, ( \frac{qs}{c_1} )\}^{\nu_2 - 1} + (\sigma_3 + \nu_3)* \{ k_3 \, ( \frac{qs}{c_1} )\}^{\nu_3 - 1} } \\ \end{align}

qs is the stationary concentration of free ligands [mM] typically between 50 and 1500

$\Lambda$ is initial concentration of ligands; can be a fit parameter; typically between 100 and 1700

$k$ is a fit parameter typically between 0.1 and 40

$\sigma $ is a fit parameter typically between 1 and 50

$\nu $ is a fit parameter typically between 1 and 15

Since the solver fails on some sets of input parameters and i do not want to fiddle around, i use scipy`s basin-hopping solver together with the method L-BFGS-B (quite slow). The parameters are bound so they make physical sense. Starting points are generated randomly within the bounded domain.

I was wondering how someone used to optimization would set this problem up ? Would you try to provide an analytical gradient ?

Edit:

I do minimize the sum of squares. I thought about scipy`s leastsq, but it does not accept bounds. From the local solvers, fmin (using L-BFGS-B) was the most robust and fastest one.

Answer 1

Some general thoughts from someone with optimization experience:

• If you can provide an analytical gradient, do it. In addition, providing an analytical Hessian is even better. Calculating a gradient with finite differences results in a loss of precision for ill-conditioned problems; you'll lose about half your significant figures in the best case, and in worse cases, more than that. Low-rank approximations to the Hessian such as BFGS and TR1 will retard convergence somewhat, but the improvement from a low-rank Hessian approximation to an analytical Hessian is not as great as the improvement in going from a finite difference gradient to an analytical gradient.

• For a bound-constrained least-squares problem, methods like Gauss-Newton and L-BFGS-B are good first choices. For nontrivial constraints, you should look at an active set method (for instance, SLSQP in SciPy is an SQP method; OpenOpt interfaces to ALGENCAN, an augmented Lagrangian approach) or an interior point method (pyipopt interfaces to IPOPT). Active set methods tend to be good for small- to medium-scale problems; interior point methods tend to be good for large-scale problems.

• If your nonlinear least squares problem is nonconvex, you should consider a globalization strategy like multistart to augment your choice of local optimization solver in the previous bullet. An actual nonconvex optimizer would be preferable, but those methods tend to be limited to small scale problems, are intrusive, and generally do not have Python interfaces.