3
$\begingroup$

How can one evaluate $-\text{log}~ \text{Det}(I_{i-j}(x))$ where $x>200$ and $I_{\nu}(x)$ is the modified Bessel function of first kind and $i,j = 1 \cdots 100$ and the determinant is take of the matix whose elements are given by $A_{ij} = I_{i-j}(x)$. I used Scipy/Python and it gives me an incorrect answer. I want to see how much better Mathematica does. Note that $I_{i-j}(200)$ is a huge number for any $i-j$.

$\endgroup$
  • $\begingroup$ What have you tried? $\endgroup$ – C. E. Nov 7 '19 at 6:58
4
$\begingroup$

This is what I get with Mathematica:

logdet[x_, prec_] := - Log @ Det @ N[Array[BesselI[#1-#2, x]&, {100, 100}], prec]

Presumably due to subtractive cancellation, a high precision is needed to get an accurate result, so:

r500 = logdet[200, 500]

-9033.801557200126884639983232381821638986208454697027595029420498495859730519\ 029590617989373466937836223259754375229723311300136058858941346024154635908053\ 762729182813641676207541530844596440371730202099559364747933131151396757544867\ 410989540289820352083687426358192116464335718811026707179015898854249615451052\ 308452992411926795485314268493234726335846547885562794131295349448973511839300\ 690475561361508703872728426715140348231043410320195012668497801821524162092336\ 9400443012392291448771312756532037106

Comparing this to a computation with a bit more precision:

r510 = logdet[200, 510]

-9033.801557200126884639983232381821638986208454697027595029420498495859730519\ 029590617989373466937836223259754375229723311300136058858941346024154635908053\ 762729182813641676207541530844596440371730202099559364747933131151396757544867\ 410989540289820352083687426358192116464335718811026707179015898854249615451052\ 308452992411926795485314268493234726335846547885562794131295349448973511839300\ 690475561361508703881209596637699666119599037581235795792730971525348781541540\ 179968386785031399656551557996358738271830794709

The difference is:

Log10[r500 - r510]

-404.0715442354949104419778748391483390025461143878688319727113963282042725909\ 57302327785726278047318

So, I would expect r500 to have an accuracy of about 404:

SetAccuracy[r500, 404]

-9033.801557200126884639983232381821638986208454697027595029420498495859730519\ 029590617989373466937836223259754375229723311300136058858941346024154635908053\ 762729182813641676207541530844596440371730202099559364747933131151396757544867\ 410989540289820352083687426358192116464335718811026707179015898854249615451052\ 308452992411926795485314268493234726335846547885562794131295349448973511839300\ 6904755613615087039

| improve this answer | |
$\endgroup$
  • $\begingroup$ Very informative answer. And the precision becomes even more important for x >> 200. Thanks. I wonder why Python can't handle this but that's not something for here. $\endgroup$ – R.G.J Nov 6 '19 at 22:57

Not the answer you're looking for? Browse other questions tagged or ask your own question.