Certified type-1 Gaussian gridding error and type-2 adjoint bound

If $N=0$, [F17] and the finite sums in [F1] give exact zero type-1 output and an empty type-2 output; [F12] assigns zero to empty norms, so both error assertions hold with zero radii. For the remaining proof assume $N>0$. Let $S_{\mathrm{\infty }}$ be the full spread matrix with entries $(S_{\mathrm{\infty }}{)}_{\mathrm{\ell }j}={\mathrm{\Phi }}_{\tau }({\mathrm{\ell }}_1\mathrm{/}{M}_1-x_{j,1},\dots ,{\mathrm{\ell }}_d\mathrm{/}{M}_d-x_{j,d})$, and let $S$ keep only the entries in each source's product neighborhood. Define the mode-by-grid matrix $Q_{k\mathrm{\ell }}=P^{-1}{e}^{-2\pi i{\sum }_r{k}_r{\mathrm{\ell }}_r\mathrm{/}{M}_r}$ and the diagonal mode matrix $D_{kk}=d_k$, with $D_{kk}^{-1}=1\mathrm{/}{d}_k$ by the positive-multiplier premise. The products below use the finite matrix multiplication of [F19]. The normalized grid sum in [F2], applied to the unit source vector $e_j$, says that every entry of $D^{-1}Q{S}_{\mathrm{\infty }}-A$ has modulus at most $\alpha$. For each column $j$, [F3] bounds the sum of the nonnegative deleted entries by $T_d$; since $\mathrm{\mid }{Q}_{k\mathrm{\ell }}\mathrm{\mid }=P^{-1}$, it follows that $\mathrm{\mid }(Q(S_{\mathrm{\infty }}-S){)}_{kj}\mathrm{\mid }\le T_d\mathrm{/}P$. By [F10], $\mathrm{\mid }{k}_r\mathrm{\mid }\le K_r$, and [F14]–[F15], $d_k\ge d_{\min }>0$, so $\mathrm{\mid }(D^{-1}Q(S_{\mathrm{\infty }}-S){)}_{kj}\mathrm{\mid }\le T_d\mathrm{/}(P{d}_{\min })$. Adding the two matrix differences gives $\mathrm{\mid }(D^{-1}QS-A{)}_{kj}\mathrm{\mid }\le {ϵ}_0$ for every $k,j$. For each row, multiplication by $c$ and the complex modulus laws [F21] give $\mathrm{\mid }((D^{-1}QS-A)c{)}_k\mathrm{\mid }\le {ϵ}_0{\sum }_j\mathrm{\mid }{c}_j\mathrm{\mid }$, hence $\mathrm{\parallel }(D^{-1}QS-A)c{\mathrm{\parallel }}_{\mathrm{\infty }}\le {ϵ}_0\mathrm{\parallel }c{\mathrm{\parallel }}_1$ by [F12]. Finally [F9] gives ${ϵ}_0\le {ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}$ for the returned outward bound.

