Specifically in the context of local solvers with and without constraints.
What is the standard way people report error bars on fitting parameters while using Optimization.jl?
In principle this can be calculated if the covariance matrix is known but not even that is returned upon solving.
Need more context. If you are fitting a curve using Optimization.jl, then it should just compute the best parameters according to your objective function; there is no error to plot (assuming the solver converged). Maybe you want to show a confidence interval for the predicted function values?
Precisely, calculating confidence intervals (what I call “error bars”) can be done straightforwardly if a covariance (or non-singular Hessian) is returned by the solver, see e.g. https://arxiv.org/pdf/1210.3781.pdf
This appears to not be the case in Optimization.jl
It’s generally not computed since most optimization would be more expensive if that was used. You can just directly calculate it.
I think it makes sense to add maybe a flag in the solve call to include confidence interval calculations in the output.
For a fit to be meaningful (or to be ruled out as meaningless) the best-fit parameters, some goodness-of-fit statistics as well as confidence intervals for the parameters are necessary.
If the computation of the cost function is efficient (it should be) then, the added computational cost in the way of calculating CIs is absolutely worth it.
Is this a feature that makes sense to add? Is it requested by people?
I imagine it would take Optimization.jl closer in utility to scipy.optimize’s functions.
Instead of a flag, maybe a post-processing function?
I don’t think most people would want this confidence interval since it’s not generally correct for nonlinear problems. It’s only correct in the linear approximation so it’s making quite a bit of assumptions in the nonlinear. case. Using a Bayesian approach like Turing.jl would generally give more accurate and interpretable confidence intervals so that would likely be recommended in most cases which are nonlinear.
Yeah, for nonlinear problems, you need to actually have found the global optimum for the estimate to make any sense at all, and if your residuals are non-Gaussian, it still makes no sense.
True, but in the case where error in user’s data is assumed as Gaussian and uncorrelated then it is a fair approximation to make. Indeed this is the assumption lots of optimization libraries make in order to approximate a confidence interval, this is not the case for Optimization.jl as of now.
My above comment applies here as well, under assumptions about the nature of the uncertainty in the data. If the errors are a-priori known to be non-Gaussian, other optimization route is of course necessary.
My main stance is that I think that Optimization.jl would strongly benefit from offering a comparable utility to general purpose fitting programs such as e.g. gnuplot, lmfit in python or even classic fitting frameworks such as Minuit without needing to resort to pycall or C++ bindings.
@ChrisRackauckas you mentioned a post-processing function to estimate confidence-intervals around a found local minimum. Do you see something like that ever making it in Optimization.jl for the rough but albeit common gaussian-error approximation?
Yes, if someone makes a PR I’d accept it. I think it’s a feature worth having and documenting well, though it has enough caveats that I personally won’t be prioritizing implementing it.
What other optimization libraries do this? Optimization.jl isn’t a statistics package.
I don’t think it would be appropriate to have a post-processing function like you recommend in Optimization.jl since it’s easily misused by people (e.g. I’ve seen to many +/- 2se type intervals for cases where it makes no sense…) / it doesn’t make sense for the majority of inputs. Wouldn’t it be better off in some other statistics package? Especially since (last I checked?) Optimization.jl doesn’t really have any post-processing functions in it.
I agree; Optimization.jl is neither a statistics package nor a plotting package and should do one thing (optimization) and do it well.
That said, I am surprised that OP’s use case isn’t covered by one of the existing nonlinear regression tools elsewhere in the Julia ecosystem. @abehersan , perhaps you can post a minimal version of your model and someone may be able to point you in the right direction.
My use case is extremely simple. Peak lineshape fitting to experimental data. Error bars on the data are known. A minimal example of a simple peak fit is shown below when working with python and lmfit:
"""
lorentzian.py
Generate noisy data that is Lorentzian in nature, fit a Lorentzian to the
data and report parameter uncertainties using LMFIT
"""
import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model
np.random.seed(777)
def lorentzian(x, a, c, s):
L = (a / np.pi) * (s / ((x - c)**2 + s**2))
return L
"""true parameters"""
a = 10.0
c = 0.0
s = 0.02
x_data = np.linspace(-0.5, 0.5, 101)
noise = np.random.normal(0, 7, len(x_data))
y_data = lorentzian(x_data, a, c, s) + noise
dy_data = np.sqrt(y_data - noise) # assumed poissonian noise
data_xye = np.array([x_data, y_data, dy_data]).T
np.savetxt("lorentzian_peak.dat", data_xye, delimiter="\t")
"""fit with lmfit"""
model = Model(lorentzian)
model.nan_policy = "omit"
params = model.make_params(a=17.0, c=0.0, s=0.04)
fit_result = model.fit(y_data, params, x=x_data, weights=1.0/dy_data)
fit_ax = fit_result.plot_fit(numpoints=701, fit_kws={"lw":3})
fit_fig = plt.gcf()
fit_fig.suptitle("Lorentzian fit with LMFIT")
fit_fig.set_size_inches(11.7, 8.3)
"""fit statistics"""
fit_stats = fit_result.fit_report()
fit_ax.text(0.05, 0.5, fit_stats, transform=fit_ax.transAxes)
plt.tight_layout()
plt.savefig("lmfit_fit.pdf", dpi=300)
plt.show()
An analogous work flow in Julia is as follows:
"""
lorentzian.jl
Fit a Lorentzian using Optimization.jl
"""
using Optimization, OptimizationNLopt
using Plots
default(size=(1080, 720))
default(margin=5Plots.mm)
default(linewidth=3.5)
function lorentzian(x::Real, a::Real, c::Real, s::Real)::Float64
(a / pi) * (s / ((x - c)^2 + s^2))
end
function chisqr_lorentzian(u, p)::Float64
x_data = p[1]
y_data = p[2]
dy_data = p[3]
y_calc = [lorentzian(x, u...) for x in x_data]
sum( ((y_data .- y_calc) ./ dy_data) .^ 2 )
end
"""load data generated with python"""
dfile = open("lorentzian_peak.dat", "r")
x_data = Float64[]
y_data = Float64[]
dy_data = Float64[]
for line in eachline(dfile)
parts = split(line)
push!(x_data, parse(Float64, parts[1]))
push!(y_data, parse(Float64, parts[2]))
push!(dy_data, parse(Float64, parts[3]))
end
u0 = Float64[17, 0, 0.04]
data_xye = [x_data, y_data, dy_data]
optfunc_lor = OptimizationFunction(chisqr_lorentzian)
optprob_lor = OptimizationProblem(optfunc_lor, u0, data_xye)
sol = solve(optprob_lor, NLopt.LN_NELDERMEAD())
s_names = ["amplitude", "center", "sigma"]
sol_str = [n*": "*string(u) for (n, u) in zip(s_names, sol)]
push!(sol_str, "FWHM: $(2*sol[3])")
p = plot(xlabel="x", ylabel="y", title="Lorentzian fit with Optimization.jl")
scatter!(p, x_data, y_data, yerr=dy_data, label="data")
plot!(p, x_data, [lorentzian(x, sol...) for x in x_data], label="fit")
annotate!(p, (-0.3, 100, join(sol_str, "\n")), annotatiohnalign=:left)
savefig(p, "optimization_fit.pdf")
The fit results naturally agree. Of course the added benefit of lmfit is the calculated CI of the fit parameters from scipy´s covariance matrix computation.
For my use case, in principle the Hessian can be calculated since for most lineshapes I work with, analytical expressions for the derivative of the model w.r.t. the fit parameters are straightforward to calculate.
I understand Optimization.jl is not a statistics package per se, but it strikes me as odd that it optimizes very well but then it does not provide confidence interval for the optimized parameters (even in the Gaussian approximation) nor any goodness-of-fit metric beyond the value of the objective function at the found minimum.
I also understand that for more complex cases pushing a “+/-” on the users’ fit results willy-nilly is dangerous. However I would contend that most use cases are not all use cases.
This exact surprise is what I felt and what prompted me to start this conversation topic in the first place. And while LsqFit is a counter-example, it’s still quite disheartening to have to give up the whole amazing machinery of Optimization.jl just because you want an error bar on your parameters post-fit.
The terms “confidence interval” and “goodness of fit” are quite foreign to the domain of optimization. Optimization solvers solve math problems of the form minimize f(x) subject to g(x) <= 0. You either find the answer or you don’t.
It would be weird if Julia returned a confidence interval on sum(x) just because it is “common” for x to be a vector of observations from a statistical distribution, right? 
If you really wanted that, you would do using StatsBase; mean_and_std(x) .* length(x) or whatever.
While there would be utility in a separate library that does nonlinear regression using Optimization.jl on the back end and returns the covariance matrices you need—indeed, Julia’s composability means it shouldn’t be terribly hard for someone to build such a library—as someone who solves optimization problems every day, I would much prefer Optimization.jl be kept at a more limited feature and dependency set.
This wouldn’t require any dependencies. And no, having such functionality would be inline with other parts of the SciML interface. DifferentialEquations has nice solution interface features that people like to use:
Etc. LinearSolve.jl, NonlinearSolve.jl, Integrals.jl, Optimization.jl, etc. should all follow suit. And I think a feature that allows one to do a post-analysis that adds uncertainties to the parameter results under some (well-documented) assumptions is precisely the kind of post-solution analysis functionality that would be good to offer, along with nice pre-defined logging functions, benchmarking tooling, etc.
Thanks for considering the utility of such post-processing!
I think the challenging part of this extension to Optimization.jl is that different solvers return a solution struct with different fields and is also not type-inferrable, see e.g. `solve` not type-inferrable · Issue #556 · SciML/Optimization.jl · GitHub
I would be glad to contribute if some pointers are hinted at since I’m not familiar with the whole of SciML suite, maybe with “good first-issue” tags in the repo’s issues(?
Follow up with @Vaibhavdixit02. I recommend joining the Julia slack and the #sciml-bridged channel where a lot of the dev chat takes place. @Vaibhavdixit02 can probably spend a bit of time today tagging some good first issues and helping get you up to speed. There’s a lot of good first issues in that repo.
I think better error throwing on things not supported by a specific algorithm would be a good one. Or another would be to contribute to the benchmark system. We just got Optimization.jl benchmarks up:
https://docs.sciml.ai/SciMLBenchmarksOutput/stable/Optimizaton/2drosenbrock/
but it’s really jank right now. Cleaning that up, but also systematizing Get BlackBoxOptimizationBenchmarking.jl updated to Optimization.jl's interface and run on SciMLBenchmarks · Issue #640 · SciML/SciMLBenchmarks.jl · GitHub and other optimization benchmarks would be good.
