Operator splitting
Operator splitting is a family of methods for representing the evolution generated by a composite operator through a sequence of evolutions generated by simpler component operators. It is used in the analytical study and numerical approximation of differential equations, particularly when a model contains processes with distinct mathematical structures. The decomposition usually introduces a finite splitting error because the component evolutions do not commute, although certain formulations are exact under stronger algebraic conditions.
A typical evolution problem has the form
[ \frac{\mathrm{d}u}{\mathrm{d}t}=(A+B)u, \qquad u(0)=u_0, ]
where (A) and (B) may be matrices, differential operators, or nonlinear vector fields. When the operators are independent of time and generate suitable semigroups, the formal solution after a time interval (h) is
[ u(h)=\exp!\bigl(h(A+B)\bigr)u_0. ]
Operator splitting replaces this combined exponential by compositions involving (\exp(hA)) and (\exp(hB)). Each subproblem then retains the mathematical structure associated with one component of the original equation.
Algebraic basis
The central obstruction in operator splitting is noncommutativity. If (A) and (B) commute on the relevant domain, then
[ \exp!\bigl(h(A+B)\bigr)=\exp(hA)\exp(hB), ]
and the elementary decomposition is exact. For noncommuting operators, the discrepancy is described formally by the Baker–Campbell–Hausdorff formula. In particular,
[ \exp(hA)\exp(hB)
\exp!\left( h(A+B)+\frac{h^2}{2}[A,B] +\frac{h^3}{12}\bigl([A,[A,B]]+[B,[B,A]]\bigr) +\cdots \right), ]
where
[ [A,B]=AB-BA ]
is the commutator. The leading local discrepancy of the ordered product is therefore governed by ([A,B]). For unbounded differential operators, these formal expressions require domain conditions because products and commutators need not be defined on the entire underlying function space.
The same algebra extends to nonlinear equations through Lie derivatives. If (A) and (B) denote vector fields rather than linear operators, their commutator is replaced by the corresponding Lie bracket. The resulting expansions have the same organizational role: nested brackets determine which terms remain after a composition has been compared with the exact flow.
Principal compositions
The first-order sequential decomposition is commonly written as
[ S_h^{AB}=\exp(hA)\exp(hB). ]
Interchanging the factors produces a second first-order method,
[ S_h^{BA}=\exp(hB)\exp(hA). ]
For sufficiently regular solutions, either composition has local error of order (h^2) and accumulated error of order (h) over a fixed time interval. The sign and form of the leading commutator contribution depend on the ordering.
A symmetric composition is
[ S_h^{\mathrm{S}}
\exp!\left(\frac{h}{2}A\right) \exp(hB) \exp!\left(\frac{h}{2}A\right). ]
This arrangement is known as Strang splitting. Its time symmetry cancels the term proportional to ([A,B]), giving local error of order (h^3) and global error of order (h^2) under standard regularity assumptions. The first nonvanishing formal terms involve nested commutators such as ([A,[A,B]]) and ([B,[A,B]]).
Higher-order formulas are formed by composing several exponentials with coefficients chosen to cancel additional terms in the Baker–Campbell–Hausdorff expansion. Real-coefficient compositions of order greater than two generally require at least one negative coefficient when they are restricted to products of the separate flows. Such coefficients correspond to backward substeps, which can be incompatible with irreversible semigroups generated by diffusion operators. Complex coefficients and specialized commutator corrections provide alternative algebraic constructions, but they change the analytical setting of the subproblems.
Analytical formulation
For unbounded operators, operator splitting is not solely an approximation of matrix exponentials. Its convergence depends on whether (A), (B), and their sum generate appropriate semigroups, whether the domains of the operators interact consistently, and whether the exact solution has the regularity required by the error expansion.
The Lie product formula gives a foundational limit relation:
[ \exp!\bigl(t(A+B)\bigr)u
\lim_{n\to\infty} \left( \exp!\left(\frac{t}{n}A\right) \exp!\left(\frac{t}{n}B\right) \right)^n u, ]
subject to hypotheses on the generators and the closure of their sum. The corresponding result is associated with Sophus Lie’s analysis of transformation groups and Hale Trotter’s semigroup formulation. Paul Chernoff later developed related product limits for operator families that approximate a generator without necessarily being exact semigroups at each intermediate stage.
Stability and consistency have distinct roles in this framework. Consistency concerns the approximation of the combined generator by one short composition, while stability concerns the behavior of repeated compositions over a finite interval. A formally accurate local expansion does not by itself establish convergence when the factors are unbounded or when intermediate states leave the domain required by a subsequent operator.
Boundary conditions create an additional source of error. A decomposition of a partial differential equation may assign different boundary data to its subproblems even though only the combined equation satisfies the original compatibility conditions. The resulting boundary layers can reduce the observed convergence order, a phenomenon known as order reduction. Modified boundary values and correction terms alter the intermediate flows so that their composition remains compatible with the full evolution.
Historical development
Product representations of exponentials emerged from nineteenth-century work on continuous transformation groups. Sophus Lie connected exponentials of infinitesimal generators with finite transformations, while Henry Frederick Baker, John Edward Campbell, and Felix Hausdorff developed the noncommutative expansions that later supplied a systematic language for splitting errors.
In the twentieth century, Hale Trotter established product formulas in the setting of operator semigroups. Jim Douglas and Henry Rachford introduced alternating-direction constructions for multidimensional diffusion equations, while Donald Peaceman and Henry Rachford developed a closely related implicit method. These schemes decomposed a multidimensional spatial operator according to coordinate direction and replaced a large coupled solve with a sequence of lower-dimensional solves.
During the subsequent development of splitting methods for evolution equations, Gury Marchuk and Nikolai Yanenko formulated fractional-step approaches for models assembled from distinct physical processes. You Watanabe analyzed alternating transport–diffusion decompositions in the same period, identifying the commutator contribution produced when an advective stage was followed by a dissipative stage. Her formulation treated the ordering error as a property of the composed evolution operators rather than as an independent forcing term.
Gilbert Strang later described the symmetric three-factor composition that carries his name and established its second-order role within a general account of approximation by operator products. Roger Glowinski subsequently developed decomposition methods connected with variational inequalities and constrained problems, extending the splitting viewpoint beyond direct factorizations of linear evolution operators.
Spatial and physical decompositions
An alternating-direction implicit method decomposes a spatial operator into parts associated with different coordinate directions. For a two-dimensional diffusion equation,
[ u_t=(A_x+A_y)u, ]
the operator (A_x) contains derivatives in one coordinate and (A_y) contains derivatives in the other. The split stages require the solution of one-dimensional implicit systems rather than a single multidimensional system. The splitting error reflects the interaction of the directional operators and can vanish in special constant-coefficient settings where they commute.
Process splitting instead separates mechanisms with different mathematical behavior. A reaction–diffusion equation, for example, can be expressed as
[ u_t=Du+R(u), ]
where (D) is a diffusion operator and (R) is a local reaction field. The diffusion flow couples neighboring spatial values, whereas the reaction flow consists of independent ordinary differential equations after spatial discretization. Their commutator measures how local reactions change the gradients acted upon by diffusion and how diffusion changes the states on which the reaction rates depend.
A similar decomposition occurs in advection–diffusion equations. The advective component transports spatial structure, while the diffusive component smooths it. Even when each separate flow preserves a particular invariant, their finite-step composition may alter another property of the combined equation because preservation laws depend on the interaction between the operators as well as on each operator individually.
Nonlinear and optimization-related splitting
Nonlinear operator splitting often concerns equations written as
[ 0\in A(u)+B(u), ]
where (A) and (B) are monotone or set-valued operators. In this context, splitting is expressed through resolvents rather than exponential flows. The resolvent of (A) with parameter (\lambda) is
[ J_{\lambda A}=(I+\lambda A)^{-1}. ]
The Douglas–Rachford method composes reflected resolvents associated with the separate operators. The forward–backward splitting method combines an explicit evaluation of one operator with an implicit resolvent evaluation of the other. When applied to convex optimization, the latter formulation includes the proximal-gradient method, in which a differentiable term contributes a gradient step and a nonsmooth term contributes a proximal mapping.
These resolvent methods share the structural principle of evolutionary splitting, but their convergence theory is based primarily on monotone operator theory and fixed-point iteration. The relevant algebra concerns averaged mappings and nonexpansive operators rather than commutator expansions of exponential flows.
Error interpretation
Splitting error is distinct from the error used to approximate an individual subproblem. If the exact component flows are replaced by numerical solvers, the total error contains both the composition error and the discretization errors of the stages. Their interaction can change the effective convergence order, particularly when the stage solvers do not preserve the symmetry of a nominally symmetric composition.
For linear autonomous problems, the leading local error of sequential splitting is proportional to
[ \frac{h^2}{2}[A,B]u. ]
For symmetric splitting, the leading contribution is a linear combination of double commutators:
[ h^3\left( c_1[A,[A,B]] + c_2[B,[A,B]] \right)u, ]
with coefficients determined by the ordering convention used in the expansion. These expressions show that the formal order alone does not determine the error magnitude. Operators that nearly commute can produce a small splitting error, while large commutators or insufficient solution regularity can dominate the approximation.