LsqFit curve fit does not consider bounds on intercepts correctly

I am trying to get some fits with LsqFit with bounds, but I am observing strange results which seem to imply that LsqFit is minimizing the correlation coefficient but not the residue of the fit. This is a minimal example, with a linear fit, in which I add a bound to the intercept coefficient:

using LsqFit
x = collect(1:10); y = copy(x);
@. model(x,p) = p[1]*x + p[2]
p0 = [ 0. , -5. ];
fit = curve_fit(model,x,y,p0); # no bounds
fit_bounded = curve_fit(model,x,y,upper=[+Inf,-5.]); # with bounds
scatter(x,y)
plot!(x,model(x,fit.param),label="no bounds")
plot!(x,model(x,fit_bounded.param),label="with bounds")

Result:

fit

Of course that is not the line with the minimal residue sum, it is just the same fit, but shifted downwards.

Am I doing or interpreting something wrong?

You’re constraining the parameter vector to not allow an intercept greater than -5,so how could you possibly get anything else? Did you mean to wrote +5?

I was expecting a greater slope, such that the fit cross the data, that would have a smaller residue relative to the data, while satisfying the bound.

Edit: which can be demonstrated if the model is defined with that parameter equal to -5 from start:

julia> @. model(x,p) = p[1]*x - 5.
model (generic function with 1 method)

julia> x = collect(1:10); y = x;

julia> fit = curve_fit(model, x, y, p0);

julia> scatter(x,y)

julia> plot!(x,model(x,fit.param))


Results in (the expected result):

fit

Ah, good point! The gradient for the slope parameter is probably zero in the point in the figure?

No, it is not. The result is independent on the initial point as well.

Also, this kind of problem occurs with any model, for example:

julia> x = collect(1:10); y = x.^2; scatter(x,y)

julia> @. model(x,p) = p[1]*x^2 + p[2]
model (generic function with 1 method)

julia> p0 = [ 5. , -20. ];

julia> fit = curve_fit(model, x, y, p0);

julia> plot!(x,model(x,fit.param),label="no bounds")

julia> fit = curve_fit(model, x, y, p0, upper=[+Inf,-20.]);

julia> plot!(x,model(x,fit.param),label="with bounds")

Result:

fit

Edit: I have opened an issue: https://github.com/JuliaNLSolvers/LsqFit.jl/issues/167

If that is the expected result for some reason, please let me know.

1 Like