On Sep 9, 5:27 pm, Robert Bradshaw <[email protected]>
wrote:
> On Wed, Sep 8, 2010 at 6:39 PM, KvS <[email protected]> wrote:
> > Thanks all, however I am not very successful so far :(.
>
> > I tried both options mentioned before:
>
> > - only optimize the loops in Cython and keep using symbolic
> > expressions/infinite precision, but this is unfortunately rather
> > slow;
>
> What do you mean by symbolic? Are you actually using maxima and the
> various calculus modules?
By symbolic/infinite precision I just mean keeping quantities like say
log(2) as they are, rather than approximating them by some finite
precision real number, sorry for the confusion. No, I am not using
Maxima or anything more advanced than the elementary operations, the
exponential function and the find_root function.
>
>
>
> > - fully optimize in Cython by turning to doubles everywhere, although
> > it gets very fast here happens what I was already afraid of: the
> > algorithm goes belly up after 15 steps :( (I would need a lot more
> > steps). In step 14 values like 2.51885860e-34 appear and in
> > combination with the numerical noise gathered this turns out to
> > produce very wrong values at the next step.
>
> > @Robert: I'm trying to implement an algorithm that computes a sequence
> > of functions v_k according to the rule v_k(x) = \max \{ K-e^x, \int
> > v_{k-1}(y) p(x+y) dy \}, v_0(x)=K-e^x, where p is a prob. density. The
> > background is an optimal stopping problem in stochastics. So at each
> > step I essentially first compute the integral and then determine where
> > the integral and x -> K-e^x intersect, this just by the basic
> > find_root function in Sage.
>
> > It's not difficult (however very annoying...) to write down a formula
> > for the integral (given the specific density p I'm using) and work out
> > formulae for all the coefficients that appear in this formula for the
> > integral in terms of those appearing in v_{k-1}. The number of
> > coefficients (and hence computations) involved in the formula for v_k
> > increases very quickly with k, that explains why the algorithm gets
> > slow in 'symbolic mode'. However 'double mode' is not good enough as
> > mentioned above.
>
> It sounds like you're trying too extreems--how about using your
> "double mode" algorithm, but instead using elements of RealField(N)
> where N > 53. This should be much faster than "symbolic" (if I
> understand you right, calling find_root and integrate) but have higher
> precision than using raw doubles. This may be enough, without getting
> your hands on MPFR directly yourself.
(I don't need to do any integration in the code, sorry if my previous
post seemed to suggest that, the integral that occurs there I have
worked out completely in terms of elementary operations.) Yes, I have
tried your suggestion of using RealField (without Cython though) and
this works good in the sense that the added precision in comparison
with floats makes the algorithm produce reliable output, thanks!
However, it is still too slow. If I could further speed it up by a few
factors that would be great. I will further try if I can figure out
how to use Cython to speed up parts of the algorithm where a lot of
looping and computing is done.
>
> You might also consider looking at how the dickman rho function is
> implementedhttp://hg.sagemath.org/sage-main/file/b5dab6864f35/sage/functions/tra...
> which sounds similar in spirit.
>
> > I will try your suggestion for using mpfr, thanks Robert and Greg!
> > (however it looks a bit scary to me ;)). I need to store the huge
> > amounts of coefficients somewhere, I guess numpy is my best guess even
> > though I have to declare the dtype as object then, no?
>
> What do you mean by huge? Would a list suffice? Or if it's a
> polynomial, what about a polynomial in RR? For high precision, the
> per-element overhead is constant, so I don't see an advantage to using
> NumPy here just as a storage divice if the dtype is object.
The function v_k is not a polynomial but a piecewise function
consisting of linear combinations of terms a*x^b*exp(c*x) (*) in each
part of its domain. 'Huge' means that the number of coefficients (i.e.
all the a, b and c's in (*)) to be computed in the k-th step of the
algorithm roughly grows like k^2. These coefficients are currently
organized in 4-dimensional Numpy-arrays: the k runs along one axis and
there are three other 'counters' (one to keep track of the part of the
domain etc.). I would prefer not to give up this way of representing
the coefficients as my notes use this form as well. Would you say I
could better use 4-d Python lists rather?
Thanks again, Kees
>
> - Robert
--
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
