BUG: Increase accuracy of log1p for small inputs #22611

mfkasim1 · 2022-11-17T16:09:26Z

Fixes #22609 to increase the accuracy of log1p for small inputs. With this PR, log1p can deserve its reason for existence for complex numbers.

import numpy as np
print(np.log1p(1e-18 + 1e-18j))  # 1e-18 + 1e-18j, as expected

seberg · 2022-11-17T18:03:57Z

Thanks for looking into a fix so quickly! Before I (or someone else dives in), maybe you have some quick tips to help. Do you have a reference to quickly make sure this is a good implementation? Does it have a speed tradeoff?

rossbar

Also FWIW the test_log1p_complex tests all pass on main, so they're not actually testing the desired behavior. I think this is expected given the precision involved.

WarrenWeckesser · 2022-11-17T18:38:31Z

Thanks @mfkasim1. This fixes the one case reported in gh-22609, but it breaks several other cases that are correctly handled by the existing code. The new problems are created by overflow or underflow in the expression xx * (2 + xx) + yy * yy. For example, try these inputs:

         z                expected value of log1p(z)
    -1 + 1e-18j        -41.44653167389282 + 1.5707963267948966j        
    -1 + 1e250j        575.6462732485114 + 1.5707963267948966j
    1e250 + 1j         575.6462732485114 + 1e-250j
    1e250 + 1e250j     575.9928468387914 + 0.7853981633974483j  
    1e-250 + 1e250j    575.6462732485114 + 1.5707963267948966j
    1e250 + 2e-250j    575.6462732485114 + 0.0j

A correct implementation will require scaling the inputs appropriately, similar to what is done for npy_hypot:

numpy/numpy/core/src/npymath/npy_math_internal.h.src

Lines 223 to 253 in 9e144f7

    
           NPY_INPLACE double npy_hypot(double x, double y) 
        
           { 
        
           #ifndef NPY_BLOCK_HYPOT 
        
               return hypot(x, y); 
        
           #else 
        
               double yx; 
        
               if (npy_isinf(x) || npy_isinf(y)) { 
        
                   return NPY_INFINITY; 
        
               } 
        
               if (npy_isnan(x) || npy_isnan(y)) { 
        
                   return NPY_NAN; 
        
               } 
        
               x = npy_fabs(x); 
        
               y = npy_fabs(y); 
        
               if (x < y) { 
        
                   double temp = x; 
        
                   x = y; 
        
                   y = temp; 
        
               } 
        
               if (x == 0.) { 
        
                   return 0.; 
        
               } 
        
               else { 
        
                   yx = y/x; 
        
                   return x*npy_sqrt(1.+yx*yx); 
        
               } 
        
           #endif 
        
           }

mfkasim1 · 2022-11-18T11:23:30Z

@rossbar Thanks for the comment. I've added the decimal arguments in the test to make sure it's checked with a lot of decimal points.

seberg · 2022-11-18T11:38:26Z

numpy/core/tests/test_umath.py

+    def test_log1p_complex(self):
+        assert_almost_equal(ncu.log1p(0.2 + 0.3j), ncu.log(1.2 + 0.3j))
+        assert_almost_equal(ncu.log1p(1e-19 + 1e-18j), 1e-19 + 1e-18j,
+                            decimal=23)


I suspect that assert_almost_equal_nulp is the more convenient way of checking things here.

@WarrenWeckesser if you could want to dig into whether the implementation looks right, that would be great. I would still feel things would be much simpler if we can "steal" an implementation from elsewhere (that has a clearly BSD compatible license).

