Add the capability to do adjoint transforms #633

mreineck · 2025-02-18T10:14:42Z

This is a first outline how I propose to add adjoint transforms; it is mainly meant as a basis for discussions and measurements. The computation of the actual adjoint is completely untested yet, but the standard functionality should still be OK, as far as my tests show.

@ahbarnett, @DiamonDinoia please let me know your thoughts on this one!

mreineck · 2025-02-18T10:35:28Z

Just to be clear: at this point, there is no interface support for the new functioality. I just want to show which kind of impact the new feature has on the existing implementation.

mreineck · 2025-02-18T16:22:18Z

I've added the C and Python interfaces, as well as basic Python unit tests.

mreineck · 2025-02-18T16:54:00Z

Test failures seem to be "near misses" in adjoint type 3 transforms. No idea why that direction should be less accurate, and why it only happens in some of the tests.
This needs more investigation, but the approach seems to work well in principle.

…type 3

mreineck · 2025-02-19T09:14:14Z

Implements #571 and #566 (on CPU).

…s. Perhaps 1e-6 is just a bit too close to the machine epsilon for single precision

mreineck · 2025-02-19T13:29:23Z

I think this is ready for "technical" review. If there is agreement that the change is desirable, I can try to

provide more language-specific interfaces
document the new function

mreineck

Added a few explanatory comments

mreineck · 2025-02-24T13:53:53Z

python/finufft/examples/guru1d1f.py

@@ -20,7 +20,7 @@
 strt = time.time()

 #plan
-plan = fp.Plan(1,(N,),dtype='single')


This change is unrelated to the PR, but without it, CI simply fails. I don't know why this hasn't caused issues so far.

It seems because we removed the conversion for real dtypes in https://github.com/flatironinstitute/finufft/pull/606/files

mreineck · 2025-02-24T13:54:12Z

python/finufft/examples/guru2d1f.py

@@ -34,7 +34,7 @@

 # instantiate the plan (note n_trans must be set here), also setting tolerance:
 t0 = time.time()
-plan = finufft.Plan(nufft_type, (N1, N2), eps=1e-4, n_trans=K, dtype='float32')


This change is unrelated to the PR, but without it, CI simply fails. I don't know why this hasn't caused issues so far.

Thanks for catching this, it seems all the *f.py examples(in python/finufft/examples) should change according to PR606? otherwise the check is_single_dtype(dtype)

finufft/python/finufft/finufft/_interfaces.py

Line 106 in d2e0ff7

is_single = is_single_dtype(dtype)

will fail.

mreineck · 2025-02-24T13:56:16Z

python/finufft/test/test_finufft_plan.py

@@ -86,6 +120,19 @@ def test_finufft3_plan(dtype, dim, n_source_pts, n_target_pts, output_arg):

    utils.verify_type3(source_pts, source_coefs, target_pts, target_coefs, 1e-6)

+    # test adjoint type 3
+    plan = Plan(3, dim, dtype=dtype, isign=-1, eps=1e-5)


I'm increasing eps from 1e-6 to 1e-5 here, because I get occasional failures with single precision otherwise. Given that 1e-6 is uncomfortably close to machine epsilon, I'm not too worried about this change.

I agree. type 3 errors are usually 2-3x bigger than type 1 or 2 at the same tolerance.

mreineck · 2025-02-24T13:57:59Z

python/finufft/test/utils.py

@@ -154,7 +154,7 @@ def verify_type1(pts, coefs, shape, sig_est, tol):

    type1_rel_err = np.linalg.norm(fk_target - fk_est) / np.linalg.norm(fk_target)

-    assert type1_rel_err < 25 * tol


Switching from assert to np.testing.assert_allclose here, because the latter will provide more information in case of failure, which speeds up debugging a lot.

mreineck · 2025-02-24T14:02:05Z

src/fft.cpp

      }
    }
  }
 #endif
 #else
-  p->fftPlan->execute(); // if thisBatchSize<batchSize it wastes some flops