Let $g=Sc$ be the exact truncated spread and $u$ the dyadic center grid computed by the algorithm. The finite image sums and symmetric remainder intervals enclose every exact kernel by [F4]; isotonicity of the outward DAG operations then gives $\mathrm{\parallel }g-u{\mathrm{\parallel }}_2\le E_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}$, with norms as in [F12]. By [F18], the counting Euclidean norm has the triangle inequality. For coordinate axis $r$, regard its $P\mathrm{/}{M}_r$ length-$M_r$ FFTs as one block-diagonal operator. Squaring and summing the linewise supplier bound [F5] gives a normwise error at most ${\rho }_r$ times the block input norm; the exact block transform is unitary by [F6]. If $x_r$ and $y_r$ are the computed and exact outputs after axes $1,\dots ,r$, the triangle inequality [F18], the local error bound and unitarity give $\mathrm{\parallel }{x}_r{\mathrm{\parallel }}_2\le (1+{\rho }_r)\mathrm{\parallel }{x}_{r-1}{\mathrm{\parallel }}_2$ and $\mathrm{\parallel }{x}_r-y_r{\mathrm{\parallel }}_2\le \mathrm{\parallel }{x}_{r-1}-y_{r-1}{\mathrm{\parallel }}_2+{\rho }_r\mathrm{\parallel }{x}_{r-1}{\mathrm{\parallel }}_2$. Starting with $x_0=y_0=u$ and iterating gives $\mathrm{\parallel }{x}_d-y_d{\mathrm{\parallel }}_2\le U{\sum }_{r=1}^d{\rho }_r{\prod }_{s<r}(1+{\rho }_s)=U({\prod }_r(1+{\rho }_r)-1)=\rho U$. For each extracted coordinate, the FFT discrepancy has modulus at most $\rho U$, since a coordinate modulus is bounded by the Euclidean norm. The exact scaling by $P^{-1\mathrm{/}2}$ and exact division by the mode multiplier therefore make its deconvolved discrepancy at most $\rho U\mathrm{/}(\sqrt{P}{d}_k)\le \rho U\mathrm{/}(\sqrt{P}{d}_{\min }^{})$, using [F15] and $0<d_{\min }^{}\le d_{\min }\le d_k$. By [F4], the computed scaling/deconvolution rectangle encloses the exact scaling/deconvolution of the computed FFT coordinate, so its distance from the returned lower-left center is at most $b_k$. The algorithm forms $R_{\mathrm{F}\mathrm{F}\mathrm{T}}$ outward from the actual input norm bound and the verified call factors; [F7] and the complex triangle inequality [F21] then give distance at most $b_k+R_{\mathrm{F}\mathrm{F}\mathrm{T}}\le r_k$ from that center to the exact deconvolution of the dyadic center grid. This justifies the full containment claimed in [F7], including the previously omitted computed FFT error. To transfer the grid-center enclosure, define $U_P$ as the exact unitary DFT on the full rectangular grid and $R_{\kappa }$ as restriction to the slots $\kappa (k)$. Since $M_r>2K_r$, two modes in $[-K_r,K_r]$ with the same residue differ by a multiple of $M_r$ of magnitude less than $M_r$, hence are equal; coordinatewise, $\kappa$ is injective. The grid transform in step 1.1 is $Q=P^{-1\mathrm{/}2}{R}_{\kappa }{U}_P$. For any finite array $z$, each $\mathrm{\mid }{z}_k{\mathrm{\mid }}^2$ is a nonnegative summand of ${\sum }_{\mathrm{\ell }}\mathrm{\mid }{z}_{\mathrm{\ell }}{\mathrm{\mid }}^2$, so [F12] and [F23] give $\mathrm{\parallel }z{\mathrm{\parallel }}_{\mathrm{\infty }}\le \mathrm{\parallel }z{\mathrm{\parallel }}_2$. Parseval and contraction under coordinate restriction in [F6] then yield $\mathrm{\parallel }Q(u-g){\mathrm{\parallel }}_{\mathrm{\infty }}\le \mathrm{\parallel }Q(u-g){\mathrm{\parallel }}_2\le E_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}\mathrm{/}\sqrt{P}$. By [F10], [F14] and [F15], each $d_k\ge d_{\min }$. Thus the additional deconvolution error is at most $E_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}\mathrm{/}(\sqrt{P}{d}_{\min })$. The algorithm defines $B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}={\max }_k{r}_k+E_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{l}}\mathrm{/}(\sqrt{P}{d}_{\min }^{})$; using $d_{\min }^{}\le d_{\min }$, taking the maximum returned $r_k$ and adding its grid-center term gives the asserted radius.

By [F19], the finite gridding matrix is $\tilde{A}=D^{-1}QS$. Step 1.1 bounds its analytic discrepancy by ${ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}\mathrm{\parallel }c{\mathrm{\parallel }}_1$, and step 2.1 bounds the finite-arithmetic discrepancy from $\tilde{A}c$ by $B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}$. The coordinatewise complex triangle inequality [F21], followed by taking the finite maximum, gives $\mathrm{\parallel }\widehat{F}-Ac{\mathrm{\parallel }}_{\mathrm{\infty }}\le {ϵ}_{\mathrm{a}\mathrm{n}\mathrm{a}}\mathrm{\parallel }c{\mathrm{\parallel }}_1+B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}$, and [F1] identifies $Ac=F$.

