Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updated Lambert W function implementation #371

Open
wants to merge 24 commits into
base: master
Choose a base branch
from

Conversation

alyst
Copy link

@alyst alyst commented Dec 31, 2021

This is a fork of #84. I did:

@codecov
Copy link

codecov bot commented Dec 31, 2021

Codecov Report

Patch coverage: 96.74% and project coverage change: -0.07 ⚠️

Comparison is base (07a890c) 94.25% compared to head (ed3c825) 94.19%.

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #371      +/-   ##
==========================================
- Coverage   94.25%   94.19%   -0.07%     
==========================================
  Files          14       15       +1     
  Lines        2908     3031     +123     
==========================================
+ Hits         2741     2855     +114     
- Misses        167      176       +9     
Flag Coverage Δ
unittests 94.19% <96.74%> (-0.07%) ⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files Coverage Δ
src/SpecialFunctions.jl 100.00% <ø> (ø)
src/lambertw.jl 96.74% <96.74%> (ø)
src/beta_inc.jl 91.46% <0.00%> (-0.75%) ⬇️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

@ViralBShah
Copy link
Member

cc @jlapeyre

@CarloLucibello
Copy link

Can this be merged? Seems well tested and well documented

@devmotion
Copy link
Member

No, I don't think it should be merged. My impression from the discussion in #84 was that there was a tendency/agreement to keep the function in a separate package.

@alyst
Copy link
Author

alyst commented Feb 28, 2022

At this point I wonder if you are interested to improve the implementation of this function at all.

@devmotion
Copy link
Member

There's a package with this function, and the general opinion (not just my personal one) was to keep the function there. So if it can be improved, I guess you should open a PR in that repo.

@alyst
Copy link
Author

alyst commented Feb 28, 2022

The "general opinion" is a bit of a stretch. But thank you for your efforts to guard the quality of Julia ecosystem.

As for me, I think trying to put this code elsewhere would be a waste of my resources, given that nobody have reviewed it here (I have noted "if it can be improved").

@ViralBShah
Copy link
Member

@alyst thanks for putting it together. It is generally true that it has been hard to get reviews here.

@devmotion
Copy link
Member

I'm sorry for you that the PR is not reviewed, but I guess many people, including myself, are not eager to review a PR in detail if there's already a longer discussion that indicates that there is no general agreement that the Lambert W function should be added to SpecialFunctions but instead shows that people seem to be in favour of keeping it in LambertW.jl. Therefore my suggestion was to make a PR to LambertW.jl instead if you would like to improve its implementation. I think you will receive more reviews and comments there.

@ViralBShah
Copy link
Member

ViralBShah commented Feb 28, 2022

It did seem that the original author @jlapeyre (IIRC, and if I have my facts straight) expressed interest in keeping it separate - and thus I feel that we should respect that opinion.

@ViralBShah
Copy link
Member

Since I was the one who triggered the original thought of getting this in - I would like to apologize for causing wasted effort here. :-(

@jlapeyre
Copy link

jlapeyre commented Mar 1, 2022

I'm also sorry if I caused wasted effort. Like many people, I'm over-extended and this is not the highest priority. I was also reluctant to put effort into a PR because I made one (for the same question) a few years ago, and it sat around till it bitrotted. I'm not blaming anyone: Most of us are not paid for this and have other obligations. But, my best estimation was that the best first step is to put LambertW under the JuliaMath org. I asked how to do that in the appropriate (I think) channel in slack a few weeks ago. But, I never got a response.

@ViralBShah
Copy link
Member

ViralBShah commented Mar 1, 2022

@jlapeyre I can help with moving LambertW.jl to the JuliaMath org. I'm inviting you as an owner of JuliaMath, and you should be able to then transfer the package.


### root finding, iterative solution

# Use Halley's root-finding method to find
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Worth noting that there's a paper from Lóczi published last year in Applied Mathematics and Computation called "Guaranteed- and high-precision evaluation of the Lambert W function" that describes algorithms that improve significantly on e.g. Halley's method. I used it to implement W-1 for real arguments in (-1/e, 0) in https://github.com/JuliaStats/Distributions.jl/blob/master/src/univariate/continuous/lindley.jl#L137-L156.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've actually tried the iteration scheme evaluated in the paper, but it looks like it's on par with Halley or the latter requires 1-2 iterations less from the same starting point (but I have not tested very extensively).
IIUC the paper doesn't claim that the Iacono-Boyd method provides faster convergence than Halley method, it just proves the bounds for the Iacono-Boyd one (because it's a simpler formula).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah good to know, it's been a bit since I read it so the details have been largely GC'd. Probably should have reread the paper before attempting to summarize it. 😅 One thing I recall was that when I was testing things out for the -1 branch in (-1/e, 0), I was comparing LambertW.jl vs. Iacono-Boyd vs. Mathematica as "ground truth" and found that Iacono-Boyd produced more accurate results near -1/e and 0 but was otherwise the same. Your comment makes me realize that I had actually neglected to check the number of iterations that LambertW.jl required for convergence. LambertW.jl explicitly keeps track of the change between iterations and will stop early based on that, whereas my naive implementation of Iacono-Boyd just iterates a set number of times. In developing it, I was playing around with the number of iterations and found that most inputs required only 2 or so iterations to converge, so I'd be interested to see how that compares with LambertW.jl.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've tried to implement Iacono-Boyd scheme here: https://github.com/alyst/SpecialFunctions.jl/tree/lambertw_iacono_iters.
Actually, the iteration scheme is just an input parameter to the root finding function, so that code could be easily adopted for more comprehensive evaluation of efficiency of both schemes over the whole LambertW(x) domain and for different numeric types.

@jlapeyre
Copy link

The separate package has been moved here https://github.com/JuliaMath/LambertW.jl

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants