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

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?

• I tried running your code and I get a lot of warnings because NMinimize failed to converge to a solution. Do you not get the same? – MarcoB Nov 13 '15 at 15:13
• There is a mistake in the plotting command. It should be f[y,a,b,c] not f[a,b,c,y]. – Anton Antonov Nov 13 '15 at 17:59
• Belisarius is correct, the optimal quadratic fit is the constant zero function, whose infinity-norm distance from $v(x)$ is 1. For a function $f(x)$ to have distance less than that, it would have to be negative at $x=0$, positive at $x=1/2$, negative at $x=\sqrt3/2$, and positive at $x=1$ (because those are the extrema of $v(x)$ on $[0,1]$: i.stack.imgur.com/zLP5g.png). This is impossible if $f(x)$ is quadratic. – Rahul Nov 13 '15 at 19:23
• Thanks to Anton for spotting the mistake! However the comment of MarcoB remains valid too -- I was also getting the warnings of "failure to converge to a solution" (independently of plotting mistake, of course). – Kass Nov 13 '15 at 21:56
• In case you want to minimize the relative error, there is a built-in function for this: MiniMaxApproximation. – Silvia Nov 15 '15 at 13:56

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}} *)


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

• And what is the result (the norm) . I can't run it on v9! .Can you compare it with mine, please? – Dr. belisarius Nov 13 '15 at 19:04
• Interesting approach. Note that you could also find the least squares approximation symbolically (and possibly faster) by minimizing the square of the area between the curves, e.g. Minimize[Integrate[(ChebyshevT[6, x] - a x^2 - b x - c)^2, {x, 0, 1}], {a, b, c}]. – MarcoB Nov 13 '15 at 20:20
• I agree with Marco, that one can do it directly with 'Minimize' of the 'NIntegrate' of the squared difference between functions and it works for all $x \in [0,1]$. But it is baffling that one cannot go exactly the same way and use NMinimize to find parameters given the infinite norm.... – Kass Nov 13 '15 at 20:50
• doh! I didn't think to check that that integral could be done symbolically. In any case its a nonlinear problem with numerous local minima. I wouldn't call it baffling that NMinimize fails to find the global. – george2079 Nov 13 '15 at 21:20
• I can use this code, which seems to provide a better solution, but with minimization of the sum of squared distance between two curves. I was looking for finding the version of the code where the distance is similar to the infinite norm: – Kass Nov 13 '15 at 21:53

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.

• isn't the fit supposed to be over 0-1? – george2079 Nov 13 '15 at 18:57
• @george2079 The range is not specified in the question. I answered what is wrong with the code, and how to get a fitting parabola close to the given function. – Anton Antonov Nov 13 '15 at 19:00
• @belisariushassettled Right, the range is specified. Still the question asked what is wrong with the original code, and I found that by using the solution I posted. – Anton Antonov Nov 13 '15 at 19:06
• @belisariushassettled I am a little surprised by your comments. I used your solution to find what is wrong with the plotting code in the original post and your solution. The original post does finish with the question "Do I do anything wrong with the code? – Anton Antonov Nov 13 '15 at 19:13
• @anton I don't think it is that important, and I don't believe it deserves downvotes. Deleting my comments! – Dr. belisarius Nov 13 '15 at 22:08