Incorrect Results for Extra Floating Functions

[julmb]

The library calculates incorrect results for at least some of the extra functions log1p, expm1, log1pexp, and log1mexp of the Floating class.

For example:

>>> import Numeric.AD
>>> log1pexp <$> [-1000, -100, -10, -1, 0, 1, 10, 100, 1000]
>>> fst . diff' log1pexp <$> [-1000, -100, -10, -1, 0, 1, 10, 100, 1000]
[0.0,3.720075976020836e-44,4.539889921686465e-5,0.31326168751822286,0.6931471805599453,1.3132616875182228,10.000045398899218,100.0,1000.0]
[0.0,0.0,4.5398899216870535e-5,0.31326168751822286,0.6931471805599453,1.3132616875182228,10.000045398899218,100.0,Infinity]

Near zero, there are only mild inaccuracies, but further away, the results are completely wrong.

I believe that this is due to missing definitions for these functions in the type class instances. This causes them to use the default implementations for those functions, which are pretty bad. For example, log1pexp x = log1p (exp x), which exhausts the precision of Double quite quickly.

I will do some further testing and then most likely start working on a pull request.

[RyanGlScott]

I believe this was addressed by #111, so I'll close this. Please re-open if there is more to do here.

[julmb]

All good, thank you very much!

[julmb]

Welp, looks like I will be working on this some more:

>>> import Numeric.AD
>>> diff (diff log1pexp) (-1000)
NaN

This is after #111 has been applied. Is there a way to reopen this issue? As far as I am concerned, the topic still applies. I will make a fresh PR though.

[julmb]

The problem here seems to be with derivatives of reciprocals of large numbers. This problem is not unique to the new implementations of the extra floating-point functions in #111. For example:

>>> import Numeric.AD
>>> diff' (recip . exp) 1000
>>> diff' (\ x -> recip $ x * x) 1e400
(0.0,NaN)
(0.0,NaN)

In both cases, the function evaluates to recip Infinity, which is just 0 and that would be fine. The derivatives also exist and are well-behaved. However, they are calculated as follows:
$$f(x) = \frac{1}{g(x)} \quad f'(x) = -\frac{g'(x)}{g(x)^2}$$
So we end up with - exp x / exp x ^ 2 and - 2 * x / x ^ 4 respectively. This evaluates to Infinity / Infinity = NaN.

This seems to be a fundamental problem for the ad library. My guess is that it does not usually come up as one would have to be operating with extreme values. However, with log1pexp being just the Softplus function, it is quite reasonable to evaluate it at values like -1000 or 1000, so it is easy to run into this issue.

[cartazio]

How many digits of excess precision would we need to calculate it correctly?

…

On Thu, Mar 14, 2024 at 8:39 AM Julian Brunner ***@***.***> wrote: The problem here seems to be with derivatives of reciprocals of large numbers. This problem is not unique to the new implementations of the extra floating-point functions in #111 <#111>. For example: >>> import Numeric.AD>>> diff' (recip . exp) 1000>>> diff' (\ x -> recip $ x * x) 1e400 (0.0,NaN) (0.0,NaN) In both cases, the function evaluates to recip Infinity, which is just 0 and that would be fine. The derivatives also exist and are well-behaved. However, they are calculated as follows: $$f(x) = \frac{1}{g(x)} \quad f'(x) = -\frac{g'(x)}{g(x)^2}$$ So we end up with - exp x / exp x ^ 2 and - 2 * x / x ^ 4 respectively. This evaluates to Infinity / Infinity = NaN. This seems to be a fundamental problem for the ad library. My guess is that it does not usually come up as one would have to be operating with extreme values. However, with log1pexp being just the Softplus function, it is quite reasonable to evaluate it at values like -1000 or 1000, so it is easy to run into this issue. — Reply to this email directly, view it on GitHub <#108 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAABBQUFVQAUO4GV4OUPDP3YYGLAPAVCNFSM6AAAAABEJCP4WGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSOJXGM3DIMJXHA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[julmb]

How many digits of excess precision would we need to calculate it correctly?

I am not sure what you mean by that, but the derivative of log1pexp includes exp in the denominator. With the largest Double being on the order of $10^{308}$, for a calculation to not overflow to Infinity, we would have to represent numbers as large as $e^{10^{308}}$.

