# Searching for a simple polynomial curve fitting example

**URL:** <https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [August 9, 2017, 12:50pm UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302 "2017-08-09T12:50:08Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![J\_T](https://avatars.discourse-cdn.com/v4/letter/j/57b2e6/32.png) [@J\_T](https://discourse.julialang.org/u/J_T)\
**Post date:** [August 9, 2017, 12:50pm UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302/1 "2017-08-09T12:50:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [August 9, 2017, 3:12pm UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302/2 "2017-08-09T15:12:49Z")

</div>

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

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

```

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

```julia
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](https://github.com/JuliaOpt/JuMP.jl/tree/master/examples) and [JuMP — Julia for Mathematical Optimization — JuMP -- Julia for Mathematical Optimization 0.18 documentation](http://www.juliaopt.org/JuMP.jl/0.18/).

---

<div class="post-metadata">

**Author:** ![J\_T](https://avatars.discourse-cdn.com/v4/letter/j/57b2e6/32.png) [@J\_T](https://discourse.julialang.org/u/J_T)\
**Post date:** [August 10, 2017, 6:53am UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302/3 "2017-08-10T06:53:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Humphrey\_Lee](https://avatars.discourse-cdn.com/v4/letter/h/a88e4f/32.png) [@Humphrey\_Lee](https://discourse.julialang.org/u/Humphrey_Lee)\
**Post date:** [April 8, 2020, 5:52am UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302/4 "2020-04-08T05:52:28Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![J\_T](https://avatars.discourse-cdn.com/v4/letter/j/57b2e6/32.png) [@J\_T](https://discourse.julialang.org/u/J_T)\
**Post date:** [March 18, 2021, 2:03pm UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302/5 "2021-03-18T14:03:55Z")

</div>

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

```julia
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.

```julia
( [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

```julia
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

```julia
( [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.]

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 18, 2021, 7:52pm UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302/6 "2021-03-18T19:52:09Z")

</div>

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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![J\_T](https://avatars.discourse-cdn.com/v4/letter/j/57b2e6/32.png) [@J\_T](https://discourse.julialang.org/u/J_T)\
**Post date:** [March 19, 2021, 5:10am UTC](https://discourse.julialang.org/t/searching-for-a-simple-polynomial-curve-fitting-example/5302/7 "2021-03-19T05:10:51Z")

</div>

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](https://discourse.julialang.org/t/fixing-vs-constraining-a-variable-in-jump/27609/6)

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

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