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$
  • $\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 '15 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$ – Anton Antonov Nov 13 '15 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$ – Rahul Nov 13 '15 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 '15 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 '15 at 13:56
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$
  • $\begingroup$ And what is the result (the norm) . I can't run it on v9! .Can you compare it with mine, please? $\endgroup$ – Dr. belisarius Nov 13 '15 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 '15 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 '15 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 '15 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 '15 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$
  • 1
    $\begingroup$ isn't the fit supposed to be over 0-1? $\endgroup$ – george2079 Nov 13 '15 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$ – Anton Antonov Nov 13 '15 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$ – Anton Antonov Nov 13 '15 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$ – Anton Antonov Nov 13 '15 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$ – Dr. belisarius Nov 13 '15 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.