This needs discussion: how do we want to deal with this situation? The trick used here is nice and simple, but it will make the adjoint call slower than the forward one.

I actually don't understand the trick or why it's slower - sorry for my stupidity. I'd also like to understand how FFTW is handled. Since your change to the fft.h adds the data pointer, can that work with FFTW too (and allow us to allocate-at-execute stage) ?

The trick is the equivalence:
fft_adjoint(x) = conj(fft(conj(x)))
(since we do not have a plan for fft_adjoint, we just mess around with the sign of the imaginary part of input and output to get te same effect as if we had switched isign).
The additional complex conjugate of input and output is of course not completely free, but it might be sub-dominant compared to the cost of the FFT itself.

And yes, we do allocate the temporary at the execute stage in this PR, also for FFTW.

ahbarnett · 2025-02-25T16:48:27Z

I just fixed py examples in 19a24c6

…

On Mon, Feb 24, 2025 at 9:30 AM Libin Lu ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In python/finufft/examples/guru2d1f.py <#633 (comment)> : > @@ -34,7 +34,7 @@ # instantiate the plan (note n_trans must be set here), also setting tolerance: t0 = time.time() -plan = finufft.Plan(nufft_type, (N1, N2), eps=1e-4, n_trans=K, dtype='float32') Thanks for catching this, it seems all the *f.py examples(in python/finufft/examples) should change according to PR606? otherwise the check is_single_dtype(dtype) https://github.com/flatironinstitute/finufft/blob/d2e0ff7cd3a7bd66384d961d647060f497fde13c/python/finufft/finufft/_interfaces.py#L106 will fail. — Reply to this email directly, view it on GitHub <#633 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZRSUEOSDPVX33Q6B3WTT2RMUHHAVCNFSM6AAAAABXLHSEQ2VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDMMZXGI4TSMZSGQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

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

ahbarnett

Hi Martin,

[This is supposed to go as a PR review but appeared as a comment]

Your work here is great. As usual you are able to slot something in with minimal disruption. It will certainly allow power users to slip in the adjoint when needed. I don't understand how the FFTW part works - there seem to be no changes which would allow the fftw_execute to change type (or handle the new data ptr in the fft.h interface) - or did I miss it?

However, I am thinking about the interface from the new user perspective, and think it is getting more confusing than needs be. In your PR a user plans a type 1, then can do its adjoint (type 2 with flipped isign) freely. Users will want to know if that's the same a planning a type 2 and doing a forward - it's not due to isign. (Extra confusion: in Europe t1 is called "adjoint" and t2 "forward" NUDFT.) I see the isign flip between going from 1 to 2 vs 1 to its adj being a major source of confusion. The motivation for your PR was users wanting a single plan for t1 and t2, or for t1 and its adj. I think we should bounce around ideas for a new interface that fits the need. An idea to reduce confusion is to have a "plan 1 and 2", and then execute_1(c,f,isign) execute_2(c,f,isign). No confusion. 1+ is adj of 2-, and 1- adj of 2+. Question is can we do it in a neat way that preserves the existing guru interface (add a type=12 to the plan?) since legacy users need to still be able to plan a type 2 explicitly and have plan execute(c,f) do that t2.
Other questions are: 4 if FFTW plan can handle the isign swtiching?
This would all be for the plan (guru) interface. It would expand that interface from 4 to six commands. Under the hood (as you do in this PR with execute_internal()), execute() could call execute_1 or _2.

I have lots of deadlines in the next month so will have to wait a bit, but I'd love to discuss this and hash out the best interface, since I don't want a confusing interface to become locked in...

Best, Alex

ahbarnett · 2025-04-04T14:38:32Z

src/fft.cpp

      }
    }
  }
 #endif
 #else
-  p->fftPlan->execute(); // if thisBatchSize<batchSize it wastes some flops


I actually don't understand the trick or why it's slower - sorry for my stupidity. I'd also like to understand how FFTW is handled. Since your change to the fft.h adds the data pointer, can that work with FFTW too (and allow us to allocate-at-execute stage) ?

