Online Spectral Deflation for State Constrained Optimal Control Problems

Abstract: Parametric PDE-constrained optimal control with pointwise state constraints requires repeated solution of restricted Schur-complement systems on parameter-dependent inactive sets. In a primal active-set method, each inactive-set system is symmetric positive definite, but the active set can change nonsmoothly with the parameter. The resulting operator may vary in dimension, sparsity pattern, and spectrum, limiting reuse of sparse factorizations, multigrid hierarchies, and Krylov information. We propose a reusable spectral-deflation strategy anchored to one full-domain reference Schur complement. Low reference eigenmodes are computed once, restricted online to each inactive set, and used as an A-DEF2 deflation basis for Jacobi-preconditioned CG. The framework also supports POD enrichment, Rayleigh-Ritz reselection, coarse-grid or analytical reference modes, and conditioning safeguards. Given the active set, the method preserves the high-fidelity inactive-set system and solves it to the prescribed CG tolerance; it accelerates the linear algebra rather than replacing the optimal-control solve with a surrogate. We explain the method through a spectral-coherence view, motivated by interlacing and perturbation arguments and assessed with principal-angle diagnostics. Across diffusion, convection-diffusion, nonlinear thermal, and conjugate-heat-transfer benchmarks, deflation reduces CG iterations by about 55 to 98 percent. GPU deployments also show wall-time gains over CPU sparse-direct and algebraic-multigrid baselines, because the reference basis is built once whereas competing solver structures are rebuilt per instance. Coarse-grid or analytical modes amortize the offline cost within a single parameter sweep; fine-grid eigensolves remain more precompute-limited. Timings isolate the inactive-set linear-solve kernel; reducing the active-set outer loop is outside the present scope.

• The paper introduces a reference-based spectral deflation method that reuses full-domain eigenmodes across changing inactive sets, reducing conjugate-gradient iterations by approximately 55–98%.
• The method combines restricted and QR-orthonormalized eigenmodes with optional Ritz stabilization, POD enrichment, coarse-grid prolongation, and analytical Laplacian modes to reduce setup costs and maintain numerical reliability.
• Experiments show up to 973× deployment speedups over CPU sparse direct solvers, 13–18× speedups over CPU AMG, and 56–85% iteration reductions in space–time problems, while preserving relative errors below 1.3×10⁻⁸.

This paper addresses a recurring computational bottleneck in parametric PDE-constrained optimal control with pointwise state constraints: after eliminating control and adjoint variables, each active-set iteration requires solving a Schur-complement system restricted to a parameter-dependent inactive set. Because the active set changes discontinuously across parameter instances, these restricted SPD systems vary in dimension, sparsity pattern, and spectrum, defeating standard reuse strategies such as sparse factorizations, AMG hierarchies, and instance-to-instance Krylov recycling. The authors propose anchoring reusable spectral information to a single full-domain reference Schur complement rather than to previously solved restricted instances (2606.17971).

Problem setting and reduction

The problem class is a distributed tracking-type optimal control problem with Tikhonov regularization $\alpha > 0$ and a pointwise state bound $y \le \psi$, both possibly parameter-dependent. After discretization (uniform finite differences on Cartesian grids throughout), elimination of the control and adjoint yields the full-domain SPD operator

$M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$

with $\lambda_{\min}(M) \ge 1$ regardless of whether $A$ is symmetric or invertible. The primal active-set iteration — a simplification of PDAS that updates the active set from primal violation alone — fixes $y_A = \psi_A$, sets the multiplier to zero on the inactive set, and produces the restricted system $M_{II} y_I = b$. This restricted solve is the kernel accelerated by the paper; all timings isolate it, and the outer active-set loop is explicitly out of scope. The authors acknowledge that this deliberate SPD route squares the state operator, giving biharmonic-like conditioning $\kappa(M) = O(1 + \alpha N^{4/d})$; they argue this is precisely what leaves room for deflation to recover, whereas instance-tuned saddle-point block preconditioners would have to be rebuilt as the inactive set changes.

Reusable spectral deflation

