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}}
?NumericQ
to the argument in your two versions offunc
... you should find the result educational. $\endgroup$ – ciao Apr 21 '16 at 23:20func
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 threefunc
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:09Clear[func2]; func2[u_]...
gives me a factor of 2 or so improvement. As for general speed issue, others note it has more to do withFindRoot
than function evaluation speed. Which means various options and possiblyMethod
suboptions will need to be tried for purposes of speed tuning. $\endgroup$ – Daniel Lichtblau Apr 22 '16 at 4:40func1
andfunc2
, nowfunc2
is 25 times faster thanfunc1
. But the speed of whole process is not that different. What didFindRoot
spend its time on? $\endgroup$ – matheorem Apr 22 '16 at 5:21