Find parameters that minimize the distance between two curves in terms of the infinite norm

Question

I am trying to implement a very simple NMinimize code that would search for parameters of a polynomial that minimize (globally) the distance between this polynomial and one given curve. The distance that I need to use is the uniform distance (i.e. $\max\limits_{x\in [0,1]} |f(x)-g(x)|$) , or infinite norm.

I do not want to use the FindFit command as it has less flexibility and also globality of the found solution is not ensured (especially for more involved problems).

So I have constructed the following code:

v[x_] := ChebyshevT[6, x] (* the given function to be approximated *)

f[x_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
 a x^2 + b x + c   (* funcn whose parameters I want to find *)

abs[a_?NumericQ, b_?NumericQ, c_?NumericQ, x_?NumericQ] := 
 Abs[v[x] - f[x, a, b, c]]   

maxabs[a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
NMaxValue[abs[a, b, c, y], {y, 0, 1}]

n = NMinimize[{maxabs[a, b, c]}, {a, b, c}, Method -> "NelderMead"]

Plot[{v @ y, f[a, b, c, y] /. n[[2]]}, {y, 0, 1}, 
PlotStyle -> {{Red}, {Dashed, Blue}}]

However this gives me a fit which is obviously suboptimal, (a straight line, while even the naked eye suggests that a simple parabola would be a way better fit)... Did I do anything wrong with the code?

Answer 1

The following is fast and suggests an almost straight and horizontal line:

p = Range[0, 1, 1/100];
v[x_] := v[x] = ChebyshevT[6, x]
f[x_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ] := a x x + b x + c
abs[x_?NumericQ, a_?NumericQ, b_?NumericQ,  c_?NumericQ] := (v[x] - f[x, a, b, c])^2
maxabs[a_?NumericQ, b_?NumericQ, c_?NumericQ] := Max[abs[#, a, b, c] & /@ p]

n = Monitor[NMinimize[maxabs[a, b, c], {a, b, c}], {a, b, c}]

(* {1.00244, {a -> 0.00419288, b -> 0.00395157, c -> -0.00424115}}*)

While this "analytical" solution gets the same result (only more horizontal and straighter):

es[a_, b_, c_, x_] = (ChebyshevT[6, x] - (a x x + b x + c))^2 //  Expand;
extrema[a_, b_, c_] := Join[{0, 1}, 
                       Select[x /. Solve[D[es[a, b, c, x], x] == 0, x] // N, 
                              Head[#] =!= Complex && 0 < # < 1 &]]

maxabs[a_?NumericQ, b_?NumericQ, c_?NumericQ]:=Max[es[a, b, c, #] & /@ extrema[a, b, c]]

sol = NMinimize[maxabs[a, b, c], {a, b, c}, Method -> "NelderMead"]

(* {1.00014, {a -> -0.00256894, b -> 0.00390044, c -> -0.00137957}} *)

Answer 2

My approach.. first do a least squares fit, which gives a global minimum, although with a different error measure:

f[x_, a_, b_, c_] := a (x)^2 + b (x) + c;
s1 = First@
  Solve[(D[ 
        Simplify[
         Total[((f[#, a, b, c] - v[#])^2 & /@ 
            Range[0, 1, .001])]] , #] & /@ {a, b, c, d}) == 0, {a, b, 
    c}]
Plot[{v[y], f[y, a, b, c] /. s1 }, {y, 0, 1}]

Then use that solution as a start point for FindMinimum:

crit[a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
     Norm[f[#, a, b, c] - v[#] & /@ Range[0, 1, .001], Infinity]
{a0, b0, c0} = {a, b, c} /. s1;
s2 = Last@FindMinimum[ crit[a, b, c] , {{a, a0}, {b, b0}, {c, c0}}]
ep = First@
  MaximalBy[Range[0, 1, .2], ({Abs[v[#] - f[#, a, b, c] /. s2]}) &]
Plot[{v[y], f[y, a, b, c] /. s1, f[y, a, b, c] /. s2 }, {y, 0, 1},
 Epilog -> Arrow[{{ep, v[ep]}, {ep, f[ep, a, b, c] /. s2 }}]]

The max error is 1.3822 @ x=1 (Which we can see see straight away is not a global minimum since the line y=0 has a max error of 1. )

Answer 3

As I mentioned in a comment the original code has a mistake in the plotting commands -- it should be used f[y,a,b,c], not f[a,b,c,x].

I found that mistake using the solution by "belisarius has settled" as a base with the following code, which is specially made to get a parabola that fits the curve.

First we select sampling points close to the characteristic points of the curve:

vps = N@Select[
    Flatten[Solve[D[ChebyshevT[6, x], x] == 0, x]][[All, 2]], 
    0 < # < 0.7 &];
p = Flatten@Map[# + Range[-1, 1, 0.1]/20 &, vps] ;

Now using those points p we define the minimization functions and then run the minimization itself:

ClearAll[v, f, abs, maxabs]
v[x_] := v[x] = ChebyshevT[6, x]
f[x_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ, d_?NumericQ] := 
  a (x - d)^2 + b (x - d) + c;
abs[x_?NumericQ, a_?NumericQ, b_?NumericQ, c_?NumericQ, 
   d_?NumericQ] := ((v[x] - f[x, a, b, c, d])/v[x])^2;
maxabs[a_?NumericQ, b_?NumericQ, c_?NumericQ, d_?NumericQ] := 
 maxabs[a, b, c, d] = Max[abs[#, a, b, c, d] & /@ p]
sol = Monitor[
  NMinimize[{maxabs[a, b, c, d], a < 0, 0.45 < d < 0.55}, {a, b, c, 
    d}], {a, b, c, d}]

(* {4.77768*10^-6, {a -> -22.2079, b -> 1.53664, c -> 0.971457, 
  d -> 0.464512}} *)

Here is the plot with the result:

Remarks:

1. The plot grid lines show the points over which the minimization is done.

2. I have provided some constraints to guide the minimization process.

3. I changed the minimization function to use relative error not absolute error, but in this case with both types we get very similar results.