ahbarnett · 2025-04-04T14:39:58Z

src/finufft_core.cpp

@@ -424,72 +425,76 @@ static void deconvolveshuffle3d(int dir, T prefac, std::vector<T> &ker1,
 // --------- batch helper functions for t1,2 exec: ---------------------------

 template<typename T>
-static int spreadinterpSortedBatch(int batchSize, FINUFFT_PLAN_T<T> *p,
-                                   std::complex<T> *fwBatch, std::complex<T> *cBatch)
+static int spreadinterpSortedBatch(int batchSize, const FINUFFT_PLAN_T<T> &p,


THis is just part of your ongoing cleanup, right?

Partly yes. But the change to const also ensures that we don't do things like

innerT2plan->ntrans = thisBatchSize;

any more (which we had in execute() before). This sort of "messing with the state of another object" can cause a lot of trouble once plans are perhaps invoked in parallel ... and that's what I wanted to make possible with this patch.
The switch from pointer to reference is cosmetics, and I should perhaps have delayed that until later ... sorry,

ahbarnett · 2025-04-04T14:41:04Z

src/finufft_core.cpp

  This is mostly a loop calling deconvolveshuffle?d for the needed dim, batchSize
  times.
  Barnett 5/21/20, simplified from Malleo 2019 (eg t3 logic won't be in here)
 */
 {
  // since deconvolveshuffle?d are single-thread, omp par seems to help here...
+  int dir = p.spopts.spread_direction;
+  if (adjoint) dir = 3 - dir;


Ok, there's the meat - please add some description of what this obscure subtraction does :) (I know, but most readers will not...) Thanks!

You are totally right, I'll do that!

ahbarnett · 2025-04-04T15:06:57Z

src/fft.cpp

 #endif
 }
-template void do_fft<float>(FINUFFT_PLAN_T<float> *p);
-template void do_fft<double>(FINUFFT_PLAN_T<double> *p);
+template void do_fft<float>(const FINUFFT_PLAN_T<float> &p, std::complex<float> *fwBatch,


Is there an FFTW version of passing the data ptr like this?

https://fftw.org/fftw3_doc/New_002darray-Execute-Functions.html#New_002darray-Execute-Functions

ahbarnett · 2025-04-04T15:07:40Z

src/finufft_core.cpp

  /* See ../docs/cguru.doc for current documentation.

   For given (stack of) weights cj or coefficients fk, performs NUFFTs with
   existing (sorted) NU pts and existing plan.
-   For type 1 and 3: cj is input, fk is output.
-   For type 2: fk is input, cj is output.
+   For adjoint == false:


these comments are super helpful

ahbarnett · 2025-04-04T15:08:58Z

src/finufft_core.cpp

+     - 0 < ntrans_actual <= batchSize: instead of doing ntrans transforms,
+         perform only ntrans_actual
+
+   scratch_size, aligned_scratch:


Is this something we should do anyway with t3? (regardless of adjoint?)

ahbarnett · 2025-04-04T15:09:15Z

src/finufft_core.cpp

+
+   Additional parameters introducd by MR, 02/2025
+
+   adjoint: if false, behave as before; if true, compute the adjoint of the


see my overall PR review discussion

I'm very happy to adjust the interface in any way that makes it easier for finufft users to get used to. This PR is (at the moment) mostly a demonstration of how things can be done. How they are presented is always up to the Finufft team.

mreineck · 2025-04-14T09:25:03Z

When a new interface is perhaps on the horizon, and there are, as @ahbarnett points out, different conventions concerning what is type 1 and type 2, perhaps using descriptive names may be a way out, like u2nu, nu2u, nu2nu. Combined with an isign argument, this would solve almost all issues with adjoint invocation os plans.
The one thing it does not solve, unfortunately, is how to tell the library to run a nu2nu transform "backwards" ...

first version; adjoint still completely untested

1180e27

mreineck added 2 commits February 18, 2025 17:13

add Python interface; fix a few bugs

e7bd659

add tests; add overlooked file

fe0dcdc

use assertions that are more helpful in cse of errors

e148935

test commit to see if tests pass with incresed tolerance for adjoint …

b7d7ce1

…type 3

mreineck added 2 commits February 19, 2025 10:42

increase overall requested plan accuracy for adjoint type 3 transform…

f54a1dd

…s. Perhaps 1e-6 is just a bit too close to the machine epsilon for single precision

be more clever about memory consumption

d50a912

mreineck added 3 commits February 24, 2025 11:48

document new parameters for execute_internal()

a4c313f

Merge remote-tracking branch 'origin/master' into add_adjoint

2692658

update CHANGELOG

eafcf3e

mreineck commented Feb 24, 2025

View reviewed changes

mreineck added 4 commits February 26, 2025 10:10

merge master

7f19185

more comments

0fa2a4e

small tweak

da3df92

merge master

ec969b5

ahbarnett mentioned this pull request Mar 20, 2025

consider whether the big array fwbatch should be malloced at plan or at execute. #480

Open

mreineck added 3 commits March 22, 2025 09:12

fix typos

5a61c86

merge master

9bc0e02

maerge master

3dbd742

ahbarnett reviewed Apr 4, 2025

View reviewed changes

mreineck added 2 commits April 14, 2025 10:00

merge master

ee21762

add explanations for obscure tricks

d17dc7a

merge master

886cdf9

ahbarnett mentioned this pull request Apr 30, 2025

towards release 2.5 to-do list #670

Open

6 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add the capability to do adjoint transforms #633

Add the capability to do adjoint transforms #633

mreineck commented Feb 18, 2025

mreineck commented Feb 18, 2025

mreineck commented Feb 18, 2025

mreineck commented Feb 18, 2025

mreineck commented Feb 19, 2025 •

edited

Loading

mreineck commented Feb 19, 2025

mreineck left a comment

mreineck Feb 24, 2025

lu1and10 Feb 24, 2025

mreineck Feb 24, 2025

lu1and10 Feb 24, 2025

mreineck Feb 24, 2025

ahbarnett Apr 4, 2025

mreineck Feb 24, 2025

ahbarnett Apr 4, 2025

mreineck Feb 24, 2025

ahbarnett Apr 4, 2025

mreineck Apr 14, 2025

mreineck Apr 14, 2025

ahbarnett commented Feb 25, 2025 via email

ahbarnett left a comment •

edited

Loading

ahbarnett Apr 4, 2025

ahbarnett Apr 4, 2025

mreineck Apr 14, 2025 •

edited

Loading

ahbarnett Apr 4, 2025

mreineck Apr 14, 2025

ahbarnett Apr 4, 2025

DiamonDinoia Apr 4, 2025

ahbarnett Apr 4, 2025

ahbarnett Apr 4, 2025

ahbarnett Apr 4, 2025

mreineck Apr 14, 2025

mreineck commented Apr 14, 2025

		@@ -154,7 +154,7 @@ def verify_type1(pts, coefs, shape, sig_est, tol):

		type1_rel_err = np.linalg.norm(fk_target - fk_est) / np.linalg.norm(fk_target)

		assert type1_rel_err < 25 * tol


		Additional parameters introducd by MR, 02/2025

		adjoint: if false, behave as before; if true, compute the adjoint of the

Add the capability to do adjoint transforms #633

Are you sure you want to change the base?

Add the capability to do adjoint transforms #633

Conversation

mreineck commented Feb 18, 2025

mreineck commented Feb 18, 2025

mreineck commented Feb 18, 2025

mreineck commented Feb 18, 2025

mreineck commented Feb 19, 2025 • edited Loading

mreineck commented Feb 19, 2025

mreineck left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ahbarnett commented Feb 25, 2025 via email

ahbarnett left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mreineck Apr 14, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mreineck commented Apr 14, 2025

mreineck commented Feb 19, 2025 •

edited

Loading

ahbarnett left a comment •

edited

Loading

mreineck Apr 14, 2025 •

edited

Loading