bugGNU Octave - Bugs: bug #65230, Poor accuracy of betaln for high...


bug #65230: Poor accuracy of betaln for high ratios of a/b

Submitter:  None
Submitted:  Tue 30 Jan 2024 04:38:42 PM UTC
Category:  Octave Function Severity:  2 - Minor
Priority:  3 - Low Item Group:  Inaccurate Result
Status:  Confirmed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Tue 13 Feb 2024 07:46:05 PM UTC, comment #8: 

I'm fine with just having the generic warning about incompatible broadcasting arguments be emitted.

Rik <rik5>
Group administrator
Tue 13 Feb 2024 10:33:37 AM UTC, comment #7: 

Should we still keep input validation to check that the input dimensions allow broadcasting? Or are the generic error messages good enough in this case?

Markus M├╝tzel <mmuetzel>
Group administrator
Sun 11 Feb 2024 02:56:21 AM UTC, comment #6: 

On point #2, you're quite right.  There was no reason to be so particular about the inputs.  It was so easy I immediately fixed that in this changeset https://hg.savannah.gnu.org/hgweb/octave/rev/4fb466fb717e.

On point #3, we don't necessarily have to follow Matlab, but their betaln function only accepts real, nonnegative arguments.  Perhaps the easiest thing to do would be to keep that restriction.  As far as I know, there haven't been bug reports or requests for extending this function to accept complex values.  The current implementation of betaln in Octave uses gammaln which accepts negative arguments, but not complex ones.  So there is no problem with nonnegative arguments.

On point #1, I think the answer is generally that more accuracy is certainly desirable as long as the trade-off in computational time isn't too excessive.  I benchmarked the current implementation with

octave:1> x = 1:1000;
octave:2> y = x';
octave:3> tic; z = betaln (x,y); toc
Elapsed time is 0.020807 seconds.

This calculates 1e6 values and took 20.8 milliseconds which is effectively instantaneous compared to human reaction times.  I would think even 100 milliseconds would be a worth trade-off.

Rik <rik5>
Group administrator
Sun 11 Feb 2024 12:37:26 AM UTC, comment #5: 

OP here again.  Thank you for the suggestions involving the symbolic package.

> The accuracy is limited not by the gammaln() function, but by the addition operation using IEEE 8-byte doubles.  eps() is 2e-16 and if the difference in scale between two arguments is greater than 1e16 then the smaller argument is lost.

The accuracy of betaln(a, b) can be poor even if the addition operation a + b is exact.  For example betaln(2^62, 2^10) returns -32768 for me, but the actual value is about -37935.2477864453572569.

I am interested in the possibility of contributing a version of betaln.m that improves accuracy at some cost to speed and simplicity.  I have a few questions:

  1. How likely is this to be accepted?  Assume for the sake of argument that accuracy is at least 40 bits for the vast majority of the input space, and that speed and simplicity are about 10 times worse than the existing implementation.
  2. The existing implementation has input validation that requires the arguments to have equal size or at least one of them to be scalar.  Is there any reason to retain this and not fully support broadcasting?
  3. Some definitions of the beta function seem to require (the real parts of) both arguments to be positive, but others relax this requirement.  What should betaln return with one or more non-positive arguments?

Any input would be appreciated.  Thanks.

Wed 07 Feb 2024 01:22:42 AM UTC, comment #4: 


pkg load symbolic

n = sym('2') ^ 128
k = sym('30')

double (log (nchoosek (n, k)))

But Octave is numerical software not mathematical software.

Fri 02 Feb 2024 10:51:03 PM UTC, comment #3: 

I've never seen input values to nchoosek like 2^128.  I think that is well outside the norm which is probably why this hasn't been noticed before.  I'm not saying your application doesn't call for these numbers, just that overwhelming number of use cases don't call for these values.

Generally nchoosek blows up very quickly such that it isn't conceivable to try the direct approach of

log (nchoosek (n, k))

Something simple like

nchoosek (1e6,4)
warning: nchoosek: possible loss of precision
warning: called from
    nchoosek at line 151 column 7

ans = 4.1666e+22

is already complaining about lost precision.

Back to betaln:

lnb = gammaln (a) + gammaln (b) - gammaln (a + b)

The accuracy is limited not by the gammaln() function, but by the addition operation using IEEE 8-byte doubles.  eps() is 2e-16 and if the difference in scale between two arguments is greater than 1e16 then the smaller argument is lost.

octave:19> format long
octave:20> 1 + 1e-16
ans = 1
octave:21> sprintf ('%.17d\n', 1 + 2e-16)
ans = 1.0000000000000002

