Filtered Spectral Projection Algorithm (FSPA)

• FSPA is a framework that uses spectral filtering to isolate or amplify targeted invariant subspaces across applications such as linear solvers, data assimilation, and quantum PCA.
• It employs diverse operators—including resolvents, Fourier projections, and power maps—to suppress unwanted spectral components and improve convergence and stability.
• FSPA’s design pattern adapts to different computational challenges, offering benefits like enhanced deflation in ill-conditioned systems, synchronization in fluid dynamics, and efficient projection in quantum algorithms.

Filtered Spectral Projection Algorithm (FSPA) denotes a class of procedures in which a spectral filter or spectral projector is used to isolate a targeted invariant subspace and then transfer that subspace information to a downstream task. In the available literature, this includes a contour-integral spectral projection preconditioner for ill-conditioned linear systems, a spectrally-filtered discrete-in-time downscaling data assimilation scheme for the two-dimensional Navier–Stokes equations, and a projection-first framework for quantum principal component analysis (Yeung et al., 2016, Celik et al., 2018, Hossain et al., 13 Mar 2026). The operators differ—resolvents $(zI-A)^{-1}$, Fourier projections $P_\lambda$, and powers of a density operator $\rho$—but all three formulations use filtering to preserve a chosen spectral region and suppress complementary components.

1. Conceptual structure

In the linear-systems setting, the central object is the classical spectral projector

$$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$

where $\Gamma$ encloses a selected subset of eigenvalues of $A$. This projector is applied to a random matrix $Y$ to form $Z=P_\Gamma Y$, and $Z$ is then used to build a deflation preconditioner (Yeung et al., 2016). In the data-assimilation setting, the key filtered projector is

$J=P_\lambda P_\sigma I_h,\qquad E=I-J,$

where $P_\lambda$0 is a low-resolution interpolant observable, $P_\lambda$1 is the Leray projector, and $P_\lambda$2 is the orthogonal projection onto low Fourier modes; the algorithm repeatedly replaces the observed low-mode content by $P_\lambda$3 at observation times (Celik et al., 2018). In the qPCA setting, the filter is the power map $P_\lambda$4, and the normalized iterate

$P_\lambda$5

amplifies overlap with the dominant eigenspace of $P_\lambda$6 without explicit eigenvalue estimation (Hossain et al., 13 Mar 2026).

These formulations suggest a common abstraction: first specify a spectral region of interest, then realize a filter that isolates or amplifies that region, and finally use the filtered subspace in a solver, an assimilation update, or a dimensionality-reduction primitive. A plausible implication is that FSPA is better understood as a design pattern centered on spectral isolation than as a single fixed algorithm.

2. Contour-integral FSPA for ill-conditioned linear systems

For linear systems

$P_\lambda$7

the 2016 formulation targets ill-conditioned problems, especially those whose condition number is large because some eigenvalues are very close to the origin. The projector $P_\lambda$8 is chosen so that $P_\lambda$9 encloses those troublesome eigenvalues, typically a small circle around $\rho$0. Applying $\rho$1 to a random matrix $\rho$2 produces

$\rho$3

whose columns almost surely span the targeted spectral subspace when $\rho$4 is suitably chosen. Because the contour integral cannot be evaluated exactly, it is approximated on a circle $\rho$5 by Legendre–Gauss quadrature,

$\rho$6

so the projector is realized as a rational filter built from shifted resolvents (Yeung et al., 2016).

The filtered subspace $\rho$7 defines the deflation matrices

$\rho$8

The preconditioned problem is the singular but consistent deflated system

$\rho$9

When $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$0 consists of eigenvectors $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$1, the spectrum of $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$2 becomes

$$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$3

so the deflated eigenvalues are replaced by zeros and no longer affect the nonsingular part on which GMRES effectively iterates. The associated solution decomposition is

$$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$4

