pytensor icon indicating copy to clipboard operation
pytensor copied to clipboard

Remove `Iv` core `Op` in favor of `Ive` and add rewrite rule for log of it

Open ricardoV94 opened this issue 1 year ago • 3 comments

Description

Similarly to how we approached Kv/Kve in #1081, we only need one core Op: Ive.

We should include the equivalent log(iv(x)) rewrite for stability, where iv is the expression based on ive returned by the helper.

The goal is to have one canonical form and reduce the number of Ops we need to test/support.

ricardoV94 avatar Nov 15 '24 21:11 ricardoV94

Hey @ricardoV94, for this issue I followed the PR for kv and I expressed iv in the form of a ive helper as iv(v,z) = ive(v,z) * exp(z), but it overflows for iv(v,z) when z>100 for v=4.5 that is around 1e36 which I don't think should be happening. Please let me know where I'm going wrong. Thanks!

AdvH039 avatar Dec 07 '24 13:12 AdvH039

Sorry for missing this thread. How long after does iv itself start overflowing? It's not completely surprising that the multiplication will be a bit less precise.

It might or not be a deal breaker.

The rewrite for the log will probably still make things more stable via the ive route

ricardoV94 avatar Jan 23 '25 21:01 ricardoV94

I'm not seeing a very early overflow?

import numpy as np
from scipy.special import iv, ive
print(
    iv(4.5, 500),
    ive(4.5, 500) * np.exp(500)
) # 2.4545477241973476e+215 2.4545477241973484e+215

It overflows a tiny bit before iv:

import numpy as np
from scipy.special import iv, ive

v = 4.5
for z in [709, 710, 711, 712, 713, 714, 715]:
    print(z, iv(v, z), ive(v, z) * np.exp(z))

709 1.2140731682798808e+306 1.2140731682798845e+306
710 3.2979337389547726e+306 inf
711 8.958584547932482e+306 inf
712 2.433533186526627e+307 inf
713 6.610518938992669e+307 inf
714 1.7957018810395855e+308 inf
715 inf inf

For now I think it's fine trade-off for simpler codebase. We can revisit later. Most applications will probably want to work on log-scale anyway, where the ive approach should win clearly

import numpy as np
from scipy.special import iv, ive

v = 4.5
for z in [709, 710, 711, 712, 713, 714, 715]:
    print(z, np.log(iv(v, z)), np.log(ive(v, z)) + z)
# 709 704.7850194174108 704.7850194174108
# 710 705.7843345888324 705.7843345888324
# 711 706.7836506961244 706.7836506961244
# 712 707.7829677367349 707.7829677367349
# 713 708.7822857081223 708.7822857081223
# 714 709.7816046077556 709.7816046077554
# 715 inf 710.7809244331132

ricardoV94 avatar Jan 31 '25 10:01 ricardoV94