Why timing of FindRoot is almost the same even if the function is speedup significantly

Question

This question is almost the same issue as this one. However, it seems that it didn't draw enough attention. Now I encountered the same problem.

I previously make this post, I will copy the code here. below code defines a function func1 and a value num,

n = 100;
block = Symmetrize[RandomReal[1., {n, n}], Symmetric[{1, 2}]];
diagF = SparseArray[Band[{1, 1}] -> Normal /@ #] &;
mat = diagF[{block, block, block, block}];
dim = Length@mat;
eigs = Transpose@Sort@Transpose@Eigensystem@mat;
eigsdensity = {eigs[[1]], Abs[eigs[[2]]]^2};
fermiCom = 
  Compile[{{e, _Real}, {u, _Real}, {t, _Real}}, 
   Module[{tmp}, tmp = Exp[(u - e)/t]; tmp = 1./(1. + tmp); 
    tmp = 1. - tmp], RuntimeOptions -> "Speed"];
Clear[fermi];
fermi[e_?NumericQ, u_?NumericQ, t_?NumericQ] := fermiCom[e, u, t];
num = Length@Select[eigs[[1]], # < 0 &] + 1;
Clear[func1];
func1[u_] := 
 Total[Sum[
   eigsdensity[[2, i]] fermi[eigs[[1, i]], u, 0.0001], {i, 1, dim}]]

now I want to find the root like this

FindRoot[func1[u] - num == 0, {u, -100,100}]// RepeatedTiming
(*{0.179, {u -> 0.0407121}}*)

it is pretty fast, however not fast enough, since it is the bottleneck of my code, I need it to be faster.

I found the func1 itself have performance issue. Timing shows

func1[1.] // RepeatedTiming
(*{0.0500, 256.}*)

I could significantly improve the speed of func1 by defining it like

Clear[func2];
func2[u_] := (fermiList = 
   Table[fermi[eigs[[1, i]], u, 0.0001], {i, 1, dim}]; 
  Total[eigsdensity[[2]]* fermiList, 2])

Now

func2[1.] // RepeatedTiming
(*{0.0021, 256.}*)

25 times speedup, but... look at FindRoot again

FindRoot[func2[u] - num == 0, {u, -100,100}]// RepeatedTiming
(*{0.117, {u -> 0.0407121}}*)

Nothing significantly speed up !! Why?? Shouldn't be 25 times faster?

Making func2[u_] into func2[u_?NumericQ] can reduce timing to 0.041, but still not the anticipated speed up.

How to correctly make it faster?

---

update

current comment and answer argue that the timing of FindRoot is much related to the internal process which could be much longer than the evaluation of my function.

I definitely can't agree on this. Here is a simple demonstration to prove this

FindRoot[Sin[1/x] == 0.5, {x, -1, 100000000}, 
 StepMonitor :> Print["Step to x = ", x]]

When evaluate the above, it shows that there are approximately 30 steps, and how much will this FindRoot takes?

FindRoot[Sin[1/x] == 0.5, {x, -1, 100000000}] // RepeatedTiming
(*{0.00060, {x -> 1.90986}}*)

This is so fast. You see there are no supposed overhead in the internal process

---

update 2

After reconsider the timing result, I think I have to ask why func1 is faster than expected.

As you can see func1[1.] takes 0.05s while FindRoot using func1 takes 0.179s, this is the time of only 3 evaluations!! However the StepMonitor clearly shows that there are much more than 3 steps.

As for FindRoot using func2, as I said, it should be faster.

Sander has mentioned the root finding algorithm. According to wikipedia, Brent method is actually a combination of the bisection method, the secant method and inverse quadratic interpolation. It is not that complex, just more condition branch which I think won't take too much time.

For my particular func, as I mentioned in another post, it has a good property, that is monotonic. So bisection is simple and robust for this task. Here I implemented a bisection, and shows a normal linear timing dependence on func.

Clear[bisection]
bisection[func_, a1_, a2_] := Module[{x1, x2, x3, f1, f2, f3},
  x1 = a1; x2 = a2;
  f1 = func[x1]; f2 = func[x2];
  If[f1*f2 < 0,
   Do[
    x3 = x1 - f1*(x2 - x1)/(f2 - f1);
    f3 = func[x3];
    If[Abs[f3] < 10.^-6,(*Print["got it","\t","step = ",i]*)Break[]]; 
    If[f1*f3 < 0, x2 = x3; f2 = f3, x1 = x3; f1 = f3];
    , {i, 1, 100}], 
   Print["two initial function values must be one positive and one \
negative"]]; x3]

define

Clear[func1new];
func1new[u_] := func1[u] - num;
Clear[func2new];
func2new[u_] := func2[u] - num;

and we have these timings

In[229]:= Table[func1new[u], {u, 1, 20}]; // AbsoluteTiming

Out[229]= {2.32441, Null}

In[230]:= Table[func2new[u], {u, 1, 20}]; // AbsoluteTiming

Out[230]= {0.102519, Null}

In[231]:= bisection[func1new, -100, 100] // AbsoluteTiming

Out[231]= {2.7258, 0.0697169}

In[232]:= bisection[func2new, -100, 100] // AbsoluteTiming

Out[232]= {0.119729, 0.0697169}

In[233]:= FindRoot[func2new[u] == 0, {u, -100, 100}] // AbsoluteTiming

Out[233]= {0.289189, {u -> 0.0697169}}

In[234]:= FindRoot[func1new[u] == 0, {u, -100, 100}] // AbsoluteTiming

Out[234]= {0.423125, {u -> 0.0697169}}

As you can see, my bisection shows linear timing dependence for func1new and func2new. While FindRoot's behavior is odd, mysteriously fast for func1new and slow for func2new

---

update 3

If we add ?_NumericQ

Clear[func1new];
func1new[u_?NumericQ] := func1[u] - num;
Clear[func2new];
func2new[u_?NumericQ] := func2[u] - num;

Then linear dependence is restored for FindRoot

In[568]:= FindRoot[func1new[u] == 0, {u, -100, 100}] // AbsoluteTiming

Out[568]= {2.50487, {u -> 0.0785305}}

In[569]:= FindRoot[func2new[u] == 0, {u, -100, 100}] // AbsoluteTiming

Out[569]= {0.112653, {u -> 0.0785305}}

Answer 1

FindRoot is slowing this down, rather than your function func.

---

UPDATE 1

If every FindRoot evaluation step takes 1 second, and every func evaluation takes 0.02 seconds, then after 20 steps, you will have spent 20.4 seconds. Improving func's performance with 90%, hence it will now take 0,002 seconds, brings down your evaluattion time to 20.04 seconds.

---

UPDATE 2

Have a look at the search algorithm that is used here at Wikipedia's page on Brent's search algorithm.

Your func is difficult and Sin roots are easy to discover. Without going into detail on the algorithm, imagine finding the highest peak when there is one mountain peak and there is a very thick fog, so you cannot see the peak (mathematical similar to having one place where the derivative is zero) which is easy ... continue taking a step towards the steepest slope and you will get there. If there are however many mountains, with different peak altitudes, this steepest ascent method does not work and you have to take a much more complex approach (as taking many more steps may simply not work or take too long). There are many places with a Derivative of zero, it is expensive to guarantee you found them all, remember there is fog, and your func is not simply predictable like Sin.

This is where your Sin example is treated fundamentally straightforward by the search algorithmm and func less so.

---

END OF UPDATES

That may explain it, but does not solve your issue. Hopefully the approach below can help you out speeding up FindRoot instead.

What I typically do is finetune the use of FindRoot by setting the Options for Method, MaxIterations, AccuracyGoal etc.

To find out about how this works out for this context, find for example the speed of conversion as follows:

{res, data} = Reap[FindRoot[func[u] - num == 0, {u, -100, 100}, EvaluationMonitor :> Sow[u]]];
ListLogPlot[data[[1]], PlotRange -> All]

Check out that Method -> Brent performs best:

TableForm[
 Table[
   Block[{e = 0},
   {method,
    FindRoot[func[u] - num == 0, {u, -100, 100}, Method -> method, 
        EvaluationMonitor :> e++],
   e}], 
   {method, {Automatic, "Newton", "Brent", "Secant"}}], 
 TableHeadings -> {{}, {"Method", "fun[u]", "Evaluations"}} ]

After playing around with other settings, you can subsequently increase performance (at a penalty of decreased accuracy in this example) and cutting the first several iterations out by improving starting point approaches an acceptable accuracy:

FindRoot[func[u] - num == 0, {u, -100, 100}] // AbsoluteTiming
FindRoot[func[u] - num == 0, {u, -100, 100}, Method -> "Brent", 
    MaxIterations -> 15] // AbsoluteTiming
FindRoot[func[u] - num == 0, {u, -1, 1}, Method -> "Brent", 
    MaxIterations -> 15] // AbsoluteTiming

{0.306293, {u -> 0.00242511}}
{0.217473, {u -> 0.00207048}}
{0.223452, {u -> 0.00242421}}