The convergence analysis explains the gain. Without deflation, for diagonalizable $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$5, the GMRES residual satisfies a standard ellipse-based bound involving $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$6. With exact-eigenvector deflation, the analogous bound replaces $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$7 by $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$8, where $$P_\Gamma=\frac{1}{2\pi \sqrt{-1}}\oint_\Gamma (zI-A)^{-1}\,dz,$$9 comes from a QR factorization of the eigenvector matrix. The remaining ellipse needs to contain only the undeflated eigenvalues, which can substantially improve the convergence factor (Yeung et al., 2016).

The numerical experiments illustrate the mechanism. For a convection–diffusion problem with $\Gamma$0, $\Gamma$1, and a contour centered at $\Gamma$2 with radius $\Gamma$3, plain GMRES required $\Gamma$4 iterations with $\Gamma$5 and $\Gamma$6. Deflated GMRES using $\Gamma$7 from spectral projection with $\Gamma$8 and $\Gamma$9 required $A$0 iterations, with $A$1 and $A$2. For the sparse benchmark matrix bcsstm27, unpreconditioned MBiCG did not converge within $A$3 iterations and after $A$4 iterations still had $A$5, whereas a spectral-projection-based deflated run with improved $A$6 required $A$7 iterations and achieved $A$8 (Yeung et al., 2016).

A distinctive feature of this formulation is that the method is problem-independent in the sense stated in the paper: it uses only the matrix $A$9 and the contour $Y$0, not problem-specific eigenvector guesses or subdomain structures. The dominant practical cost is computing $Y$1, which requires solving $Y$2 shifted systems and motivates parallelization, multigrid acceleration, and numerical rank refinement through CGE-based column selection when $Y$3 is rank-deficient or nearly rank-deficient.

3. Spectrally filtered projection in discrete-in-time data assimilation

In the 2D incompressible Navier–Stokes setting on the periodic box $Y$4, the reference solution $Y$5 solves

$Y$6

with $Y$7 the Stokes operator and $Y$8. Observations are available only at discrete times

$Y$9

through a low-resolution interpolant observable $Z=P_\Gamma Y$0. The filtered observation operator is defined by

$Z=P_\Gamma Y$1

where $Z=P_\Gamma Y$2 is the orthogonal projection onto Fourier modes with $Z=P_\Gamma Y$3. The algorithm initializes with

$Z=P_\Gamma Y$4

updates by

$Z=P_\Gamma Y$5

and evolves continuously between updates via

$Z=P_\Gamma Y$6

The paper describes this as inserting new observational data directly into the dynamical model as it is being evolved over time, rather than nudging (Celik et al., 2018).

The filtering role of $Z=P_\Gamma Y$7 is mathematically explicit. Because $Z=P_\Gamma Y$8, the operator replaces only low Fourier modes, and the complementary operator $Z=P_\Gamma Y$9 retains the part not directly overwritten by observations. The paper proves norm estimates for $Z$0 using the approximation-of-the-identity properties of the interpolant. For type-I interpolants,

$Z$1

and for type-II interpolants,

$Z$2

These bounds allow the filtered interpolant $Z$3 to behave like a near-identity on the relevant low-mode content (Celik et al., 2018).

The main analytical result is global synchronization: for every fixed $Z$4, there exist $Z$5 and $Z$6, depending only on the interpolant constant $Z$7, on $Z$8, on $Z$9, and on the domain, such that

$J=P_\lambda P_\sigma I_h,\qquad E=I-J,$0

For type-I interpolants, the contraction is established in the $J=P_\lambda P_\sigma I_h,\qquad E=I-J,$1-norm through a discrete inequality of the form

$J=P_\lambda P_\sigma I_h,\qquad E=I-J,$2

and then converted into exponential continuous-time decay. For type-II interpolants, the argument is carried out in vorticity form and yields exponential convergence of vorticity, which implies exponential convergence of velocities in $J=P_\lambda P_\sigma I_h,\qquad E=I-J,$3 on the periodic domain (Celik et al., 2018).

