So, we want to find a polynomial P with integer coefficients and a fixed degree
that minimizes max(abs(P(point)-value)
where the max is for all (point,value)
in mydata
.
With such a formulation, we see that there is a possibility with integer programming.
But the abs()
function is not linear. So we add a error
variable (to be minimized) such that both linear inequalities P(point)-value <= error
and value-P(point) <= error
hold.
Putting all together, we get:
sage: degree = 3
sage: mydata = [(0,0),(1,5),(2,1),(4,6)]
sage: p = MixedIntegerLinearProgram(maximization=False)
sage: a = p.new_variable(integer=True)
sage: error = p.new_variable(real=True)
sage: p.set_objective(error[0])
sage: for point,value in mydata:
....: p.add_constraint(sum(a[d]*point^d for d in range(degree+1)) - value <= error[0])
....: p.add_constraint(value - sum(a[d]*point^d for d in range(degree+1)) <= error[0])
sage: error = p.solve()
sage: dico = p.get_values(a)
sage: R = PolynomialRing(ZZ, 'x')
sage: P = sum(R(v*x^d) for d,v in dico.items())
sage: P
x + 2
sage: error
3.0
If you want to find the best polynomial with real entries, you just have to replace the two lines
sage: a = p.new_variable(integer=True)
sage: R = PolynomialRing(ZZ, 'x')
by
sage: a = p.new_variable(real=True)
sage: R = PolynomialRing(RR, 'x')
In our example, you will get :
sage: P
0.500000000000000*x + 2.25000000000000
sage: error
2.25
For huge data, it should also be possible to approximate the best solution with lattice reduction algorithms.