The method has three layers. First, the $r_{\mathrm{pool}}$ smallest eigenpairs of a fixed reference Schur complement $$M_{\mathrm{ref}} = \alpha A(\theta_{\mathrm{ref}})^\top A(\theta_{\mathrm{ref}}) + I$$ are computed once via shift-invert Lanczos (ARPACK). Second, at each query the modes are restricted to the current inactive index set and QR-orthonormalized, forming the deflation basis $y \le \psi$0 used in an A-DEF2 deflated CG scheme: a coarse-correction predictor $y \le \psi$1 with $y \le \psi$2, followed by Jacobi-preconditioned CG on the projected operator $y \le \psi$3, with convergence declared on the unprojected residual. Third, optional stabilization includes Rayleigh–Ritz reselection from an overcomplete pool, online POD enrichment from previously converged states merged under a conditioning guard ($y \le \psi$4 during construction, $y \le \psi$5 at solve time with fallback to plain CG), coarse-grid eigenmode prolongation, and closed-form analytical sine-product eigenmodes for tensor-product Laplacian references.

Two cost-reduction devices matter for amortization. Coarse-grid prolongation computes eigenmodes on a grid coarsened by factor $y \le \psi$6 and lifts them by multilinear interpolation; measured eigensolve speedups exceed the nominal $y \le \psi$7 DOF reduction because the shift-invert sparse factorization grows superlinearly (e.g., $y \le \psi$8–$y \le \psi$9 at $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$0 in 3D versus $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$1). For tensor-product Laplacian references, analytical Kronecker-product sine modes reduce the 500-mode precompute at $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$2 from roughly 1,344 s to 0.73 s — a $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$3 speedup — verified against ARPACK to principal angles of $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$4 at $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$5.

Spectral-coherence rationale

The paper motivates the reuse through two classical results, carefully labeled as scaffolding rather than proof. Cauchy interlacing bounds the eigenvalues of the principal submatrix $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$6 in windows whose width depends on the active fraction, not on the symmetric-difference distance between active sets; Davis–Kahan bounds eigenspace rotation under operator drift $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$7, which vanishes identically when $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$8 is parameter-independent. Directional coherence after restriction is not covered by either result and is assessed empirically: on three diagnostic problems (40K DOF, 60 instances), the principal angle at the deflation cutoff stays below $M(\theta) = \alpha A(\theta)^\top A(\theta) + I,$9, $\lambda_{\min}(M) \ge 1$0, and $\lambda_{\min}(M) \ge 1$1 respectively, even though worst-case angles over all modes approach $\lambda_{\min}(M) \ge 1$2. This two-regime structure — stable leading subspace, randomized trailing modes — is the empirical basis for choosing deflation ranks within the coherent regime. The authors also derive a heuristic conditioning wall, $\lambda_{\min}(M) \ge 1$3 in 2D versus $\lambda_{\min}(M) \ge 1$4 in 3D under Weyl asymptotics, which correctly predicts that raw-eigenmode deflation fails at moderate-to-high rank only in 2D; they present this as interpretation organizing the measurements, not as a rigorous bound.

Numerical results

Across 14 steady configurations plus one space–time family, deflation reduces CG iterations relative to cold Jacobi-CG by roughly 55–98%. Representative iteration reductions:

Configuration class	Cold CG	Reduction range
3D Laplacian / CHT ($\lambda_{\min}(M) \ge 1$5, $\lambda_{\min}(M) \ge 1$6)	390–1,604	60–83%
3D CHT ($\lambda_{\min}(M) \ge 1$7)	up to 1,604	up to 92%
2D ($\lambda_{\min}(M) \ge 1$8, Ritz-stabilized, $\lambda_{\min}(M) \ge 1$9)	2,995–29,092	84–95%
2D combined QR(500,20), solve-guarded	—	92–98%

Several findings deserve emphasis. Convection-dominated CHT problems are simultaneously the hardest (1,400–1,600 cold iterations) and among the most deflation-friendly, reaching 92% reduction; the authors note this may reflect convection-induced spectral separation but flag that a dedicated gap analysis is needed to confirm the mechanism. Warm-starting from the previous solution yields no measurable benefit (within ±1 iteration) across all configurations, because abrupt active-set changes make solution continuation useless — a clear negative result supporting spectral over state recycling. Raw eigenmodes hit the 2D conditioning wall (one of thirty instances falls back at $A$0 on 2d_asym, all at $A$1); Ritz reselection eliminates all divergences but is strongly detrimental on 3D convection-dominated problems ($A$2 to $A$3 at moderate $A$4), so the paper treats it as a regime-dependent rescue, not a universal enhancement. Notably, condition monitoring cannot detect the harmful-Ritz regime, since a well-conditioned basis can still be deflation-suboptimal; the policy to disable Ritz there is a static configuration choice based on the regime study.

