Revision history [back]

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.