docs tut inv1d2 basic version done, two matvec versions

docs/tutorial/inv1d2.rst

@@ -16,9 +16,8 @@ Note that there are many other methods to fit smooth functions from
 nonequispaced samples, eg high-order local Lagrange interpolation.
 However, often your model for the function is a global Fourier series,
 in which case, what we describe below is a good starting point.
-We first assume that the samples are noise-free and no regularization is
-needed, then proceed to the noisy regularized case. We also illustrate
-well- and ill-conditioned cases.
+We assume that the samples are noise-free and no regularization is
+needed. We illustrate a well-conditioned then an ill-conditioned case.
 
 Here's code to set up a random complex-valued
 test problem of size 600000 by 300000 (way too large to solve directly):
@@ -30,12 +29,12 @@ test problem of size 600000 by 300000 (way too large to solve directly):
   M = 2*N;                             % overdetermined by a factor 2
   x = 2*pi*((0:M-1)' + 2*rand(M,1))/M; % jittered-from-uniform points on the periodic domain
   ftrue = randn(N,1) + 1i*randn(N,1);  % choose known Fourier coeffs at k inds
-  ftrue = ftrue/sqrt(N);               % make signal f(x) variance=1, Re or Im part
+  ftrue = ftrue/sqrt(N);               % scale signal f(x) to variance=1, Re or Im part
   y = finufft1d2(x,+1,1e-12,ftrue);    % eval noiseless data at high accuracy
 
 
-Case of no noise and no regularization
---------------------------------------
+Well-conditioned example
+------------------------
 
 The linear system to solve is
 
@@ -104,8 +103,10 @@ estimate on a fine grid:
   yg = finufft1d2(xg,+1,1e-12,f);          % eval recovered series there
   fprintf('abs max err: %.3g\n', norm(yg-ytrueg,inf))
 
-This returns ``abs max err: 0.00146`` indicating that the conditioning
-is not 
+This returns ``abs max err: 0.00146``, 3 digits worse than the $\ell_2$
+coefficient error, indicating that there are some locations for the problem
+which are not entirely well-conditioned. Yet, almost everywhere we see excellent matching of the recovered to the true function to 5-6 digits, for instance
+in this zoom in of the first 0.03% of the periodic domain:
 
 .. image:: ../pics/inv1d2err_wellcond.png
    :width: 90%
@@ -129,61 +130,82 @@ unit strengths:
   v = finufft1d1(x, ones(size(x)), -1, tol, 2*N-1);  % Toep vec, inds -(N-1):(N+1)
   vhat = fft([v;0]);                                 % pad to length 2N
 
-We now use a pair of padded FFTs to apply the discrete convolution to
-any vector $f$.
+We now use a pair of padded FFTs in a function (see ``tutorial/applyToep.m`` for the documented self-testing version) which applies the discrete convolution to any vector $f$:
 
-sensible padding
+.. code-block:: matlab
+
+  function Tf = applyToep(f,vhat)
+    N = numel(f);
+    fpadhat = fft(f(:),2*N);               % first zero-pads out to size of vhat
+    Tf = ifft(fpadhat .* vhat(:));         % do periodic convolution
+    Tf = Tf(N:end-1);                      % extract correct chunk of padded output
+  end
 
+.. note::
 
+   Since FFTs are periodic, the minimum length that padded FFTs can be to correctly compute the central $N$ entries of the nonperiodic convolution of a length $N$ vector with a length $2N-1$ vector is $2N-1$. However, for $N=3\times 10^5$, $2N-1=599999$ is prime! Its FFT is several times slower than one of length $2N$. Thus we choose $2N$ as the padded length; a more optimized code might pad to the next 5-smooth even number above $2N-1$, using, eg, `next235even <https://github.com/flatironinstitute/finufft/blob/9a1fae7ab1c2f6b1e51c8907b4d6483d5b55f716/src/utils_precindep.cpp#L15>`_.
 
+The resulting iteration count is identical to that for the NUFFT-based matvec,
+but the CPU time is now 0.65 seconds, ie, 2.5x faster.
+As a reminder, this is because the spreading/interpolation operations in the NUFFTs are avoided (the FFT sizes in the NUFFTs being similar to those in this Toeplitz matvec). The errors and plots are very similar to before.
 
-
-  This reaches ``relres<1e-6`` in 1461 iterations
-(a large count indicating poor conditioning),
-taking about 100 seconds on an 8-core laptop.
-The relative residual for the desired system $A{\bf f}={\bf y}$
-is ``2.7e-05``, indicating that *the linear system was solved
-reasonably accurately*,
-but the relative coefficient error is a much larger
-``2.4e-02``. Their ratio places a lower bound on the condition
-number $\kappa(A)$ of about 900, explaining the large iteration count
-for the normal equations.
-Note that 0.0001% residual error in the normal equations resulted
-in 2.4% coefficient error.
+Ill-conditioned example
+-----------------------
+
+The conditioning of the inverse NUFFT problem is set by the nonuniform (sample) point distribution. To illustrate, we now switch to iid random points:
+
+.. code-block:: matlab
 
-The error in the signal $f(x)$ is in fact very unequally distributed
+  x = 2*pi*rand(M,1);
+                
+We use the Toeplitz FFT matvec (method 2 above), and now find that CG
+reaches the requested ``relres<1e-6`` in a huge 1461 iterations
+(the large count implying poor conditioning), taking 35 seconds. The above error
+diagnosis code lines now print::
+
+  rel l2 resid of Af=y: 2.62e-05
+  rel l2 coeff err: 0.0236
+  abs max err: 2.4
+
+Here the residual shows that *the linear system was solved reasonably accurately*, but that the coefficient error is a much worse.
+This is typical behavior for an ill-conditioned problem.
+Their ratio (both in $\ell_2$ norm) places a lower bound on the condition
+number $\kappa(A)$ of about 900. This explains the large iteration count
+for the normal equations, since their condition number is $\kappa(A^*A) = \kappa(A)^2$.
+
+The error in the signal $f(x)$ turns out to be very unequally distributed
 for this problem: it is correct to 4-5 digits almost everywhere,
 including at almost all the data points,
-but errors are ${\cal O}(1)$ in the very largest gaps
-between the (iid random) sample points. Here is such a gap:
+while errors are ${\cal O}(1)$ only near the very largest gaps
+between the (iid random) sample points. Here is a picture near such a gap:
 
 .. image:: ../pics/inv1d2err.png
    :width: 90%
 
-Notice the large error around 0.9212. However, the problem of
+The gap near $x\approx 0.9212$ has size 0.00009, which is about
+two wavelengths at the Nyquist frequency $N/2$.
+The observed ill-conditioning is a feature of the *problem*, and cannot
+be changed by a different solution method.
+Indeed, it can be shown mathematically that the problem of
 interpolating a band-limited function
 is exponentially ill-conditioned with respect to the length of
-any node-free gap measured in wavelengths. The gap near 0.9212 is
-about 0.00009, ie, two wavelengths at the frequency $N/2$.
-A sampling point distribution without large gaps would improve the conditioning
-and make the reconstruction error in $f$ uniformly closer to the residual
-error.
-
-
-CG-Toep relres 9.97e-07 done in 1465 iters, 35 s
+any node-free gap measured in Nyquist wavelengths. This
+partially explains the ill-conditioning observed above.
 
+..  note::
+
+  The coefficient and $f(x)$ reconstruction error could be reduced in the above demo (without changing the conditioning) by reducing the residual (ie, setting a smaller CG stopping criterion); however, the ``1e-6`` relative residual stopping value used already presumes that the data were at least 6-digit accurate (0.0001% measurement error), which is already a stretch in any practical problem. In short, it is not meaningful to demand residuals much lower than the data measurement error.
 
-The solution and plot is essentially identical to that from the
-NUFFT-pair method.
-	rel l2 resid of Ax=y: 2.63e-05
-	rel l2 coeff err: 0.0236
+Two ways to change the problem to reduce ill-conditioning include:
+1) changing the sampling point distribution to avoid large gaps, and
+2) changing the problem by introducing a regularization term.
 
