nufft precision

[PascalCarrivain]

Dear all,

I am trying to implement a Python version of the algorithm from the paper A nonuniform fast Fourier transform based on low rank approximation.
I think nufft.jl from FastTransforms.jl uses the same algorithm.

I first did a comparison between nufft and its naive implementation for a given x and uniform samples samples.
As expected I get the same result for k=1.

julia> x = complex([-0.1, 0.2, 0.3, 0.01, -1.0, 4.0])
6-element Vector{ComplexF64}:
 -0.1 + 0.0im
  0.2 + 0.0im
  0.3 + 0.0im
 0.01 + 0.0im
 -1.0 + 0.0im
  4.0 + 0.0im

julia> samples = collect(0:6-1) / 6
6-element Vector{Float64}:
 0.0
 0.16666666666666666
 0.3333333333333333
 0.5
 0.6666666666666666
 0.8333333333333334

julia> sum(x .* exp.((2.0 * im * pi * k) .* samples)) / sqrt(6)
0.9553009996854398 - 0.8838834764831844im

julia> f = getindex(nufft2(x, samples, eps()) / sqrt(6), k + 1)
0.9553009996854397 + 0.8838834764831844im

However, for nonuniform samples I get an error:

julia> samples = [0.05, 0.051, 0.3, 0.47, 0.65, 0.9]
6-element Vector{Float64}:
 0.05
 0.051
 0.3
 0.47
 0.65
 0.9

julia> sum(x .* exp.((2.0 * im * pi * k) .* samples)) / sqrt(6)
1.557891331279937 - 0.49922137697035407im

julia> f = getindex(nufft2(x, samples, eps()) / sqrt(6), k + 1)
-0.03047518500012294 - 1.3431665848287448im

Do I miss something or the error is expected ?
Do you have a way to estimate the error for each entry of the output ?

Thank you.

Best,
Pascal

[MikaelSlevinsky]

If you look at equation (1.2) in the paper, you will find that your sum above is not looking at the k'th entry in the output of the matrix-vector product. For a comparison, you should see correctness like this:

julia> k = 1
1

julia> x = complex([-0.1, 0.2, 0.3, 0.01, -1.0, 4.0])
6-element Vector{ComplexF64}:
 -0.1 + 0.0im
  0.2 + 0.0im
  0.3 + 0.0im
 0.01 + 0.0im
 -1.0 + 0.0im
  4.0 + 0.0im

julia> samples = collect(0:6-1) / 6
6-element Vector{Float64}:
 0.0
 0.16666666666666666
 0.3333333333333333
 0.5
 0.6666666666666666
 0.8333333333333334

julia> sum(exp.(-(2.0 * im * pi * (0:5)) .* samples[k+1]) .* x)
2.339999999999998 + 2.165063509461098im

julia> nufft2(x, samples, eps())[k+1]
2.3400000000000003 + 2.1650635094610964im

julia> samples = [0.05, 0.051, 0.3, 0.47, 0.65, 0.9]
6-element Vector{Float64}:
 0.05
 0.051
 0.3
 0.47
 0.65
 0.9

julia> sum(exp.(-(2.0 * im * pi * (0:5)) .* samples[k+1]) .* x)
-0.07464865306722146 - 3.2900727723871217im

julia> nufft2(x, samples, eps())[k+1]
-0.0746486530672209 - 3.2900727723871217im

I think you saw a fluke in the uniform test because samples and frequencies are transposes (apart from a scaling).