While I am familiar to using solvers and Levenberg-Marquardt, I find it a pity that LMSolver, more precisely LMSolver::Callback, is so poorly documented.

Especially the part output Jacobian: J_ij = d(err_i)/d(param_j)
is d(err_i) really intended, or should it be d_i(f) ?

Since most of the time, users do not really derivate the function but compute an estimation (f(x+eps)-f(x-eps))/(2eps), would it be possible that such an estimated Jacobian could be the default, only to be overwritten if derivation is affordable ?