bugGNU Octave - Bugs: bug #34363, [octave forge] (statistics) Poor...


bug #34363: [octave forge] (statistics) Poor performance for hypergeometric distributions

Submitter:  Rik <rik5>
Submitted:  Fri 23 Sep 2011 04:50:54 AM UTC
Category:  Octave Forge Package Severity:  2 - Minor
Priority:  3 - Low Item Group:  Performance
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Sat 06 Aug 2022 11:00:05 AM UTC, comment #19: 

Thanks Nicholas for the patch. Git commit 2f33c7c push at https://github.com/gnu-octave/statistics/commit/2f33c7cfdab5b9e31530d542438543e372ce7953
This bug can be closed.

Andreas Bertsatos <pr0m1th3as>
Fri 05 Aug 2022 08:57:16 PM UTC, comment #18: 

can only add 4 files at once. patch attached here.

(file #53530)

Nicholas Jankowski <nrjank>
Group Member
Fri 05 Aug 2022 08:55:59 PM UTC, comment #17: 

ok. functions vectorized. fairly straightforward. I've attached a git patch as well as the individual m files to this report.

the fixes rely on an optional "vectorexpand" parameter i added to hygepdf. hygecdf, hygeinv, hygernd all made the same use of a vector v = 0:n(i), which was the challenge identified in comment #14. that was simply worked around by generating a single v = 0:max(n) vector and then after calculating the hygepdf for each t,m,n set against that. then it was easy enough to only work with the portion needed for each n(i). the hygepdf option returns a 2D array where each pdf(i,j) is hygepdf(v(j), t(i),m(i),n(i)).

after that it was just a matter of re-implementing the equivalent of lookup to operate row-wise on the array (that would be a nice feature add to the current lookup), and selecting the correct item from the list.

if adding the optional hygepdf parameter is undesired, it's likely possible to incorporate those changes in each of the three files, or as an extra helper file, but I think incorporating it in hygepdf is cleaner.

BISTs were added to hygepdf for the new option, vectorized versions of the other functions still pass all tests. I didn't see any other tests needed on those.

quick speed tests show 3-10x improvement depending a lot on what's expanded.

old looped functions

tic; for idx = 1:100, hygecdf(50,1000,400,[1:200]);endfor, toc
Elapsed time is 10.6302 seconds.

tic; for idx = 1:100, hygeinv(0.5,1000,[200:400],100);endfor, toc
Elapsed time is 11.9475 seconds.

tic; for idx = 1:100, hygernd(1000,[200:400],100);endfor, toc
Elapsed time is 9.30343 seconds.

new functions

tic; for idx = 1:100, hygecdf(50,1000,400,[1:200]);endfor, toc
Elapsed time is 1.48031 seconds.

tic; for idx = 1:100, hygeinv(0.5,1000,[200:400],100);endfor, toc
Elapsed time is 1.26353 seconds.

tic; for idx = 1:100, hygernd(1000,[200:400],100);endfor, toc
Elapsed time is 1.3026 seconds.

(file #53526, file #53527, file #53528, file #53529)

Nicholas Jankowski <nrjank>
Group Member
Wed 03 Aug 2022 11:20:39 PM UTC, comment #16: 

Yes, it would be easier to make a patch or a pull request to github directly. Thanks

Andreas Bertsatos <pr0m1th3as>
Wed 03 Aug 2022 06:50:06 PM UTC, comment #15: 

after the other bug about hygepdf i've taken a look at vectorizing the last three functions here (hygernd, hygecdf, hygeinv). hygernd (with a vectorizing mod to hygpdf) wasn't too difficult, but the other two call it through discrete_cdf and discrete_inv (core octave functions) which may be a bit trickier. will see how it goes then upload changes. (can post as a git patch if that's easier for the new repo.)

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Feb 2016 04:08:51 AM UTC, comment #14: 

A brief look at the hygeXXX case convinced me that they won't be easy to vectorise efficiently.  It requires simultaneously working on PDFs of different lengths.  If each is made a column of a 2D array, then there can be excessive padding of the short columns.  If they are concatenated to a single long array, there can be catastrophic loss of precision when taking cumsum.

One approach would be to code these entirely in C++.  A less drastic change would be to implement a "cumsum with reset points" in C++ and use that.

All supplied patches have been submitted, so I'm changing the status back to "none", and making the summary specific to hypergeometric distributions.  I'm also tentatively lowering the severity and priority.

Lachlan Andrew <lachlan>
Mon 12 Oct 2015 04:25:58 AM UTC, comment #13: 

I reviewed and checked in binoinv.m (http://hg.savannah.gnu.org/hgweb/octave/rev/9d2023d1a63c) and nbininv.m (http://hg.savannah.gnu.org/hgweb/octave/rev/d6d04088ac9e).  I had to update some of the code to follow Octave coding conventions.  In particular, using a space after a function and before the opening parenthesis as in "sin (x)" rather than "sin(x)".  The new BIST tests also didn't run at first and required some debugging.

The speed-up for binoinv was 28X and an incredible 85X for nbininv.

For gaminv and betainv, I re-wrote these to use do-until loops with a maximum loop count of 40, rather than 10,000.  If the results ever do fail to converge within 40 cycles a warning will be issued.  But at least these won't "hang" anymore while looping from 40 to 10,000.

I'll leave this bug report open for the hygeXXX.m functions which remain.
As an example,

octave:13> x = rand (100, 1);
octave:14> tic; y = hygeinv (x, 100*ones (size (x)), 30, 20); toc
Elapsed time is 0.839279 seconds.
octave:15> tic; y2 = hygeinv (x, 100, 30, 20); toc
Elapsed time is 0.00971079 seconds.

which shows that the scalar case is much, much faster than the vector case.

Rik <rik5>
Group administrator
Tue 15 Sep 2015 11:28:47 AM UTC, comment #12: 

Here are drafts for new  binoinv.m  and nbininv.m.

In version 4.0.0,  binoinv.m  didn't have a while(1) loop and so was "safe".  It would run at most n iterations.  However, the current version is about 10 times faster, on my tests.

These are both probably overkill given how often they are used, but it was fun writing them :)

I chose not to use Gile's algorithm for these.  The quick read I had of the paper Rik found suggested that it is only useful for n>500.  That paper was by a summer student -- presumably an undergrad getting a feel for research -- and so I think it was a case of "Look! We can apply this technique more broadly", rather than aiming for something practical.

Regarding the other distributions, I think betainv.m and gaminv.m are doing a Newton search, which is pretty efficient.  It can't be vectorised because each iteration relies on the result of the previous.  I recommend we keep them as they are.  (They may have been updated since the original bug report.)

I don't understand the hygeXXX.m code well enough to comment.

(file #34889, file #34890)

Lachlan Andrew <lachlan>
Tue 18 Aug 2015 07:32:41 PM UTC, comment #11: 

The two remaining functions which are critical to upgrade are


Both of these use while (1) loops which may not converge. 

Interestingly, there seems to have been recent progress on this exact issue.  There is an academic project from 2013 which calculates new ways of calculating the inverse binomial function (http://www.lancaster.ac.uk/postgrad/moorhead/Work/SummerProject13.pdf).  What is doubly interesting, is that the research is based on Mike Gile's approximation to the inverse Poission distribution which Lachlan just implemented.  I didn't have time to read and digest the whole paper, but Appendix B contains Matlab code for the calculation already.  Presumably since this is an academic publication we might ask the author for permission to use and then just change the coding style to match Octave conventions.

Appendix C contains Matlab code for the general calculation of the inverse of a CDF.  If this is accurate, we might be able to extend this to all of the other missing inverse functions.

Rik <rik5>
Group administrator
Tue 18 Aug 2015 06:40:45 PM UTC, comment #10: 

@Lachlan: Thanks, the new version works perfectly.  Where the old and new implementations disagree, I prefer the new implementation  as it is closer to the expected value.

The test script I uploaded also had benchmarking code in it.  I find that for scalar lambda the speed-up is dramatic, ~6X the old implementation or -80% of the run time.  For lambda as an array the speed-up is smaller because there is still a for loop involved.  Nevertheless, the improvement is still -33%.

I checked in the changes here http://hg.savannah.gnu.org/hgweb/octave/rev/b3f56ed8d1f9.

Rik <rik5>
Group administrator
Tue 18 Aug 2015 09:21:13 AM UTC, comment #9: 

Lachlan: I looked through your code and just want
to congratulate you on a job well-done.

And, of course, is there any chance that you would
be willing, and have the time, to work on another
of the "needy" functions that Rik listed?

By the way, the documentation of this function in the
Manual is very brief. Feel free to add some detail if
you think that it would help.

In any case, thanks for your help.


Michael Godfrey <godfrey>
Group Member
Tue 18 Aug 2015 06:46:59 AM UTC, comment #8: 

Thanks for the extra info, Rik.

The bug is fixed.  The problem was that I was assuming  X  to be a vector, not matrix.  I used  if (any (k))  instead of  if (any (k(:))).

Previously, I only tested it with vectors.  I tried your script with a 100x100x100 case, and it gave eight out-by-one errors.  I've also tried it when the mean is a matrix and X is a scalar, and both are matrices, and both are three dimensional.

(file #34672)

Lachlan Andrew <lachlan>
Mon 17 Aug 2015 03:01:33 AM UTC, comment #7: 

@Lachlan: Here is a sample debugging session showing an off-by-2 error.  What's very interesting is that if poissinv is called standalone on the failing example it works.  This seems to indicate that there is some state involved.

octave:1> tst_poissinv
stopped in /home/rik/wip/Projects_Mine/octave-dev/scripts/statistics/distributions/tst_poissinv.m at line 16
16:     keyboard;
Failed at index 8
debug> find (y1 - y2)
ans =  1632
debug> y1(1632)
ans =  22
debug> y2(1632)
ans =  20
debug> format long
debug> x(1632)
ans =  0.999984299566334
debug> poissinv3 (x(1632), 8)
ans =  22

Rik <rik5>
Group administrator
Sun 16 Aug 2015 07:45:36 AM UTC, comment #6: 

Thanks, Rik.

I had already spotted discrepancies by 1.  I should have commented on them.  The are an off-by-epsilon error, as follows.

The Poisson distribution is discrete.  The answer is always an integer, and so if the answer is wrong at all, it will be wrong by an integer.

These functions are finding  n = argmin (F(n) < x).  The error occurs when  F(n)=x+epsilon  or  F(n)=x-epsilon, and the  n  we return is off by one.

This "epsilon" isn't the machine  eps.  It is the error associated with the asympmtotic approximation used by the algorithm, which will be a bit larger.

I'll check out the cases with a discrepancy of 2 when I get a chance.

Lachlan Andrew <lachlan>
Sun 16 Aug 2015 03:32:54 AM UTC, comment #5: 

I reviewed and made some minor changes to follow Octave coding guidelines.  Appendix C and D of the documentation have some of the guidelines we use.  Things I changed

1) Use '##' for comments that standalone on a line
2) Use '#' for comments that are appended to a line of code
3) Use parentheses around the test clause in an if statement
4) Use space between function name and opening parenthesis: 'find (x)' rather than 'find(x)'

I've uploaded the latest version as poissinv3.m.  I've also uploaded a file I was using to test the new implementation against the old one.  I found that they occasionally disagree by an integer (1 or 2) and not just a small number on the order of eps.  Could you investigate?  I'm pretty sure it's something small having to do with the cutover between the limit and the analytic approximation, but I can't check the new implementation in until I'm certain it is producing correct results.

(file #34654, file #34655)

Rik <rik5>
Group administrator
Sat 15 Aug 2015 04:10:11 AM UTC, comment #4: 

Thanks, Rik.  There's no rush to review the patch, and thanks for checking Google Scholar; I'm trying to help, not make work :)

Lachlan Andrew <lachlan>
Tue 11 Aug 2015 06:58:50 PM UTC, comment #3: 

Thank you Lachlan for the patch.  I hope I can review it within the next two weeks (awfully busy and all of Octave is volunteer-run).

Since Octave is used widely we generally try to write basic good code that will perform well under a range of circumstances.  This includes doing things like eliminating for loops which is good practice regardless of the particular value of lambda.  A direct answer to your question is that I don't know what value of lamda is more common.  I tried doing a search on scholar.google.com to sample a number of academic papers and just see what values of lamda are common in the literature, but it didn't work very well.

Rik <rik5>
Group administrator
Mon 10 Aug 2015 02:16:57 AM UTC, comment #2: 

Update: The file  poissinv_newnew.m  is an improvement on  poissinv_new.m with the following changes:
- For scalar  lambda,  it avoids the loop for small cases
- it is optimised for small lambda instead of large lambda, by doing the explicit summation before doing the asymptotic approximation.

It is about twice as fast as poissinv_new.m when the answer is less than 20, but slower when the answer is greater than 20 (but is still several orders of magnitude faster than poissinv.m in that case).

Do we have any usage data to say whether small or large  lambda  is more common?

Do we have any usage data to say which of these functions is used most often?

(file #34618)

Lachlan Andrew <lachlan>
Sat 08 Aug 2015 02:13:50 PM UTC, comment #1: 

The file  poissinv_new.m  is a replacement for  poissinv.m  that is much faster for large  lambda,  and has no infinite loop.

(file #34596)

Lachlan Andrew <lachlan>
Fri 23 Sep 2011 04:50:54 AM UTC, original submission:  

The following functions in scripts/statistics/distributions are implemented with loops which are slow in Octave.




The first three are problematic because they use a 'while (1)' construction and the stopping condition is not always met. For example,


will create an infinite loop that can only be broken with Ctrl-C.

At the very least, the 'while (1)' construction should be replaced with a for loop which will eventually terminate.

The best possible outcome would be to recode these entirely to get rid of the loops.

Rik <rik5>
Group administrator


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

Attach Files:

Attached Files
file #53530:  hygexxx-function-vectorization-bug-34363.patch added by nrjank (11KiB - application/octet-stream - git patch implementing hyge vectorizations)
file #53529:  hygernd.m added by nrjank (5KiB - text/plain - vectorized hygexxx functions)
file #53528:  hygeinv.m added by nrjank (5KiB - text/plain - vectorized hygexxx functions)
file #53526:  hygepdf.m added by nrjank (7KiB - text/plain - vectorized hygexxx functions)
file #53527:  hygecdf.m added by nrjank (4KiB - text/plain - vectorized hygexxx functions)
file #34889:  binoinv.m added by lachlan (7KiB - text/x-objcsrc)
file #34890:  nbininv.m added by lachlan (6KiB - text/x-objcsrc)
file #34672:  poissinv4.m added by lachlan (6KiB - text/x-objcsrc - fixed bug when X is a matrix)
file #34654:  poissinv3.m added by rik5 (6KiB - text/x-objcsrc)
file #34655:  tst_poissinv.m added by rik5 (241B - text/x-objcsrc)
file #34618:  poissinv_newnew.m added by lachlan (6KiB - text/x-objcsrc)
file #34596:  poissinv_new.m added by lachlan (6KiB - text/x-objcsrc)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by pr0m1th3as (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by nrjank
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by lachlan (Updated the item)
  • -email is unavailable- added by rik5 (Submitted the item)

    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 project members can vote.


    Follow 25 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-08-06 nrjank StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2022-08-05 nrjank Attached File- Added hygexxx-function-vectorization-bug-34363.patch, #53530
        StatusNone Patch Submitted
        Assigned tonrjank None
    2022-08-05 nrjank Attached File- Added hygepdf.m, #53526
        Attached File- Added hygecdf.m, #53527
        Attached File- Added hygeinv.m, #53528
        Attached File- Added hygernd.m, #53529
    2022-08-03 nrjank Assigned toNone nrjank
        Carbon-Copy- Added pr0m1th3as
    2019-06-04 mtmiller SummaryPoor performance for hypergeometric distributions [octave forge] (statistics) Poor performance for hypergeometric distributions
    2019-06-01 rik5 CategoryLibraries Octave Forge Package
    2016-02-08 lachlan Severity3 - Normal 2 - Minor
        Priority5 - Normal 3 - Low
        StatusPatch Submitted None
        Operating SystemGNU/Linux Any
        SummaryPoor performance for certain statistical distributions Poor performance for hypergeometric distributions
    2015-09-15 lachlan Attached File- Added binoinv.m, #34889
        Attached File- Added nbininv.m, #34890
    2015-08-18 lachlan Attached File- Added poissinv4.m, #34672
    2015-08-16 rik5 Attached File- Added poissinv3.m, #34654
        Attached File- Added tst_poissinv.m, #34655
    2015-08-11 rik5 StatusNone Patch Submitted
    2015-08-10 lachlan Attached File- Added poissinv_newnew.m, #34618

    Back to the top

    Powered by Savane 3.12