Searching for a simple polynomial curve fitting example

Specific Domains Optimization (Mathematical)
1

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 ?

thanks in advance,
JT

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

2

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 JuMP — Julia for Mathematical Optimization — JuMP -- Julia for Mathematical Optimization 0.18 documentation.

1 Like
3

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

4

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 Gadfly
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
	optimize!(m) # added new 2020-04-08
	#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

Gadfly.plot(df, x="X", y="Y",
	Geom.point,
)
5

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
using Gadfly

# simple domain
dom = ( 0.75, 1.5 )

xvals = collect(range( dom[1], dom[2]; 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[1] == 0 );
	@constraint(m, alpha[2] == 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[1]), 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[1] == 0 );
	@constraint(m, alpha[2] == 1 );

	@constraint(m, alpha[3] <= 0 );
	@constraint(m, 0 <= alpha[4] );
	@constraint(m, alpha[5] <= 0 );
	@constraint(m, 0 <= alpha[6] );
 	@constraint(m, alpha[7] <= 0 );
	@constraint(m, 0 <= alpha[8] );
	@constraint(m, alpha[9] <= 0 );
	@constraint(m, 0 <= alpha[10] );

	@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[1]), 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. ]

6

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[1], 0.0)
fix(alpha[2], 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
7

Thank you!

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

Edit: I found the answer in another thread: fixing-vs-constraining-a-variable-in-jump

P.S. For others reference, I ended up adding

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).