Just to centralise ideas (pytorch/pytorch#89214 (comment)), I agree. Now, it's not clear where is this function implemented for complex numbers really...

mfkasim1 · 2022-11-19T15:58:24Z

I found a good implementation from scipy: https://github.com/scipy/scipy/blob/5b8490a2b80e73aec59af9ac93bca2ebd178fee1/scipy/special/_cunity.pxd#L25-L66
The idea is similar to handle it specially if the input is small. However, they also added a conditional if xx ** 2 + yy ** 2 + 2 * xx is close to zero, they use double-double precision.

WarrenWeckesser · 2022-11-19T19:12:06Z

@mfkasim1, I'm testing a reformulation of the calculation that seems to work well and shouldn't require higher precision for special cases.

Update: The reformulation below doesn't work as well as I originally thought. If I check the relative errors of the real and imaginary parts separately, I see that there are cases where the real parts can have errors that are large. For some inputs, the expression np.log1p(x) + 0.5*np.log1p((y/(x + 1))**2) in the code below will suffer loss of precision. A similar phenomenon is noted in the comment of the SciPy code that @mfkasim1 linked to.

The idea is to use this reformulation of $(x+1)^2 + y^2$:

$$ (x+1)^2 + y^2 = (x+1)^2 \left(1 + \left(\frac{y}{x + 1}\right)^2\right) $$

so

$$ \begin{split} \log\left((x+1)^2 + y^2\right)^{1/2} & = \frac{1}{2} \log\left((x+1)^2 + y^2\right) \\ & = \frac{1}{2} \log\left((x+1)^2 \left(1 + \left(\frac{y}{x+1}\right)^2\right)\right) \\ & = \log(x + 1) + \frac{1}{2}\log\left(1 + \left(\frac{y}{x+1}\right)^2\right) \\ & = \textrm{log1p}(x) + \frac{1}{2}\textrm{log1p}\left(\left(\frac{y}{x + 1}\right)^2\right) \end{split} $$

In the code, this formulation would be used when $x$ is small (e.g. $|x| < 1/2$) and when $|y| < x + 1$ (so $(y/(x+1))^2$ does not cause problems with overflow).

Here's an implementation in Python that I've been using to test this:

def log1p(z):
    x, y = z.real, z.imag
    theta = np.arctan2(y, x + 1)
    if abs(x) < 0.5 and abs(y) < x + 1:
        lnr = np.log1p(x) + 0.5*np.log1p((y/(x + 1))**2)
    else:
        lnr = np.log(np.hypot(x + 1, y))
    return lnr + 1j*theta

So far, quite a bit of fuzz testing shows that it works well. I use mpmath with high precision to compute the "true" values for tests. E.g.:

In [469]: from mpmath import mp

In [470]: mp.dps = 125

In [471]: z = -2e-17 + 2e-8j

Compute the expected value of log1p(z) with mp.log1p. Here is the 125 digit result

In [472]: mp.log1p(z)
Out[472]: mpc(real='0.00000000000000017999999999999997473817585095895590238446612672721306825082019861479410189297518927535665999267474842690586183589671360222800567',
imag='0.000000019999999999999998151784549935903144445373327991386300560662953457748900829013377241693442896709281068311236443103345569367960809')

And this is the expected double precision result:

In [473]: complex(mp.log1p(z))
Out[473]: (1.7999999999999997e-16+1.9999999999999997e-08j)

For this z, the function log1p(z) defined above gets the expected result:

In [474]: log1p(z)
Out[474]: (1.7999999999999997e-16+1.9999999999999997e-08j)

Fuzz testing shows that the relative error is generally less than 5e-16.

I'm still checking this, so don't update your code--it could definitely use another set of eyes. If anyone finds a case where this doesn't work, add a comment here.

mfkasim1 · 2022-11-21T15:54:25Z

@WarrenWeckesser Thanks for testing the reformulation! Could you give us an example where the formulation above produces large error?

seberg · 2022-11-22T17:58:11Z

CircleCI is doing some weird thing here, I hope with the next ~~stuff~~ push it will clear out, but this PR shouldn't affect documentation, so just ignore it if not.

mfkasim1 · 2022-11-23T11:06:47Z

It seems now the best option we have is to follow scipy's implementation. However, I can't see numpy have dependency on cephes's double double (required to handle the special case of cancellation terms) or any other cephes files. Would having the dependency on cephes of interest to numpy? Otherwise, should we just sacrifice the accuracy in that very special case (i.e. similar to calculate 1 + 1e-30 - 1 == 0.

WarrenWeckesser · 2022-11-25T14:01:25Z

@mfkasim1, for the function log1p that I showed above, any point in the complex plane on or sufficiently close to the circle $|z + 1| = 1$ (i.e. $(x + 1)^2 + y^2 = 1$, or $x^2 + 2x + y^2 = 0$, which is the same condition as in the SciPy code) will suffer loss of precision in the real part of the result. For example, z = -0.1 + 0.43588989435j is close to that circle:

In [311]: z = -0.1 + 0.43588989435j

In [312]: abs(z + 1)
Out[312]: 0.9999999999982271

The function log1p that I defined earlier gives

In [313]: log1p(z)
Out[313]: (-1.772956781387336e-12+0.4510268117926018j)

The true value (as computed with mpmath) is

In [314]: from mpmath import mp

In [315]: mp.dps = 100

In [316]: complex(mp.log1p(z))
Out[316]: (-1.7729356247190155e-12+0.4510268117926018j)

That shows that the relative error in the real part computed by log1p(z) is on the order of 1e-5.

However, if we look at the closeness of the complex numbers (rather than the relative error of just the real part), we get a relative error of 4.7e-17:

In [319]: w = log1p(z)

In [320]: w_true = complex(mp.log1p(z))

In [321]: abs(w - w_true)/abs(w_true)
Out[321]: 4.6907784121193136e-17

It seems now the best option we have is to follow scipy's implementation. However, I can't see numpy have dependency on cephes's double double (required to handle the special case of cancellation terms)

Even if there was interest in having a double-double library to use as a last resort for maintaining precision in special cases, I don't think the SciPy double-double code (in its current state) is a good source. Some of the transcendental functions in that code are buggy (see scipy/scipy#16227).

seberg · 2022-12-01T09:45:59Z

As much as I like getting things fully right. @WarrenWeckesser do you have a recommendation for the best improvement that we can easily get right now?
It sounds like we have a options that we can be confident about being better at little or no cost otherwise, and considering that cephes falls back to higher precision operations, it sounds like a fully precise option likely doesn't even exist.

seberg · 2023-02-17T12:02:32Z

Ping once more since I remembered this as better than nothing. I think we should do something (maybe just copy whatever pytorch did).

I am happy to base this on gut feeling and confidence that the new method is better rather than ideal in this case. But I do trust your gut feeling much more than mine @WarrenWeckesser :).