This formulation extends earlier orthogonal-projection-based discrete-in-time schemes because $J=P_\lambda P_\sigma I_h,\qquad E=I-J,$4 need not itself be an orthogonal projection. The spectral filter $J=P_\lambda P_\sigma I_h,\qquad E=I-J,$5 transforms rough or non-orthogonal observables—such as local averages or nodal measurements—into a smooth low-mode correction that is compatible with the dissipative structure of the PDE. A plausible implication is that, in this setting, FSPA functions as a regularizing interface between irregular measurement operators and the continuous Navier–Stokes dynamics.

4. Projection-first FSPA for quantum principal component analysis

In the 2026 qPCA formulation, the objective is not explicit eigenvalue recovery but projection onto the dominant spectral subspace of a covariance-encoded density operator. The starting point is a Hermitian positive semidefinite operator

$J=P_\lambda P_\sigma I_h,\qquad E=I-J,$6

For amplitude-encoded centered classical data, the paper states that the ensemble density matrix

$J=P_\lambda P_\sigma I_h,\qquad E=I-J,$7

coincides with the covariance matrix. For uncentered data, $J=P_\lambda P_\sigma I_h,\qquad E=I-J,$8 corresponds to PCA without centering, and the paper derives eigenvalue interlacing bounds quantifying the deviation from standard centered PCA. It also states that ensembles of quantum states admit an equivalent centered covariance interpretation (Hossain et al., 13 Mar 2026).

The algorithm is an adaptive normalized power iteration. Given an initial state

$J=P_\lambda P_\sigma I_h,\qquad E=I-J,$9

FSPA repeatedly applies $P_\lambda$00 and normalizes after each application. In the version described in the paper, the input is a Hermitian operator $P_\lambda$01 with $P_\lambda$02, an initial warm-start state $P_\lambda$03, and a number of rounds $P_\lambda$04. One first normalizes $P_\lambda$05, sets an amplification parameter $P_\lambda$06, and for each round applies $P_\lambda$07 exactly $P_\lambda$08 times, normalizing after every application, then doubles $P_\lambda$09. Conceptually, the effective iterate is

$P_\lambda$10

The filter is therefore the spectral map $P_\lambda$11, which amplifies components corresponding to larger eigenvalues relative to the rest (Hossain et al., 13 Mar 2026).

A central structural property is eigenvalue magnitude invariance: replacing $P_\lambda$12 by $P_\lambda$13 with $P_\lambda$14 does not change the normalized iterates. This removes the magnitude sensitivity associated with estimation-first qPCA pipelines based on phase estimation. The paper highlights two claimed drawbacks of the estimation-first approach: precision overhead or magnitude sensitivity, and overkill for projection tasks. FSPA is introduced specifically as a projection-first framework that bypasses explicit eigenvalue estimation while preserving the essential spectral structure (Hossain et al., 13 Mar 2026).

Warm-start overlap is essential. In the nondegenerate case $P_\lambda$15, if $P_\lambda$16, the fidelity with the top eigenvector after $P_\lambda$17 effective power applications is

$P_\lambda$18

If the dominant eigenvalue is degenerate, $P_\lambda$19, the algorithm does not force a basis choice inside the dominant eigenspace; instead it converges to the normalized projection of the initial state onto

$P_\lambda$20

This subspace-focused behavior is one of the paper’s explicit design goals in near-degenerate regimes (Hossain et al., 13 Mar 2026).

The implementation viewpoint is compatible with block-encoding and QSVT-style polynomial transformations. Because normalization is non-unitary, the paper interprets it through post-selection or amplitude amplification rather than as a standalone unitary step.

5. Convergence, stability, and empirical behavior

The three formulations use different convergence mechanisms. In contour-integral deflation for linear systems, the improvement is spectral and geometric: once the targeted eigenvalues are deflated, GMRES effectively sees a smaller nonsingular problem whose residual bound depends on $P_\lambda$21 rather than $P_\lambda$22, and the enclosing ellipse need only contain the remaining eigenvalues (Yeung et al., 2016). In the Navier–Stokes assimilation scheme, convergence is a synchronization result: suitable choices of $P_\lambda$23, $P_\lambda$24, and $P_\lambda$25 produce a discrete contraction that yields exponential decay of the state error in time (Celik et al., 2018). In qPCA, the guarantee is overlap amplification: to obtain

