bugGNU Octave - Bugs: bug #34363, Poor performance for certain...

 
 

bug #34363: Poor performance for certain statistical distributions

Submitted by:  Rik <rik5>
Submitted on:  Fri 23 Sep 2011 04:50:54 AM UTC  
 
Category: LibrariesSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Performance
Status: Patch SubmittedAssigned to: None
Originator Name: Open/Closed: Open
Release: devOperating System: GNU/Linux

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Tue 18 Aug 2015 07:32:41 PM UTC, comment #11:

The two remaining functions which are critical to upgrade are

binoinv
nbininv

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>
Project 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>
Project 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

Michael Godfrey <godfrey>
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.

Rik <rik5>
Project 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>
Project 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>
Project 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.

binoinv
nbininv
poissinv

betainv
gaminv

hygecdf
hygeinv
hygernd

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

poissinv(1-eps,4)

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>
Project Administrator

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #34672:  poissinv4.m added by lachlan (6kB - text/x-objcsrc - fixed bug when X is a matrix)
file #34654:  poissinv3.m added by rik5 (6kB - text/x-objcsrc)
file #34655:  tst_poissinv.m added by rik5 (241B - text/x-objcsrc)
file #34618:  poissinv_newnew.m added by lachlan (6kB - text/x-objcsrc)
file #34596:  poissinv_new.m added by lachlan (6kB - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 6 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue 18 Aug 2015 06:46:59 AM UTClachlanAttached File-=>Added poissinv4.m, #34672
    Sun 16 Aug 2015 03:32:54 AM UTCrik5Attached File-=>Added poissinv3.m, #34654
      Attached File-=>Added tst_poissinv.m, #34655
    Tue 11 Aug 2015 06:58:50 PM UTCrik5StatusNone=>Patch Submitted
    Mon 10 Aug 2015 02:16:57 AM UTClachlanAttached File-=>Added poissinv_newnew.m, #34618
    Sat 08 Aug 2015 02:13:50 PM UTClachlanAttached File-=>Added poissinv_new.m, #34596

    Back to the top


    Powered by Savane 3.1-cleanup