ADI Splitting Method Overview

• The ADI splitting method is a numerical scheme that decomposes high-dimensional PDEs into one-dimensional or structured subproblems for efficient computation.
• It provides unconditional stability and second-order convergence for parabolic and elliptic PDEs while reducing computational cost and memory usage.
• Recent extensions like the GADI framework, high-order IMEX schemes, and mixed-precision implementations enhance its applicability to large-scale linear systems and control problems.

The alternating direction implicit (ADI) splitting method is a class of operator-splitting numerical schemes for high-dimensional evolutionary partial differential equations (PDEs) and large-scale algebraic systems. By decomposing high-dimensional problems into a series of one-dimensional or structured subsystems, ADI methods provide computationally efficient, unconditionally stable, and often easily parallelizable solvers for parabolic, elliptic, and certain nonlinear or fractional-derivative problems. The method originated as a tool for multidimensional parabolic PDEs but has since been unified, extended, and optimized in general frameworks, including the generalized ADI (GADI), high-order IMEX-Runge–Kutta, and low-rank matrix equation classes. ADI splitting also underpins some of the fastest iterative schemes for large sparse linear or matrix equations through careful choice of operator splitting, parameter tuning, preconditioning, and mixed-precision implementation.

1. Operator Decomposition and Canonical ADI Schemes

The underlying principle of ADI is a splitting of the problem operator into lower-dimensional or structured sub-operators. For a prototypical linear parabolic PDE in $d$ spatial dimensions (such as $\partial_t u = \sum_{j=1}^d L_j u$), ADI methods alternate, within each time step, between implicit solves in each spatial direction while treating other directions explicitly. This yields schemes of the form:

• Douglas (θ) scheme:

$$\text{Stage}~1:~ (I - \theta \Delta t L_1) U^{(1)} = \ldots$$

Continue with $L_2,\ldots,L_d$ in succession per time step, with θ a parameter (commonly $1/2$ for second order). Each substep matrix is typically tridiagonal or block-banded.

• Peaceman–Rachford (PR) scheme:

$$\text{Half step in direction 1, half in direction 2},~ \text{alternating each time step}$$

Crucially, such splitting reduces large multidimensional implicit solves to a sequence of 1D-like subproblems, usually with linear or near-linear computational cost per step and significant memory reduction relative to fully coupled implicit schemes. This separation is at the core of modern ADI applications for multidimensional time-dependent PDEs, including mean curvature flow (Zhou et al., 2023), image osmosis (Calatroni et al., 2017), and option pricing PDEs (Mashayekhi et al., 2023).

2. Generalized Alternating Direction Implicit (GADI) Framework

The GADI framework systematically extends classical ADI methods to arbitrary large sparse linear systems and matrix equations via flexible operator splitting and the introduction of tunable shift ($\alpha$) and relaxation ($\omega$) parameters. Given $A = M + N$, the GADI iteration for $Ax=b$ reads:

$$\begin{cases} (M + \alpha I) x^{(k+\frac{1}{2})} = (\alpha I - N) x^{(k)} + b \ (N + \alpha I) x^{(k+1)} = [N - (1-\omega) \alpha I] x^{(k)} + (2-\omega)\alpha x^{(k+\frac{1}{2})} \end{cases}$$

Sub-problems $(M+\alpha I)$ and $(N+\alpha I)$ are chosen for efficient direct or iterative inversion (e.g., tridiagonal, block, or Kronecker structure). This framework unifies and generalizes Peaceman–Rachford, Douglas–Rachford, Hermitian/skew-Hermitian splitting (HSS), and many matrix-equation-specific ADIs for Sylvester/Lyapunov equations (Ge et al., 24 Dec 2025, Zhang et al., 2024, Jiang et al., 2021, Zhang et al., 2024).

Key properties:

• Convergence holds under broad definiteness and splitting conditions for any $\alpha > 0$ and $0 \leq \omega < 2$.
• The spectral radius and actual convergence are highly sensitive to the choice of splitting and parameters, motivating data-driven parameter prediction (see Section 4 and (Ge et al., 24 Dec 2025, Jiang et al., 2021)).

3. Practical Stability, Convergence, and Complexity

ADI methods are renowned for stability properties markedly superior to explicit or even standard implicit schemes in high dimensions:

• Unconditional stability for linear parabolic or symmetric problems (no CFL time step restriction).
• Mild constraints for certain nonlinear and fractional-order problems, with small explicit terms or mixed derivatives sometimes requiring moderate parameter tuning ($\theta$, $\epsilon$, $\Delta t$).
• Convergence: Second-order in time for Douglas and Peaceman–Rachford; first-order or optimal for high-order IMEX/ADI-GLM/ADI-GARK schemes as long as IO and stage order meet the required coupling conditions (González-Pinto et al., 2021, Sarshar et al., 2019).
• Efficiency: For $d$-dimensional spatial grids of size $N$, the cost to reach fixed terminal time is $O(N^d)$ (per step cost $O(N^{d-1})$), orders of magnitude lower than full implicit solvers, especially in 3D and above (Zhou et al., 2023, Mashayekhi et al., 2023, Los et al., 2019).

4. Parameter Selection, Mixed Precision, and Model Selection

Performance of ADI/GADI schemes is critically determined by shift and relaxation parameters.

