2

I know discrete convolutions are supposed to be associative. So if I have some array, x, then x * (x * (x * x)) should equal (x * x) * (x * x). But there are some situations where that doesn't happen.

Here is code that exercises that formula with four examples:

[1, 1]  # Works
[-54298238, 5998425]  # Works
[-95.3720828588315, 52.6296402253613]  # Works
[-94.25133703348938, 90.41999267567854]  # Broken

import numpy as np

def main():

    int_ok = np.array([1, 1], dtype=np.int)
    int_larger = np.array(
        [-54298238,5998425], dtype=np.int
    )

    float_ok = np.array(
        [-95.3720828588315, 52.6296402253613], dtype=np.float
    )
    float_broken = np.array(
        [-94.25133703348938, 90.41999267567854], dtype=np.float
    )

    fixtures = [int_ok, int_larger, float_ok, float_broken]

    for fixture in fixtures:
        reference = np.convolve(
            fixture,
            np.convolve(
                fixture,
                np.convolve(
                    fixture,
                    fixture
                )
            )
        )

        tmp = np.convolve(fixture, fixture)
        test = np.convolve(tmp, tmp)

        print('input', fixture)
        print('reference output', reference)
        print('output', test)
        all_equal = all(reference == test)
        print('all equal', all_equal)
        if not all_equal:
            print('error', np.abs(reference - test))


if __name__ == '__main__':
    main()

I assume this is due to some kind of numerical instability, but I can't quite figure out what's going on.

Does anyone have any ideas?

Thanks.

John P
  • 65
  • 7
  • 2
    Even plain addition (or multiplication) isn't associative with binary floating-point. I don't think there's any reason to expect something more complicated like convolution to be perfectly associative. – Mark Dickinson Sep 14 '19 at 15:45
  • 1
    Possible duplicate of [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) –  Sep 14 '19 at 15:47
  • @MarkDickinson That's a good point about floating point not being associative, but I would expect the results to at least be _close_, even if not exact. That's why I have that second check. I was thinking the differences were significant, but after calling numpy.allclose on the results, it returned True. So I guess they're not significantly different, numerically. – John P Sep 14 '19 at 21:40
  • 1
    @JohnP They're not just close; they're *extremely* close. At least on my machine, in your failing case, `reference` is out from `test` by exactly one unit in the last place (ulp) in three of the five entries, and exactly equal in the other two. In other words, you're getting a relative error of better than 2e-16 in all entries. Given the limited precision of floating-point, you couldn't really ask for anything better here. – Mark Dickinson Sep 15 '19 at 08:46

1 Answers1

0

Inspired by this comment Why does numpy.convolve not behave associatively? by Mark Dickinson, I reworked my original code so that I'm no longer assuming the results are strictly the same, but rather using numpy.allclose.

Here is the modified code:

import numpy as np


def main():

    int_ok = np.array([1, 1], dtype=np.int)
    int_larger = np.array(
        [-54298238,5998425], dtype=np.int
    )

    float_ok = np.array(
        [-95.3720828588315, 52.6296402253613], dtype=np.float
    )
    float_broken = np.array(
        [-94.25133703348938, 90.41999267567854], dtype=np.float
    )

    fixtures = [int_ok, int_larger, float_ok, float_broken]

    for fixture in fixtures:
        reference = np.convolve(
            fixture,
            np.convolve(
                fixture,
                np.convolve(
                    fixture,
                    fixture
                )
            )
        )

        tmp = np.convolve(fixture, fixture)
        test = np.convolve(tmp, tmp)

        if np.allclose(reference, test):
            print('OK', fixture)
        else:
            print('FAIL', fixture)


if __name__ == '__main__':
    main()

numpy.allclose is returning True, so I think this is the best solution I'm going to get.

John P
  • 65
  • 7