Re: Time, Inverse, Simulation simple dynamical system, speed issue
- To: mathgroup at smc.vnet.net
- Subject: [mg77438] Re: [mg77380] Time, Inverse, Simulation simple dynamical system, speed issue
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Fri, 8 Jun 2007 05:41:34 -0400 (EDT)
- References: <21255603.1181208515635.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
1) Some of that is gibberish (= E2 = 89 =, for instance); copy as plain
text, THEN paste into e-mail.
2) Since you're setting RelWel={} every iteration AND Clearing it ever y
iteration, why append anything to it or use it at all? You're actually
appending to Value[t]. But see point (3).
3) AppendTo is a slow way to gather results over hundreds or thousands of
iterations. Look into Sow and Reap, instead. For instance,
AppendTo[ThetaTable, ThetaMatrix];
would become
Sow[ThetaMatrix,"thetaTable"]
AppendTo[RelWel, Rho[[1, 1]]*Wel[[1]]/Total[Rho.Wel]] would become
Sow[Rho[[1, 1]]*Wel[[1]]/Total[Rho.Wel],"Value"]
and Do[AppendTo[Value[t], RelWel[[t]]], {t, 1, 1000}] would just go away.
You'd start the loop with something like results=Reap at Timing... or
Timing@Reap@For... and "results" would include a tagged list for each Sow
in the loop. You'd no longer use ThetaTable, RelWel, or Value[t].
It's possible WRI has greatly improved AppendTo at some point and I didn't
hear about it, but it used to be, at least, well known that appending
elements to a long list was prohibitively slow. I suspect AppendTo is
using far more time than your calculations.
4) N[Inverse[IdentityMatrix[2] - ThetaMatrix.LamRho, Method ->
DivisionFreeRowReduction].ThetaMatrix.divi, 1000] is NOT faster than the
same thing without N. But using Wel in a later statement IS faster
(probably), since Wel is approximate. Otherwise, the rational results are
likely to grow huge numerators and denominators.
It's possible you'd do just as well replacing N with Simplify, but
probably not.
5) Do[AppendTo[Value[t], RelWel[[t]]], {t, 1, 1000}] is MUCH slower than
the equivalent
Value[t] = Join[Value[t],RelWel] (* if RelWel has 1000 elements *)
or
Value[t] = Join[Value[t],Take[RelWel,1000]] (* if RelWel has more than
1000 elements, but you only want 1000 *)
But, in view of (3) above, RealWel and Value should go away entirely, and
so does this Do loop.
6) If the theta function will be called thousands of times, it's much more
efficient to define it this way:
theta[i_, k_] :=
theta[i, k] =
Lam[[i, k]]*Rho[[i, i]]*
Wel[[i]]/Sum[Lam[[j, k]]*Rho[[j, j]]*Wel[[j]], {j, 1, 2}];
That doesn't actually work, since Wel changes during the run, but you can
do this instead:
Clear[a, b, temp]
Wel = {a, b};
Lam = {{1/2, 2/3}, {1/2, 1/3}};
Rho = {{1/2, 0}, {0, 648/1000}};
temp[i_, k_] :=
Lam[[i, k]]*Rho[[i, i]]*
Wel[[i]]/Sum[Lam[[j, k]]*Rho[[j, j]]*Wel[[j]], {j, 1, 2}];
ThetaMatrix := Evaluate@Table[Simplify@temp[i, k],
{i, 1, 2}, {k, 1, 2}];
ThetaMatrix
{{(125 a)/(125 a + 162 b), (125 a)/(125 a + 81 b)}, {(162 b)/(
125 a + 162 b), (81 b)/(125 a + 81 b)}}
ThetaMatrix stores those results, and it always uses the current value of
a and b. The theta function is not needed and neither is Clear[Wel],
Clear[ThetaMatrix], or ThetaMatrix = Table[theta[i, k], {i, 1, 2}, {k,1,
2}]. (temp isn't needed again, either... unless you want to use a new Lam
or Rho.)
Replace every mention of Wel with {a,b}. That is,
{a,b} = {5,5} instead of Wel = {5,5}, a instead of Wel[[1]], Rho.{a,b}
rather than Rho.Wel, and of course
{a,b} = N[Inverse[IdentityMatrix[2] - ThetaMatrix.LamRho,
Method -> DivisionFreeRowReduction].ThetaMatrix.divi, 1000];
7) There's a missing semicolon after ListPlot, before Clear. (Clear wasn't
needed anyway.) In version 6 you'd need Print@ListPlot to see the charts .
The ListPlot would be awfully dull, anyway... doesn't RealWel have only
two points? You're supposedly appending to it each iteration, but you're
also reinitializing it to {} each iteration, so isn't it just...
ListPlot[Rho[[1, 1]]*Wel[[1]]/Total[Rho.Wel], PlotStyle -> {Hue[.8]},
PlotRange -> {0, 1},
ImageSize -> 350]
I hope some of that helps.
Bobby
On Thu, 07 Jun 2007 03:03:33 -0500, kristoph <kristophs.post at web.de> wrote:
> Attached you find a simulation of 5 dynamical systems each consisting
> of 1000 Periods. For each period t = 1,...,1000 an inverse of a 2x2m
> matrix needs to be computed.
> It takes about 6 seconds to simulate the 5 dynamical systems. I would
> be grateful for any hint in getting the simulation done much quicker.
> It seems that calculating the inverse
>
> Wel = N[Inverse[
> IdentityMatrix[
> 2] - ThetaMatrix.LamRho, Method -> \
> DivisionFreeRowReduction].ThetaMatrix.divi, 1000]
>
> takes most of the time, therefore I tried to solve the problem
> numerically (see above). This is a lot faster then using just
> Inverse[...]. I tried LinearSolve, but it is slower then the above.
> Decreasing the precision can result in a significant error, depending
> on the parameter constellation.
>
> As you might have guessed I have to simulate not only 5 dynamical
> systems for different parameter constellations. I would also like to
> increase the order of the matrix that needs to be inverted.
>
> Here is the code. Thanks in advance. Kristoph
>
>
> << Statistics`ContinuousDistributions`;(*packages and starting values
> needed for the simulation*)
> << Graphics`MultipleListPlot`;
>
> Wel = {5, 5};
> Lam = {{1/2, 2/3}, {1/2, 1/3}};
> Rho = {{1/2, 0}, {0, 648/1000}};
> RelDiv = {{2/3, 1/3}, {1/3, 2/3}, {1, 0}};
> theta[i_, k_] := Lam[[i, k]]*Rho[[i, i]]*Wel[[i]]/Sum[Lam[[j,
> k]]*Rho[[j, j]]*Wel[[j]], {j, 1, 2}];
>
> Timing[For[n = 1, n =E2=89=A4 5, n++,(*begin loop for the 5 simulations
> of the
> dynamical systeme*)
> Wel = {5, 5};(*starting values*)
> RelWel = {};(*needed for the plots, see below*)
> ThetaTable = {};(*values needed for recursive calculations, saves
> time*)
>
> For[t = 0, t =E2=89=A4 1000, t++,(*begin loop for the dynamical
> system
> consiting of t = 1000 Periods*)
> Clear[ThetaMatrix];
> ThetaMatrix = Table[theta[i, k], {i, 1, 2}, {k, 1, 2}];
> AppendTo[ThetaTable, ThetaMatrix];
> AppendTo[RelWel, Rho[[1, 1]]*Wel[[1]]/Total[Rho.Wel]];
> Clear[Wel];
> x[t] = Random[];(*random pertubations are drawn from the above
> matrix*)
> If[x[t] =E2=89=A4 1/3, divi = Div[[1, All]], If[1/3 < x[t] =E2=89=
> =A4 2/3,
> divi = Div[[2, All]], divi = Div[[3, All]]]];
>
> Wel = N[Inverse[IdentityMatrix[2] -
> ThetaMatrix.LamRho,
> Method -> DivisionFreeRowReduction].ThetaMatrix.divi,
> 1000];(*the inverse needed for t + 1*)
> ];(*end loop dynamical system*)
>
> Do[AppendTo[Value[t], RelWel[[t]]], {t, 1, 1000}];
> ListPlot[RelWel, PlotStyle -> {Hue[.8]}, PlotRange -> {0, 1},
> ImageSize -> 350]
> Clear[RelWel, ThetaTable];
>
> ] (*end loop simulations*)
> ]
>
>
>
--
DrMajorBob at bigfoot.com