WarrenWeckesser · 2023-07-27T14:57:30Z

@mfkasim1 and @seberg, sorry for not getting back to this for so long. I'm running more fuzz tests this week to compare this PR, my suggested formulation, and the version that is used in the Julia implementation at https://github.com/JuliaLang/julia/blob/c43e5a10be27b7f93b5368875aa1d2596b4d4947/base/complex.jl#L745. I'll report back in a day or so with the results.

WarrenWeckesser · 2023-08-03T02:08:49Z

FYI: I ended up exploring a few rabbit holes:

fixed a Julia bug: Fix log1p(z) for extremely small |z|. JuliaLang/julia#50769;
added a two implementations to ufunclab: log1p_theorem4 and log1p_doubledouble;
found the Jax/XLA version (source code; it has a bug, but I haven't created an issue for it yet);
ran many tests on the variety of ways that log1p(z) can be computed.

I'm still working on summarizing and coming up with a recommendation for NumPy.

WarrenWeckesser · 2023-08-07T14:16:19Z

After lots of experimentation with several versions of the complex log1p, here's the version that I recommend, written as a Python function that accepts scalars:

def log1p_theorem4(z):
    z = complex(z)
    u = z + 1
    if u == 1:
        return z
    elif u - 1 == z:
        return np.log(u)
    else:
        return np.log(u)*(z/(u - 1))

(The final version should probably have a check for either z.real or z.imag being nan, and return complex(nan, nan) in that case.)

It is based on the formula for log1p(x) (real argument) given in Theorem 4 of Goldberg's paper "What every computer scientist should know about floating-point arithmetic". The only tweak that I've made is the extra check for u - 1 == z, where log(u) is returned. Computing the complex division z/(u - 1) in that case can only introduce an additional (small) error.

I fuzz-tested with large samples clustered around sets such as: the circle |z + 1| = 1; specific points on that circle including 0+0j, -2+0j, -1+1j, -1-1j; the point -1+0j; and a few other scatter points. I also checked on various rays emanating from -1+0j through and beyond the circle |z + 1| = 1. I computed the relative error of the complex result w = log1p_theorem4(z) as abs(w - w_true)/abs(w_true), where w_true is computed to full precision by using mpath, and I also looked at the relative error of the real part, abs(w.real - w_true.real)/abs(w_true.real). In all cases, the relative error of the complex result was less than 6e-16.

The worst case relative error of the real part of the result depends the relative error of the real part of the underlying clog(z) function provided by the platform/compiler. On my Linux machine, this relative error was also very small, generally not exceeding 8e-16. That is because the complex log function clog(z) provided by the standard library is apparently very accurate, even on the troublesome unit circle.

On my Mac (m2), using clang, the relative error of the real part of the result could be quite high. That's because the relative error of the real part of clang's clog(z) on the unit circle can be high. For example,

In [36]: from mpmath import mp

In [37]: mp.dps = 200

In [38]: za = -0.0026626976230049726-0.07292671175487339j

This point is on the circle |z + 1| = 1:

In [39]: abs(za + 1)
Out[39]: 1.0

In [40]: w_true = complex(mp.log1p(za))

In [41]: w = log1p_theorem4(za)

In [42]: print(w_true)
(1.5951249881265712e-22-0.07299150803396787j)

In [43]: print(w)
(4.336808689942018e-19-0.07299150803396787j)

The only thing correct about the real part is the sign. The values are so small, however, that the overall relative error of the complex result is also very small:

In [44]: abs(w - w_true)/abs(w_true)
Out[44]: np.float64(5.93933963240823e-18)

The only way I could get the relative error of the real part reliably small on the Mac was to use an approach similar to that of SciPy: when the input z is close to the circle |z + 1| = 1, compute the real part of the result using double-double numbers. (This version is now implemented in ufunclab as log1p_doubledouble.) Using double-double can be an effective way to maintain precision when you really need it, but it is not clear to me that it is important for us to ensure that the real part has very low relative error; just ensuring that the overall relative error of the complex result is small seems good enough. Does that sound reasonable? I don't think we have an explicit policy about the acceptable relative errors for our functions.

seberg · 2023-08-07T14:26:26Z

I would say go for it. Anything seems like an improvement here and you covered a lot of bases it seems (enough so to fix a few bugs elsewhere!).
One thing I am just now wondering is whether the system cmath library may have an implementation available that we could use (at least sometimes)?

WarrenWeckesser · 2023-08-07T14:44:45Z

One thing I am just now wondering is whether the system cmath library may have an implementation available that we could use (at least sometimes)?

I've only seen the function for real arguments (e.g. C99 has log1pf/log1p/log1pl).

WarrenWeckesser · 2023-08-10T17:07:40Z

Well, here's another interesting wrinkle: I expected the version based on Theorem 4 to be faster than the double-double implementation. But in fact, when I measure the performance of the two functions implemented in ufunclab, the double-double version is often faster!

That needs more investigation, but for now, I think the version based on Theorem 4 is very good. We could implement that while I spend more time than is reasonable on figuring out if a follow-up based on the double-double implementation is worthwhile. 😄

@mfkasim1, are you interested in updating this pull request to implement what I called log1p_theorem4(z) above? FWIW, my C version is here. If not, I can start a new PR with the proposed function.

mfkasim1 · 2023-08-11T00:16:58Z

@WarrenWeckesser please go ahead with the new PR. I'm currently on holiday.

lezcano · 2023-08-11T09:21:10Z

Rather interesting that, after so much coming and going, the final formula is simply to compute log(1+z) with a "minor" tweak :)

