TuringLang / NestedSamplers.jl

Implementations of single and multi-ellipsoid nested sampling
https://turinglang.org/NestedSamplers.jl/stable/
MIT License
42 stars 8 forks source link

use relative eigen threshold #81

Open aplavin opened 2 years ago

aplavin commented 2 years ago

upd: this post describes background and the first commit of this PR

It's probably more reasonable? But still not a complete fix, see below.

I encountered an error arising from the A matrix here having (kinda) negative eigenvalues: https://github.com/TuringLang/NestedSamplers.jl/blob/eb6999e0da004524a5848eb2c99b1f53e8a57f1d/src/bounds/ellipsoid.jl#L108 This matrix is passed to the Ellipse constructor, and fails in sqrt at https://github.com/TuringLang/NestedSamplers.jl/blob/eb6999e0da004524a5848eb2c99b1f53e8a57f1d/src/bounds/ellipsoid.jl#L49 Actually, eigenvalues of A are not negative as-is, but become negative in Symmetric(A). Digging into this error, I found that the https://github.com/TuringLang/NestedSamplers.jl/blob/eb6999e0da004524a5848eb2c99b1f53e8a57f1d/src/bounds/ellipsoid.jl#L168-L178 function applies a huge correction to the cov matrix in my case at https://github.com/TuringLang/NestedSamplers.jl/blob/eb6999e0da004524a5848eb2c99b1f53e8a57f1d/src/bounds/ellipsoid.jl#L105 Not sure why it is so, and I don't understand the make_eigvals_positive! function, what exactly should it do and why? Maybe, adapting the function https://github.com/joshspeagle/dynesty/blob/2c5f1bbe5a0745c6625876f23ec6aa710c845fd4/py/dynesty/bounding.py#L1230 from dynesty would be better?

This PR fixes the error with the particular matrix (below) because it has the condition number < 1e10. But I don't think it's a proper solution.

To reproduce my original error in fit():

x = [0.08191723085341904 0.13545261362441807 0.15338757993573843 0.08837059081152063 0.08728351269161129 0.08976317416597561 0.08761130920214773 0.08923971980362749 0.10418815131134627 0.08992877418862986 0.09147367077861655 0.09316510988605262 0.09563856390024207 0.08817882872618149 0.08642233931372963 0.07606003971830035 0.08600136293862232 0.09001296711562165 0.08306752534109188 0.08743590459596759; 0.612801232286909 0.714270171203782 0.7432347258487771 0.6171628893505378 0.6146034675622163 0.618426961051714 0.6153771179867582 0.618237345987351 0.6433971878217596 0.6167304187461977 0.6230179763674526 0.6241807458756284 0.6289414291480476 0.6155424785652647 0.6135085496409891 0.6092673979927283 0.6118844405903802 0.6206484791282253 0.6152761860541198 0.616198978017134; 0.7383544297901901 0.7354965419762426 0.7428745987891984 0.7327667778877038 0.7339120753059504 0.7387398098330907 0.7354058136867483 0.7381832299207463 0.7413389235114898 0.742224279223494 0.7348810278237009 0.7436972298155209 0.7364770089790889 0.7383629710796894 0.7413569492327434 0.7362735912601804 0.7336810546823396 0.7352224953280587 0.7410045201586068 0.7399545331614852; 0.1612003103669555 0.17894146736390518 0.17029641990856784 0.16778130063246074 0.1638315675787379 0.16416912251875704 0.1654360888624566 0.1676542140972558 0.17099877889821316 0.16798666146800578 0.16456076899222052 0.16939681645739854 0.1645159868982537 0.16184076817446968 0.16408355461735 0.1678911971274453 0.16279537950489986 0.16982106496996155 0.16446940610576574 0.16910125322597902; 0.3015394069935417 0.30026425694930864 0.2994049623941189 0.3010710663986933 0.301071146149588 0.3010033751883286 0.3010944773251966 0.30106316849967435 0.30083348313915326 0.30099797021522834 0.30111934488216846 0.30081150727006567 0.3009656556536787 0.3010426343916504 0.30126176295770074 0.30065104712971746 0.3010888432022384 0.3009836404276063 0.3015359771562088 0.3011239551672445; 0.5549266001517701 0.5588638010499112 0.5640153323753487 0.5491925108916669 0.550612494463565 0.5519747776187889 0.5504569238284182 0.5509517066765967 0.552547020571869 0.5503630125806462 0.5515650101461752 0.551034543383568 0.5504353956960134 0.5524011219734949 0.5517674292812791 0.5484217147622831 0.5500281378739291 0.5492799977097876 0.5545623488063589 0.5495472914747409; 0.5008162500177611 0.5005795652908372 0.5005407968494694 0.5007867594557819 0.5007961883413762 0.5007973556753138 0.5007956461849565 0.5007938403571855 0.5007215310985919 0.500789895035134 0.5007714776575859 0.500786916916394 0.5007457448675278 0.5007975364393731 0.5008094040718447 0.5009042631669718 0.5007997994974079 0.5007858508689621 0.500813250005897 0.5008008595761184; 0.49908223425184267 0.49941058059052995 0.49938214318498764 0.4991085755098695 0.49908581311112277 0.4990930025933883 0.49909752581993566 0.49910926460063043 0.49918939209707686 0.4991116032909387 0.49913072391024976 0.4991138899112875 0.49913976610253785 0.49908712776257685 0.4991476982112339 0.4989846710998443 0.49909690598505785 0.49911684472308426 0.4991024180081014 0.4991380573498765];
using NestedSamplers
Bounds.fit(Bounds.Ellipsoid, x, pointvol=1e-13)
codecov[bot] commented 2 years ago

Codecov Report

Merging #81 (3e6b6ea) into main (eb6999e) will increase coverage by 0.09%. The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##             main      #81      +/-   ##
==========================================
+ Coverage   88.62%   88.72%   +0.09%     
==========================================
  Files          13       13              
  Lines         554      550       -4     
==========================================
- Hits          491      488       -3     
+ Misses         63       62       -1     
Impacted Files Coverage Δ
src/bounds/ellipsoid.jl 85.88% <100.00%> (+0.48%) :arrow_up:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update eb6999e...3e6b6ea. Read the comment docs.

aplavin commented 2 years ago

Pushed the second commit that tries to make a proper fix.

mileslucas commented 2 years ago

Hey @aplavin thanks for the work on this, sorry I haven't had time to respond to any of this yet, I'm swamped under my thesis work nowadays. I've been wanting to rewrite the bounding space algorithms for a while now, I think they're the slowest part of the NS code (not sure the impact on total performance, but the fits have tons of allocations and use mutable structs and are pretty haphazard). I don't really like the interface, either, and I think it can all be rewritten to utilize Distributions.jl since these "bounding spaces" are just statistical distributions.

So, with all that being said, I just want you to know that I've seen this, and that I have ideas for making this part of the interface better! I really appreciate what you've contributed already 😃

aplavin commented 2 years ago

Yeah, no real hurry. I've used my modifications (this PR) on some problems already, and sampling basically never fails - as in "throws", while it did before. Note that the PR is independent from the exact interface, it's about computations themselves.