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

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

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Aug 24 '15 at 4:29
• I'm not sure if this is answerable without more information about your function – Dr. belisarius Aug 24 '15 at 4:50
• My function is the phase difference between the transfer functions of two digital filters. Since phase is found by taking Arg[transferFunction[w]], it is generally not possible to do useful symbolic computation on the output. I want to produce a pair of quadrature filters, so the phase difference should be pi/2 in the region of the frequency domain for which I am doing the minimisation. I minimise the quantity (pi/2 - phase difference)^2, so when I am near the desired values for a and b, the function is nearly 0 everywhere in the region of integration. – sirhans Aug 24 '15 at 6:26
• @user3551064 Could you forge a simple example of your scenario? I'm sure it will be much easier for those trying to help you to start with concrete functions – Dr. belisarius Aug 24 '15 at 7:33
• Thanks for trying. I posted code in the question. I found an effective way to speed up the calculation by estimating the integral before calling NIntegrate. It's not as nice as actually handling the exception from NIntegrate but it will probably be enough to get the job done. I'd still like to know if there's a less kluge way to do it though. – sirhans Aug 24 '15 at 7:58

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.

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