This PR improves the implementation of `torch.log1p` for complex inputs as mentioned in issue #107022. The new implementation is based on the insights provided in numpy/numpy#22611 (comment). Pull Request resolved: #107100 Approved by: https://github.com/lezcano

This PR improves the implementation of `torch.log1p` for complex inputs as mentioned in issue pytorch#107022. The new implementation is based on the insights provided in numpy/numpy#22611 (comment). Pull Request resolved: pytorch#107100 Approved by: https://github.com/lezcano

BUG: Increase accuracy of log1p for small inputs

3c06dd4

github-actions bot added the 00 - Bug label Nov 17, 2022

rossbar reviewed Nov 17, 2022

View reviewed changes

mfkasim1 mentioned this pull request Nov 18, 2022

Added log1p for complex in c10 pytorch/pytorch#89214

Closed

mfkasim1 added 4 commits November 18, 2022 11:06

BUG: Handle overflow and underflow for log1p complex

c57b5a2

STY: split lines that are too long for log1p complex test

59ab5bd

TST: Added infj test for log1p complex

efd5c82

TST: Added decimal arg in testing small numbers

d559bca

seberg reviewed Nov 18, 2022

View reviewed changes

mfkasim1 added 2 commits November 18, 2022 12:27

TST: Using assert_array_almost_equal_nulp for small numbers

c1729cf

STY: Put a space after a comma

3ad9ec3

mfkasim1 mentioned this pull request Nov 28, 2022

Log1p for complex in CPU pytorch/pytorch#89691

Closed

seberg requested a review from WarrenWeckesser December 8, 2022 12:04

lezcano mentioned this pull request Aug 11, 2023

Improve log1p(complex) speed/accuracy pytorch/pytorch#107022

Closed

tringwald mentioned this pull request Aug 13, 2023

Improved log1p implementation for complex inputs pytorch/pytorch#107100

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

BUG: Increase accuracy of log1p for small inputs #22611

BUG: Increase accuracy of log1p for small inputs #22611

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

BUG: Increase accuracy of log1p for small inputs #22611

Are you sure you want to change the base?

BUG: Increase accuracy of log1p for small inputs #22611

Uh oh!

Conversation

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!