5
$\begingroup$

As far as I understand, the merit function measures how far the data is from the fitted function.

I am calulcating the second partial derivative of a function, and then numerically finding its maxima and minima with FindRoot.

It always outputs an error message along the lines of:

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than 30.` digits of working precision to meet these tolerances.

Which makes me think: can I just change the settings of the merit function? I.e., lower the threshold for what is considered a sufficient change? And by the way, what is the default setting?


CODE:

k = 2 π;
k1 = {k, 0};
k2 = FullSimplify[RotationMatrix[(2 π)/5]. k1];
k3 = FullSimplify[RotationMatrix[(2 π)/5]. k2];
k4 = FullSimplify[RotationMatrix[(2 π)/5] .k3];
k5 = FullSimplify[RotationMatrix[(2 π)/5] .k4];
E1[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k1.{x, y} - α1]
E2[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k2.{x, y} - α1]
E3[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k3.{x, y} - α1]
E4[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k4.{x, y} - α1]
E5[x_, y_, I1_, α1_] := Sqrt[I1] Exp[I k5.{x, y} - α1]
Eges[x_, y_] := 
 E1[x, y, 1, 0] + E2[x, y, 1, 0] + E3[x, y, 1, 0] + E4[x, y, 1, 0] + 
  E5[x, y, 1, 0]
IIges[x_, y_] := 
  Conjugate[Eges[x, y]] Eges[x, y] // ComplexExpand // 
   Simplify; (* can't differentiate Conjugate, had to change this *)
{dx[x_, y_], dy[x_, y_]} = 
  D[IIges[x, y], {{x, y}}]; (* calculating partial derivatives *)
hes[x_, y_] = 
  D[IIges[x, y], {{x, y}, 
    2}]; (* Hessian matrix: Second derivative matrix *)

(* Find critical points *)
crit = Cases[
   Normal[ContourPlot[dx[x, y] == 0, {x, -xm, xm}, {y, -xm, xm}, 
     ContourStyle -> None, Mesh -> {{0}}, 
     MeshFunctions -> Function[{x, y, z}, dy[x, y]]]], 
   Point[{x0_, y0_}] :> ({\[FormalX], \[FormalY]} /. 
      FindRoot[{dx[\[FormalX], \[FormalY]], 
        dy[\[FormalX], \[FormalY]]}, {{\[FormalX], x0}, {\[FormalY], 
         y0}}, WorkingPrecision -> 30, AccuracyGoal -> 20, 
       PrecisionGoal -> 20, MaxIterations -> 1000]), ∞];

(* Identify points as minima, maxima or saddle points *)
hl = hes @@@ crit;
mnp = PositiveDefiniteMatrixQ /@ hl; (*pick minima*)
mxp = PositiveDefiniteMatrixQ /@ (-hl); (*pick maxima*)
sdp = Thread[mnp ⊽ mxp];(*saddle points are leftovers*)

(* Finding coordinates of critical points *)
mini = Pick[crit, mnp];
maxi = Pick[crit, mxp];
sadl = Pick[crit, sdp];

(* Plot stuff!
   I changed the Image Size because it was too small *)
{Legended[
  ContourPlot[IIges[x, y], {x, -xm, xm}, {y, -xm, xm}, 
   ImageSize -> Scaled[0.3], ColorFunction -> "DarkTerrain", 
   Contours -> 15, 
   Epilog -> {AbsolutePointSize[6], {Cyan, Point[mini]}, {Orange, Point[maxi]}}], 
  PointLegend[{Cyan, Orange}, {"Minima", "Maxima"}]], 
 Show[Plot3D[IIges[x, y], {x, -xm, xm}, {y, -xm, xm}, 
   BoundaryStyle -> None, Boxed -> False, 
   ColorFunction -> "DarkTerrain", Mesh -> 15, 
   MeshFunctions -> {#3 &}], 
  Graphics3D[{{Cyan, Sphere[mini, 1/20]}, {Orange, 
      Sphere[maxi, 1/20]}} /. {x_?NumericQ, y_?NumericQ} :> {x, y, IIges[x, y]}],
  ImageSize -> Scaled[0.43], PlotRange -> All, Frame -> False]}
$\endgroup$
  • $\begingroup$ It would help to have an example to play with. $\endgroup$ – Michael E2 Nov 12 '15 at 12:28
  • 1
    $\begingroup$ I don't have time to write a fully fleshed-out answer at the moment, but did you notice that Conjugate[Eges[x, y]] == Eges[-x, -y]? $\endgroup$ – J. M. is away Nov 12 '15 at 13:04
  • $\begingroup$ that fixed it, thanks. But why? $\endgroup$ – SuperCiocia Nov 12 '15 at 14:38
  • $\begingroup$ actually sorry no, it didn't fix it $\endgroup$ – SuperCiocia Nov 12 '15 at 19:18
5
$\begingroup$

Here is the merit function:

Both step control methods [line search and trust region] were developed originally with minimization in mind. However, they apply well to finding roots for nonlinear equations when used with a merit function. In the Wolfram Language, the 2-norm merit function r(x).r(x) is used. -- Introduction to Step Control

The "sufficient decrease" criterion is controlled by AccuracyGoal. We can check by solving for one of the OP's troublesome points and determining how far from zero the merit function is. Its square-root should be closer to zero than 10^-ag, where ag is the setting for AccuracyGoal.

sol = Block[{x0 = 1.5107119013044863`, y0 = -1.6225965143185352`}, 
  FindRoot[{dx[\[FormalX], \[FormalY]], dy[\[FormalX], \[FormalY]]},
    {{\[FormalX], x0}, {\[FormalY], y0}},
    WorkingPrecision -> 30, AccuracyGoal -> 20, PrecisionGoal -> 20, MaxIterations -> 100]]

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than 30.` digits of working precision to meet these tolerances. >>

(*
  {\[FormalX] ->  1.51066343431975599208008475284,
   \[FormalY] -> -1.62258704517760281788841182199}
*)

{dx[\[FormalX], \[FormalY]], dy[\[FormalX], \[FormalY]]} . 
 {dx[\[FormalX], \[FormalY]], dy[\[FormalX], \[FormalY]]} /. sol  (* merit function *)
ag = -Log10[%]/2                                                  (* AccuracyGoal *)
(*
  0.065933684103633462001724717
  0.59044632849264420129427203
*)

Block[{x0 = 1.5107119013044863`, y0 = -1.6225965143185352`}, 
 FindRoot[{dx[\[FormalX], \[FormalY]], dy[\[FormalX], \[FormalY]]},
  {{\[FormalX], x0}, {\[FormalY], y0}}, 
  WorkingPrecision -> 30, AccuracyGoal -> ag, PrecisionGoal -> 20, 
  MaxIterations -> 100]]
(*
  {\[FormalX] ->  1.51066407892396160887133952303,
   \[FormalY] -> -1.62258717078894941887987024328}
*)

A slight higher setting, AccuracyGoal -> acc + 1*^-6, results in the FindRoot::lstol message.

Of course, this change does not make the solution returned any closer to a root. It just hides the fact that the objective function is not very close to zero. Quiet has the same effect. :-)

$\endgroup$
  • $\begingroup$ Very nice answer indeed. This kind of confusion is so frequent that I believe this Q/A should be upvoted to death. Good work $\endgroup$ – Dr. belisarius Nov 12 '15 at 15:35

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.