8
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – MarcoB
    Nov 13, 2015 at 15:13
  • 1
    $\begingroup$ There is a mistake in the plotting command. It should be f[y,a,b,c] not f[a,b,c,y]. $\endgroup$ Nov 13, 2015 at 17:59
  • 1
    $\begingroup$ 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. $\endgroup$
    – user484
    Nov 13, 2015 at 19:23
  • 1
    $\begingroup$ 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). $\endgroup$
    – Kass
    Nov 13, 2015 at 21:56
  • 1
    $\begingroup$ In case you want to minimize the relative error, there is a built-in function for this: MiniMaxApproximation. $\endgroup$
    – Silvia
    Nov 15, 2015 at 13:56

3 Answers 3

9
$\begingroup$

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}} *)
$\endgroup$
4
$\begingroup$

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}]

enter image description here

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 }}]]

enter image description here

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

$\endgroup$
6
  • $\begingroup$ And what is the result (the norm) . I can't run it on v9! .Can you compare it with mine, please? $\endgroup$ Nov 13, 2015 at 19:04
  • 1
    $\begingroup$ 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}]. $\endgroup$
    – MarcoB
    Nov 13, 2015 at 20:20
  • $\begingroup$ 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.... $\endgroup$
    – Kass
    Nov 13, 2015 at 20:50
  • $\begingroup$ 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. $\endgroup$
    – george2079
    Nov 13, 2015 at 21:20
  • $\begingroup$ 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: $\endgroup$
    – Kass
    Nov 13, 2015 at 21:53
1
$\begingroup$

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:

enter image description here

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.

$\endgroup$
6
  • 1
    $\begingroup$ isn't the fit supposed to be over 0-1? $\endgroup$
    – george2079
    Nov 13, 2015 at 18:57
  • $\begingroup$ @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. $\endgroup$ Nov 13, 2015 at 19:00
  • $\begingroup$ @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. $\endgroup$ Nov 13, 2015 at 19:06
  • $\begingroup$ @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? $\endgroup$ Nov 13, 2015 at 19:13
  • $\begingroup$ @anton I don't think it is that important, and I don't believe it deserves downvotes. Deleting my comments! $\endgroup$ Nov 13, 2015 at 22:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.