# Searching for a simple polynomial curve fitting example

Im new to mathematical optimisation, so I’m looking for a simple example I can learn from …

My (initial) problem looks like
1000 cartesian points in Float16
I want to find a polynomial p(x) (to be evaluated in Float16) to fit these points.

Conceptually
my constraints are |p(x) - y| <= \text{ulp}(y) and
the objective is \text{Min} of \text{norm(|p(x) - y|, Inf)}

But Im not seeing how to declare it within JuMP and I haven’t found a suitable example to learn from.

Any pointers ?

JT

Edit: I corrected the constraints refer to \text{ulp}(y) instead of x.

You can use linear programming if you transform your model as follows:

Minimize t
s.t.
-t <= p(x) - y <= t for each (x,y)


It’s an LP in the parameters of p(x).

using JuMP
using GLPKMathProgInterface
#using Clp

function fit_poly_inf_norm(x,y,n)
m = Model(solver=GLPKSolverLP())
#m = Model(solver=ClpSolver())
@variable(m, alpha[1:n+1])
@variable(m, t >= 0.)
@constraint(m, [i=1:length(y)], sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] <= t);
@constraint(m, [i=1:length(y)], -t <= sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i]);
@objective(m, Min, t)
solve(m)
return getvalue(alpha)
end

x = rand(100)
y = rand(100)
n = 4 #order
coeffs = fit_poly_inf_norm(x,y,n)


There are plenty of examples in https://github.com/JuliaOpt/JuMP.jl/tree/master/examples and http://www.juliaopt.org/JuMP.jl/0.18/.

1 Like

Thanks Mohamed.

I had found (many of) the examples, but failed to recognise my problem in the Julia implementation.
But I think I can easily extend this to solve my problem.

thanks again,
JT

I’ve updated the code to make it work with Julia ver. 1.4, JuMP ver. 0.21 and GLPK ver. 0.12. I always found good example but need to be updated with updated Julia and its packages.

using JuMP
#using GLPKMathProgInterface # original code
using GLPK # added new 2020-04-08
#using Clp

using Plots
using DataFrames

function fit_poly_inf_norm(x,y,n)
#m = Model(solver=GLPKSolverLP()) # original code
m = Model(GLPK.Optimizer) # added new 2020-04-08
#m = Model(solver=ClpSolver())
@variable(m, alpha[1:n+1])
@variable(m, t >= 0.)
@constraint(m, [i=1:length(y)], sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] <= t);
@constraint(m, [i=1:length(y)], -t <= sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i]);
@objective(m, Min, t)
#solve(m) #original code
#return getvalue(alpha) #original code
return JuMP.value.(alpha)
end

x = rand(100)
y = rand(100)
n = 4 #order
coeffs = fit_poly_inf_norm(x,y,n)

println(coeffs)

df = DataFrame()
df[!, :Y] = y
df[!, :X] = x

Geom.point,
)

Here is a modified example to fit a curve, specifically log on the domain [ 3/4, 3/2 ].

This will use a degree 9 polynomial where the first two coefficients are explicitly set to { 0, 1 } respectively and the simple transform f=x-1 is used because of the root at x=1.

using JuMP, GLPK

# simple domain
dom = ( 0.75, 1.5 )

xvals = collect(range( dom, dom; length=101 ));

yvals = log.( xvals );

function fit_poly_inf_norm(x,y,n)
m = Model(GLPK.Optimizer)
@variable(m, alpha[1:n+1])
@variable(m, t >= 0.)
@constraint(m, alpha == 0 );
@constraint(m, alpha == 1 );
@constraint(m, [i=1:length(y)], sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] <= t );
@constraint(m, [i=1:length(y)], -t <= sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] );
@objective(m, Min, t)
optimize!(m)
coeffs = JuMP.value.(alpha)
err = JuMP.value.(t)
return (coeffs, err)
end

degree = 9
# root at x=1 so approximate with f=x-1
(coeffs,err) = fit_poly_inf_norm( xvals .- one(dom), yvals, degree )


which gives results.

( [0.0, 1.0, -0.5000011437016926, 0.33333600206898395, -0.2499042550870806, 0.19972612897284364, -0.16863298896795648, 0.1506453512060692, -0.1197307463685523, 0.05204078783867752], -2.1019126300592106e-9 )


Note that the coefficients (after the zero) have alternating sign.
Its my understanding that adding constraints to force this should be redundant, but the following

function fit_poly_inf_norm2(x,y,n)
m = Model(GLPK.Optimizer)
@variable(m, alpha[1:n+1])
@variable(m, t >= 0.)
@constraint(m, alpha == 0 );
@constraint(m, alpha == 1 );

@constraint(m, alpha <= 0 );
@constraint(m, 0 <= alpha );
@constraint(m, alpha <= 0 );
@constraint(m, 0 <= alpha );
@constraint(m, alpha <= 0 );
@constraint(m, 0 <= alpha );
@constraint(m, alpha <= 0 );
@constraint(m, 0 <= alpha );

@constraint(m, [i=1:length(y)], sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] <= t );
@constraint(m, [i=1:length(y)], -t <= sum(alpha[j]*x[i]^(j-1) for j in 1:n+1) - y[i] );
@objective(m, Min, t)
optimize!(m)
coeffs = JuMP.value.(alpha)
err = JuMP.value.(t)
return (coeffs, err)
end

(coeffs2,err2) = fit_poly_inf_norm2( xvals .- one(dom), yvals, degree )


give the different results

( [0.0, 1.0, -0.5000012569584985, 0.33333324995783115, -0.24988908708604318, 0.19981205548969036, -0.16904719456430753, 0.15022553801306932, -0.11640673034878432, 0.04871565717295707], 1.2662657797313237e-9 )


Can anybody explain why the results change when adding redundant constraints ?
[ Julia 1.5.3 and 1.6.0-rc3 give same results. ]

The results are the same to within the tolerance of the solver. (~3.5e9, where as the default tolerance is ~1e-8).

Your second formulation adds extra linear constraints. You may get a better solution by modifying the variable bounds instead:

fix(alpha, 0.0)
fix(alpha, 1.0)
for a = 3:10
if isodd(a)
set_upper_bound(alpha[a], 0.0)
else
set_lower_bound(alpha[a], 0.0)
end
end


Thank you!

Can I ask for learning purposes what is the practical difference between fixing coefficients (or setting bounds) and defining constraints ?

set_optimizer_attribute(m, "tol_bnd",  1e-11 )

before the results stabilised and the absolute errors looked suitably like Chebychev equi-oscillation (with obvious root at x=1).