The type-2 exact finite matrix computed by the algorithm is ${\tilde{A}}^{\ast }$. By [F19], $(\tilde{A}{)}_{kj}=d_k^{-1}{\sum }_{\mathrm{\ell }}{Q}_{k\mathrm{\ell }}{S}_{\mathrm{\ell }j}$. Since $d_k$ and $S_{\mathrm{\ell }j}$ are real by [F14]–[F15] and [F20], conjugating this entry gives $d_k^{-1}{\sum }_{\mathrm{\ell }}\overline{Q_{k\mathrm{\ell }}}{S}_{\mathrm{\ell }j}$; this is exactly the $(j,k)$ entry of $S^{\ast }{Q}^{\ast }{D}^{-1}$ by the same finite-product formula. In the slot notation of step 2.1, $Q^{\ast }=P^{-1\mathrm{/}2}{U}_P^{\ast }{R}_{\kappa }^{\ast }$, where $R_{\kappa }^{\ast }$ embeds a mode vector into the grid slots $\kappa (k)$ and zero-fills every other slot; this is the inverse-FFT construction in step 4. Thus ${\tilde{A}}^{\ast }=S^{\ast }{Q}^{\ast }{D}^{-1}$. Step 1.1 gives $\mathrm{\mid }(\tilde{A}-A{)}_{kj}\mathrm{\mid }\le {ϵ}_0$ for each entry, so the complex modulus laws [F21] give $\mathrm{\mid }(({\tilde{A}}^{\ast }-A^{\ast })f{)}_j\mathrm{\mid }\le {\sum }_{k\in \mathcal{K}}\mathrm{\mid }(\tilde{A}-A{)}_{kj}\mathrm{\mid }\mathrm{\mid }{f}_k\mathrm{\mid }\le {ϵ}_0\mathrm{\parallel }f{\mathrm{\parallel }}_1\mathrm{.}$ It remains to bound the finite-arithmetic discrepancy. Let $v=D^{-1}f$ be the exact deconvolved mode vector and $v^c$ the dyadic centers of the component rectangles computed in step 4 of the algorithm. That step defines $E_f,V$ with $\mathrm{\parallel }v-v^c{\mathrm{\parallel }}_2\le E_f$ and $\mathrm{\parallel }{v}^c{\mathrm{\parallel }}_2\le V$ by [F12]–[F13]. The positive-sign computed unitary FFT has axis error factors ${\rho }_r$ and total $\rho ={\prod }_r(1+{\rho }_r)-1$; its bound is the same as the type-1 calculation in step 2.1 because [F5] gives the same factor for the conjugated inverse routine. Thus its computed center differs from the exact transform of the embedded $v^c$ by at most $\rho V$, since zero-fill embedding preserves the Euclidean norm. The exact unitary transform preserves the embedded input uncertainty norm by [F6], so [F18] bounds the computed FFT center's distance from the exact transform of embedded $v$ by $E_f+\rho V$. After scaling by $P^{-1\mathrm{/}2}$, each used grid coordinate is therefore within $(E_f+\rho V)\mathrm{/}\sqrt{P}$ of the exact $Q^{\ast }v$ coordinate. The algorithm encloses this discrepancy, evaluates every real spread-window coefficient with a finite image interval containing its exact value by [F4] and [F20], and applies only sound outward rectangle operations. Hence the final rectangle at source $j$ contains $(S^{\ast }{Q}^{\ast }{D}^{-1}f{)}_j=({\tilde{A}}^{\ast }f{)}_j$; its center has error at most its radius $r_j$. Taking the maximum over the finite output set gives $\mathrm{\parallel }\widehat{g}-{\tilde{A}}^{\ast }f{\mathrm{\parallel }}_{\mathrm{\infty }}\le B_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}^{\ast }$. Combining this finite bound with the displayed analytic entry estimate proves the type-2 claim, using [F12] for the finite norms.

There are $P\mathrm{/}{M}_r$ length-$M_r$ transforms on axis $r$, so [F8] gives $O(P{\log }_2{M}_r)$ work on that axis. Since $1\le M_r\le P$ and $2>1$, [F11] gives $\log 2>\log 1=0$. Formula [F16] therefore shows that ${\log }_2$ is increasing on positive arguments, hence ${\log }_2{M}_r\le {\log }_{2}P$; fixed $d$ gives $O(P\log P)$ over all axes. By [F24], forming the coordinate neighborhoods takes $O({\sum }_r(2w_r+1))$ index operations per source and produces at most ${\prod }_r(2w_r+1)$ retained product cells. For fixed $d$ and $L_r\ge 1$, the support-generation cost is bounded by the displayed per-source grid term. Evaluating its coordinatewise image sums costs at most a constant multiple of ${\sum }_r(2L_r+1)$ per retained product cell, and the mode deconvolution pass costs $O(\mathrm{\mid }\mathcal{K}\mathrm{\mid })$; interval-series evaluator internals are excluded as stated. These counts yield the displayed bound. This analysis covers types 1 and 2 only; the type-3 frequency ranges require separate scaling and work estimates.