# Is there a C implementation of \`erfcinv\` (or \`erfinv\`) that we could add to Pytensor to allow loop fusion? Can we write one?

**URL:** <https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646>\
**Category:** v5\
**Tags:** pytensor\
**Created:** [August 3, 2023, 7:22am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646 "2023-08-03T07:22:48Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![jonsedar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jonsedar/32/2590_2.png) [@jonsedar](https://discourse.pymc.io/u/jonsedar)\
**Post date:** [August 3, 2023, 7:22am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/1 "2023-08-03T07:22:48Z")

</div>

I’m using an [`erfcinv`](https://github.com/pymc-devs/pytensor/blob/4cc13bcea9d49f165d68774f6aab41fb3e46b4e0/pytensor/scalar/math.py#L229) in a model and get a warning from [`initialize_fuseable_mappings`](https://github.com/pymc-devs/pytensor/blob/4cc13bcea9d49f165d68774f6aab41fb3e46b4e0/pytensor/tensor/rewriting/elemwise.py#L691) to the effect that we can’t do loop fusion because there’s no C implementation available / used.

I’m out of my depth here, but a quick search suggests this function simply doesn’t exist in the standard C library and developers have to create their own. [algorithm - Computing the inverse of the complementary error function, erfcinv(), in C - Stack Overflow](https://stackoverflow.com/questions/60472139/computing-the-inverse-of-the-complementary-error-function-erfcinv-in-c). Relatedly, I see there’s a TODO in the docstring of the related [`Erfinv`](https://github.com/pymc-devs/pytensor/blob/4cc13bcea9d49f165d68774f6aab41fb3e46b4e0/pytensor/scalar/math.py#L183) to go find such a C implementation.

**My question: Is it worth us (me?) trying to implement such a function in Pytensor directly? Or is this liable to just lead to implementation pain?**

I think [I’ve found an C implementation of `erfinv` here](https://github.com/lakshayg/erfinv) based on the implementation in [go](https://github.com/golang/go/blob/master/src/math/erfinv.go) - which could get us halfway there…

Additional reading:

- [https://people.maths.ox.ac.uk/gilesm/files/gems\_erfinv.pdf](https://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf)
- [Inverse complementary error function - MATLAB erfcinv - MathWorks United Kingdom](https://uk.mathworks.com/help/matlab/ref/erfcinv.html)
- [Inverse Error Function in C - Stack Overflow](https://stackoverflow.com/questions/27229371/inverse-error-function-in-c)
- [Error function - Wikipedia](https://en.wikipedia.org/wiki/Error_function)

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [August 3, 2023, 7:45am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/2 "2023-08-03T07:45:53Z")

</div>

It’s a bit rough but would be a great addition. Perhaps the c code for the gamma functions is a good template (and how it’s called from the related gamma Ops)? [https://github.com/pymc-devs/pytensor/blob/main/pytensor/scalar/c\_code/gamma.c](https://github.com/pymc-devs/pytensor/blob/main/pytensor/scalar/c_code/gamma.c)

---

<div class="post-metadata">

**Author:** ![jonsedar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jonsedar/32/2590_2.png) [@jonsedar](https://discourse.pymc.io/u/jonsedar)\
**Post date:** [August 3, 2023, 7:59am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/3 "2023-08-03T07:59:00Z")

</div>

Okay cool - I’ll look into it a little more then 🙂

Do you have any feel for the kind of speedup one might get from this? I know that’s fairly impossible to guess, but are we talking fractional percentage gains in speed, or multiples?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [August 3, 2023, 8:19am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/4 "2023-08-03T08:19:02Z")

</div>

It can provide nice speedups but it’s very context dependent.

You can benchmark your d+logp function in JAX to get an idea of a lower bound (they also do loop fusion and are generally faster than PyTensor C)

You could introduce a rewrite that replaces erfcinv by the less stable `erfinv(1-x)` if that has C code to get a more direct answer.

Having said that I wouldn’t suggest doing a C implementation as an attempt to speedup your model as a user. Too much work for that.

---

<div class="post-metadata">

**Author:** ![jonsedar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jonsedar/32/2590_2.png) [@jonsedar](https://discourse.pymc.io/u/jonsedar)\
**Post date:** [August 3, 2023, 8:25am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/5 "2023-08-03T08:25:29Z")

</div>

Hmmm… food for thought, thanks!

I’ve not used JAX before and TBH am a little apprehensive of probably days of work to find & fix numerical issues etc (unfounded?)

But I do see JAX has an [`erfinv`](https://jax.readthedocs.io/en/latest/_autosummary/jax.scipy.special.erfinv.html) - and (without knowing anything about how it works), might it be possible in theory to JAXiify `scipy.special.erfcinv` too?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [August 3, 2023, 8:40am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/6 "2023-08-03T08:40:15Z")

</div>

Those functions are automatically used by Pytensor jaxify machinery based on the names, there should be no need for you to do anything.

I mentioned it more as a benchmark to assess the max speedup you could hope for by doing a C impl.

Using numpyro/nutpie for sampling is likely to shift sampling speed much more than adding C code for this one Op, but that requires more care so I didn’t suggest that directly.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [August 3, 2023, 9:32am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/7 "2023-08-03T09:32:49Z")

</div>

Scipy already has a C implementation of erfcinv [here](https://github.com/scipy/scipy/blob/a4ed9450b6d5c1f36495b29a0471f03b471ca8ae/scipy/special/cephes/erfinv.c#L57) I think? I’m not sure what would be required to bring it into pytensor.

It seems like it should be possible to write a numba overload for it though, and use `nutpie` for a speedup. I’m surprised it’s not already overloaded in `numba-scipy` like most of the `scipy.special` module.

---

<div class="post-metadata">

**Author:** ![jonsedar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jonsedar/32/2590_2.png) [@jonsedar](https://discourse.pymc.io/u/jonsedar)\
**Post date:** [August 3, 2023, 11:10am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/8 "2023-08-03T11:10:02Z")

</div>

Indeed it looks like [pytensor’s erfcinv](https://github.com/pymc-devs/pytensor/blob/4cc13bcea9d49f165d68774f6aab41fb3e46b4e0/pytensor/scalar/math.py#L229) already falls back onto the scipy implementation…

I’m all ears about a nutpie / JAX implementation that overloads the scipy one… but that’s a little above my pay grade!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [August 3, 2023, 11:25am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/9 "2023-08-03T11:25:51Z")

</div>

> Indeed it looks like [pytensor’s erfcinv](https://github.com/pymc-devs/pytensor/blob/4cc13bcea9d49f165d68774f6aab41fb3e46b4e0/pytensor/scalar/math.py#L229) already falls back onto the scipy implementation…

Yes but that will run in Python mode and indeed break the loop fusion as the warning tells you.

_Edited_

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [August 3, 2023, 11:27am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/10 "2023-08-03T11:27:42Z")

</div>

> Scipy already has a C implementation of erfcinv [here](https://github.com/scipy/scipy/blob/a4ed9450b6d5c1f36495b29a0471f03b471ca8ae/scipy/special/cephes/erfinv.c#L57) I think? I’m not sure what would be required to bring it into pytensor.

Just copy the files and make sure the Op loads the headers / references the function in the `c_code` method. It’s not too bad, but requires some trial and error.

Importantly that would be good for users in general, but I doubt it will fix any bottleneck in a specific model (probabilistically)

---

<div class="post-metadata">

**Author:** ![jonsedar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jonsedar/32/2590_2.png) [@jonsedar](https://discourse.pymc.io/u/jonsedar)\
**Post date:** [August 3, 2023, 11:45am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/11 "2023-08-03T11:45:23Z")

</div>

> [@ricardoV94](#):
>
> Just copy the files and make sure the Op loads the headers / references the function in the `c_code` method. It’s not too bad, but requires some trial and error.

Okay, that sounds like something a rank C/C++ amateur (like yours truly) could at least attempt 😉

I’ll look into it in a week or so

---

<div class="post-metadata">

**Author:** ![julien](https://avatars.discourse-cdn.com/v4/letter/j/71e660/32.png) [@julien](https://discourse.pymc.io/u/julien)\
**Post date:** [August 9, 2023, 3:49pm UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/12 "2023-08-09T15:49:08Z")

</div>

In terms of standard C libraries, “erfcinv” is actually in the GSL, it’s just slightly hidden: one has to use the inverse CDF gsl\_cdf\_gaussian\_Pinv.

> **[The Gaussian Distribution - GNU Scientific Library -- Reference Manual](http://gnu.ist.utl.pt/software/gsl/manual/html_node/The-Gaussian-Distribution.html)**
>
> GNU Scientific Library -- Reference Manual

Otherwise, less standard but still common, you can find it directly in Boost for a C++ code:  
[https://live.boost.org/doc/libs/1\_67\_0/libs/math/doc/html/math\_toolkit/sf\_erf/error\_inv.html](https://live.boost.org/doc/libs/1_67_0/libs/math/doc/html/math_toolkit/sf_erf/error_inv.html)

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [August 10, 2023, 7:48am UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/13 "2023-08-10T07:48:48Z")

</div>

Thanks for sharing. I don’t know how easy it is to include GSL stuff for our Windows users, so my inclination would be to go with the scipy/cephes implementation.

---

<div class="post-metadata">

**Author:** ![twiecki](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/twiecki/32/6930_2.png) [@twiecki](https://discourse.pymc.io/u/twiecki)\
**Post date:** [August 11, 2023, 7:04pm UTC](https://discourse.pymc.io/t/is-there-a-c-implementation-of-erfcinv-or-erfinv-that-we-could-add-to-pytensor-to-allow-loop-fusion-can-we-write-one/12646/14 "2023-08-11T19:04:54Z")

</div>

Moreover, GSL is GPL licensed which is incompatible with BSD/Apache v2.