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

[Rik]

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

Agree with most of this.

Specifically, Octave is for numerical computation, not pure math.  If F(x) and
F(x+delta) yield the same answer there is no point distinguishing between
them.  This concept still trips up a lot of users who expect either infinite
precision (pure math) or exact Matlab-equivalency.

Second, I agree that algorithm should be agnostic about singles/doubles and
perform all operations with the desired precision.  This is how other Octave
functions that support mixed inputs behave.

I do think df_new calculation can be simplified/eliminated.  The calculation


df_new = (x + step) - x;


will be 0 if step is less than eps (x), but otherwise will be equal to step. 
In this case, x is restricted to the range (0, 1) so using a tolerance of eps
(class (x)) as a check on step is a sufficient stopping condition.

The existing code is


    step = -F (x(todo), a(todo), b(todo), y(todo)) ./ ...
       JF (x(todo), a(todo), b(todo), Bln(todo));
    df_new = (x(todo) + step) - x(todo);
    x(todo) += step;
    ind = df_new .* df > 0;
    todo = todo(ind);
    df = df_new(ind);


and I propose


    step = -F (x(todo), a(todo), b(todo), y(todo)) ./ ...
       JF (x(todo), a(todo), b(todo), Bln(todo));
    x(todo) += step;
    ind = (step > tol);
    todo = todo(ind);


which gets rid of both the df and df_new variables.  It also seems to have a
performance benefit for certain values.  Here is a run of the current code
against a particularly hard value.


octave:60> badx
badx = 5.717842473940138e-07
octave:61> newbeta (badx, 2, 4)
iter: 1, # todo: 1, x: 0.0759261969940135, max (df): 0.174073803005986
iter: 2, # todo: 1, x: 0.0347141003450014, max (df): 0.0412120966490121
iter: 3, # todo: 1, x: 0.0167224581112431, max (df): 0.0179916422337583
iter: 4, # todo: 1, x: 0.00821961591239634, max (df): 0.00850284219884681
iter: 5, # todo: 1, x: 0.0040791709326122, max (df): 0.00414044497978414
iter: 6, # todo: 1, x: 0.00203830950551201, max (df): 0.00204086142710019
iter: 7, # todo: 1, x: 0.00103118309152024, max (df): 0.00100712641399177
iter: 8, # todo: 1, x: 0.000542869669144377, max (df): 0.000488313422375863
iter: 9, # todo: 1, x: 0.000324036344981909, max (df): 0.000218833324162468
iter: 10, # todo: 1, x: 0.000250279883120058, max (df): 7.3756461861851e-05
iter: 11, # todo: 1, x: 0.00023942338681535, max (df): 1.08564963047085e-05
iter: 12, # todo: 1, x: 0.000239177433992639, max (df): 2.45952822710551e-07
iter: 13, # todo: 1, x: 0.000239177307623453, max (df): 1.26369186572188e-10
iter: 14, # todo: 1, x: 0.000239177307623419, max (df): 3.31765864780564e-17
iter: 15, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 16, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 17, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 18, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 19, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 20, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 21, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 22, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 23, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 24, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 25, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 26, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 27, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
iter: 28, # todo: 1, x: 0.000239177307623419, max (df): 2.71050543121376e-20
ans = 2.391773076234194e-04


The algorithm spends iterations 14 - 28 changing the guess by an amount less
than eps ("double").  When I use the criteria I suggested the algorithm stops
after the 13th iteration.


octave:62> biinv3 (badx, 2, 4)
iter: 1, # todo: 1, x: 0.0759261969940135, max (df): 0.174073803005986
iter: 2, # todo: 1, x: 0.0347141003450014, max (df): 0.0412120966490121
iter: 3, # todo: 1, x: 0.0167224581112431, max (df): 0.0179916422337583
iter: 4, # todo: 1, x: 0.00821961591239634, max (df): 0.00850284219884681
iter: 5, # todo: 1, x: 0.0040791709326122, max (df): 0.00414044497978414
iter: 6, # todo: 1, x: 0.00203830950551201, max (df): 0.00204086142710019
iter: 7, # todo: 1, x: 0.00103118309152024, max (df): 0.00100712641399177
iter: 8, # todo: 1, x: 0.000542869669144377, max (df): 0.000488313422375863
iter: 9, # todo: 1, x: 0.000324036344981909, max (df): 0.000218833324162468
iter: 10, # todo: 1, x: 0.000250279883120058, max (df): 7.3756461861851e-05
iter: 11, # todo: 1, x: 0.00023942338681535, max (df): 1.08564963047085e-05
iter: 12, # todo: 1, x: 0.000239177433992639, max (df): 2.45952822710538e-07
iter: 13, # todo: 1, x: 0.000239177307623453, max (df): 1.26369186577854e-10
ans = 2.391773076234193e-04



    _______________________________________________________

Reply to this item at:

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

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