Group Abstract Group Abstract

Message Boards Message Boards

Speed up this code computing a fixed distribution from a large data amount?

Alessio Lapolla
Posted 9 years ago

Dear community, I have a program that should give me a fixed distribution through a recurrence recursion. Because it is something related to the statistic I need to get the result for the program using a large amount of data (around 10^4 at least), this is the code:

RecurrenceDistributionRandomGraph[J_, H_, \[Beta]_, y_, Pf_, P0_, n_, 
  m_] := Block[  {x0, q, x, z},
  x0 = RandomVariate[P0, n];
  (*x0=Table[3,n];*)
  R = {x0};
  Do[
   x = Table[0, n];
   Do[
     q = Round[RandomVariate[Pf]];
    z = RandomSample[R[[i]], q];
    x[[j]] = N[y[z, J, H , \[Beta]]]
    , {j, 1, n}];
   R = Append[R, x];
   , {i, 1, m}]
  ]

and the inputs are, for example (they can be modified):

P0 = PoissonDistribution[5];
Pf = ProbabilityDistribution[q*(E^-3 3^q/q!)/3, {q, 0, +\[Infinity]}];
y[z_, J_, H_, \[Beta]_] := (
 E^(\[Beta] (J + H)) Product[i, {i, z}] + E^(-\[Beta] (J + H)))/(
 E^(\[Beta] (-J + H)) Product[i, {i, z}] + E^(\[Beta] (J - H))).

y is the recursion relation. Unfortunately my code seems to work properly but it really takes an huge quantity of time. For example even with few data:

RecurrenceDistributionRandomGraph[1., 1., 0.2, y, Pf, P0, 100, 
  20] // Timing
{305.652, Null}.

So how can I speed up the algorithm or maybe there is no solution? Unfortunately I have to store every set of data in the list (matrix) R, so I cannot make a quicker program in that sense. Thanks you for the answers
POSTED BY: Alessio Lapolla
4 Replies
Sort By:

rather than 

R = {x0}

you can:

R = ConstantArray[x0, m];

and then change:

AppendTo[R, x];

to

R[[i]] = x;

which is 5% faster because you don't need to append every time...
POSTED BY: Sander Huisman

Great, thanks you very much! I did not know that creations of random number was very expensive.
POSTED BY: Alessio Lapolla

Actually by doing some more 'refinements' you can speed it up even more:

ClearAll[RecurrenceDistributionRandomGraph,y]
RecurrenceDistributionRandomGraph[J_,H_,\[Beta]_,y_,Pf_,P0_,n_,m_]:=
    Block[{x0,q,qs,x,z},
            x0=RandomVariate[P0,n];
            R={x0};
            qs=Round[RandomVariate[Pf,{m,n}]];
            Do[
                x=ConstantArray[0.0,n];
                Do[
                q=qs[[i,j]];
                z=RandomSample[R[[i]],q];
                x[[j]]=y[z,J,H,\[Beta]]
                ,
                {j,n}
                ];
                AppendTo[R,x];
                ,
                {i,m}
            ]
        ]
P0=PoissonDistribution[5];
Pf=ProbabilityDistribution[q*(E^-3 3^q/q!)/3,{q,0,+\[Infinity]}];
y[z_,J_,H_,\[Beta]_]:=With[{prod=Times@@z},(E^(\[Beta] (J+H)) prod+E^(-\[Beta] (J+H)))/(E^(\[Beta] (-J+H)) prod+E^(\[Beta] (J-H)))]

which runs in 0.1 seconds on my machine. So more than 1000x faster than the original...
POSTED BY: Sander Huisman

Hi Alessio,

The main problem is that you create q everytime (in two loops). If you move q outside the inner Do loop it goes 100x faster:

ClearAll[RecurrenceDistributionRandomGraph,y]
RecurrenceDistributionRandomGraph[J_,H_,\[Beta]_,y_,Pf_,P0_,n_,m_]:=
    Block[{x0,q,qs,x,z},
        x0=RandomVariate[P0,n];
        R={x0};
        Do[
            x=ConstantArray[0.0,n];
            qs=Round[RandomVariate[Pf,n]];
            Do[
                (*q=Round[RandomVariate[Pf]];*)
                q=qs[[j]];
                z=RandomSample[R[[i]],q];
                x[[j]]=N[y[z,J,H,\[Beta]]]
                ,
                {j,1,n}
            ];
            R=Append[R,x];
            ,
            {i,1,m}
        ]
    ]
P0=PoissonDistribution[5];
Pf=ProbabilityDistribution[q*(E^-3 3^q/q!)/3,{q,0,+\[Infinity]}];
y[z_,J_,H_,\[Beta]_]:=(E^(\[Beta] (J+H)) Product[i,{i,z}]+E^(-\[Beta] (J+H)))/(E^(\[Beta] (-J+H)) Product[i,{i,z}]+E^(\[Beta] (J-H)))

executing:

RecurrenceDistributionRandomGraph[1., 1., 0.2, y, Pf, P0, 100, 20] // AbsoluteTiming

runs in 1.12 seconds on my laptop (down from 108 seconds). Note that I also used ConstantArray rather than Table, and I put 0.0 in there so the data-type does not change during evaluation (from integer to doubles).

POSTED BY: Sander Huisman
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard