Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: how to speed up lenghty formulas??

  • To: mathgroup at smc.vnet.net
  • Subject: [mg84145] Re: [mg84104] how to speed up lenghty formulas??
  • From: DrMajorBob <drmajorbob at bigfoot.com>
  • Date: Tue, 11 Dec 2007 15:28:46 -0500 (EST)
  • References: <18972263.1197362311113.JavaMail.root@m35>
  • Reply-to: drmajorbob at bigfoot.com

The obvious construct is FoldList, and the recursion formula has to be

f = a0 + a1 #2 + b1 #1 &;

The only question is what the initial value must be and how much of the 
list I'll get this way.

I start with a small and very general X matrix, I use just v rather than 
Variance[X], and I name the initial value i.

Here's h:

Clear[h]
h[0] := v;
h[1] := a0 + h[0];
h[t_] := h[t] = a0 + a1*X[[t - 1]]^2 + b1*h[t - 1]
Array[h, 6]

{a0 + v, a0 + b1 (a0 + v) + a1 x[1]^2,
  a0 + b1 (a0 + b1 (a0 + v) + a1 x[1]^2) + a1 x[2]^2,
  a0 + b1 (a0 + b1 (a0 + b1 (a0 + v) + a1 x[1]^2) + a1 x[2]^2) +
   a1 x[3]^2,
  a0 + b1 (a0 +
      b1 (a0 + b1 (a0 + b1 (a0 + v) + a1 x[1]^2) + a1 x[2]^2) +
      a1 x[3]^2) + a1 x[4]^2,
  a0 + b1 (a0 +
      b1 (a0 +
         b1 (a0 + b1 (a0 + b1 (a0 + v) + a1 x[1]^2) + a1 x[2]^2) +
         a1 x[3]^2) + a1 x[4]^2) + a1 x[5]^2}

and here's the FoldList:

f = a0 + a1 #2^2 + b1 #1 &;
FoldList[f, i, X]

{i, a0 + b1 i + a1 x[1]^2,
  a0 + b1 (a0 + b1 i + a1 x[1]^2) + a1 x[2]^2,
  a0 + b1 (a0 + b1 (a0 + b1 i + a1 x[1]^2) + a1 x[2]^2) + a1 x[3]^2,
  a0 + b1 (a0 + b1 (a0 + b1 (a0 + b1 i + a1 x[1]^2) + a1 x[2]^2) +
      a1 x[3]^2) + a1 x[4]^2,
  a0 + b1 (a0 +
      b1 (a0 + b1 (a0 + b1 (a0 + b1 i + a1 x[1]^2) + a1 x[2]^2) +
         a1 x[3]^2) + a1 x[4]^2) + a1 x[5]^2}

A little thought shows that i must be a0 + v, and we can quickly check  
that:

Array[h, 6] == FoldList[f, a0 + v, X]

True

and your H list is

Table[h@i, {i, 2, Length@X}] == Most@Rest@FoldList[f, a0 + v, X]

True

Finally we have

