S4 convolution-filter arithmetic and numerical cost

By [F5], each fixed matrix $G^S,G^T,G^L$ admits equivalent DPLR coordinates of bounded rank. When such a matrix is used as a constant state matrix and satisfies the stated invertibility hypotheses, it is an admissible input here. This does not identify the original time-dependent LegS ODE with a constant-matrix model. The operation argument applies to any supplied DPLR matrices of the same dimensions. By [F2], each required root value has the stated Cauchy-product expression, including the direct value at $-1$. Under the noncollision and correction-invertibility assumptions, the filter procedure in [F1] therefore computes all $M$ generating values exactly.

For $z\mathrm{\ne }-1$, $R=\beta D$. Substitution in [F2] gives $H_M(z)=x-\beta u(I_r+\beta T{)}^{-1}v$ and $W=I_r+\beta T$. Also $\mathrm{\mid }\beta \mathrm{\mid }\le 1$. The pole mapping in [F2] expresses each entry of $D$ as $-2\mathrm{/}(a+{\lambda }_j)(z-{\mu }_j{)}^{-1}$ for ${\lambda }_j\mathrm{\ne }-a$, and as the constant $\mathrm{\Delta }\mathrm{/}2$ otherwise. Thus the weights in [F6] are exactly $-2\overline{f_j}{g}_j\mathrm{/}(a+{\lambda }_j)$ for the applicable vector pair $f,g$. Their preparation and the exceptional constant sums cost $O(N)$ for fixed $r$. At $-1$, retain the direct value $(\mathrm{\Delta }\mathrm{/}2){\tilde{C}}_M^{\ast }B$ from [F2], without a resolvent division.

Each scalar product has at most $N$ poles and is evaluated on $M$ roots. If it has at least one pole, [F3] bounds its cost by $O((N+M){\log }^2(2+N+M))$ operations and $O(N+M)$ space; if it has no poles, its constant value is formed by a single finite sum in $O(N+M)$ operations and space. There are at most $(r+1{)}^2$ such products by [F2], a constant number because $r$ is fixed. The per-root correction matrix products and inversions cost $O(M(r+1{)}^3)$ operations by finite-dimensional elimination, which is $O(M)$ for fixed $r$.

The generating values are the DFT of the $M$ coefficients $C^{\ast }{\bar{A}}^j\bar{B}$ for $0\le j<M$. By [F4], the inverse radix-two transform recovers these coefficients exactly; [F1] then returns their first $L$ entries, which are precisely the length-$L$ filter.

Let $\delta =\epsilon \mathrm{/}[32r{w}^2(b+1{)}^2]$ and first consider exact arithmetic on truncated Cauchy sums having scalar entry errors at most $\delta$, with finite constant sums exact. Denote the resulting blocks by $\widehat{x},\widehat{u},\widehat{v},\widehat{T}$. Each of their infinity-norm errors is at most $e=r\delta$, since each row has at most $r$ entries and $r\ge 1$. Hence $\mathrm{\parallel }\widehat{W}-W{\mathrm{\parallel }}_{\mathrm{\infty }}\le e$, $we<1\mathrm{/}2$ and $e\le 1$. Writing $E=W^{-1}(\widehat{W}-W)$, the geometric series ${\sum }_{j\ge 0}(-E{)}^j$ converges in infinity norm because $\mathrm{\parallel }E{\mathrm{\parallel }}_{\mathrm{\infty }}\le we<1\mathrm{/}2$, is the inverse of $I+E$ by finite geometric identities and passage to the norm limit, and has norm at most $2$. Therefore $\mathrm{\parallel }{\widehat{W}}^{-1}{\mathrm{\parallel }}_{\mathrm{\infty }}\le 2w$. The identity ${\widehat{W}}^{-1}-W^{-1}=-W^{-1}(\widehat{W}-W){\widehat{W}}^{-1}$ gives $\mathrm{\parallel }{\widehat{W}}^{-1}-W^{-1}{\mathrm{\parallel }}_{\mathrm{\infty }}\le 2w^{2}e$.

