Strange discrepancy in correction factor calculations

[mreineck]

Lines

finufft/src/finufft_core.cpp

Line 198 in 4565d4b

int q = (int)(2 + 3.0 * J2); // not sure why so large? cannot exceed MAX_NQUAD

and

finufft/src/finufft_core.cpp

Line 258 in 4565d4b

int q = (int)(2 + 2.0 * J2); // > pi/2 ratio. cannot exceed MAX_NQUAD

read slightly differently, but I don't really understand why they should be different at all. The first version agrees with the FINUFFT paper (see the line below eq 3.10), but the reduced accuracy caused by the second version doesn't seem to cause any test failures either ... strange.

Given the comment in the source, I'm apparently not the first one wondering about this.

[ahbarnett]

Hi Martin, Well, source comments are by me, I expect. I imagine, after paper was published, I wanted to speed up type-3 by lowering it until errors started growing. But for type1,2 there was little incentive to lower it because the winding formula makes it very fast. If you have new evidence that both should be tweaked to improve error and/or speed, make a PR :) Thanks! Alex

…

On Mon, Jan 6, 2025 at 8:07 AM mreineck ***@***.***> wrote: Lines https://github.com/flatironinstitute/finufft/blob/4565d4b60a3b69b598cc5ff89712986a34d43fbc/src/finufft_core.cpp#L198 and https://github.com/flatironinstitute/finufft/blob/4565d4b60a3b69b598cc5ff89712986a34d43fbc/src/finufft_core.cpp#L258 read slightly differently, but I don't really understand why they should be different at all. The first version agrees with the FINUFFT paper (see the line below eq 3.10), but the reduced accuracy caused by the second version doesn't seem to cause any test failures either ... strange. Given the comment in the source, I'm apparently not the first one wondering about this. — Reply to this email directly, view it on GitHub <#601>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZRSVC3SKTUGRT3CALIML2JJ5ZLAVCNFSM6AAAAABUVSMNV2VHI2DSMVQWIX3LMV43ASLTON2WKOZSG43TANRSGUZTCNI> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- *-------------------------------------------------------------------~^`^~._.~' |\ Alex Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942

[mreineck]

Hi Alex, thanks a lot for the explanation! You are right, there is no need for speeding up the correction factor calculation for types 1 amd 2, where it should require negligible time (except perhaps in 1D).
It's true, I'm thinking about ways to make the precalculation faster for type 3, since its cost is not too far away from the spreading/interpolatiion cost itself, of course with the difference that it is only required during plan generation and not furing every execution.
One thing I have been wondering about: in the range where we need the correction function, it is pretty smooth and well-behaved, right? So perhaps we could just use the (piecewise?) polynomial approximation trick there as well? This could give about an order of magnitude acceleration, but of course requires lots of testing to gather sufficient confidence.

[SepandKashani]

One thing I have been wondering about: in the range where we need the correction function, it is pretty smooth and well-behaved, right? So perhaps we could just use the (piecewise?) polynomial approximation trick there as well? This could give about an order of magnitude acceleration, but of course requires lots of testing to gather sufficient confidence.

I have tried using the ppoly approximator to compute the deconvolution factors in fourier_toolkit, and I could not see any significant rel-error difference compared to the analytic form. I did not use it for this task in the end to avoid having yet another knob in the system, but I confirm it is doable.

You can test this easily by replacing this line https://github.com/SepandKashani/fourier_toolkit/blob/15dfbdb260cd7959c64088de21d0a67bd5229c77/src/fourier_toolkit/nu2nu.py#L1127

phiF = ftk_kernel.KaiserBesselF(beta)

with

phiF = ftk_kernel.PPoly.from_kernel(
    ftk_kernel.KaiserBesselF(beta),
    B=B_ppoly,
    N=N_ppoly,
    sym=True,
)

[mreineck]

Thank you everyone! I'm convinced now that this is not a problem. Experimenting with potential speedups can be done as a separate task in the future.