norm="ortho" does not take into account the actual transformed axes

[vincefn]

Using norm="ortho" the array norm should be preserved but it is not:

import numpy as np
import mkl_fft

s = 16,18
ax = 0,

a = np.random.uniform(0,1,s) + 1j * np.random.uniform(0,1,s)
a /= np.sqrt((abs(a)**2).sum()) # normalise

a_np = np.fft.fftn(a, axes=ax, norm="ortho")
a_mkl = mkl_fft.fftn(a, axes=ax, norm="ortho")

print(f"  Array norm: {(abs(a)**2).sum():.6f}")
print(f" np-fft norm: {(abs(a_np)**2).sum():.6f}")
print(f"mkl_fft norm: {(abs(a_mkl)**2).sum():.6f}")

Gives:

  Array norm: 1.000000
 np-fft norm: 1.000000
mkl_fft norm: 0.055556

The ratio of the norms here is 18, the size of the non-transformed axis. So norm="ortho" does not take into account the actual list of transformed axes.

Tested with numpy 2.3.5 and mkl_fft 2.2.1, all from conda-forge.

[vincefn]

I made more tests - the same behaviour is seen with fftn, ifftn and rfftn. However irfftn does not have the issue - the norm is correctly computed.

Here are snippets to compare numpy and mkl normalisation, for even and odd sizes (only matters for r2c):

for c2c:

import numpy as np
import mkl_fft

for s in [(16, 18), (15,17)]:
    ax = -1,

    a = np.random.uniform(0,1,s)
    a /= np.sqrt((abs(a)**2).sum()) # normalise
    a0 = a.copy()

    a_np = np.fft.fftn(a, axes=ax, norm="ortho")
    n_np = (abs(a_np) ** 2).sum()

    a_mkl = mkl_fft.fftn(a, axes=ax, norm="ortho")
    n_mkl = (abs(a_mkl) ** 2).sum()

    a_mkl = mkl_fft.ifftn(a_np, axes=ax, norm="ortho")
    n_mkli = (abs(a_mkl) ** 2).sum()

    a_np = np.fft.ifftn(a_np, axes=ax, norm="ortho")
    n_npi = (abs(a_np) ** 2).sum()


    print(f"   Array norm: {(abs(a)**2).sum():.6f}")
    print(f"  np-fft norm: {n_np:.6f}")
    print(f" mkl_fft norm: {n_mkl:.6f}")
    print(f" np-ifft norm: {n_npi:.6f}")
    print(f"mkl_ifft norm: {n_mkli:.6f}")

Gives:

   Array norm: 1.000000
  np-fft norm: 1.000000
 mkl_fft norm: 0.062500
 np-ifft norm: 1.000000
mkl_ifft norm: 16.000000
   Array norm: 1.000000
  np-fft norm: 1.000000
 mkl_fft norm: 0.066667
 np-ifft norm: 1.000000
mkl_ifft norm: 15.000000

for r2c

import numpy as np
import mkl_fft

for s in [(16, 18), (15,17)]:
    ax = -1,

    a = np.random.uniform(0,1,s)
    a /= np.sqrt((abs(a)**2).sum()) # normalise
    a0 = a.copy()

    a_np = np.fft.rfftn(a, axes=ax, norm="ortho")
    print(s, "->", a_np.shape)
    if s[-1] % 2 != 0:
        n_np = (abs(a_np) ** 2).sum() + (abs(a_np[..., 1:]) ** 2).sum()
    else:
        n_np = (abs(a_np) ** 2).sum() + (abs(a_np[..., 1:-1]) ** 2).sum()

    a_mkl = mkl_fft.rfftn(a, axes=ax, norm="ortho")
    if s[-1] % 2 != 0:
        n_mkl = (abs(a_mkl) ** 2).sum() + (abs(a_mkl[..., 1:]) ** 2).sum()
    else:
        n_mkl = (abs(a_mkl) ** 2).sum() + (abs(a_mkl[..., 1:-1]) ** 2).sum()

    a_mkl = mkl_fft.irfftn(a_np, axes=ax, norm="ortho", s=[s[i] for i in ax])
    n_mkli = (abs(a_mkl) ** 2).sum()

    a_np = np.fft.irfftn(a_np, axes=ax, norm="ortho", s=[s[i] for i in ax])
    n_npi = (abs(a_np) ** 2).sum()
    print(a_np.shape)


    print(f"   Array norm: {(abs(a)**2).sum():.6f}")
    print(f"  np-rfft norm: {n_np:.6f}")
    print(f" mkl_rfft norm: {n_mkl:.6f}")
    print(f" np-irfft norm: {n_npi:.6f}")
    print(f"mkl_irfft norm: {n_mkli:.6f}")

gives:

(16, 18) -> (16, 10)
(16, 18)
   Array norm: 1.000000
  np-rfft norm: 1.000000
 mkl_rfft norm: 0.062500
 np-irfft norm: 1.000000
mkl_irfft norm: 1.000000
(15, 17) -> (15, 9)
(15, 17)
   Array norm: 1.000000
  np-rfft norm: 1.000000
 mkl_rfft norm: 0.066667
 np-irfft norm: 1.000000
mkl_irfft norm: 1.000000