Why does Maximize return a wrong answer?

Question

I want to find the maximum of the following function, but get a wrong result. What am I doing wrong?

Clear[x, y, xx, yy]
peaks[x_, y_] := 
  3*(1 - x)^2.*Exp[-(x^2) - (y + 1)^2] - 
   10*(x/5 - x^3 - y^5) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2];
N[Maximize[peaks[xx, yy], {xx, yy}]]
NMaximize[{peaks[x, y], x >= -2, x <= 2, y >= -2, y <= 3}, {x, y}]
peaks[0, 2]
peaks[1.28568, -0.00484756]

When you look at these graphs, you realize that the results are strange:

peaks[x_, y_] := 3*(1 - x)^2*Exp[-(x^2) - (y + 1)^2] - 
                 10*(x/5 - x^3 - y^5) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2];
Plot3D[peaks[x, y], {x, -2, 2}, {y, -2, 3}, 
    ColorFunction -> "DarkRainbow", 
    AxesLabel -> Automatic]
ContourPlot[peaks[x, y], {x, -2, 2}, {y, -2, 3},  
    ColorFunction -> "DarkRainbow", 
    ContourLines -> False,  
    Contours -> 25]

Mathematica graphics

More on the related problems you'll find examining these two posts: Prove (or check) an expression is positive given constraints on variables? and about the Maximize issue: How do I determine the maximum value for a polynomial, given a range of x values? — Artes, Jun 14 '14 at 18:26

score 7 · Answer 1 · answered Jun 14 '14 at 18:01

From the documentation for Maximize: "If f and cons are linear or polynomial, Maximize will always find a global maximum." and "If Maximize is given an expression containing approximate numbers, it automatically calls NMaximize."

From documentation for NMaximize: "If f and cons are linear, NMaximize can always find global maxima, over both real and integer values." and "Otherwise, NMaximize may sometimes find only a local maximum."

Your function is neither linear nor polynomial and you get just a local maximum.

Clear[x, y, xx, yy]
peaks[x_, y_] := 
  3*(1 - x)^2*Exp[-(x^2) - (y + 1)^2] - 
   10*(x/5 - x^3 - y^5) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2];

Try different methods for NMaximize

NMaximize[{peaks[x, y], -2 <= x <= 2, -2 <= y <= 3}, {x, y}, 
   Method -> #] & /@ {"NelderMead", "DifferentialEvolution", 
  "SimulatedAnnealing", "RandomSearch"}

{{3.59249, {x -> 1.28568, y -> -0.00484756}}, {8.10621, {x -> -0.00931759, y -> 1.58137}}, {8.10621, {x -> -0.0093176, y -> 1.58137}}, {8.10621, {x -> -0.00931757, y -> 1.58137}}}

eldo · Accepted Answer · 2014-06-14T18:02:42.723

3

You should use FindMaximum:

FindMaximum[peaks[xx, yy], {xx, yy}]

{8.106, {xx -> -0.009, yy -> 1.581}}

plot = Plot3D[peaks[x, y], {x, -2, 2}, {y, -2, 2}, PlotRange -> All];

point = Graphics3D[{PointSize@0.05, Point[{-0.009, 1.581, 8.106}]}];

Show[plot, point]

enter image description here

On the other hand, Maximize returns in this case the last local maximum it can find:

enter image description here

edited Jun 14 '14 at 18:02

answered Jun 14 '14 at 17:49

eldo

67,911
5
60
168

ok, but... Why maxime gives bad results? – Mika Ike Jun 14 '14 at 17:52
When I dlouldn´t use Maximize? – Mika Ike Jun 14 '14 at 17:54
@Mika Ike - If you have to use NMaximize you can go with Bob Hanlon's answer. – eldo Jun 14 '14 at 18:22

score 2 · Answer 3 · edited Jun 14 '14 at 23:52

Your function peaks[x_, y_] := 3*(1 - x)^2.*Exp has a decimal point after the 2, which makes the 2 a Real and not an Integer. Remove it and Maximize[peaks[xx, yy], {xx, yy}] just returns the input, while NMaximize[peaks[xx, yy], {xx, yy}] gives {3.59249, {xx -> 1.28568, yy -> -0.00484755}}.

The function peaksis just to complex to allow a symbolic solution, I guess.

Added 19:52 CET

NMaximize[peaks[xx, yy], {xx, yy}, Method -> RandomSearch]`

gives

{8.10621, {xx -> -0.00931759, yy -> 1.58137}}`

See the help file for other options. Finding a global maximum numerically can always be made to fail; just add a very sharp spike, say

1000*Exp[-((x - u)^2 + (y - v)^2)/2/w^2]

and let w get extremely small. No numeric algorithm can ever locate that spike.

I tried, but the results are a little strange. See the edited post — Mika Ike, Jun 14 '14 at 17:50

Why does Maximize return a wrong answer?

3 Answers3

Linked