f = Expand[a0 + a1 #2^2 + b1 #1] &;
X = RandomReal[NormalDistribution[0, 5], {10}];
Timing[H = Most@Rest@FoldList[f, a0 + Variance[X], X];
  logLike = Total[-Log[H] - (Drop[X, 1]^2/H)];
  NMinimize[{-logLike,
    a0 > 0 && 1 > a1 > 0 && 1 > b1 > 0 && a1 + b1 < 1}, {a0, a1, b1}]
  ]

{0.422, {38.2035, {a0 -> 0.310425, a1 -> 0., b1 -> 1.}}}

or

f = Expand[a0 + a1 #2 + b1 #1] &;
X = RandomReal[NormalDistribution[0, 5], {10}];
Timing[H = Most@Rest@FoldList[f, a0 + Variance[X], X^2];
  logLike = Total[-Log[H] - (Drop[X, 1]^2/H)];
  NMinimize[{-logLike,
    a0 > 0 && 1 > a1 > 0 && 1 > b1 > 0 && a1 + b1 < 1}, {a0, a1, b1}]
  ]

{0.282, {33.527, {a0 -> 4.17668, a1 -> -1.49002*10^-30,
    b1 -> 0.699952}}}

That's still pretty slow, however.

I'm not sure the following will work, but we may as well try it:

logLikely[x_List] :=
  Block[{a0, a1, b1, sq = x^2, f = Expand[a0 + a1 #2 + b1 #1] &, h},
   h = Most@Rest@FoldList[f, a0 + Variance[x], sq];
   Total[-Log[h] - Rest@sq/h]
   ]
nMin[x_List] := Block[{a0, a1, b1, logLike = logLikely@x},
   {a0, a1, b1} /.
    Last@NMinimize[{-logLike,
       a0 > 0 && a1 > 0 && b1 > 0 && a1 + b1 < 1}, {a0, a1, b1}]
   ]

Now I solve with n = 100 _and_ I also solve 5 lists of length 20.

X = RandomReal[NormalDistribution[0, 5], {100}];
Timing@nMin@X
Timing[nMin /@ Partition[X, 20]]

{43.938, {21.6771, 0.0205958, 5.5216*10^-8}}

{7.797, {{12.8, 0.210692, 0.373676}, {0.568711, 5.86338*10^-13,
    1.}, {8.11917, 0.734091, 1.28712*10^-8}, {17.6268, 9.57188*10^-9,
    9.69176*10^-8}, {19.6105, 3.46363*10^-8, 3.15567*10^-7}}}

It's much faster to solve several small problems, but the results are not  
very close to one another!!

Still, it's worth a try to average and use FindMinimum:

Clear[findMin]
findMin[x_List]/;Length@x>10&&Divisible[Length@x,10]:=Block[{a0,a1,b1,logLike,avg,start},
Print[Row@{"average: ",avg=Mean[nMin/@Partition[x,10]]}];
start=Transpose@{{a0,a1,b1},avg};
logLike=logLikely@x;
FindMinimum[{-logLike,a0>=0&&a1>=0&&b1>=0&&a1+b1<=1},start]
]

Timing[{t2,r3}=findMin@X]
s3={a0,a1,b1}/.r3;

average: {6.25198,0.12138,0.531861}
{10.968,{419.771,{a0->23.8944,a1->0.063187,b1->0.00296808}}}

Column@{s1,s3}

{23.8769,0.0630664,0.00374499}
{23.8944,0.063187,0.00296808}

The last answer (gotten in 11 seconds) is pretty close to the first answer  
(which took 46 seconds). That's great.

But notice... the AVERAGE was nowhere close, so what good did solving the  
smaller problems do??

So... the best strategy may be to find a point in the feasible region and  
use FindMinimum, never mind all the sub-problems.

Bobby

On Mon, 10 Dec 2007 19:35:13 -0600, <Arkadiusz.Majka at gmail.com> wrote:

> DearAll,
>
> I have a general problem with lenghty formula
>
> E.g for  a data set X
>
> X=Table[Random[NormalDistribution[0,5],{100}]];
>
> I define function h:
>
>
>  h[0] := Variance[X];
> h[1] := a0+ h[0];
> h[t_] := h[t] = a0 + a1*X[[t - 1]]^2 + b1*h[t - 1]
>
> and H which is the collection of h in all times up to the length of
> data set X
>
> H = Table[h[i], {i, 2, Length[X]}];
>
> Next I define a function:
>
> logLike = Total[-Log[H] - (Drop[X, 1]^2/H)];
>
> which I want to minimize in order to find a0,a1,b1
>
> NMinimize[{-logLike, a0> 0 && 1 >a1 > 0 && 1 > b1> 0 && a1+ b1 <1},
> {a0, a1, b1}]
>
>
> It works but takes a lot of time and does not scale at all (i.e. when
> I add an element to X the computational time increases very much)
>
>
> The form of logLike and NMinimize function is only an example of how
> to deal with very long function but with only a few arguments (here 3)
> to obtain fast results. Instead of NMinimize you can try to
> contourplot logLike function or do whatever you want - its time
> consuming.
>
> I know that there existing a solution ,i.e to speed up such
> computation. The example above is a parameter estimation of GARCH
> process which in TimeSeries packet runs very fast. But how did
> developers do so???
>
> Thanks for any help,
>
> Arek
>
>



-- 

DrMajorBob at bigfoot.com


  • Prev by Date: Re: FullSimplify with Pi
  • Next by Date: Re: efficient programming problem
  • Previous by thread: Re: how to speed up lenghty formulas??
  • Next by thread: FullSimplify with Pi