use relative eigen threshold

[aplavin]

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:

NestedSamplers.jl/src/bounds/ellipsoid.jl

Line 108 in eb6999e

A = pinv(cov)

This matrix is passed to the Ellipse constructor, and fails in sqrt at

NestedSamplers.jl/src/bounds/ellipsoid.jl

Line 49 in eb6999e

axlens = @. 1 / sqrt(E.values)

Actually, eigenvalues of A are not negative as-is, but become negative in Symmetric(A). Digging into this error, I found that the

NestedSamplers.jl/src/bounds/ellipsoid.jl

Lines 168 to 178 in eb6999e

	function make_eigvals_positive!(cov::AbstractMatrix, targetprod)
	E = eigen(cov)
	mask = E.values .< 1e-10
	if any(mask)
	nzprod = prod(E.values[.!mask])
	nzeros = count(mask)
	E.values[mask] .= (targetprod / nzprod)^(1 / nzeros)
	cov .= E.vectors * Diagonal(E.values) / E.vectors
	end
	return cov
	end

function applies a huge correction to the cov matrix in my case at

NestedSamplers.jl/src/bounds/ellipsoid.jl

Line 105 in eb6999e

make_eigvals_positive!(cov, targetprod)

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]

Codecov Report

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

@@            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%)	⬆️

---

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]

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

[mileslucas]

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]

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.