Composite Polynomial Filtering

• Composite polynomial filtering is a technique that synthesizes complex spectral responses by combining simpler polynomial blocks, enabling precise control over approximation and modularity.
• It underpins diverse applications such as Hermitian eigensolvers, graph learning, and digital filter design by strategically composing inner and outer polynomial functions.
• Its design leverages functional composition, bandwise summation, and iterative products to balance selectivity, locality, and computational efficiency in various domains.

Composite polynomial filtering denotes a family of constructions in which a target spectral transformation, transfer function, or projection is realized through a structured combination of polynomial components rather than by a single undifferentiated polynomial. In the cited literature, that combination appears as outer–inner composition $F(G(z))$, as a sum of bandwise spectral polynomials, as a multivariate polynomial $h({\bf S}_1,\ldots,{\bf S}_d)$ in commuting shifts, or as an iterationwise product $P^{(s)}P^{(s-1)}\cdots P^{(1)}$ accumulated by repeated filtering. The concept therefore spans Hermitian interior eigensolvers, graph filtering and collaborative filtering, modular DSP design, multiband filter banks, and factorization-free approximation of matrix functions (Li et al., 2015, Lingam et al., 2021, Demirtas et al., 2015, Kang et al., 12 Jul 2025).

1. Algebraic forms and scope

The cited literature does not use a single canonical formalization. One recurrent form is functional composition, in which a subfilter $G$ is transformed by an outer polynomial

$$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$

A second form is multivariate shift composition,

