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): $$P_L := P P_R^{-1} = (\beta L_1)^{1/2}\otimes I_N + (\beta L_1)^{-1/2}\otimes T$$.

The two-sided preconditioned system to solve is

$$P_L^{-1} A P_R^{-1} \widehat U = P_L^{-1} F, \quad U = P_R^{-1} \widehat U,$$

where $F$ is the discretized right-hand side.

2. Derivation and Splitting Strategy

Starting from $P = \beta L_1 \otimes I_N + I_J \otimes T$, the requirement that $P_R$ contains only a fractional power of $L_1$ and no $T$ leads uniquely to the given split. Algebraically,

$$P_L = (\beta L_1)^{1/2}\otimes I_N + (\beta L_1)^{-1/2} \otimes T,$$

by factoring $P = P_L P_R$ with $P_R = (\beta L_1)^{1/2}\otimes I_N$.

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 $L_1$ can be diagonalized by a fast sine transform (FST);

$L_1 = Q \Lambda Q^T,$

with $Q = S_{m_1} \otimes ... \otimes S_{m_d}$ ($S_m$ the $m\times m$ sine-transform matrix) and $$\Lambda = \operatorname{diag}(\lambda_1, ..., \lambda_J)$$.

• The right preconditioner inverse is

$$P_R^{-1} = (Q\otimes I_N)\left[(\beta \Lambda)^{-1/2} \otimes I_N\right](Q\otimes I_N).$$

• The left preconditioner, under $Q\otimes I_N$, block-diagonalizes into $J$ independent $N\times N$ blocks,

$$T_i := (\beta \lambda_i)^{1/2} I_N + (\beta \lambda_i)^{-1/2} T, \quad i=1, ..., J,$$

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

$$P_L^{-1} = (Q\otimes I_N) \operatorname{diag}(T_1^{-1}, ..., T_J^{-1}) (Q \otimes I_N).$$

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

To solve the preconditioned system $P_L^{-1} A P_R^{-1} \widehat{U} = P_L^{-1} F$, any standard Krylov subspace iterative algorithm (e.g., GMRES, normalized CG for symmetric systems) may be used. The core matrix-vector product $w = (P_L^{-1}AP_R^{-1})v$ is implemented in three phases:

1. Pre-apply right preconditioner: Transform $v$ by $P_R^{-1}$ using joint FST, scaling, and inverse FST: cost $O(JN\log J)$.
2. Apply system matrix $A$: Multiply by $L_a\otimes I_N$ and $I_J\otimes T$ efficiently. Cost for spatial multiplication is $O(N\,\mathrm{nnz}(L_a))=O(NJ)$ (if $L_a$ is sparse); temporal multiplication via Toeplitz structure is $O(JN\log N)$.
3. Pre-apply left preconditioner: Transform by FST, invert the block Toeplitz systems $T_i$ in parallel for each spatial index (cost $O(JN\log N)$), and transform back. All $J$ temporal solves are independent and can be distributed for ideal PinT scaling.

The total cost per matrix-vector product is $O(JN(\log J + \log N))$, with storage dominated by the $JN$ all-at-once vector and $JN$ entries for first columns of $T_i^{-1}$. Only $O(J)$ storage for $\Lambda$ is also needed. FFT-based routines enable in-place application of $Q$.

5. Theoretical Properties and Robustness

Under the assumption

$$a_{\min} I_J \preceq L_1^{-1/2} L_a L_1^{-1/2} \preceq a_{\max} I_J,$$

the two-sided preconditioned matrix

$M := P_L^{-1} A P_R^{-1}$

has condition number

$\kappa_2(M) \leq \frac{a_{\max}}{a_{\min}},$

independent of $J$ and $N$. 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 $O(JN\log (JN))$. If $P$ processors are available, distributing the $J$ 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 $a(x)$ varies slowly, it suffices to recompute scaled $T_i$ blocks as $\lambda_i$ changes, requiring at most $O(JN\log N)$ per update. In the special case $a(x) \equiv \mathrm{const}$, the matrix $A$ coincides with $P$, and the two-sided preconditioner becomes a direct solver with complexity $O(JN\log(JN))$.

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).