This is neither possible nor necessary, as log1pexp is nearly constant and linear at the far ends, respectively. The instances for Float and Double use this property to define log1pexp piecewise. However, we do not have that option when defining an instance like Floating a => Floating (Forward a) as there is no Ord a constraint in scope. I have not yet found a solution to this problem.

[julmb]

I might be overlooking something, but none of this would even be an issue in the first place if log1pexp and log1mexp were just functions of type (Ord a, Floating a) => a -> a rather than members of the Floating type class.

1. Their default implementations only yield correct (or even any!) results for a tiny fraction of the domains of these functions
2. The only two instances that actually implement them are Float and Double and they both use the same (Ord a, Floating a) => a -> a implementation that is partly duplicated
3. We would get correct derivatives for free from these implementations, without maintenance burden or dealing with issues like https://gitlab.haskell.org/ghc/ghc/-/issues/17125

I understand why log1p and expm1 are class members, they are translated to primops during compilation. But I feel like log1pexp and log1mexp should just be ordinary functions.

[julmb]

Yesterday I was wondering why log1pexp and log1mexp were members of the Floating class. So today I took a deep dive and read through https://gitlab.haskell.org/haskell/prime/-/wikis/libraries/proposals/expand-floating and https://mail.haskell.org/pipermail/libraries/2014-April/022667.html to understand the ideas and motivations behind the design.

It turns out that one of the motivations was to support instances of vector spaces (https://mail.haskell.org/pipermail/libraries/2014-April/022741.html). In the mean time, I did overlook the fact that types like Complex a could not use the Ord-based implementations. So there are good reasons for log1pexp and log1mexp to be members of the Floating class. However, base does not seem to make good use of this capability, as instances like Complex a still use the default implementations:

ghci> log1pexp (1000 :+ 0)
NaN :+ NaN

It also looks like @ekmett, who orignally made the proposal, already had ad in mind when coming up with it. However, it seems like he never got around to implementing the new functions in the AD mode instances.

[julmb]

As for the actual issue at hand here, I have been considering several ideas:

1. Implement log1*exp as lift1 log1*exp d in such a way that d and all the higher-order derivatives of it are well-behaved throughout the entire domain. I do not know how to do this.
2. Implement log1*exp piece-wise to exploit behavior near the extremes of the domain. I do not know how to do this without having access to an Ord instance.
3. Add an Ord a constraint to the Floating instances of AD modes. This would mean we can no longer have Floating (Forward (Complex Double)).
4. Add a warning about higher-order derivatives of log1*expOrd to the documentation and resign to using out-of-class implementations like log1*expOrd :: (Ord a, Floating a) => a -> a for applications that require AD. This would make functions that use AD needlessly less general since there are good implementations for types like Complex a despite the lack of an Ord instance. The whole point of having log1*exp in the Floating class is to enable abstraction over these different implementations.

It sounds like this should not be such a hard problem, but I am pretty lost and do not really see a good solution here. I feel like the issue boils down to this:

• Concrete instances like Floating Double are easy since they have access to all the tools of the concrete types
• Newtype-like instances like Floating a => Floating (Down a) are easy since they just pass through the values
• The instance RealFloat a => Floating (Complex a) requires the base type to be RealFloat in order to get access to the necessary tools
• The instance Floating a => Floating (Forward a) has to squeeze blood from a stone. It has to implement nontrivial functions using nothing but a Floating instance on the base type. And we cannot really strengthen the constraints on the base type since users really do expect to be able to use instances like Floating (Forward (Complex Double)). In a sense, we have to come up with a single implementation for the derivative of log1*expOrd that works accurately for both Double and Complex Double simultaneously.

Sorry if I am rambling a bit. Just trying to put my thoughts down before the weekend since I will not have time to work on this until next week. I am curious if anyone has either more ideas or thoughts on the ideas listed here already.

[ekmett]

I think one part of the answer might be to provide a custom Complex type that has the corrected instance if base can't be brought to fix it and then use it to shame base into correct behavior.

I also vaguely recall that log1mexp might have had wrongly defined behavior in base, which is yet another issue with getting these fixed up to a fully usable state.