1
$\begingroup$

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}}
$\endgroup$
  • 6
    $\begingroup$ What on earth makes you think there's a linear relationship between function speed and finding the roots? In many (most) cases, the time is spent in the surrounding internal machinations. Add step/evaluation monitors, and put a dummy counter in your function to count how many calls it gets - you may be surprised. Then, when you've observed that, add ?NumericQ to the argument in your two versions of func... you should find the result educational. $\endgroup$ – ciao Apr 21 '16 at 23:20
  • $\begingroup$ @ciao Hi, ciao. I don't quite understand your claim. The three func versions produced the same result. Why will there be any difference in internal steps? As I tested, the internal processes are exactly the same. FindRoot[func[u] - num == 0, {u, -100, 100}, StepMonitor :> Print["Step to u = ", u]] for all three func version. On the other hand, the most time cosuming part in finding root is calculate the function values at different point, Why shouldn't we expect a linear relationship?? $\endgroup$ – matheorem Apr 22 '16 at 3:09
  • $\begingroup$ Clear[func2]; func2[u_]... gives me a factor of 2 or so improvement. As for general speed issue, others note it has more to do with FindRoot than function evaluation speed. Which means various options and possibly Method suboptions will need to be tried for purposes of speed tuning. $\endgroup$ – Daniel Lichtblau Apr 22 '16 at 4:40
  • $\begingroup$ @DanielLichtblau Hi, DanielLichtblau. I see many upvotes on ciao's comment and also Sander's, but I just still don't get the point. The stepmonitor shows the same for func1 and func2, now func2 is 25 times faster than func1. But the speed of whole process is not that different. What did FindRoot spend its time on? $\endgroup$ – matheorem Apr 22 '16 at 5:21
  • $\begingroup$ I just noticed a similar discussion here as one below my answer. I hope the following clarifies it better. 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 evaluation time to 20.04 seconds. I updated my answer. $\endgroup$ – Sander Apr 22 '16 at 8:13
6
$\begingroup$

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]

enter image description here

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

enter image description here

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}}
$\endgroup$
  • $\begingroup$ Hi, Sander, thank you very much for the answer. But my point in my post is that "why speed up func has no effect on speed of FindRoot even though the internal steps are exactly the same?" It seems that you havn't explained this in your answer $\endgroup$ – matheorem Apr 22 '16 at 4:16
  • 1
    $\begingroup$ I guess it is implicit. 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. I updated my answer for completeness. $\endgroup$ – Sander Apr 22 '16 at 8:10
  • $\begingroup$ I update my post again, have a look again : ) $\endgroup$ – matheorem Apr 22 '16 at 15:33
  • $\begingroup$ Nice research, will look into it ..... but it will be tomorrow :) $\endgroup$ – Sander Apr 22 '16 at 17:00

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.