Improved spreading kernel, tolerance matching and testing, for v2.5.0 #798
Replies: 18 comments 64 replies
|
@ahbarnett I did a quick experiment, without tweaking upsampling 1.25: upsampling 2: It might be worth optimizing the c(beta) formula or tol fudge factor for pswf and including pswf in future releases. Updates: upsampling 1.25: upsampling 2.0: While a little bit tweak on pswf c for upsample fact 1.25 and low ns gives: It still might be worth optimizing the c(beta) formula or tol fudge factor for pswf and including pswf in future releases. |
|
I didn't find the pswf in std lib, I plugged in the esp code's pswf in my local machine, I can clean a bit and push to the github though. I use your c_beatty for KB, only tweak the upsampling 1.25, ns < 6, turning beta a bit larger. Not too much improvement with KB Beatty's parameters. Maybe there's some math that can derive some optimal parameters for pswf, which I don't know. The ekw prints And |
@ahbarnett I pushed the pswf to a test branch as kernel type 7 with minimal changes, adding a separate large file pswf.cpp...
I didn't commit any changes to the logic of ker ns and beta chosen. In the branch, pswf uses the same parameters as KB. Another thing I noticed is a slight error spike in pswf ns 15, as in the following plot, to remove the spike, I have to change back to I didn't commit to the banch back to ns - 1, the branch still has ns - 4, if you play with pswf, watch out for this. 
wsweepkerrcomp.m diff to reproduce the plot diff --git a/matlab/test/wsweepkerrcomp.m b/matlab/test/wsweepkerrcomp.m
index 58d530db..822a5c27 100644
--- a/matlab/test/wsweepkerrcomp.m
+++ b/matlab/test/wsweepkerrcomp.m
@@ -19,12 +19,12 @@ dim = 1; Ntot = 300; % which dimensionality to test, tot #modes (not too small)
%dim = 3; Ntot = 1e3;
ntr = 20; % #transforms to average error over at each tol (was 10)
isign = +1;
-sigma = 2.0; %1.25;
+sigma = 1.25;
tolsperdecade = 10;
tolstep = 10 ^ (-1 / tolsperdecade); % multiplicative step in tol, < 1
% following names must match kernel.hpp:
-kfnam = {"ES legacy", "ES Beatty", "KB Beatty", "cont-KB Beatty"};
-kfs = [1 3]; %1:4; % kernel formulae to test
+kfnam = {"ES legacy", "ES Beatty", "KB Beatty", "cont-KB Beatty", "5", "6", "pswf"};
+kfs = [1 3 7]; %1:4; % kernel formulae to test
o.upsampfac = sigma;
%o.debug = 0;
@@ -87,7 +87,7 @@ figure('name','mean rel err vs tol, comparing kernels, for 3 NUFFT types',...
for y=1:3 % types
subplot(1,3,y);
legs = {};
- symb = '+.';
+ symb = '+.o';
tt = tols(toloks); % plot only the non-warning tol domain
for i=1:nkf % kernels
loglog(tt, squeeze(errs(i,y,toloks)), symb(i)); % 'markersize',10); |
|
got it - I brought in your branch as PR and see the ns=15 spike... upsetting. It happens when nc_fit=18 (trim to nc=12). |
I'm not sure, it seems in the coefficient filtering, the first coeff that is less than cutoff * tol is found, then the loop breaks and all the remaining coeffs filtered, is it possible later coeffs may still larger than cutoff * tol and being filtered mistakenly?( Lines 577 to 582 in c6125b1 In the esp code, I loop over all the coeffs(no middle break) and find the last(highest degree) coeff that is larger than cutoff * tol, mark that as cutoff for filtering? I was hoping there might be similar analysis for pswf as Beatty's beta formula for KB, i.e, find the first crossing zero of prolate Fourier transform. May be that's not so easy, if we don't have any good formula, we may just use Marco's parameter sweep + regression fit. |
|
I chedked, it seems no problem with the coeffs filtering, should filter to nc=12, not sure why the poly coeffs are small while the error is spike... |
|
I now summary improvements in accuracy testing (originally brought in #780, but fully explored in MATLAB in #787). [Updated 1/24/26 for PWSF kernel] First I discuss one improvement in the ns (width w) choice logic that prevents accuracy loss at small requested tolerance (eps) in single precision, due to the dynamic range (r_dyn) being large at small sigma (upsampfac). The way to visualize all this is with the new 
The pink circles are runs where the 
In 3D note that even in double-prec sigma=1.25, in type 3 the r_dyn-related catastrophic cancellation causes the floor to have to be raised to 3e-8 (this is still deemed acceptable for this low sigma): 
Note that the type 1 and type 2 errors are misleadingly good at small tolerances in the upper right plot above - this is because the numbers of modes in each dimension is small (around 13, chosen for rapid direct evaluations), so that the find grid size restriction nf>=2w leads to better than typical accuracies. Normally in 3D you would expect similar accuracies to the above 2D and 1D plots, which show a close match to tol itself. These dotted lines are exactly the error criteria used in the C++ CI test |
|
I used @lu1and10 pswf and tweaked beta like this: const double s = sigma;
const double n = static_cast<double>(ns);
const double s2 = s * s;
const double A = -0.19638654 + 2.31685991 * s - 0.53110991 * s2;
const double B = 2.29051829 - 2.82937718 * s + 0.91381927 * s2;
const double C = -0.61525503;
spopts.beta = A * n + B + C / n;I did a sweep over beta,sigma and ns and then tried to fit a something in the (separable) form: A * n + B + C/nThen rewrote this as a LS problem and solved it. Hence, the coeffs there. The results (my matlab plots dark IDK why). |
|
Sorry, ignore my test from last night of Marco's. I had faulty logic that caused kf=7 always to be used 👎 And I will redo the coefficient drop-off detection logic - Libin is correct that Marco's logic would fail when a coefficient is small followed by larger later ones (odd ns triggers this by symmetry reasons - odd poly coeffs vanish - amazing we didn't catch this before...). |
I double checked last night, I think I was wrong, coeffs are ordered from high degree to low degree. |
|
And I don't know the reason why drop off tol * 0.05 coeffs gives the spike. In the esp code, I was using coeff/max_coeff the relative ratio to compare with cut off tolerance which I think it's too strict. |
double t = (double)ns * (1.0 - 1.0 / (2.0 * sigma));
spopts.beta = ((-0.00149087 * t + 0.0218459) * t + 3.06269) * t - 0.0365245;Using t this fitting seems to work well. |
Thanks, very nice. I was hoping beta could be linear function of t. |
|
Please note I have updated the top two Discussion boxes for the PSWF kernel, dated today (1/24/26). We don't have the bandwidth for any more kernel logic changes before v2.5.0-rc1. (Apart from if CI fails on some platforms I may have to loosen the tolsweep slack & floor params). We are just going to do stability tests. |
|
Just reporting from some downstream projects (nifty-ls and jax-finufft): everything seems okay with v2.5.0-rc1! In nifty-ls, there seemed to be a minor decrease in accuracy in some tests (a fraction of a digit), so I just relaxed the test tolerance. In terms of performance, everything seemed to be about the same. In one latency-sensitive CPU test there was a degradation from 70 ms to 100 ms, but for most tests there was no difference. I didn't do any performance tests for jax-finufft, but the test suite passed without any issues (CPU and GPU). |
|
Hi @lgarrison , thanks for testing out. It will be great to know in what cases the accuracy decreases. If it's on the CPU end, there is a new developing option called |
|
The particular test case that's failing is one that compares a periodogram computed with two different code paths. One runs mostly in Python and uses the finufft Python API, and one runs mostly in C++ and uses the finufft C API. The periodograms have plenty of other operations done to them besides the finufft step that could amplify tiny differences. So I shouldn't have said that the accuracy "decreased", just that now the two paths are producing results that are slightly more different than before. I did pull out a minimal example from the Python API path, which shows that a particular transform with eps=1e-9 differs by more than that amount (about 1e-7) in 2.4.1 and 2.5.0rc1. Perhaps that's expected, or maybe my definition of a "difference" is different from what finufft uses. But maybe this is a good place to start. Here's the example: Example
# /// script
# requires-python = ">=3.11"
# dependencies = [
# "finufft",
# "numpy",
# ]
# ///
"""
We saved off the exact inputs and outputs of a particular finufft call that nifty-ls is making.
This script verifies the inputs are identical, then prints the max difference in the outputs.
Finally, it reproduces the finufft call the code is making to verify that the outputs are reproducible.
Recommended usage:
$ uv run compare.py
$ uv run --with "finufft==2.5.0rc1" compare.py
"""
import finufft
import numpy as np
input24 = np.load("debug_finufft_inputs-v2.4.1.npz")
input25 = np.load("debug_finufft_inputs-v2.5.0rc1.npz")
output24 = np.load("debug_finufft_outputs-v2.4.1.npz")
output25 = np.load("debug_finufft_outputs-v2.5.0rc1.npz")
def max_rel_diff(a, b):
r = np.abs(a - b) / np.abs(a)
i = np.argmax(r)
m = r.flat[i]
if m == 0:
return m
return (m, (a.flat[i], b.flat[i]))
print("Max relative differences between v2.4.1 and v2.5.0-rc1:")
# Inputs are all identical
print("t1:", max_rel_diff(input24["t1"], input25["t1"]))
print("yw_w:", max_rel_diff(input24["yw_w"], input25["yw_w"]))
print("t2:", max_rel_diff(input24["t2"], input25["t2"]))
print("w2:", max_rel_diff(input24["w2"], input25["w2"]))
print("f1:", max_rel_diff(output24["f1"], output25["f1"]))
# print('fw:', max_rel_diff(output24['fw'], output25['fw']))
print("f2:", max_rel_diff(output24["f2"], output25["f2"]))
# now reproduce
plan_solo = finufft.Plan(
nufft_type=1,
n_modes_or_dim=(1000,),
n_trans=4,
dtype=np.complex128,
nthreads=1,
eps=1e-9,
upsampfac=1.25,
)
plan_solo.setpts(input24["t1"])
f1_fw = plan_solo.execute(input24["yw_w"])
f1 = f1_fw[:4]
# fw = f1_fw[4:]
print('Matches 2.4.1 exactly: ', np.all(f1 == output24["f1"]))
print('Matches 2.5.0-rc1 exactly: ', np.all(f1 == output25["f1"]))The data files are just a few hundred KB, so I'll upload them here for provenance. |
|
Thanks Lehman, and thanks Libin.
I will add that our accuracy is much closer to user-requested tol, as
documented in
#798 (comment)
Look at the "error vs eps" plots, eg for 1D type 3.
Also see the new error-matching plots here:
#798 (comment)
It is possible that you were sitting at an eps for which it became less
accurate (eg for upsampfac=2, if you chose eps=0.03 then v2.4.1 chose w=3
and got around err=0.01), even though it should be uniformly no worse than
about 2*eps. (ie, now for eps=0.03 v2.5.0 would choose w=2 because of the
better kernel - hence be faster - but maybe the error is 0.04). For the
old kernel/logic (as Libin says, kerformula=1), there were eps values for
which the error was higher than 2*eps.
We'd be interested in the latency-intensive example for 70ms going up to
100ms. I wonder if it's the move to allocate-upon-execute?
Cheers, Alex
…On Wed, Feb 11, 2026 at 7:16 PM Libin Lu ***@***.***> wrote:
Hi @lgarrison
<https://urldefense.com/v3/__https://github.com/lgarrison__;!!DSb-azq1wVFtOg!TX2ckAQuI6fLBuCYNrkxoGd8SpJH2YpvyOluIIFSX8tDO_X3P1j6kOFkDy9dLmeQb-2Zoh6PR5kmnzDLUY_ZwuGvfwl2oX9s$>
, thanks for testing out. It will be great to know in what cases the
accuracy decreases. If it's on the CPU end, there is a new developing
option called spread_kerformula, currently, if nothing is set, it will
use spread_kerformula=8, and setting spread_kerformula=1 should use the
spread kernel setup in v2.4.1. You might want to try this option to see if
it restores the old accuracy behavior. And it's also worth trying
spread_kerformula=9 and spread_kerformula=3. If you can provide a minimal
example, that will be great too. This will help us decide the default
spread kernel formula in the next stable release. Thanks!
—
Reply to this email directly, view it on GitHub
<https://urldefense.com/v3/__https://github.com/flatironinstitute/finufft/discussions/798*discussioncomment-15775524__;Iw!!DSb-azq1wVFtOg!TX2ckAQuI6fLBuCYNrkxoGd8SpJH2YpvyOluIIFSX8tDO_X3P1j6kOFkDy9dLmeQb-2Zoh6PR5kmnzDLUY_ZwuGvfyDRPqR-$>,
or unsubscribe
<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACNZRSWZXJYZ7IWR3MEHRL34LPA6FAVCNFSM6AAAAACSLGX6FGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKNZXGU2TENA__;!!DSb-azq1wVFtOg!TX2ckAQuI6fLBuCYNrkxoGd8SpJH2YpvyOluIIFSX8tDO_X3P1j6kOFkDy9dLmeQb-2Zoh6PR5kmnzDLUY_ZwuGvf49H-kup$>
.
You are receiving this because you were mentioned.Message ID:
***@***.***
com>
--
*-------------------------------------------------------------------~^`^~._.~'
|\ Alex Barnett Center for Computational Mathematics, Flatiron Institute
| \ http://users.flatironinstitute.org/~ahb 646-876-5942
|
|
I extracted an example from the latency-sensitive case. I'm doing a 1D type 1 transform with 3000 NU points and 1000 modes. With 2.4.1, I get: With 2.5.0rc1, I get: In my application, I'm executing a few hundred of these, so the 250 us difference becomes a few dozen ms. I'm running on a Rome node. It's pretty easy to verify that the plan creation is main driver of the difference. While I'd like to reuse plans, the transforms in this particular application have different lengths (nifty-ls periodograms of heterogeneous time series). I would also expect this to matter for jax-finufft, since we're not exposing the plan interface directly there. (Within a single jax-finufft call, transforms are aggressively grouped and stacked into as few plans as possible, but if the user is making many calls they'll end up creating many plans.)
|
|
Hi @lgarrison, @ahbarnett, I run some tests and it seems that in master the new hotspots that cause the regression are prolate0_eval, FFTW planning takes longer for some reason and memory allocation overhead is higher in master. Changing to spread_kerformula makes the code a bit faster but not by much. Since the code is deleting the plan caching things to be reused in this case seems not possible. I don't want a global state to store these. My suggestion would be to make finufft_makeplan callable multiple times and then based on what changed re-use previous computation as much as possible. This requires some thought and work tough so it might not happen soon. What do you all think? Cheers, |
|
FYI I committed this benchmark to master so that we can run it in case we make changes. |
|
That sounds fine to me, I don't think this needs to be a blocker for the 2.5 release. Thanks! |
|
I tried just make plan and destroy plan, the following is v2.4.1 vs master. On the same es kernel, master has about 40us more than v2.4.1. In master, the result is: |
|
comment out |
|
So we are still a bit slower due to the memalloc. I wonder how much faster ducc would be in these cases. |
|
From what I see 1,2, 5,6 are vectorized by the compiler the others are not. |
Yes, if poly coeffs computation is the most expensive part, we could reuse coeffs if ns, upsampfact are the same for different make plan calls. I guess this is a future release task not for v2.5.0. |
|
Thanks for collecting this data and investigating. I agree that 300us plan
time is acceptable, and any user who wants high throughput for small
problems should plan once, then use combos of setpts and execute to do
things.
There is test/testmanysmallprobs.cpp that I wrote to compare these two
cases.
In the old days before we had a plan interface, for FFTW the plan look-up
time was 100us *per thread*, giving multi-milliseconds per call, which we
deemed unacceptable.
It is conceivable that in v2.6+ we might want a "plan for new number of
modes without changing the kernel" which would avoid the 300us in the
FFT=DUCC case. But I suggest leaving for now.
…On Thu, Feb 19, 2026 at 4:04 PM Lehman Garrison ***@***.***> wrote:
That sounds fine to me, I don't think this needs to be a blocker for the
2.5 release. Thanks!
—
Reply to this email directly, view it on GitHub
<https://urldefense.com/v3/__https://github.com/flatironinstitute/finufft/discussions/798*discussioncomment-15864477__;Iw!!DSb-azq1wVFtOg!Qk3ZJjI0muwXM-VAQBql6pRGWSnKtV2Z4NxnqU3Z2ZXmzjA9KcjE9Ru2mYldEgXJE-1LKQ6AJpWOHRqPbQGkkWgq35sGqRrx$>,
or unsubscribe
<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACNZRSRWJDHCCFEN325N72T4MYQO7AVCNFSM6AAAAACSLGX6FGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKOBWGQ2DONY__;!!DSb-azq1wVFtOg!Qk3ZJjI0muwXM-VAQBql6pRGWSnKtV2Z4NxnqU3Z2ZXmzjA9KcjE9Ru2mYldEgXJE-1LKQ6AJpWOHRqPbQGkkWgq38LRRno3$>
.
You are receiving this because you were mentioned.Message ID:
***@***.***
com>
--
*-------------------------------------------------------------------~^`^~._.~'
|\ Alex Barnett Center for Computational Mathematics, Flatiron Institute
| \ http://users.flatironinstitute.org/~ahb 646-876-5942
|
|
I'm curious: how many evaluations of the PSWF does finufft need for, say, a support of 8 and a polynomial degree of 12? If I remember correctly, the polynomial is obtained by evaluating the PSWF at points in the complex plane and solving a Vandermonde system. Would an approximation using just the function values at Chebyshev points be quicker? That should require 4 (instead of 8 because of kernel symmetry) x 12 function evaluations. |
|
Adding my thoughts on this. Using symmetry and making the code vectorizable would give a good speedup on top of adding caching to it. I am not sure it is worth the effort. |
|
For PSWF, vectorization will be an immense task, I think. |
|
Actually there may even be another factor 2 in the PSWF evaluation, since at the moment all function values are normalized by repeated calls to Lines 1050 to 1054 in 3a96a20 Caching |
Thanks for catching. With all your suggestions, we should be able to reduce the plan for kb/pswf from 260 us to around 150 us level, which is quite similar to ES. It's helpful for small problems with multiple planning calls. |
|
We can also propagate tot to the Legendre eval instead of going for machine precision. That can be done to the bessel functions too. |
|
Marco observed above that FFTW planning has become somewhat slower. That may be because we now create two plans: forward and adjoint. |
|
Not sure if this is of interest, but I have a PSWF implementation that is pretty compact and also seems faster than the one in It's not terribly polished, but you should be able to play around with the /*
The code in this file was converted from specfun.f from scipy, which
is in turn based on the work by Shanjie Zhang and Jianming Jin in
their book "Computation of Special Functions", 1996, John Wiley &
Sons, Inc.
*/
#include <cmath>
#include <vector>
#include <iostream>
double sdp_pswf_cipow(double base, int exp)
{
double result = 1;
// This is never called with negative exponents. Avoids potential recursion.
if (exp==1) return base;
while (exp)
{
if ((exp & 1) != 0) result *= base;
exp >>= 1;
base *= base;
}
return result;
}
static inline void pswf_sdmn(
const int m,
const int n,
const double c,
const double cv,
const int kd,
double* df
)
{
/* ===================================================== */
/* Purpose: Compute the expansion coefficients of the */
/* prolate and oblate spheroidal functions, dk */
/* Input : m --- Mode parameter */
/* n --- Mode parameter */
/* c --- Spheroidal parameter */
/* cv --- Characteristic value */
/* KD --- Function code */
/* KD=1 for prolate; KD=-1 for oblate */
/* Output: DF(k) --- Expansion coefficients dk; */
/* DF(1), DF(2), ... correspond to */
/* d0, d2, ... for even n-m and d1, */
/* d3, ... for odd n-m */
/* ===================================================== */
int nm = int ((n-m) * .5 + c) + 25;
if (c < 1e-10)
{
for (int i=0; i<nm; ++i)
df[i] = 0.;
df[(n-m) / 2] = 1.;
return;
}
double cs = c * c * kd;
int ip = 1;
int k = 0;
if (n-m == (n-m)/2<<1)
ip = 0;
std::vector<double> a(nm+2), d_(nm+2), g(nm+2);
for (int i=0; i<nm+2; ++i)
{
if (ip==0)
k = i << 1;
if (ip==1)
k = (i<<1) + 1;
double dk0 = double (m + k);
double dk1 = double (m + k + 1);
double dk2 = double ((m + k) << 1);
double d2k = double ((m << 1) + k);
a[i] = (d2k+2.) * (d2k+1.) / ((dk2+3.) * (dk2+5.)) * cs;
d_[i] = dk0 * dk1 + (dk0*2.*dk1 - m*2.*m - 1.) / ((dk2-1.) * (dk2+3.)) * cs;
g[i] = k * (k-1.) / ((dk2-3.) * (dk2-1.)) * cs;
}
double fs = 1.;
double f1 = 0.;
double f0 = 1e-100;
int kb = 0;
df[nm] = 0.;
double fl = 0.;
for (k=nm; k>=1; --k)
{
double f = -((d_[k] - cv) * f0 + a[k] * f1) / g[k];
if (std::abs(f) > std::abs(df[k]))
{
df[k-1] = f;
f1 = f0;
f0 = f;
if (std::abs(f) > 1e100)
{
for (int k1=k; k1<=nm; ++k1)
df[k1-1] *= 1e-100;
f1 *= 1e-100;
f0 *= 1e-100;
}
}
else
{
kb = k;
fl = df[k];
f1 = 1e-100;
double f2 = -(d_[0] - cv) / a[0] * f1;
df[0] = f1;
if (kb==1)
fs = f2;
else if (kb==2)
{
df[1] = f2;
fs = -((d_[1] - cv) * f2 + g[1] * f1) / a[1];
}
else
{
df[1] = f2;
for (int j=3; j<=kb+1; ++j)
{
double f = -((d_[j-2] - cv) * f2 + g[j-2] * f1) / a[j - 2];
if (j <= kb)
df[j-1] = f;
if (std::abs(f) > 1e100)
{
for (int k1=0; k1<j; ++k1)
df[k1] *= 1e-100;
f *= 1e-100;
f2 *= 1e-100;
}
f1 = f2;
f2 = f;
}
fs = f;
}
break;
}
}
double r1 = 1.;
for (int j=m+ip+1; j<=(m+ip)<<1; ++j)
r1 *= j;
double su1 = df[0] * r1;
for (k=2; k<=kb; ++k)
{
r1 = -r1 * (k + m + ip - 1.5) / (k-1.);
su1 += r1 * df[k-1];
}
double su2 = 0.;
double sw = 0.;
for (k = kb + 1; k <= nm; ++k)
{
if (k!=1)
r1 = -r1 * (k + m + ip - 1.5) / (k-1.);
su2 += r1 * df[k-1];
if (std::abs(sw-su2) < std::abs(su2) * 1e-14)
break;
sw = su2;
}
double r3 = 1.;
for (int j=1; j <= (m+n+ip) / 2; ++j)
r3 *= j + (n+m+ip) * .5;
double r4 = 1.;
for (int j=1; j <= (n-m-ip) / 2; ++j)
r4 *= -4.*j;
double s0 = r3 / (fl * (su1 / fs) + su2) / r4;
for (int k=0; k<kb; ++k)
df[k] *= fl/fs*s0;
for (k=kb; k<nm; ++k)
df[k] *= s0;
}
static inline void pswf_sckb(
int m,
int n,
double c,
const double* df,
double* ck
)
{
/* ====================================================== */
/* Purpose: Compute the expansion coefficients of the */
/* prolate and oblate spheroidal functions */
/* Input : m --- Mode parameter */
/* n --- Mode parameter */
/* c --- Spheroidal parameter */
/* DF(k) --- Expansion coefficients dk */
/* Output: CK(k) --- Expansion coefficients ck; */
/* CK(1), CK(2), ... correspond to */
/* c0, c2, ... */
/* ====================================================== */
// df[nm+1], ck[nm]
if (c <= 1e-10) c = 1e-10;
int nm = int ((n-m) * .5 + c) + 25;
int ip = 1;
if (n-m == (n-m)/2 << 1)
ip = 0;
double reg = 1.;
if (m + nm > 80)
reg = 1e-200;
double fac = -sdp_pswf_cipow(.5, m);
double sw = 0.;
for (int k = 0; k <= nm-1; ++k)
{
fac = -fac;
int i1 = (k<<1) + ip + 1;
double r = reg;
for (int i=i1; i<=i1 + (m<<1) - 1; ++i)
r *= i;
int i2 = k + m + ip;
for (int i=i2; i<=i2+k-1; ++i)
r *= i + .5;
double sum = r * df[k];
for (int i=k+1; i<=nm; ++i)
{
double d1 = i * 2. + ip;
double d2 = m * 2. + d1;
double d3 = i + m + ip - .5;
r = r * d2 * (d2-1.) * i * (d3+k) / (d1 * (d1-1.) * (i-k) * d3);
sum += r * df[i];
if (std::abs(sw-sum) < std::abs(sum) * 1e-14)
break;
sw = sum;
}
double r1 = reg;
for (int i=2; i<= m+k; ++i)
r1 *= i;
ck[k] = fac*sum/r1;
}
}
static inline void pswf_segv(
const int m,
const int n,
const double c,
const int kd,
double* cv,
double* eg
)
{
std::vector<double> a(300), b(100), d_(300), e(300), f(300), g(300), h_(100);
std::vector<double> cv0(100);
/* ========================================================= */
/* Purpose: Compute the characteristic values of spheroidal */
/* wave functions */
/* Input : m --- Mode parameter */
/* n --- Mode parameter */
/* c --- Spheroidal parameter */
/* KD --- Function code */
/* KD=1 for Prolate; KD=-1 for Oblate */
/* Output: CV --- Characteristic value for given m, n and c */
/* EG(L) --- Characteristic value for mode m and n' */
/* ( L = n' - m + 1 ) */
/* ========================================================= */
if (c < 1e-10)
{
for (int i=0; i<n-m+1; ++i)
eg[i] = (i+m+1) * (i+m);
*cv = eg[n-m];
return;
}
int icm = (n-m+2)/2;
int nm = int((n-m) * .5 + c) + 10;
double cs = c*c*kd;
int k = 0;
for (int l=0; l<=1; ++l)
{
for (int i=1; i<=nm; ++i)
{
if (l==0)
k = (i-1) << 1;
if (l==1)
k = (i<<1) - 1;
double dk0 = double (m+k);
double dk1 = double (m+k+1);
double dk2 = double ((m+k) << 1);
double d2k = double ((m<<1) + k);
a[i-1] = (d2k+2.) * (d2k+1.) / ((dk2+3.) * (dk2+5.)) * cs;
d_[i-1] = dk0*dk1 + (dk0*2.*dk1 - m*2.*m - 1.) / ((dk2-1.) * (dk2+3.)) * cs;
g[i-1] = k * (k-1.) / ((dk2-3.) * (dk2-1.)) * cs;
}
for (k=1; k<nm; ++k)
{
e[k] = sqrt(a[k-1] * g[k]);
f[k] = e[k] * e[k];
}
f[0] = e[0] = 0.;
double xa = d_[nm-1] + std::abs(e[nm-1]);
double xb = d_[nm-1] - std::abs(e[nm-1]);
int nm1 = nm-1;
for (int i=0; i<nm1; ++i)
{
double t = std::abs(e[i]) + std::abs(e[i+1]);
xa = std::max(xa, d_[i]+t);
xb = std::min(xb, d_[i]-t);
}
for (int i=0; i<icm; ++i)
{
b[i] = xa;
h_[i] = xb;
}
for (k=1; k<=icm; ++k)
{
for (int k1=k; k1<=icm; ++k1)
if (b[k1-1] < b[k-1])
{
b[k-1] = b[k1-1];
break;
}
if (k!=1)
if (h_[k-1] < h_[k-2])
h_[k-1] = h_[k-2];
double x1;
while(true)
{
x1 = (b[k-1] + h_[k-1]) / 2.;
cv0[k-1] = x1;
if (std::abs((b[k-1] - h_[k-1]) / x1) < 1e-14)
break;
int j = 0;
double s = 1.;
for (int i=1; i<=nm; ++i)
{
if (s==0.)
s += 1e-30;
s = d_[i-1] - f[i-1]/s - x1;
if (s<0.)
++j;
}
if (j<k)
h_[k-1] = x1;
else
{
b[k-1] = x1;
if (j>=icm)
b[icm-1] = x1;
else
{
h_[j] = std::max(h_[j], x1);
b[j-1] = std::min(b[j-1], x1);
}
}
}
cv0[k-1] = x1;
if (l==0)
eg[k*2-2] = cv0[k-1];
if (l==1)
eg[k*2-1] = cv0[k-1];
}
}
*cv = eg[n-m];
}
/**
* @brief Evaluate PSWF at a specific point.
*
* Compute the prolate and oblate spheroidal angular functions of the first
* kind and their derivatives.
*
* This function has been heavily specialised, as it it basically the
* inner loop of PSWF generation.
*
* @param m Mode parameter, m = 0, 1, 2, ...
* @param n Mode parameter, n = m, m + 1, ...
* @param c Spheroidal parameter.
* @param ck Expansion coefficients; CK(1), CK(2) ... correspond to c0, c2 ...
* @param x Argument of angular function, |x| < 1.0
* @return Angular function of the first kind.
*/
double sdp_pswf_aswfa(int m, int n, double c, const double* ck, double x)
{
const int nm = (int) ((n - m) / 2 + c) + 40;
const int nm2 = nm / 2 - 2;
const double x1 = 1.0 - x * x;
const double a0 = (m == 0 && x1 == 0.0) ? 1.0 : pow(x1, m * 0.5);
double su1 = ck[0];
double x1p = x1;
for (int k=1; k<=nm2; ++k, x1p*=x1)
{
const double r = ck[k] * x1p;
su1 += r;
const double t = r/su1;
// if (k >= 10 && abs(t) < 1e-14) break;
}
return ((n - m) % 2 == 0) ? (a0 * su1) : (a0 * x * su1);
}
class pswf_ducc
{
private:
double c=0;
int m=0, n=0, kd=1;
vector<double> coeff;
vector<double> df;
public:
pswf_ducc(double c_) : c(c_), coeff(2000), df(2000)
{
double cv = 0.0, eg[200] = {0.0, 0.0};
pswf_segv(m, n, c, kd, &cv, eg);
pswf_sdmn(m, n, c, cv, kd, df.data());
pswf_sckb(m, n, c, df.data(), coeff.data());
}
vector<double> eval(const vector<double> &x) const
{
vector <double> res;
res.reserve(x.size());
double v0 = sdp_pswf_aswfa(m, n, c, coeff.data(), 0);
for (auto v:x)
res.push_back(sdp_pswf_aswfa(m, n, c, coeff.data(), v));
return res;
}
}; |
|
Ah, you needn't have worried about licensing there, as far as I understand! Please note however that the version in the |
|
Great,if license is not a problem, we could add back and try out in next release. And thanks for the c value reminder. It's still in the commit history c7a56a0 I didn't cache the work array coeff that time, I only played around with it. |
|
Nice, so let's give this a try in the next release! This also gives me some time to brush up the code a bit more. |
|
After looking closely at my code, I think it's probably safer to stick to your current version (and maybe accelerate that a bit). My version starts to do really dodgy things beyond |
|
Yeah, c goes to ~42. |
|
Hello all, a couple of comments. Yes, I don't think this is worth worrying
about.
The timings were us not ns, ie 150us-250us. I don't think this is a major
problem.
Possibly throwing in an openmp pragma into the fitting loop would be an
easy fix.
Normalizing does turn out to be important (consider c=42 in 3D single-prec:
the kernel gets produced by dimension, so would numerically overflow.
I have dealt with that already. Hence normalization.
There are plenty of more important things to do to FINUFFT, like: GPU catch
up to CPU, decent upsampfac heuristics for CPU, updating benchmarking page
(CPU and GPU), etc.
Thanks!
…On Tue, Feb 24, 2026 at 10:14 AM Marco Barbone ***@***.***> wrote:
Yeah, c goes to ~42.
—
Reply to this email directly, view it on GitHub
<https://urldefense.com/v3/__https://github.com/flatironinstitute/finufft/discussions/798*discussioncomment-15912269__;Iw!!DSb-azq1wVFtOg!QT2vJ_aX8mXSw0A_uFuNN_GGRL6DLjn6GaWpUrpR5a_TgZraknoOvuv71jligOoTMQq6dko2AMIWxhCpSRW8TZo96B9SI08m$>,
or unsubscribe
<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACNZRSVD2Y3BGH2S6TXJQDL4NRTG3AVCNFSM6AAAAACSLGX6FGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKOJRGIZDMOI__;!!DSb-azq1wVFtOg!QT2vJ_aX8mXSw0A_uFuNN_GGRL6DLjn6GaWpUrpR5a_TgZraknoOvuv71jligOoTMQq6dko2AMIWxhCpSRW8TZo96INo0xpM$>
.
You are receiving this because you were mentioned.Message ID:
***@***.***
com>
--
*-------------------------------------------------------------------~^`^~._.~'
|\ Alex Barnett Center for Computational Mathematics, Flatiron Institute
| \ http://users.flatironinstitute.org/~ahb 646-876-5942
|
|
Thanks. Yes, I agree that we should focus on more important things. Considering the CPU and GPU code will/is having major refactoring, we probably need some performance benchmark regression tests to run after every major refactoring. I feel Alex tol sweep is good enough to catch accuracy changes, but the performance part seems we are not running much before releases. Minor comment on normalization. The current pswf doesn't seem to need normalization. |
|
ok, so you're proposing to drop the division by pswf(0) in the kernel
definition, fo planning speed? Just checking.
I'm happy with that if it really is O(1) normalized.
…On Tue, Feb 24, 2026 at 12:07 PM Libin Lu ***@***.***> wrote:
Thanks. Yes, I agree that we should focus on more important things.
Considering the CPU and GPU code will/is having major refactoring, we
probably need some performance benchmark regression tests to run after
every major refactoring. I feel Alex tol sweep is good enough to catch
accuracy changes, but the performance part seems we are not running much
before releases.
Minor comment on normalization. The current pswf doesn't seem to need
normalization.
Vladimir's pswf(c=43.0, 0) returns about 1.91917 before normalization.
Other kernels do need normalization.
Tested with
https://github.com/lu1and10/finufft/actions/runs/22360596864/job/64713043311
<https://urldefense.com/v3/__https://github.com/lu1and10/finufft/actions/runs/22360596864/job/64713043311__;!!DSb-azq1wVFtOg!WHT-mFwADHVDaek-1kg99YhUz3lGx6akyv1DSak-mtzeV7sD-Ba-tp3A5bhd5Q0YkXFdVO-xYmrB3-h64JaXLeFxfttI4sfy$>;
seems to pass all single- and double-sweep tests, unless we missed some
test cases.
—
Reply to this email directly, view it on GitHub
<https://urldefense.com/v3/__https://github.com/flatironinstitute/finufft/discussions/798*discussioncomment-15913393__;Iw!!DSb-azq1wVFtOg!WHT-mFwADHVDaek-1kg99YhUz3lGx6akyv1DSak-mtzeV7sD-Ba-tp3A5bhd5Q0YkXFdVO-xYmrB3-h64JaXLeFxfk6DQkth$>,
or unsubscribe
<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACNZRSWECETO5JG5RG56OPD4NSAMLAVCNFSM6AAAAACSLGX6FGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTKOJRGMZTSMY__;!!DSb-azq1wVFtOg!WHT-mFwADHVDaek-1kg99YhUz3lGx6akyv1DSak-mtzeV7sD-Ba-tp3A5bhd5Q0YkXFdVO-xYmrB3-h64JaXLeFxftBMADUh$>
.
You are receiving this because you were mentioned.Message ID:
***@***.***
com>
--
*-------------------------------------------------------------------~^`^~._.~'
|\ Alex Barnett Center for Computational Mathematics, Flatiron Institute
| \ http://users.flatironinstitute.org/~ahb 646-876-5942
|
|
For spherical harmonic transforms it is definitely a huge deal - it allowed performance improvements of roughly 30% for the standard case. |
|
Nevermind ... I was probably overthinking this way too much. [Edit:] ... aargh, perhaps not ... looking at the actual coefficients of the Legendre polynomials, evaluation via Horner looks like nightmarish cancellation errors :-( |
|
Yes, also Horner tends to blow up after degree 40. I would not go above 32 to be safe. Switching Legendre to Chebyshev is not worth either as it is just a different basis for the same stuff. |
|
I just tried it ... works like a charm until about |
|
I used the Ishioka paper to go to a recurrence that only involves the even polynomial orders. Now I manage roughly 3e8 PSWF evaluations per core per second - should be good enough for the moment :-) |
|
Hi Kacper,
Thanks for the contribution. In my experience the current way we fit the
coefficients is roughly as accurate as Remez. I expect Remez to be 1/2
digit more accurate but not much more than that. Though, feel free to give
it a try! If significantly more accurate and If you can implement the whole
remez for runtime fitting here it would be great!
…On Wed, Apr 29, 2026 at 2:43 PM Kacper Rutkowski ***@***.***> wrote:
I used the Ishioka paper to go to a recurrence that only involves the even
polynomial orders. Now I manage roughly 3e8 PSWF evaluations per core per
second - should be good enough for the moment :-)
I am new in this topi, but I have yet another stupid idea here - what if
we decided to use Remez algorithm
<https://urldefense.com/v3/__https://mathworld.wolfram.com/RemezAlgorithm.html__;!!DSb-azq1wVFtOg!W8xCSwdA-WuNnH_gGWdZAUrc46c5NZEXJS-gMG4hU1ApFUt7cQYFAmfQeZTjbwYq86gkcU57LJ_iciKj6bZbg2COSm-48oVm$>
to define the 'optimal' polynomial for absolute error value which should be
more or less fine for the double precision using even or odd powers only.
That should also be significantly lower order.
Then the resulting degree should be significantly lower than the one
required by a more elegant, arbitrary precision algorithm (my preliminary
tests have shown sth around degree 14 for double precision and degree 8 for
single, although I am not sure if my code isn't buggy). Then we could have
attempted using Horner on these low degree polynomials, possibly with double-double
arithmetic
<https://urldefense.com/v3/__https://csclub.uwaterloo.ca/*pbarfuss/dekker1971.pdf__;fg!!DSb-azq1wVFtOg!W8xCSwdA-WuNnH_gGWdZAUrc46c5NZEXJS-gMG4hU1ApFUt7cQYFAmfQeZTjbwYq86gkcU57LJ_iciKj6bZbg2COSgCMCcbX$>
if we want the results to be near-exact. Would an approach like that make
sense?
I guess it should be much faster to evaluate, than the arbitrary precision
methods, although it would also be a really mathematically "ugly" approach
as we would have to just put look-up tables with coefficients (the Remez
algorithm gets pretty computationally expensive to estimate the
coefficients, even for low order polynomials).
I'll try to implement that and provide some code in the near future,
although I believe that it might be much faster algorithm in this case.
—
Reply to this email directly, view it on GitHub
<https://urldefense.com/v3/__https://github.com/flatironinstitute/finufft/discussions/798*discussioncomment-16759710__;Iw!!DSb-azq1wVFtOg!W8xCSwdA-WuNnH_gGWdZAUrc46c5NZEXJS-gMG4hU1ApFUt7cQYFAmfQeZTjbwYq86gkcU57LJ_iciKj6bZbg2COSqNoLTCo$>,
or unsubscribe
<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACGKNQPRSUP4BQA7BTYG4LL4YJEMHAVCNFSM6AAAAACSLGX6FGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTMNZVHE3TCMA__;!!DSb-azq1wVFtOg!W8xCSwdA-WuNnH_gGWdZAUrc46c5NZEXJS-gMG4hU1ApFUt7cQYFAmfQeZTjbwYq86gkcU57LJ_iciKj6bZbg2COSuH5B2ZB$>
.
Triage notifications on the go with GitHub Mobile for iOS
<https://urldefense.com/v3/__https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675__;!!DSb-azq1wVFtOg!W8xCSwdA-WuNnH_gGWdZAUrc46c5NZEXJS-gMG4hU1ApFUt7cQYFAmfQeZTjbwYq86gkcU57LJ_iciKj6bZbg2COSmpuT1mI$>
or Android
<https://urldefense.com/v3/__https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign*3Dnotification-email*26utm_medium*3Demail*26utm_source*3Dgithub__;JSUlJSU!!DSb-azq1wVFtOg!W8xCSwdA-WuNnH_gGWdZAUrc46c5NZEXJS-gMG4hU1ApFUt7cQYFAmfQeZTjbwYq86gkcU57LJ_iciKj6bZbg2COSgihIXgs$>.
You are receiving this because you commented.Message ID:
***@***.***
com>
|
|
Thanks, I've already implemented an exetremely simplified version it, and I guess you are right. Either Remez multiplied by Gaussian or anything like that would give significantly lower order (I've ended with 4-term polynomial for single precision, but the algorithm was really problematic - even for single Remez did get lost in some cases and it took literally minutes to get correct results with Python's mpmath, so anything like compile time (pre)computation is not really feasible - also Remez gets lost quite often and I am not really able to find a starting point at which it wouldn't. So I guess in theory it might have been 2-3x faster, but also it would require tens of many different look-up tables and weighted Remez (required for correcting the Gaussian 'initial guess') gets lost quite easily. So I guess it was really just a stupid idea, which might be fast, but is exetremely difficult to implement correctly (that's why I deleted my original coment a moment after posting and sorry for that as well). Also thank you very much for yours answer, I agree that it probably is not a good idea, although I might experiment with modified Remez a little bit more in the future. I doubt adding it to the library would be practical in any way though. |
|
I would also be happy to see the results from the Remez algorithm, but I'm pretty sure that (for the values of |
|
Sorry for not being precise once again. What I did was generally multiplying exp(-cx^2/2) initial approximation by a polynomial resulting from a weighted Remez algorithm to approximate the function, not just plugging-in the polynomial. That significantly reduces the number of terms required (well… The exponential is also largely a polynomial under the hood, but still pretty low order), which ended up working pretty well as long, as it did work. The Remez algorithm gets lost more often than not though + the weighted implementation gets lots more often, than it works for degrees higher, than 6. Generally I think using exp(-cx^2/2) * poly approximation would be a good idea, although I don't really have an idea how such polynomial coefficients could be approximated quickly and efficiently… Because as I already mentioned the weighted Remez seems to struggle quite a lot. EDIT: And also I am an idiot. I've set substitution z=x^2 and I've divided number of terms by two instead of multiplying it because I can't read. I'll probably send my code here if I manage to make it work even in most cases, although currently it requires 9 terms for single precision and 20 for double, so the results for Gaussian product are not impressive at all. Thank you very much for putting me right on that issue. |
|
So… Here are results of my stupid idea after I kinda made it work. Pros: Horner's scheme polynomial evaluation works really well bc Gaussian base dampens the high-degree oscillations and evaluation seems to be pretty fast (on my laptop with Ryzen 7945hx it reaches 3e9 evaluations per second for single precision for eps=1e-7 and 8e8 evaluations for double with eps=1e-16, so I guess about 2x more, than the solution based on Ishioka paper?). Cons: Polynomial estimation is exetremely slow, generating a single look-up table can take many minutes and the modified Remez algorithm often fails to converge (original one didn't want to converge at all). Below is the code I've used, along with the polynomial generator (using eigen3, mpfr and mpreal libraries), along with a minimalistic standalone NuFFT implementation using generated polynomials (and benchmarking the PSWF generation). Also I've gotten quite a bit better accuracy ceiling when compared to Exponential of Semicircle kernel using standard library, which has been a little bit surprising to me here. Below are the results for multiple NuFFT-1s executed with 200 input points and 256 target frequencies (average errors and execution times for both kernels) - all on a single thread of course (Ryzen 7945hx once again); Double precision:
Single precision:
Well… The numerical benchmark results look pretty promising to me, although I sadly don't really see any feasible way of including a runtime coefficient generation in FiNUFFT given how computationally expensive it is (even if MPFR was replaced with simple and 'fast' Double-Double arithmetic for that step). And I don't think littering the source code with tens of LUTs for different values of epsilon is justified, although that might be potentially worth consideration. Thank you very much for everything, I guess that approach somehow works, although it's much less feasible to implement, than I originally thought (mostly due to Remez algorithm's exetreme instability when estimating psfw of 0th order. |
|
Still this is interesting thanks for coding this up! |
|
Hello all,
I haven't been following this chain much, but I would say:
* Piecewise Horner is the only eval method for spread/interp fast enough,
especially if we want robust across compilers/arch's. We can't have any
exp() in there (so, ignore my claims in the original 2019 FINUFFT paper),
and PSWF is the best kernel (established Jan 2026). We don't need to
consider any kernels other than PSWF.
* the only method fast enough for on-the-fly is linear poly fitting. The
plan stage has to be sub-100us, otherwise small transforms start to suffer.
* Remez is not going to be much better than poly-fit at Cheby nodes,
because (relative to each fitting interval) the domain of analyticity is
large - the convergence rate is high. This is covered in Trefethen's ATAP
book - please read if you want to understand this stuff. Often the error
ratio from Cheby-interp to optimal-in-infty-norm is bounded by a factor of
2. Also, if you want to avoid the python library difficulties, just use
MATLAB Chebfun which has Remez etc.
There are plenty of other things to improve in FINUFFT (bench, binsort,
GPU, ...). I think the kernel function approx & eval stuff is done already.
Ok, looking forward to getting the summary at some point. Thanks & best,
Alex
…On Thu, Apr 30, 2026 at 2:06 PM Marco Barbone ***@***.***> wrote:
Still this is interesting thanks for coding this up!
—
Reply to this email directly, view it on GitHub
<https://urldefense.com/v3/__https://github.com/flatironinstitute/finufft/discussions/798*discussioncomment-16772787__;Iw!!DSb-azq1wVFtOg!RV7g4GQp2dWp8VAOePGO3urEbNk6jT3xVv6tpVqX98r9-cTnxzx--NMIMwQrppZ-dscj-3hVbvl0u7gK91DjFQj_Q48hKAdz$>,
or unsubscribe
<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACNZRSQK4RAIXD7FDQAKFOL4YOI2FAVCNFSM6AAAAACSLGX6FGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTMNZXGI3TQNY__;!!DSb-azq1wVFtOg!RV7g4GQp2dWp8VAOePGO3urEbNk6jT3xVv6tpVqX98r9-cTnxzx--NMIMwQrppZ-dscj-3hVbvl0u7gK91DjFQj_Q1Fg6zgw$>
.
Triage notifications on the go with GitHub Mobile for iOS
<https://urldefense.com/v3/__https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675__;!!DSb-azq1wVFtOg!RV7g4GQp2dWp8VAOePGO3urEbNk6jT3xVv6tpVqX98r9-cTnxzx--NMIMwQrppZ-dscj-3hVbvl0u7gK91DjFQj_Q2F6Nmfl$>
or Android
<https://urldefense.com/v3/__https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign*3Dnotification-email*26utm_medium*3Demail*26utm_source*3Dgithub__;JSUlJSU!!DSb-azq1wVFtOg!RV7g4GQp2dWp8VAOePGO3urEbNk6jT3xVv6tpVqX98r9-cTnxzx--NMIMwQrppZ-dscj-3hVbvl0u7gK91DjFQj_Q5NZi3mx$>.
You are receiving this because you were mentioned.Message ID:
***@***.***
com>
--
*-------------------------------------------------------------------~^`^~._.~'
|\ Alex Barnett Center for Computational Mathematics, Flatiron Institute
| \ http://users.flatironinstitute.org/~ahb 646-876-5942
|
Uh oh!
There was an error while loading. Please reload this page.
Uh oh!
There was an error while loading. Please reload this page.
Dear users, [updated 1/24/26 for PSWF kernel; please ignore outdated comment/discussion below]
We will shift to a new spreading kernel function and kernel shape parameter choice for v2.5.0. Comparison was made much easier by Marco's recent introduction of on-the-fly piecewise polynomial approximation of arbitrary kernels (at arbitrary upsampling factors sigma). We compared many kernels and shape-parameter function choices, using our empirical relative l2-error measuring code
wsweepkerrcomp.m. The winner (perhaps unsurprisingly) is the prolate spheroidal wavefunction (PSWF) of order zero. A choice of shape parameter beta = beta_cutoff - 0.05 was found to equal (within the noise) all other optimizations, including more elaborate least-squares fits. (The shape parameter for the PSWF is often denoted by c and called "frequency".) Here beta_cutoff := pi.w.(1 - 2/sigma), where sigma is the upsampling factor, and w the kernel width in fine grid points. The subscript "cutoff" derives from matching the theoretical cutoff of the Fourier transform of the PSWF (namely beta) to the first aliased copy of the desired range of output modes; see (4.5) in the original [FIN] paper, setting gamma=1.The 2.5.0 release candidate will be tagged in the next couple of weeks. Here's a preview of the improvements. In this post we discuss the kernel, its shape choice, and the new logic for setting w to match the user's tolerance request.
One way to visualize the improvement is to compare the errors at each w (kernel width in fine gridpoints, a.k.a. ns or nspread). Here's this in double-precision 1D, for all 3 types. Notice the error improvement from blue (kf=1, the current ES kernel) to green (kf=8, the proposed new PSWF):
The above is for standard upsampling. The error improvement is typically a factor 2-4 (ie, 0.3 to 0.6 digits), with even more improvement at w=4. Now below we pick a low upsampling choice. There is again between 0.2 and 0.6 digits improvement, depending on overall tolerance:
This translates to more accuracy for the same (spread/interp) effort. We have explored many other kernel and shape-parameter choices (PSWF with other beta choices; Kaiser-Bessel with beta chosen using the formula of Beatty et al 2005; ES with Beatty params; continuous KB; cosh-type; and continuous cosh-type).
These plots are generated by the new code
matlab/test/wsweepkerrcomp.m, in the kerformdef branch that will shortly merge to master, then 2.5.0-rc1.This improvement in error performance at each w enables a better choice of w for each user-requested tolerance (eps), that is now a unified formula for all sigma. Here are my current proposals (again blue "+" are the current ES, while green dots are the proposed PSWF). Notice that the PSWF "steps" are quite a lot lower than the old ES ones, and the worst-case errors are much closer to the target tolerance (the red line):
Here, only the domain of tolerances that do not report a "tol too small" warning (code 1) are shown. This plot is generated by the same MATLAB code as above. We see here only a slight improvement in how far above this line the worst-case, mostly at small tol. For sigma=1.25 we see much more improvement here (mostly because the ES legacy rule for w at sigma=2.0 was different from that at other sigma's):
The choice of ns (=w) given tol (=eps) has a new form that matches the worst-case rel l2 error to be around 2*tol, uniformizing it across sigmas, types, and dimensions. This (along with the above-mentioned beta formula) is in
src/common/kernel.cpp:This combines an overall tolerance factor (tolfac) with an overall width offset (nsoff) to match the errors more uniformly to tolerance across sigmas, types, and dimensions. The type-3 specific tolfac keeps the type-3 errors more in-line with types 1,2, although it also breaks up the simple stepping (because it does not apply to the "inner type 2" of the type 3). The parameters in the above ns choice have been hand-tuned, rather than heavily optimized. The 2D and 3D error improvements are similar to the above 1D results.
See the next posting for the effect on the testing set-up.
Consequences for users: the w (ns) choices made by FINUFFT may differ from those in v2.4.1, which may affect speed. Most likely at small tolerances this will be in the good direction (smaller w), but it may be for crude tol=0.1 that ns increases from 2 to 3 for, say, 3D type-3. However the resulting errors will be more in line with the requested tol. You can always adjust tol to your liking.
One example of the improvement (lowering and uniformization) of the errors is shown by
matlab/test/fullmathtest.mcontrasting the old kernel (ES legacy):with the new PSWF kernel and parameter choices:
This is for two tolerances only, but notice how much smaller the errors are. In particular the bad errors of up to 0.005 (when requesting tol=0.001) have been reduced to at most 0.0017. Note that in 3D (for tol=1e-3) w=6 is used instead of w=5 for the lower dimensions; this explains the errors around 1.3e-4 for types 1,2.
When requesting 1e-9, all the errors lie within 4e-10 and 9e-10, which is below tol (before we saw errors up to 2.5 times tol).
Thoughts are welcome; also see PR #787.
Best wishes, Alex, Marco & Libin 1/24/26.
All reactions