Fitting two curves with linear/non-linear regression

Question

I need to fit two curves(which both should belong to cubic functions) into a set of points with JuMP.

I've done fitting one curve, but I'm struggling at fitting 2 curves into same dataset.

I thought that if I can distribute points to curves - so if each point can only be used once - I can do it like below, but it didn't work. (I know that I can use much more complicated things, I want to keep it simple.)

This is a part of my current code:

# cubicFunc is a two dimensional array which accepts cubicFunc[x,degree]

@variable(m, mult1[1:4]) // 0:3 because it's cubic
@variable(m, mult2[1:4]) // 0:3 because it's cubic

@variable(m, 0 <= includeIn1[1:numOfPoints] <= 1, Int)
@variable(m, 0 <= includeIn2[1:numOfPoints] <= 1, Int)

# some kind of hack to force one of them to 0 and other one to 1
@constraint(m, loop[i in 1:numOfPoints], includeIn1[i] + includeIn2[i] == 1)

@objective(m, Min, sum( (yPoints - cubicFunc*mult1).*includeIn1 .^2 ) + sum( (yPoints - cubicFunc*mult2).*includeIn2 .^2 ))

But it gives various errors depending on what I'm trying; *includeIn1 and, .*includeIn1 doesn't work, I've tried to do it via @NLobjective but it gave me whooping ~50 lines of errors etc.

Is my idea realistic? Can I make it into the code?

Any help will be highly appreciated. Thank you very much.

Answer 1

You can write down the problem e.g. like this:

using JuMP, Ipopt

m = Model(with_optimizer(Ipopt.Optimizer))

@variable(m, mult1[1:4])
@variable(m, mult2[1:4])
@variable(m, 0 <= includeIn1[1:numOfPoints] <= 1)
@variable(m, 0 <= includeIn2[1:numOfPoints] <= 1)

@NLconstraint(m, loop[i in 1:numOfPoints], includeIn1[i] + includeIn2[i] == 1)

@NLobjective(m, Min, sum(includeIn1[i] * (yPoints[i] - sum(cubicFunc[i,j]*mult1[j] for j in 1:4)) ^2 for i in 1:numOfPoints) +
                     sum(includeIn2[i] * (yPoints[i] - sum(cubicFunc[i,j]*mult2[j] for j in 1:4)) ^2 for i in 1:numOfPoints))

optimize!(m)

Given the constraints includeIn1 and includeIn2 will be 1 or 0 in optimum (if they are not this means that it does not matter to which group you assign the point), so we do not have to constrain them to be binary. Also I use non-linear solver as the problem does not not seem to be possible to reformulate as linear or quadratic optimization task.

However, I give the above code only as an example how you can write it down. The task you have formulated does not have a unique local minimum (that is a global one then), but several local minima. Therefore using standard non-linear convex solvers that JuMP supports will only find one local optimum (not necessarily a global one). In order to look for global optima you need to switch to global solvers like e.g. https://github.com/robertfeldt/BlackBoxOptim.jl.