• Shift selection ($\alpha$): Poorly chosen $\alpha$ degrades convergence or can destabilize low-precision subsystem solves. Practical recipes include theoretical minimizers (e.g., geometric mean extremal eigenvalues for Hermitian parts) or adaptive, data-driven strategies.
• GPR parameter prediction: A key innovation, particularly for model classes encountered in high-dimensional PDEs and large sparse algebraic systems, is the use of Gaussian Process Regression to predict nearly optimal $\alpha$ as a function of problem features. This allows for effective "one-shot" parameter selection, removing the need for hand-tuning or exhaustive grid search even at massive scale (Ge et al., 24 Dec 2025, Jiang et al., 2021).
• Mixed precision: Recent advances exploit lower arithmetic precision for structured subsystem solves (e.g., BF16/FP32 for $M+\alpha I,N+\alpha I$) while maintaining high precision for residuals and solution updates, thus achieving both throughput and memory savings without sacrificing convergence guarantees. Sharp a priori error analysis in this context ensures that only subsystem convergence rates, not final error floors, depend on low-precision conditioning (Ge et al., 24 Dec 2025).

5. Applications Across PDEs, Linear Algebra, and Control

The ADI paradigm permeates a wide spectrum of computational mathematics:

• Geometric evolution and mean curvature flow: ADI splitting underpins high-dimensional flows by graph decomposition, tangential redistribution, and local 1D evolution, yielding unconditionally stable and mesh quality-preserving solvers (Zhou et al., 2023).
• Imaging and transport: Regularized image denoising, inpainting, and osmosis models leverage ADI to reduce 2D/3D nonlinear or fourth-order diffusions to sequential banded solves (Calatroni et al., 2017, Calatroni et al., 2013).
• Fractional differential equations: For non-local operators, S-ADI and spectral-ADI schemes with careful splitting of tensor Riesz and remainder parts yield FFT-accelerated Toeplitz solvers, retaining unconditional stability and high-order convergence (Sun et al., 2023, Liu et al., 2018).
• Financial mathematics: Multi-factor option pricing, especially with stochastic volatility and interest rate (HCIR-type models), is efficiently handled by ADI, avoiding the curse of dimensionality (Mashayekhi et al., 2023).
• Matrix equations (Lyapunov, Riccati, Sylvester): Low-rank ADI/GADI schemes dominate large-scale stable control and filtering problems, especially with Newton–ADI extensions. Storage and per-iteration cost are dominated by the growing low-rank factor dimensions, but remain order(s) of magnitude cheaper than full matrix methods (Zhang et al., 2024, Zhang et al., 2024).
• Sparse linear systems: GADI and its mixed-precision and data-driven variants represent the fastest iterative solvers for large convection-diffusion–reaction systems, with speedups of $2\times-10\times$ or more over double-precision implementations (Ge et al., 24 Dec 2025).

6. Extensions: High-Order, Unified Formulations, and Modern Frameworks

Recent research generalizes ADI and its variants to arbitrary order, nonlinear equations, and flexible splitting strategies:

• General Linear and GARK frameworks: ADI schemes are now understood as members of IMEX–GARK and partitioned GLM classes, enabling rigorous derivation of coupling conditions for high-order schemes and clear understanding of order reduction phenomena (González-Pinto et al., 2021, Sarshar et al., 2019). High stage order and proper block Butcher tableau structure eliminate splitting-induced order loss.
• High-order IMEX/ADI-GLM: By embedding operator splits into high-stage-order GLM, one achieves the nominal global order in stiff, multistep, or time-dependent-inhomogeneous problems, outperforming standard DIRK-based ADI approaches.
• Isogeometric and fast direct solvers: In spatial discretizations with Kronecker-product structure (isogeometric analysis, tensor grids), ADI splitting supports $O(N)$ complexity per time step via 1D factorizations. This is crucial for hyperbolic wave propagation and elasticity, where unconditional stability is also maintained (Los et al., 2019, Los et al., 2021).

7. Implementation Considerations and Algorithmic Outline

A generic ADI iteration proceeds as follows (for $d$-dimensional spatial problems or a GADI split system):

1. Preprocessing: Identify or assemble sub-operators (e.g., $L_1,\ldots,L_d$ or $(M,N)$) with structure suitable for fast subsystems (tridiagonal/Toeplitz/banded/low-rank).
2. For each time step or outer iteration:
   • Loop over directions or blocks, alternately solving implicit equations along each direction while treating other terms explicitly or with extrapolated data.
   • For GADI, solve two shifted subsystems using fast direct or low-precision iterative algorithms.
   • In matrix equations (Lyapunov/Riccati), update low-rank factors and check residuals for convergence.
3. Update and proceed until terminal time or prescribed accuracy is reached.

Unconditional stability and optimal convergence rates require ensuring that splitting-induced commutators and explicit nonlinear terms remain controlled, with time step choices and regularization aligned to established analysis.

---

For a comprehensive technical overview, see the foundational developments in (Zhou et al., 2023, Ge et al., 24 Dec 2025, Jiang et al., 2021, Sun et al., 2023, Calatroni et al., 2017, Liu et al., 2018, Zhang et al., 2024, Zhang et al., 2024, Calatroni et al., 2013, Los et al., 2019, Los et al., 2021, González-Pinto et al., 2021), and (Sarshar et al., 2019).