Since $M$ is the least power of two at least $L$, $M<2L$ for every integer $L\ge 1$. Combining the Cauchy products, per-root corrections, and inverse transform gives $O((N+L){\log }^2(2+N+L))$ exact complex arithmetic operations. Processing the Cauchy products sequentially uses $O(N)$ product-tree workspace; their $(r+1{)}^2$ root-value arrays, the correction matrices, and the inverse transform use $O(M(r+1{)}^2)$ additional space. Fixed $r$ and $M<2L$ give total $O(N+L)$ working space. The proof counts the supplied ${\tilde{C}}_M$ as input and does not count its formation from $C$.

Expanding the difference of the two triple products as $(\widehat{u}-u){\widehat{W}}^{-1}\widehat{v}+u({\widehat{W}}^{-1}-W^{-1})\widehat{v}+u{W}^{-1}(\widehat{v}-v)$ and using step 2.3 bounds its modulus by $2we(b+1)+2b{w}^{2}e(b+1)+bwe$. Adding the error in $x$ and using $\mathrm{\mid }\beta \mathrm{\mid }\le 1$ bounds each generating-value error by $e[1+2w(b+1)+2b{w}^2(b+1)+bw]\le 9w^2(b+1{)}^{2}e\le 9\epsilon \mathrm{/}32$. The direct $-1$ value has no truncation error. The exact inverse DFT averages $M$ unit-modulus multiples in each output coordinate, so its infinity operator norm is at most one; consequently the exact inverse transform of these approximations has coefficient error at most $9\epsilon \mathrm{/}32$.

Apply [F6] with absolute tolerance $\delta$ to each scalar sum; its exact truncated circuit has error at most $\delta$, as required above. Although [F6] is stated on ${\zeta }_j=e^{2\pi ij\mathrm{/}M}$, the filter grid is $z_k=e^{-2\pi ik\mathrm{/}M}={\zeta }_{(-k<mrow><mi\; mathvariant="normal">m</mi><mi\; mathvariant="normal">o</mi><mi\; mathvariant="normal">d</mi></mrow>M)}$; run each numerical Cauchy product on the positive grid and permute its outputs to the filter order, in $O(M)$ operations and words, preserving the maximum error and pole/root separation bounds. It uses $p=O(\log (s\mathrm{/}\delta ))=O(\log (\kappa \mathrm{/}\epsilon ))$ terms. There are at most $(r+1{)}^2$ sums, processed sequentially; their arrays and the fixed-size blocks use $O(N+M)$ words. The per-root inverses can be computed by the adjugate/determinant formula, which is constant work for fixed $r$ and has nonzero determinants by step 2.3. Add the radix-two inverse transform and finite sums: the scalar count is $O((N+M)\log (2+N+M)\log (\kappa \mathrm{/}\epsilon ))$. The replacement $M<2L$ gives the stated count and space.

To include every rounding error, regard the entire selected computation, from supplied DPLR parameters, transformed vector and root/arc constants through weight formation, Cauchy sums, finite constants, adjugate inverses and inverse FFT, as one finite rational circuit. Include conjugations as exact sign changes. Require certified $T_0,V_0,D_0$ with $T_0$ an upper bound on its gates, $V_0\ge 1$ an upper bound on exact intermediate moduli, and $0<D_0\le 1$ a lower bound on every exact divisor modulus. Such bounds exist for each input because step 2.3 proves all truncated correction inverses exist and all other divisions are nonzero by the hypotheses. Put $K_0=64(V_0+1{)}^2\mathrm{/}{D}_0^2$ and select input and gate accuracy $u\le \min (1,D_0\mathrm{/}2,\epsilon \mathrm{/}2)\mathrm{/}(2K_0^{T_0+1})$, also sufficient for the certified overlapping routing tests of [F6]. The gate recurrence proved in [F6] then gives final additional error at most $\epsilon \mathrm{/}2$, including errors in supplied parameters, roots, weights, $\beta$, correction inverses and twiddle factors. Together with step 3.2 this is less than $\epsilon$. No precision selection changes the scalar count or scalar-word count, although their bit cost can grow. Bounded $\kappa$ gives $\log (\kappa \mathrm{/}\epsilon )=O(\log (1\mathrm{/}\epsilon ))$; arbitrary conditioning retains the explicit $\kappa$ term.