[Octave-bug-tracker] [bug #60539] Slow performance of betaincinv.m

[Michael Leitner]

Follow-up Comment #12, bug #60539 (project octave):

betaincinv.m, both in the new and old versions, just inverts betainc.m, that
is, it finds an x for given y so that the given implementation of betainc in
octave return y for given x. Thus, for judging whether betaincinv is working
well, you should not compare to the previous implementation or to Matlab or
Wolfram alpha, but test what betainc returns on betaincinv's results. 

What I see here is that if I say 


x=betaincinv (single (1-1e-6), 2, 4)


I get x=0.97864, and this is a single, as required. If I then say 


betainc(x,2,4)==single(1-1e-6)


I get logical true, meaning that the betaincinv really returns the best result
it could give, a single value for which betainc returns exactly the y (in
single precision) that was initially provided -- it is as accurate as in
principle can be expected (at least for this specific input). Note that this
is not hard in this case, as with the high beta of four betainc is very flat
around 1, and you have to go quite far away from 0.97864 so that the result of
betainc changes -- if you take Matlab's value


betainc(single(0.9786913),2,4)==single(1-1e-6)


you get true as well. If beta was smaller than 1, this would be the other way
round, then it would be possible that there is just no x so that the given y
can be returned by betainc (of course, this pertains both to single and double
precision). In any case, that's the test you would have to do: give any
initial y, test whether the resulting x fulfills betainc(x)==y, and if not,
see whether if you add or subtract eps(x) to x you get a closer y. If my
reasoning in implementation was correct, this can only happen if betainc is
very inaccurate (specifically non-monotonous).

I don't know whether the exact agreement between Matlab and the old algorithm
was just a fluke, whether both were modelled after a common ancestor, or
whether (God forbid!) one actually copied from the other with their
incompatible licences (we are forbidden from looking into Matlab's source code
to decide that). 

Of course, you can get more accurate single results if you do the computation
in double and at the end convert to single. I actually intended to put this up
for discussion when I posted the new algorithm, but: that is something the
user can always do, just give double instead of single arguments to
betaincinv. On the other hand, there is a reason for computing in single -- it
uses half of the (temporary) memory. That's why I thought it was intended as
it was (namely to compute in single if any input argument was single). In
particular, also the old algorithm did the computation in single precision in
this case. You could force everything to double and convert at the end, but as
I said, the users can do that themselves, and they lose the option to do the
computation in memory-tight situations.

As to df_new: yes, in pen-and-paper arithmetics df_new is equal to step. In
finite precision arithmetics it is not necessarily (in fact, it is quite
unlikely that it would be equal). The point is that the break condition is if
a given point does not get updated in a given step (then it won't be updated
in the next and all following, so you have to break there), or if the step
direction reverses (this can only happen due to numerical inaccuracies of
betainc). So for the former case, it can be that the computed step is smaller
than eps(x), in which case x will stay the same and df_new is zero, but step
is not. That's why I do it like that. Indeed, there is perhaps a possibility
for more efficiency here, in that x+step is computed twice. One could say the
following


df_new = -x(todo);
x(todo) += step;
df_new += x(todo);


which may or may not be faster (depending on cache misses and so on).

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60539>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/