$P_\lambda$26

the required number of applications of $P_\lambda$27 is

$P_\lambda$28

so complexity depends on spectral ratios and initial overlap rather than on absolute eigenvalue scale (Hossain et al., 13 Mar 2026).

The empirical behavior also differs by domain. In ill-conditioned linear systems, the filtered spectral projection preconditioner reduces iteration counts substantially and can improve accuracy in the norm of the distance to the exact solution. The convection–diffusion example and the sparse benchmark matrices bcsstm27 and mahindas show that spectral-projection-based deflation can turn essentially stagnant Krylov behavior into practical convergence, although the quality of $P_\lambda$29, the conditioning of $P_\lambda$30, and column compression of $P_\lambda$31 materially affect performance (Yeung et al., 2016). In qPCA, the reported numerical demonstrations show three features: eigenvector-level instability can coexist with stable dominant subspaces; Lloyd-style qPCA collapses below a resolution threshold under global spectral downscaling whereas FSPA remains stable; and FSPA degrades smoothly rather than abruptly as the spectral gap shrinks. On benchmark datasets including Breast Cancer Wisconsin and handwritten Digits, downstream performance remains stable whenever projection quality is preserved (Hossain et al., 13 Mar 2026).

Taken together, these results suggest that FSPA methods are usually justified less by a uniform asymptotic advantage than by a task-appropriate spectral representation. In each domain, the filter is tailored to the spectral object that actually controls computational difficulty: small eigenvalues near the origin in linear solvers, unresolved high Fourier content in data assimilation, and dominant principal subspaces in qPCA.

6. Relations to neighboring methods and interpretive issues

FSPA sits near several established research lines. The linear-systems formulation is rooted in deflation methods and projectors for Krylov solvers, and its contour-integral construction is closely related to quadrature-based spectral filtering and eigenspace extraction methods such as Sakurai–Sugiura and FEAST; the paper explicitly characterizes it as FEAST-style contour filtering applied not to eigenvalue problems per se, but to preconditioning linear systems via deflated Krylov methods (Yeung et al., 2016). The data-assimilation formulation extends earlier discrete-in-time schemes in which the observation operator itself was an orthogonal projection onto low Fourier modes; the key generalization is that the observable need only satisfy an approximation-of-the-identity property, while the added spectral projector $P_\lambda$32 restores the low-mode structure needed for analysis (Celik et al., 2018). The qPCA formulation is presented against Lloyd-style, phase-estimation-based qPCA and argues that for many applications spectral projection is the essential primitive, while explicit eigenvalue estimation is often unnecessary (Hossain et al., 13 Mar 2026).

A recurring interpretive issue is that “spectral projection” does not mean the same implementation in every setting. In the linear-systems paper, it is a contour-integral projector approximated by rational quadrature. In the Navier–Stokes assimilation paper, it is a low-mode Fourier projection composed with an interpolant observable. In the qPCA paper, it is a normalized power filter that approaches the dominant invariant subspace without numerically estimating eigenvalues. The literature therefore suggests that FSPA is not a single canonical routine but a family of filter-plus-projection constructions adapted to distinct operator classes and computational objectives.

Another common misunderstanding is that FSPA must return individual eigenvectors. The available formulations do not support that as a universal requirement. The linear-systems version needs a deflation subspace rather than explicit eigenpairs; the data-assimilation version inserts filtered low-mode observations rather than computing eigenvectors at all; and the qPCA version is explicitly subspace-oriented, with convergence to the dominant invariant subspace in degenerate cases and no artificial symmetry breaking in the absence of bias. This suggests that the defining feature of FSPA is not eigenvalue estimation, but controlled isolation of a spectral subspace and its direct use in the target algorithm.