In NMinimize, how can I efficiently handle NIntegrate errors for non-integrable functions?

Question

I am using NMinimize to find parameters that minimise the integral of a function in the least squares sense:

NMinimize[{NIntegrate[f[a,b,x]^2,{x,0,1}], -1 < a < 0, -1 < b < 0},{a,b}]

For some values of the parameters, the function f is not continuous but I know that the minimum value never occurs with parameters for which the function has a discontinuity. I want to do something like the following pseudocode:

For given parameters a and b,

if the function f[a,b,x] is discontinuous on {x,0,1} stop searching for local minima near a and b.

else

proceed with minimization of NIntegrate[f[a,b,x],{x,0,1}] as usual.

I can make this work by using a non-adaptive integration method and reducing the accuracy. But if I leave the accuracy at the default setting NIntegrate takes such a long time to determine that the function is discontinuous that the minimisation becomes computationally infeasible. If I reduce the accuracy then I can get a result in under five minutes but, of course, it's not very accurate.

I've pasted code for a simple example below. In this code I'm estimating the numerical integral before doing NIntegrate by evaluating the function at 100 points. If the estimated integral is not near zero then I use the estimated value instead of proceeding to call NIntegrate. This makes a huge improvement in the speed of the calculation. But I'd still like to know if there's a more elegant solution.

---

Here are some function definitions that aren't directly related to the question but are necessary to run the example code.

NonzeroLength[X_] := Length[X] - LengthWhile[Reverse[X], # == 0 &]

TransferFunctionNumOrDen[B_, z_] := 
  Total[Table[B[[n + 1]]*(z^(-n)), {n, 0, NonzeroLength[B] - 1}]]

TransferFunctionForPlotting[A_, B_, \[Omega]_] := 
  Module[{highestDegree}, 
   highestDegree = Max[NonzeroLength[A], NonzeroLength[B]] - 1;
   (TransferFunctionNumOrDen[B, E^(-I \[Omega])]*
      E^(-I \[Omega]*highestDegree))/(TransferFunctionNumOrDen[A, 
      E^(-I \[Omega])]*E^(-I \[Omega]*highestDegree))];

secondOrderDelayedAllPass[a_, \[Omega]_] := 
  TransferFunctionForPlotting[{1, 0, a}, {0, a, 0, 1}, \[Omega]];

secondOrderAllPass[a_, \[Omega]_] := 
  TransferFunctionForPlotting[{1, 0, a}, {a, 0, 1}, \[Omega]];

---

Here's where we start setting up the numerical minimisation

f[\[Gamma]_, \[Beta]_, lb_, ub_] := 
 Module[{phaseDifference, estimate, n = 100},
  phaseDifference[\[Omega]_] := Arg[
    (Times @@ (secondOrderDelayedAllPass[#, \[Omega]] & /@ \[Gamma])) *
    (Times @@ (secondOrderAllPass[#, \[Omega]] & /@ \[Beta]))^(-1)
 ];

  (* estimate the integral with a table *)
  estimate = Total[
     Table[(\[Pi]/2 - phaseDifference[\[Omega]])^2, {\[Omega], lb, 
       ub, (ub - lb)/n}]]/(n + 1);

  If[
   estimate > .3,
   estimate,
   NIntegrate[(\[Pi]/2 - phaseDifference[\[Omega]])^2, {\[Omega], lb, ub}]
   ]
  ]

fmin = 15/22050;
mn = -1 ;
mx = 1;

NMinimize[{f[{a}, {b}, fmin \[Pi], (1 - fmin) \[Pi]], mn  < a < mx, 
  mn < b < mx }, {a, b}, Method -> "DifferentialEvolution"]

Answer 1

The use of NumericQ as mentioned by MarcoB and Guess who it is. in the comments seems to be important. Also, estimating the integral using a Total[Table[]] as in the code I posted in the second revision of the question makes this method computational feasible enough to solve my problem.

Answer 2

You can find and exclude the part of the domain where the integrand is discontinuous.

To find where the integrand is discontinuous, look for where the argument of Arg is real and negative, since Arg has a branch cut discontinuity there:

z[\[Gamma]_,\[Beta]_,\[Omega]_]:=(Times@@(secondOrderDelayedAllPass[#,\[Omega]]&/@\[Gamma]))*
                                 (Times@@(secondOrderAllPass[#,\[Omega]]&/@\[Beta]))^(-1)
real=ComplexExpand[Re[z[{a},{b},\[Omega]]]]//FullSimplify;
imag=ComplexExpand[Im[z[{a},{b},\[Omega]]]]//FullSimplify;
\[Omega]p[a_,b_]=\[Omega]/.Solve[
Reduce[{real<0,imag==0,mn<a<mx,mn<b<mx,\[Omega] \[Element] Reals},\[Omega],Reals],\[Omega]];

which can be excluded from NIntegrate if desired. But to follow your pseudocode, I instead tried constraining the minimization to the portion of the domain where the integrand is continuous between the integration limits fmin \[Pi] and (1 - fmin) \[Pi]:

region=And @@ MapIndexed[#1/.C[1]->C[First[#2]]&, Not/@ 
       (fmin \[Pi]<#[[1]]<(1-fmin) \[Pi]&&#[[2]]& /@ \[Omega]p[a,b])]

It would have been nice if I could directly feed region to NMinimize to constrain the minimization domain. However, Mathematica is ill-equipped to deal with irregular domains, as you would find if you simply tried NMinimize[{1, region},{a,b}].

Instead, I indirectly constrained the domain by having the cost function return a large value, which I chose as the upper bound of the integral because $\infty$ doesn't work, whenever it's called outside the domain, so the minimizer would search elsewhere:

Clear[f]
discontinuousQ[a_,b_,lb_,ub_]:=AnyTrue[Simplify[#,Assumptions->C[1]\[Element] Integers]& /@
              (lb<#<ub& /@ Cases[\[Omega]p[a,b], Except[Undefined]]),#=!=False&]
f[a_?NumberQ,b_?NumberQ,lb_,ub_]:=Module[{phaseDifference},
  phaseDifference[\[Omega]_]:=Arg[z[{a},{b},\[Omega]]];
  If[discontinuousQ[a,b,lb,ub],
     (ub-lb)(3\[Pi]/2)^2,NIntegrate[(\[Pi]/2-phaseDifference[\[Omega]])^2,{\[Omega],lb,ub}]]]
AbsoluteTiming[min = NMinimize[{f[a, b, fmin \[Pi], (1 - fmin) \[Pi]], 
                     mn < a < mx, mn < b < mx}, {a, b}, Method -> "DifferentialEvolution"]]

{264.747,{0.135989,{a->-0.894316,b->-0.429818}}}

Here's a plot of f showing the minimum:

Show[ContourPlot[f[a,b,fmin \[Pi],(1-fmin) \[Pi]],{a,mn,mx},{b,mn,mx},Contours->{1,10}],
     Graphics[Point[{a,b}]/.min[[2]]]]