I have come up with the following functions, and I need to calculate a value, in order to have a condition valid. In my problem, all data are given, and I need to calculate variable l. I am using findroot, and I give an initial value of 25 (plausible assumption) in order to get what I need. Although my calculation is correct, it is really time consuming (I am using an i7 -3720 qm) and it takes almost 8-9 seconds. My question is, how can I speed up evaluation speed.
f1[x_, y_] :=
PDF[NormalDistribution[l, 3.975], x]*
PDF[NormalDistribution[53.559 - l, 3.975], y];
p = (1 - Sum[f1[i, i], {i, 0, 50}])/(1 -
Sum[f1[i, i], {i, 0, 50}]*1.4317884120211488);
f2[x_, y_] := If[x == y, f1[x, y]*1.4317884120211488`,
f1[x, y]/p];
A = SessionTime[];
sexpNEW =
FindRoot[Sum[f2[i, j], {i, 0, 50}, {j, 0, i - 6.5}] ==
0.4605263157894738, {l, 25}]
B = SessionTime[];
My results:
{l->29.8104}
{ speed:,8.6508650,seconds}
consider the following example, which is based on my previous, but this time some variables are unknown, lets call them u1, u2, u3
f1[x_,y_]:=PDF[NormalDistribution[u1,2.1],x]*PDF[NormalDistribution[u2,2.05],y];
p=(1-Sum[f1[i,i],{i,0,20}])/(1-Sum[f1[i,i],{i,0,50}]*u3);
f2[x_,y_]:=If[x==y,f1[x,y]*u3,f1[x,y]/p];
limit1=20.5;
FindRoot[
Sum[f2[i, j], {i, 0, limit1}, {j, 0, limit1 - i}] ==
0.4285714285714286`
&& Sum[f2[i, j], {i, 1, 20}, {j, 0, i - 1}] == 0.47417677642980943`
&& Sum[f2[i, j], {j, 1, 20}, {i, 0, j - 1}] ==
0.2558058925476603`, {u1, 12}, {u2, 12}, {u3,
1}] // AbsoluteTiming
{8.18311, {u1 -> 11.0489, u2 -> 10.0626, u3 -> 2.10168}}
How would I use compile as you did @xzczd, in order to use findroot for a 3*3 system?
edit: (26/2/2016)
I am studying compile
and I have a question regarding its applications. Could I use it in recursive functions too? I am using memoization, but even though I gain in speed, I was thinking that compile
could improve my speeds further. My question is, how could I define the following function with compile
? I tried to do so, but when I call my function, it keeps running forever.
This is my original function
f[t_, d_] :=
f[t, d] =
If[t > 80, 0,
If[d == A && t == B, 1,
p7h*f[t + 1, d - 7] + p5h*f[t + 1, d - 5] +
p3h*f[t + 1, d - 3] + p7a*f[t + 1, d + 7] +
p5a*f[t + 1, d + 5] + p3a*f[t + 1, d + 3] +
pns*f[t + 1, d]]]
and this is how i tried to rewrite it using Compile
fC=
Compile[{t, d},
If[t > 80, 0,
If[d == A && t == B, 1,
p7h*fC[t + 1, d - 7] + p5h*fC[t + 1, d - 5] +
p3h*fC[t + 1, d - 3] + p7a*fC[t + 1, d + 7] +
p5a*fC[t + 1, d + 5] + p3a*fC[t + 1, d + 3] +
pns*fC[t + 1, d]]]
]
t,d are intenger variables, and so are A and B. Thank you in advance.
Answer
You need Compile
, with option "EvaluateSymbolically" -> False
:
cf = Compile[{l}, #, RuntimeOptions -> "EvaluateSymbolically" -> False] &@
Sum[f2[i, j], {i, 0, 50}, {j, 0, i - 6.5}];
FindRoot[cf@l == 0.4605263157894738, {l, 25}] // AbsoluteTiming
{0.015624, {l -> 29.8104}}
Update:
The method I showed above is fully applicable to the new added sample:
cfgenerator =
Compile[{u1, u2, u3}, #, RuntimeOptions -> "EvaluateSymbolically" -> False] &;
{cf1, cf2, cf3} = cfgenerator /@ {Sum[f2[i, j], {i, 0, limit1}, {j, 0, limit1 - i}],
Sum[f2[i, j], {i, 1, 20}, {j, 0, i - 1}],
Sum[f2[i, j], {j, 1, 20}, {i, 0, j - 1}]};
FindRoot[cf1[u1, u2, u3] == 0.4285714285714286` &&
cf2[u1, u2, u3] == 0.47417677642980943` &&
cf3[u1, u2, u3] == 0.2558058925476603`, {u1, 12}, {u2, 12}, {u3, 1}] // AbsoluteTiming
{0.052929, {u1 -> 11.0489, u2 -> 10.0626, u3 -> 2.10168}}
Then I'd like to talk a little about why this method works. Simply speaking, your original code is slow because those Sum
s are very big, and Compile
is a way to speed up the evaluation of big expressions that are formed by relatively low level functions.
Well, to be honest, before answering, I hesitate for a while, because Compile
isn't really a easy-to-use function (see here for more information) and personaly I don't recommond not-that-experienced Mathematica user jumping into the world of compiling. However, your Sum
s can (luckily) evaluate to compilable algebraic expressions and are so suitable for Compile
that I can't help posting this answer. If you want to learn more about Compile
(as mentioned above, currently I don't recommend you to! ) , you can consider starting from here.
Comments
Post a Comment