Two-Sided Preconditioner Updates

• The paper introduces a two-sided preconditioning technique that splits a constant-coefficient approximation into left and right factors, resulting in a condition number independent of system size.
• The method exploits the diagonalizability of the spatial stiffness matrix via fast sine transforms and the Toeplitz structure of the temporal convolution for efficient inversion.
• The approach achieves near-optimal O(JN log(JN)) computational cost and scalable parallel-in-time updates for variable-coefficient non-local evolutionary equations.

Two-sided preconditioner updates refer to a class of preconditioning strategies for large, structured linear systems, especially those arising from the discretization of non-local evolutionary equations involving variable-coefficient elliptic operators and weakly singular temporal convolution kernels. In the context of parallel-in-time (PinT) algorithms for all-at-once systems, such preconditioners are constructed by splitting a suitably chosen constant-coefficient approximation of the system matrix into left and right factors, each of which can be efficiently inverted by exploiting the block-diagonal structure in space (via diagonalization) and the Toeplitz structure in time. The resulting two-sided preconditioning yields a preconditioned system with condition number provably independent of the system size, enables near-optimal parallel implementation, and admits efficient iterative and direct solvers, particularly for variable-coefficient problems where standard approaches are often suboptimal.

1. Construction of Two-Sided Preconditioners

Given the all-at-once matrix $A\in \mathbb{R}^{JN\times JN}$ arising from the discretization of a variable-coefficient Laplacian $L_a \in \mathbb{R}^{J\times J}$ (in space) and a Toeplitz temporal convolution matrix $T\in \mathbb{R}^{N\times N}$, the two-sided preconditioning technique begins by selecting a constant $\beta$ such that $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$ and $\beta = \sqrt{a_{\min} a_{\max}}$. The constant-coefficient approximation is then defined as

$P := \beta L_1 \otimes I_N + I_J \otimes T,$

where $L_1$ is the stiffness matrix for $-\Delta$ with unit diffusion.

The split is performed as follows:

• Right preconditioner (acts last): $P_R := (\beta L_1)^{1/2}\otimes I_N$.
• Left preconditioner (acts first): $L_a \in \mathbb{R}^{J\times J}$0.

The two-sided preconditioned system to solve is

$L_a \in \mathbb{R}^{J\times J}$1

where $L_a \in \mathbb{R}^{J\times J}$2 is the discretized right-hand side.

2. Derivation and Splitting Strategy

Starting from $L_a \in \mathbb{R}^{J\times J}$3, the requirement that $L_a \in \mathbb{R}^{J\times J}$4 contains only a fractional power of $L_a \in \mathbb{R}^{J\times J}$5 and no $L_a \in \mathbb{R}^{J\times J}$6 leads uniquely to the given split. Algebraically,

$L_a \in \mathbb{R}^{J\times J}$7

by factoring $L_a \in \mathbb{R}^{J\times J}$8 with $L_a \in \mathbb{R}^{J\times J}$9.

This split isolates the main spatial stiffness (constant-coefficient Laplacian) on the right, facilitating diagonalization, and collects the remaining convolution and compensation for variable coefficients on the left.

3. Explicit Inversion Using Diagonalization and Toeplitz Structure

For efficient inversion and matrix-vector multiplication, the following structure is exploited:

• The stiffness matrix $T\in \mathbb{R}^{N\times N}$0 can be diagonalized by a fast sine transform (FST);

$T\in \mathbb{R}^{N\times N}$1

with $T\in \mathbb{R}^{N\times N}$2 ($T\in \mathbb{R}^{N\times N}$3 the $T\in \mathbb{R}^{N\times N}$4 sine-transform matrix) and $T\in \mathbb{R}^{N\times N}$5.

• The right preconditioner inverse is

$T\in \mathbb{R}^{N\times N}$6

• The left preconditioner, under $T\in \mathbb{R}^{N\times N}$7, block-diagonalizes into $T\in \mathbb{R}^{N\times N}$8 independent $T\in \mathbb{R}^{N\times N}$9 blocks,

$\beta$0

each being a sum of diagonal and lower-triangular Toeplitz matrices, and thus invertible (with $\beta$1 also lower-triangular Toeplitz). The inverse is efficiently given by

$\beta$2

4. Krylov Iterative Algorithm and Parallel-in-Time Updates

