search

EMS-TU-Ilmenau / fastmat

A library to build up lazily evaluated expressions of linear transforms for efficient scientific computing.
https://fastmat.readthedocs.io
Apache License 2.0
24 stars 8 forks source link

Wrong results for large Fourier with Bluestein algorithm #115

Closed ChristophWWagner closed 1 year ago

ChristophWWagner commented 1 year ago

There seems to be a serious bug with the Bluestein implementation of the Fourier class.

MWE:

>>> import fastmat as fm, numpy as np
>>> N, ee = 51187, 100 # prime factorization is [17, 3011], limit to some elements
>>> np.linalg.norm(fm.Fourier(N)[:ee, :ee] - fm.Fourier(N, optimize=False)[:ee, :ee])
97.38713873069028

The problem also occurs for even Fourier orders:

>>> N, ee = 51188, 100 # prime factorization is [2, 2, 19, 67, 191], limit to some elements
>>> np.linalg.norm(fm.Fourier(N)[:ee, :ee] - fm.Fourier(N, optimize=False)[:ee, :ee])
97.3871896159737

For "simpler" prime factorizations of N, the problem does not persist:

>>> N, ee = 2 ** 8 * 7 * 5 * 11, 100
>>> np.linalg.norm(fm.Fourier(N)[:ee, :ee] - fm.Fourier(N, optimize=False)[:ee, :ee])
0.0

Checking whether Bluestein is invoked yields:

>>> fm.Fourier(51187)._numL
110592
>>> fm.Fourier(51188)._numL
110592
>>> fm.Fourier(98560)._numL
0

In conclusion, the problem occurs only when Bluestein is involved.

Next steps

  1. A test should be added that successfully catches this condition
  2. Investigate why this happens
  3. Fix it
  4. Immediate release of new subversion
ChristophWWagner commented 1 year ago

Update: Bluestein alone is not fully responsible for the trouble, since for small orders it seems to still work:

>>> fm.Fourier(127)[...] - fm.Fourier(127, optimize=False)[...]
array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       ...,
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]])
>>> fm.Fourier(127, optimize=False)._numL
0
>>> fm.Fourier(127, optimize=True)._numL
256
>>> N = 127; np.linalg.norm(fm.Fourier(N)[...] - fm.Fourier(N, optimize=False)[...])
0.0
>>> N = 191; np.linalg.norm(fm.Fourier(N)[...] - fm.Fourier(N, optimize=False)[...])
0.0
>>> N = 3011; np.linalg.norm(fm.Fourier(N)[...] - fm.Fourier(N, optimize=False)[...])
0.0
>>> N = 127*7; np.linalg.norm(fm.Fourier(N)[...] - fm.Fourier(N, optimize=False)[...])
0.0
ChristophWWagner commented 1 year ago

The previous update is misleading as it falls victim to the implicit overriding of the [...] operator by the internal override in ._getArray() (see https://github.com/EMS-TU-Ilmenau/fastmat/blob/f0e3fe95a619366c6f61e9af1ad360683b6b9398/fastmat/Fourier.pyx#L165) It seems to actually fail also for the smaller values (127, 191, 3011 included, 127*7 does not trigger Bluestein)

ChristophWWagner commented 1 year ago

Good news, everyone! Looks like reverting #d22813f0 resolves the error. The suggested fix is the following: 1) make sure the tests actually cover a Bluestein case 2) Fix the underlying issue of #d22813f0. Upon first glance for a test case of 127 there seems to be a mismatch of one or multiples factors-of-two