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