To solve the preconditioned system $\beta$3, any standard Krylov subspace iterative algorithm (e.g., GMRES, normalized CG for symmetric systems) may be used. The core matrix-vector product $\beta$4 is implemented in three phases:

1. Pre-apply right preconditioner: Transform $\beta$5 by $\beta$6 using joint FST, scaling, and inverse FST: cost $\beta$7.
2. Apply system matrix $\beta$8: Multiply by $\beta$9 and $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$0 efficiently. Cost for spatial multiplication is $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$1 (if $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$2 is sparse); temporal multiplication via Toeplitz structure is $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$3.
3. Pre-apply left preconditioner: Transform by FST, invert the block Toeplitz systems $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$4 in parallel for each spatial index (cost $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$5), and transform back. All $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$6 temporal solves are independent and can be distributed for ideal PinT scaling.

The total cost per matrix-vector product is $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$7, with storage dominated by the $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$8 all-at-once vector and $0< a_{\min} \leq a(x) \leq a_{\max} < \infty$9 entries for first columns of $\beta = \sqrt{a_{\min} a_{\max}}$0. Only $\beta = \sqrt{a_{\min} a_{\max}}$1 storage for $\beta = \sqrt{a_{\min} a_{\max}}$2 is also needed. FFT-based routines enable in-place application of $\beta = \sqrt{a_{\min} a_{\max}}$3.

5. Theoretical Properties and Robustness

Under the assumption

$\beta = \sqrt{a_{\min} a_{\max}}$4

the two-sided preconditioned matrix

$\beta = \sqrt{a_{\min} a_{\max}}$5

has condition number

$\beta = \sqrt{a_{\min} a_{\max}}$6

independent of $\beta = \sqrt{a_{\min} a_{\max}}$7 and $\beta = \sqrt{a_{\min} a_{\max}}$8. This establishes robustness with respect to both problem size and discretization. The proof proceeds via bounding Rayleigh quotients using block-Kronecker inequalities and the positive-definiteness of the temporal matrices.

6. Implementation Complexity and Practical Considerations

Each application of the two-sided preconditioners has computational cost $\beta = \sqrt{a_{\min} a_{\max}}$9. If $P := \beta L_1 \otimes I_N + I_J \otimes T,$0 processors are available, distributing the $P := \beta L_1 \otimes I_N + I_J \otimes T,$1 independent temporal solves in Phase 3 achieves near-ideal PinT scaling. Storage requirements are minimal since only the diagonals of the block systems and the first columns of the Toeplitz inverses are needed.

If the underlying spatial coefficient $P := \beta L_1 \otimes I_N + I_J \otimes T,$2 varies slowly, it suffices to recompute scaled $P := \beta L_1 \otimes I_N + I_J \otimes T,$3 blocks as $P := \beta L_1 \otimes I_N + I_J \otimes T,$4 changes, requiring at most $P := \beta L_1 \otimes I_N + I_J \otimes T,$5 per update. In the special case $P := \beta L_1 \otimes I_N + I_J \otimes T,$6, the matrix $P := \beta L_1 \otimes I_N + I_J \otimes T,$7 coincides with $P := \beta L_1 \otimes I_N + I_J \otimes T,$8, and the two-sided preconditioner becomes a direct solver with complexity $P := \beta L_1 \otimes I_N + I_J \otimes T,$9.

The method exploits the precise algebraic structure of the discretized system (Toeplitz in time, diagonalizable in space) without further approximation and is matrix-size robust in terms of conditioning. These properties distinguish it from classical one-sided and approximate preconditioners for similar equations.

7. Applications and Extensions

Two-sided preconditioners in the manner of Lin–Ng–Zhi (2020) are applicable in the fast iterative or direct solution of discretized non-local evolutionary equations with variable coefficients and weakly singular temporal kernels. They enable efficient PinT solvers for large-scale time-dependent PDEs with guaranteed convergence rates and parallel scalability. The approach is well suited to cases where temporal non-locality and spatially varying diffusion preclude straightforward block-diagonalization or classical (block) circulant preconditioners. A plausible implication is that similar splitting strategies may generalize to other hybrid-structured systems where distinct algebraic properties in different dimensions can be simultaneously exploited for parallelism and conditioning control (Lin et al., 2021).