+The 2nd idea here also fits into the iterative NUFFT or Toeplitz frameworks,
+and we plan to present it in another tutorial shortly.
 
+For the complete code for the above examples and plots, see ``tutorial/inv1d2.m``.
 
 
-
-
-
 Further reading
 ---------------
 
@@ -199,3 +221,5 @@ for CG on the normal equations, in the MRI settings, see:
 * J A Fessler et al,  Toeplitz-Based Iterative Image
   Reconstruction for MRI With Correction for Magnetic Field Inhomogeneity.
   IEEE Trans. Sig. Proc. 53(9) 3393 (2005).
+
+For background on condition number of a problem, see Chapters 12-15 of Trefethen and Bau, *Numerical Linear Algebra* (SIAM, 1997).

tutorial/applyAHA.m

@@ -0,0 +1,4 @@
+function AHAf = applyAHA(f,x,tol)         % use pair of NUFFTs to apply A^* A
+  Af = finufft1d2(x,+1,tol,f);         % apply A
+  AHAf = finufft1d1(x,Af,-1,tol,length(f));    % apply A^*
+end

tutorial/inv1d2.m

@@ -0,0 +1,68 @@
+% demo of inversion of the 1D type 2 NUFFT, ie fitting a Fourier series to
+% function values at scattered points.
+% Barnett 12/6/23
+clear; close all; addpath utils; addpath ../matlab
+
+N=3e5;   % num unknown coeffs (ie twice the max frequency)
+ks = -floor(N/2) + (0:N-1);      % row vec of the frequency indices
+M=2*N;           % number of scattered points on the periodic domain
+wellcond=0;
+toep = 1;
+rng(0);                     % fix seed for reproducibility
+if wellcond, x = 2*pi*((0:M-1)' + 2*rand(M,1))/M;   % jittered pts (will be well conditioned)
+else, x = 2*pi*rand(M,1);             % iid random (will be ill conditioned)
+end
+
+% choose known (complex) coeffs, corresponding to above freqs indices
+ftrue = (randn(N,1) + 1i*randn(N,1))/sqrt(N);
+
+tol = 1e-12;
+y = finufft1d2(x,+1,tol,ftrue);       % data = eval this Fourier series
+
+%y = y + 1e-6*(randn(M,1) + 1i*randn(M,1));   % add noise
+
+if N*M<1e7          % expensive dense direct solve, to check what it's doing
+  A = exp(1i*x(:)*ks);            % outer prod
+  fdirect = A\y;    % ouch!
+  fprintf('rel l2 coeff err of dense solve: %.3g\n', norm(fdirect-ftrue)/norm(ftrue))
+end
+
+rhs = finufft1d1(x,y,-1,tol,N);      % compute A^* y
+
+if ~toep   % iterative solve of normal eqns, each iteration a pair of NUFFTs
+  tic
+  [f,flag,relres,iter] = pcg(@(f) applyAHA(f,x,1e-6), rhs, 1e-6, N);
+  fprintf('CG-NUFFT relres %.3g done in %d iters, %.3g s\n', relres,iter,toc)
+else        % iterative solve of normal eqns, each iteration a padded FFT
+  v = finufft1d1(x, ones(size(x)), -1, tol, 2*N-1);  % Toep vec, inds -(N-1):(N+1)
+  vhat = fft([v;0]);
+  tic
+  [f,flag,relres,iter] = pcg(@(f) applyToep(f,vhat), rhs, 1e-6, N);
+  fprintf('CG-Toep relres %.3g done in %d iters, %.3g s\n', relres,iter,toc)
+end
+
+yrecon = finufft1d2(x,+1,tol,f);
+fprintf('\trel l2 resid of Af=y: %.3g\n', norm(yrecon-y)/norm(y))
+fprintf('\trel l2 coeff err: %.3g\n', norm(f-ftrue)/norm(ftrue))
+ng = 10*N; xg = 2*pi*(0:ng)/ng;       % fine plot grid
+ytrueg = finufft1d2(xg,+1,1e-12,ftrue);  % eval true series on plot grid
+yg = finufft1d2(xg,+1,1e-12,f);      % eval the recon series on plot grid
+fprintf('\tabs max err: %.3g\n', norm(yg-ytrueg,inf))
+
+figure();
+subplot(2,1,1); plot(xg,real(ytrueg),'r-'); hold on;
+plot(xg,real(yg),'b-');
+plot(x,real(y),'k.','markersize',5); legend('true f(x)', 'recon series f(x)', 'data points (x_j,y_j)')
+%axis([0 200/N -3*sqrt(N) 3*sqrt(N)]);   % show a few wiggles
+dx = 2e-3;
+if wellcond, x0=0; else, x0 = 0.920; end            % view domain start
+axis([x0 x0+dx -3 3]);   % show a few wiggles
+subplot(2,1,2); semilogy(xg,abs(ytrueg-yg),'b-'); hold on;
+plot(x,abs(yrecon-y),'k.','markersize',5); legend('error in f(x)', 'error at data points')
+axis tight; v=axis; axis([x0 x0+dx 1e-7 1])
+set(gcf,'paperposition',[0 0 8 6]);
+if wellcond, print -dpng ../docs/pics/inv1d2err_wellcond.png
+else, print -dpng ../docs/pics/inv1d2err.png
+end
+
+% [other methods (CG adj nor eqns, PCG, not as good in my tests)...?]