I would go back to the base definition of nchoosek as

c = n! / k!(n-k)!

Take the logarithm of both sides and use the properties of logarithms to write

log (c) = log (n!) - log (k!) - log ((n-k)!)

Then using the relationship

gamma (n + 1) = n!

I can write

log (c) = gammaln (n+1) - gammaln (k+1) - gammaln (n-k+1)

A quick anonymous function for this is

bignk = @(n,k) gammaln (n+1) - gammaln (k+1) - gammaln (n-k+1)

This still has a problem in that the expression "N-K+1" can't be calculated without loss of precision with IEEE doubles.  So, I tried switching to variable precision numbers which are part of the symbolic package.  You will need to have python3-sympy installed on your computer.  After that, install and load the symbolic package with

pkg install -forge symbolic
pkg load symbolic

Finally, create some numbers and try it out.

octave:47> n = vpa (2, 50)
n = (sym) 2.0000000000000000000000000000000000000000000000000
octave:48> n = n^128
n = (sym) 340282366920938463463374607431768211456.00000000000
octave:49> k = vpa (2, 50)
k = (sym) 2.0000000000000000000000000000000000000000000000000
octave:50> k = k^32
k = (sym) 4294967296.0000000000000000000000000000000000000000
octave:51> logc = bignk (n,k)
logc = (sym) 290091236578.66962544142734259366989135742187500000

Maybe this is what you are looking for?

Rik <rik5>
Group administrator
Thu 01 Feb 2024 03:14:24 AM UTC, comment #2: 

OP here.  Thank you for confirming, and for reporting the behavior of Matlab.  I see that the Matlab documentation also doesn't explicitly mention accuracy, but the "Algorithms" section does seem to indicate that the same implementation based on gammaln is used.

I don't really know what the logarithm of the beta function is commonly used for by others, but betaln has returned 0 for every real case I've tried.  For example, I needed to compute log(nchoosek(n, k)) where (n, k) was something like (2^128, 2^30) or (2^256, 2^40).  I tried

function [retval] = log_nchoosek (n, k)
  retval = -log1p(n) - betaln(k + 1, n - k + 1);

which doesn't work for the sorts of values I was interested in.  Replacing betaln with gsl_sf_lnbeta seemed to do the job.

Regarding this being a corner case, those are mighty big corners!  Of the input space where both arguments are greater than 1, betaln seems to return 0 for something like 88%, NaN for another 2%, and arguably has poor accuracy for most of the rest.  Still, I suppose I know what you mean.  This issue appears to have existed for over 17 years without much complaint.

If the current implementation is retained, I think it would be helpful for the documentation to disclose the implementation and mention the accuracy implications.

Wed 31 Jan 2024 11:22:07 PM UTC, comment #1: 

I can confirm that at high ratios of a/b the betaln() function always returns zero.  The cause is straightforward.  In betaln.m the output is calculated as

  lnb = gammaln (a) + gammaln (b) - gammaln (a + b);

When a and b are different enough in size the finite precision of IEEE floating point causes a problem and the sum "a + b" is just equal to "a".

However, it's not clear to me that this is a frequent occurrence.  For example, do constraints from real-world problems limit the sizes of a, b, and their ratio?  I would also observe that Matlab exhibits exactly the same behavior as Octave.  Maybe they are unaware of this corner case, or maybe they consider it unlikely.  I also note that Mathematica doesn't even offer a log beta function.  The help discussions simply point to writing such a function in terms of gammaln in exactly the way Octave has done.

Rik <rik5>
Group administrator
Tue 30 Jan 2024 04:38:42 PM UTC, original submission:  

The betaln(a, b) function is currently implemented as gammaln (a) + gammaln (b) - gammaln (a + b).  This seems to have reasonable accuracy when the ratio a/b is near 1, but all accuracy is lost due to catastrophic cancellation as the ratio becomes more extreme.  For most inputs, such as betaln(2^62, 2^6), 0 is returned.

My current workaround is to use the gsl_sf_lnbeta function from the gsl package.  This function can provide an error estimate in its second output argument.  The error estimates seem reasonably small for most inputs I've tried.



(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.


    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-02-03 rik5 Priority5 - Normal 3 - Low
    2024-01-31 rik5 Severity3 - Normal 2 - Minor
        StatusNone Confirmed
        SummaryPoor accuracy of betaln for most inputs due to catastrophic cancellation Poor accuracy of betaln for high ratios of a/b

    Back to the top

    Powered by Savane 3.12.
    Corresponding source code