Wall-time results are deployment comparisons (GPU deflated CG versus CPU baselines) and the paper labels them as such consistently. At $A$5 (125K DOF), per-instance speedups over CPU sparse direct reach $A$6–$A$7, driven by empirical exponents $A$8 for CPU direct versus $A$9 for GPU deflated CG over the tested range. Against CPU BoomerAMG — which achieves higher iteration reduction (91–98%) but must rebuild its hierarchy for every one of thirty instances — GPU deflated CG is $y_A = \psi_A$0–$y_A = \psi_A$1 faster per instance. With analytical precompute, amortized speedups over 30 instances reach $y_A = \psi_A$2–$y_A = \psi_A$3 with breakeven at the first instance; with coarse-grid prolongation, breakeven against AMG-RS falls to 15–32 instances at $y_A = \psi_A$4 and 10–19 at $y_A = \psi_A$5. An accuracy audit over roughly 12,000 per-instance solves shows all deflated variants match the sparse-direct reference to solver tolerance (max relative error $y_A = \psi_A$6), confirming the speedups do not trade accuracy.

In the space–time extension (48 parabolic CHT configurations, up to 540K DOF), constant-in-time spatial bases give only 1–19% reduction, while the Kronecker temporal expansion $y_A = \psi_A$7 achieves 56–85% (median 75%) — evidence that missing temporal degrees of freedom, not spatial content, dominate the residual error. Truncated cosine bases degrade sharply once $y_A = \psi_A$8 (median 27% at $y_A = \psi_A$9), showing temporal compression is safe only when the basis spans the full relevant time subspace. GPU Kronecker-deflated CG is $M_{II} y_I = b$0–$M_{II} y_I = b$1 faster per instance than CPU AMG-RS, with the gap widening with size.

Two further negative results are reported honestly: a shift-invert randomized eigensolve scales worse than ARPACK ($M_{II} y_I = b$2 vs $M_{II} y_I = b$3) with an 11–52% iteration penalty, and compressed-sensing recovery of inter-instance state updates fails because $M_{II} y_I = b$4 is not sparse in the eigenmode basis and the hardest instances are precisely where updates are densest.

Limitations and open questions

The paper is explicit about scope. Spectral coherence is empirically validated, not formally proved; strongly nonlinear regimes with substantial operator drift remain open, and the nonlinear thermal case is a verification (drift $M_{II} y_I = b$5) rather than an extension. Reference-operator choice outside the parameter-invariant setting is ad hoc (a diffusion-only baseline for CHT). The online policy — particularly when to enable Ritz reselection — requires problem-class knowledge. The largest grids tested are $M_{II} y_I = b$6 and $M_{II} y_I = b$7, limited by explicit Schur-complement assembly and ARPACK; scaling past $M_{II} y_I = b$8 DOF would require matrix-free Krylov–Schur eigensolvers, argued feasible but not demonstrated. CPU baselines use 8 MPI ranks, so larger deployments could narrow the GPU wall-time advantage, though iteration reductions are hardware-independent. Finally, the study is a kernel-level benchmark of the inner inactive-set solve with active fraction calibrated near 20%; end-to-end outer-loop costs, control constraints, and mixed formulations are excluded.

Conclusion

The paper demonstrates that a deflation basis anchored once to a full-domain reference Schur complement survives restriction to discontinuously changing inactive sets well enough to remove 55–98% of CG iterations across diffusion, convection–diffusion, nonlinear thermal, CHT, and space–time benchmarks, without replacing the high-fidelity constrained optimality system by a surrogate. Its practical value depends on cheap offline construction: analytical tensor-product modes remove the precompute barrier entirely, coarse-grid prolongation suffices otherwise, and a fine-grid eigensolve leaves gains precompute-dominated. The main unresolved issues are a formal account of directional coherence after restriction, principled reference selection under genuine operator drift, automatic online policy selection, and matrix-free eigensolves beyond the demonstrated scale.