$${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$

which requires commuting graph shifts. A third form is iterationwise composition, in which the effective filter after repeated subspace filtering is

$\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$

These formulas already show that “composite” may refer to composition in the operator argument, in multiple commuting shifts, or in algorithmic time (Demirtas et al., 2015, Emirov et al., 2020, Xu et al., 1 Apr 2026).

A separate but related interpretation appears when a global spectral response is assembled from localized polynomial pieces. In graph learning this is expressed as a sum of polynomials supported on different spectral bands, while in multiband signal processing it appears as a product of elementary polynomial matrix factors. A plausible implication is that composite polynomial filtering is best understood as a structural principle: a complicated response is synthesized from simpler polynomial blocks whose locality, degree, or placement is easier to control than that of a single global approximant.

2. Hermitian eigenvalue problems and spectrum slicing

In Hermitian interior eigensolvers, composite polynomial filtering is used to turn interior eigenvalues into dominant ones for Krylov methods. A representative construction first maps the spectrum of $A$ to $[-1,1]$ by

$$\widehat A=\frac{A-cI}{d}, \qquad c=\frac{\lambda_1+\lambda_n}{2}, \quad d=\frac{\lambda_1-\lambda_n}{2},$$

then builds a Chebyshev filter $h({\bf S}_1,\ldots,{\bf S}_d)$0 so that $h({\bf S}_1,\ldots,{\bf S}_d)$1 on a target interval $h({\bf S}_1,\ldots,{\bf S}_d)$2 and is small elsewhere. The 2015 Thick-Restart Lanczos scheme constructs $h({\bf S}_1,\ldots,{\bf S}_d)$3 from a Chebyshev approximation to a Dirac delta centered at $h({\bf S}_1,\ldots,{\bf S}_d)$4, uses damping by Jackson or Lanczos $h({\bf S}_1,\ldots,{\bf S}_d)$5-factors, chooses the degree $h({\bf S}_1,\ldots,{\bf S}_d)$6 until $h({\bf S}_1,\ldots,{\bf S}_d)$7 and $h({\bf S}_1,\ldots,{\bf S}_d)$8, and then balances the center so that $h({\bf S}_1,\ldots,{\bf S}_d)$9. The filtered operator $P^{(s)}P^{(s-1)}\cdots P^{(1)}$0 is handled by Thick-Restart Lanczos with deflation, while spectrum slicing assigns a distinct filter to each subinterval, yielding the “multi-filter / multi-slice” behavior that the paper explicitly identifies as a broad sense of composite polynomial filtering (Li et al., 2015).

The 2026 adaptive variant replaces the Dirac-centered design by a Chebyshev expansion of the step function

$P^{(s)}P^{(s-1)}\cdots P^{(1)}$1

and uses damped partial sums

$P^{(s)}P^{(s-1)}\cdots P^{(1)}$2

Here the degree $P^{(s)}P^{(s-1)}\cdots P^{(1)}$3 changes from one iteration to the next according to Ritz-value information, so the cumulative filter becomes the product polynomial $P^{(s)}P^{(s-1)}\cdots P^{(1)}$4. The same work makes the composite aspect explicit in time: each iteration contributes a different polynomial, and the total attenuation on each eigencomponent is the product of the per-iteration attenuation factors. At the implementation level it uses MaSpMM for the filtering step, and its reported average speedups are about $P^{(s)}P^{(s-1)}\cdots P^{(1)}$5 over EVSL and $P^{(s)}P^{(s-1)}\cdots P^{(1)}$6 over CJ-FEAST (Xu et al., 1 Apr 2026).

These eigensolver formulations also clarify a common misconception. Composite polynomial filtering in this setting does not require an explicit product such as $P^{(s)}P^{(s-1)}\cdots P^{(1)}$7 at the design stage. It can arise from a family of slice-specific filters, from repeated application of the same filter inside Lanczos or subspace iteration, or from an adaptive sequence of distinct filters chosen online.

3. Graph spectral learning and collaborative filtering

In graph learning, composite polynomial filtering is often bandwise rather than purely compositional in depth. PP-GNN defines a graph filter not as one global polynomial but as a sum of several polynomials, each acting on a different spectral band, together with a global GPR-style polynomial: $P^{(s)}P^{(s-1)}\cdots P^{(1)}$8 This enlarges the space of realizable filters and of adapted graphs, and the paper proves that the best approximation obtained by the “global + local” family is never worse than that of global polynomials alone. Empirically, the method reports performance gains of up to $P^{(s)}P^{(s-1)}\cdots P^{(1)}$9 over state-of-the-art models and attributes much of the benefit to the ability to learn different low- and high-frequency responses on heterophilic graphs (Lingam et al., 2021).

A second graph formulation uses Jacobi polynomial bases to construct complementary band-stop and band-pass filters on the normalized adjacency. JGCF defines

$G$0

uses this as a band-stop component that emphasizes low and high frequencies, and forms a complementary band-pass component through $G$1, followed by a $G$2 nonlinearity. The final representation concatenates the two outputs. The paper’s frequency-domain analysis argues that low- and high-frequency components correlate more strongly with future interactions than mid-frequency components, and its ablations report that removing either the band-stop or the band-pass branch degrades performance. It also reports gains of at most $G$3 on Alibaba-iFashion (Guo et al., 2023).

A third line of work uses Krylov subspaces to approximate polynomial graph filters without explicitly storing the graph operator. Mem-GF considers filters

$G$4

but replaces the full item-similarity matrix by a user-specific Lanczos basis $G$5 and tridiagonal projection $G$6, yielding

$G$7

When the polynomial degree $G$8 satisfies $G$9, the representation is exact in exact arithmetic. This turns polynomial graph filtering into a composite process in which the Krylov basis is built once and several filters, such as Mem-GF-1, Mem-GF-2, Mem-GF-3, and Mem-GF-5, are evaluated in the reduced space. The method reports up to $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$0 lower memory usage and $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$1 speedup in runtime while scaling to datasets with tens of millions of interactions (Park et al., 19 Jun 2026).

4. Distributed and spatio-temporal polynomial filtering

For distributed graph signal processing, composite polynomial filtering is frequently multishift and iterative. The multi-shift framework assumes commuting graph shifts $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$2 and defines

$$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$3

A direct inverse $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$4 usually has full geodesic-width even when $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$5 has small width, so the paper proposes to approximate the inverse by another polynomial filter $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$6 and then apply the Neumann-type iteration

$$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$7

Two concrete choices are given: IOPA, based on optimal polynomial approximation on the joint spectrum, and ICPA, based on multivariate Chebyshev approximation on a bounding box for the joint spectrum. Both are designed so that $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$8, which yields exponential convergence while preserving a distributed implementation because every application of $$F(x)=\sum_{k=0}^{K} f_k x^k, \qquad F(G(z))=\sum_{k=0}^{K} f_k G^k(z).$$9 and $${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$0 remains a local multi-shift polynomial filtering step (Emirov et al., 2020).

A spatio-temporal variant appears in the study of stationary graph signals and Kalman filtering. There, a stationary signal is generated by transmitting white noise through a polynomial graph channel

$${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$1

and state and observation models are polynomial in the same symmetric shift,

$${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$2

with $${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$3 and $${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$4. The paper shows that the Kalman gains $${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$5 and error covariances $${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$6 are themselves polynomials of the graph shift, and that Kalman filtering preserves stationarity of both estimates and errors. In this setting, composite polynomial filtering becomes a closed spatio-temporal algebra: repeated prediction, observation, and update operations remain inside the polynomial functional calculus of $${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$7 (Chen et al., 16 Sep 2025).

An older but closely related example is distributed consensus acceleration. There the filter is a polynomial $${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$8 of the consensus matrix, and one filtered update is equivalent to

$${\bf H}=h({\bf S}_1,\ldots,{\bf S}_d) =\sum_{l_1=0}^{L_1}\cdots\sum_{l_d=0}^{L_d} h_{l_1,\ldots,l_d}\,{\bf S}_1^{l_1}\cdots{\bf S}_d^{l_d},$$9

The coefficients are obtained by minimizing the spectral radius of the filtered operator through a semidefinite program, or by a Hermite-interpolation heuristic when detailed spectral information is unavailable. The paper shows that the method accelerates consensus on static networks and can also improve convergence in dynamic topologies, although high degree can become brittle when the topology changes too frequently (0802.3992).

5. Filter sharpening, modular design, and multiband filter banks

In classical DSP, composite polynomial filtering appears most explicitly as filter sharpening. Given a subfilter $\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$0, one chooses a polynomial

$\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$1

and forms the sharpened filter

$\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$2

The design problem is to minimize the infinity norm of the approximation error between the desired response and $\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$3. The paper formulates this as a generalized minimax problem and solves it with the First Algorithm of Remez, thereby extending sharpening beyond even-symmetric FIR filters to complex-valued FIR subfilters, causal IIR filters, and continuous-time filters. For causal IIR subfilters $\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$4, the composite response

$\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$5

introduces no new pole locations; it changes only multiplicities, so stability is inherited from the subfilter. The same framework also studies approximate functional decomposition, $\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$6, as a route to modular implementation (Demirtas et al., 2015).

A different DSP interpretation is given by multiband filter banks. There the global analysis/synthesis chain is encoded by an $\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$7 matrix $\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$8 with polynomial entries, and the entire multiband process is factorized into elementary upper- and lower-triangular polynomial matrices. In the $\Pi_s(A):=P^{(s)}P^{(s-1)}\cdots P^{(1)}.$9 case one has a factorization of the form

$A$0

and the paper emphasizes that finiteness follows from the Euclidean algorithm for polynomials. In this language, composite polynomial filtering is a cascade of elementary lifting or band-mixing steps, with down-sampling and up-sampling absorbed into the same matrix representation (Jorgensen et al., 2014).

These two DSP strands illuminate a common structural feature. Composite polynomial filtering may be additive, as in tapped cascades $A$1, or multiplicative, as in lifting factorizations. In both cases, the goal is modularity: a complex target response is realized through repeated use of a smaller collection of polynomial primitives.

6. Matrix-function approximation and approximation-theoretic perspective

A recent matrix-analytic example treats orthogonal projection onto the positive semidefinite cone as a spectral filtering problem. For a symmetric matrix $A$2, the exact projection is

$A$3

The factorization-free method approximates the scalar nonlinearity $A$4 by a composite polynomial

$A$5

where each $A$6 is low degree, and reconstructs the matrix approximation through a sign-based identity. After scaling $A$7, the algorithm applies the stages $A$8 and returns

$A$9

The reported half-precision implementation reaches a consistent relative error of $[-1,1]$0 with only $[-1,1]$1 matrix-matrix multiplications and roughly a $[-1,1]$2 speed-up over NVIDIA’s cuSOLVER routines; for a $[-1,1]$3 dense symmetric matrix, the paper reports approximately $[-1,1]$4 ms in half precision and $[-1,1]$5 ms in single precision on B200 GPUs (Kang et al., 12 Jul 2025).

An approximation-theoretic counterpart studies the error of the best polynomial approximation of composite functions in weighted spaces. For $[-1,1]$6, $[-1,1]$7, and Jacobi weight $[-1,1]$8, it proves that if $[-1,1]$9 and each component $$\widehat A=\frac{A-cI}{d}, \qquad c=\frac{\lambda_1+\lambda_n}{2}, \quad d=\frac{\lambda_1-\lambda_n}{2},$$0, then

$$\widehat A=\frac{A-cI}{d}, \qquad c=\frac{\lambda_1+\lambda_n}{2}, \quad d=\frac{\lambda_1-\lambda_n}{2},$$1

with $$\widehat A=\frac{A-cI}{d}, \qquad c=\frac{\lambda_1+\lambda_n}{2}, \quad d=\frac{\lambda_1-\lambda_n}{2},$$2. This is not a filtering algorithm, but it provides a direct characterization of how the smoothness of an outer map and an inner transformation controls the error of best polynomial approximation of the composite response (Fermo et al., 2023).

Taken together, these results indicate that composite polynomial filtering is less a single algorithm than a recurrent design pattern. The cited literature repeatedly replaces a monolithic high-degree or nonpolynomial target by structured combinations of lower-degree, band-limited, or stage-wise polynomials. The resulting trade-offs are domain specific but consistent: spectrum slicing reduces the number of desired eigenpairs per filter at the cost of higher degree per slice; piece-wise graph filters enlarge the filter family while introducing band-selection choices; Krylov compression lowers memory while tying exactness to polynomial degree; and low-degree matrix-function compositions exchange eigendecomposition for a fixed number of GEMMs. In that sense, composite polynomial filtering functions as an organizing principle for balancing selectivity, locality, modularity, and computational structure across several branches of numerical linear algebra, graph signal processing, and digital filter design.