Lecture 28: Linear Operator Splitting

A first example of problem-specific algorithms is afforded by consideration of the steady-state diffusion problem

where $u\left(𝒙\right),𝒙\in \Omega \subseteq {ℝ}^d$ can be understood to denote the temperature in an infinitesimal volume positioned at $𝒙$ in $d$-dimensional space, and $f$ is a local rate of heat generation. The above arises from setting the time derivative zero in ${\partial }_{t}u-▵u=f$. Though often stated in this thermal language, the Poisson equation (1) is generally valid for unresolved transport by Brownian motion. The mathematical concept of an “infinitesimal volume” is interpreted as setting some length scale $\ell$ much smaller than the length scale $L$ characterizing the size of the domain $\Omega$. As an example, for heat conduction in a column of water of length $L=1$m the length scale of a quasi-infinitesimal volume can be considered as, say, $\ell =1\mu$m. There is no physical significance to the mathematical limit process $\ell \to 0$ due to the discrete structure of matter, and for all practical purposes setting $\ell =1\mu$m is an acceptable threshold to delimit a phenomenon of interest to behavior that can be neglected. In this case the phenomenon of interest is the average transport of thermal energy in volumes of size greater than $𝒪\left({\ell }^d\right)$. The detailed Brownian motion of the water molecules that occurs on length scales $s\approx 1\mathrm{nm}\ll \ell$ can be neglected and is said to be unresolved. The only observable effect of this motion is that temperature gradients lead to a heat flux as described by $𝒇\left(u\right)=-\alpha \nabla u$ (Fourier's law). The same equation (1) arises in epidemiology when $\ell =10$m, an average separation between an infected and susceptible individual, $L=10\mathrm{km}$, the size of a large city, and $u$ is reinterpreted as the percentage of infected individuals in the population. Again, the detailed Brownian steps of size $s\approx 10\mathrm{cm}\ll \ell$ taken by individuals can be neglected.

1.Poisson equation discretization

Understanding the underlying unresolved Brownian motion is useful in constructing numerical solutions. For $f=0$, (1) becomes ${\nabla }^{2}u=0,$which states that there is no net heat flux in an infinitesimal volume, $div\cdot (\alpha \nabla u)=0$, colloquially: “what flows in on one side goes out on the other”. A function that satisfies Laplace's equation ${\nabla }^{2}u=0$ is said to be harmonic. For $d=1$, the ordinary differential equation ${\partial }_x^{2}u=0$ is obtained with solution $u\left(x\right)=a+bx$, and boundary conditions $u_0$ at $x=0$ and $u_1$ at $x=1$ gives $u\left(x\right)=u_0+(u_1-u_0)x$. A temperature difference $u_0\ne u_1$ at the boundaries induces diffusive fluxes that lead to the mean value $u(1/2)=(u_0+u_1)/2$ at the midpoint. An analogous statement is made for $d>1$ starting from Green's formula on ball $B$ with spherical boundary $S$ of radius $R$

for scalar functions $u,v$. For $d=3$ and $v=1/r$ the mean value theorem

is obtained, which states that the value of a harmonic function is the average of the values on a surrounding sphere. For $d=2$ the analogous statement is

Midpoint quadrature of (3) over four subintervals gives

Apply the above on a grid covering some arbitrary domain $\Omega$ with boundary $\Sigma =\partial \Omega$ to obtain

Complications arise for general boundaries $\Sigma$, but for a square domain $\Omega =[0,1]\times [0,1]$ grid points $(x_i=ih,y_j=jh)$ align with the boundary, $h=1/(n+1)$. Instead of two indices, one for each spatial direction, organize the grid values $u$ through a single index $k=(j-1)n+i$, with $u_k$ denoting the value at the interior points $(x_i,y_j)$. The vector of interior grid values $𝒖=\left[\begin{array}{lll}{u}_1& \dots & u_m\end{array}\right]$ has $m=n^2$ components, and the mean value theorem leads to the linear system

∴ 

n=6; m=n^2; A=zeros(m,m);

∴ 

for k=1:m
  A[k,k] = -4
  if k>1 A[k,k-1]=1 end
  if k<m A[k,k+1]=1 end
  if k>n A[k,k-n]=1 end
  if k<m-n+1 A[k,k+n]=1 end
end;

∴ 

spy(A);

∴ 

where $𝑨$ has a regular sparsity pattern induced by the uniform spacing of the grid

The vector $𝒃$ contains the boundary values, for example $b_1=u_{1,0}+u_{0,1}$.

Figure 1. Left: Mean value theorem leads to $u_0=(u_1+u_2+u_3+u_4)/4$. Middle: Five-point finite difference stencil for Laplace operator. Right: Structure of $𝑨$ matrix resulting from discretization of Laplace operator.

The system (6) was obtained by discretization of the mean-value integral (3). The same linear system is also obtained by discretization of the differential equation

where indices denote differentiation. The minus sign arises from compatibility with the unsteady form of the heat equation ${\partial }_{t}u-{\nabla }^{2}u=f$. A centered finite difference approximation of the derivatives on the uniform grid leads to

and (6) is recovered. For the Poisson equation $f\ne 0$ the right hand side changes to

with $f_k=f(x_i,y_j)$. It is often the case that the same discrete system arises from both the differential and the integral formulation of a conservation law on a uniform grid.

2.Matrix splitting iteration

The above discussion of the underlying physics of the Poisson equation productively guides construction of numerical solution procedures. Solving the linear system (7) by general factorizations such as $𝑨=𝑳𝑼$ or $𝑨=𝑸𝑹$ is costly in terms of memory usage since the sparisty pattern is lost. For the uniform grid and square domain considered above the matrix $𝑨$ need not be explicitly stored at all since $a_{i,i}=4$, and $a_{i,j}=-1$ when $|i-j|=1$ or $|i-j|=n$. The discrete mean value theorem (5) suggests that some approximation ${𝒖}^{\left(l\right)}$ can be improved by the iteration

The above is known as Jacobi iteration, and can be stated in matrix form by expressing $𝑨$ as

with $𝑳,𝑼$ containing non-zero components of $𝑨$ below, above the diagonal, and $𝑫$ containing the diagonal of $𝑨$. In contrast to the multiplicative decompositions considered up to now, the $QR$, $LU$, SVD, eigen or Schur decompositions, the decomposition (9) is now additive. Note that in (9) $𝑳,𝑼$ are strictly lower, upper diagonal matrices with zeros on the diagonal in contrast to the notation for the standard $LU$ factorization algorithm. Recall that the utility of matrix multiplication was associated with the representation of linear mapping composition. Additive decompositions such as (9) generally are useful when separating different aspects of a physical process, and are a simple example of operator splitting. For the discrete Poisson system Jacobi iteration can be expressed as

Several variants of the idea can be pursued. The matrix splitting (9) is useful in theoretical convergence analysis, but implementations directly use (8) within loops over the $(i,j)$ indices. Updated values can be immediately utilized leading to either of the following iterations

depending on loop organization. These are known as Gauss-Seidel iterations. Convergence might accelerated by extrapolation,

where the new iteration is continued by factor $\omega$ in the direction of the Gauss-Seidel update. When $\omega >1$ this is known as successive over-relaxation (SOR) and goes further in the Gauss-Seidel direction. Choosing $0<\omega <1$ leads to successive under-relaxation.

3.Convergence analysis

Turning now from algorithm construction to analysis of its behavior, simplify notation by letting ${𝒖}_k$ denote the current iterate. The previous notation ${𝒖}^{\left(l\right)}$ was convenient since individual components were referenced as in $u_{i,j}^{\left(l\right)}$, but convergence analysis is determined by the properties of the operator splitting and not of the current iterate. Introduce the error ${𝜹}_k$ at iteration $k$ as the difference between the exact solution $𝒖$ and the current iterate ${𝒖}_k$, ${𝜹}_k=𝒖-{𝒖}_k$. Also introduce the residual ${𝒓}_k=𝒄-𝑨{𝒖}_k=𝑨{𝜹}_k$, and the correction to the current iterate ${𝒆}_k={𝒖}_{k+1}-{𝒖}_k$

The above methods can be formulated as a residual correction algorithm through the steps:

1. residual computation, ${𝒓}_k=𝒄-𝑨{𝒖}_k$

2. correction computation, ${𝒆}_k=𝑩{𝒓}_k$

3. approximation update, ${𝒖}_{k+1}={𝒖}_k+{𝒆}_k$

When $𝑩={𝑨}^{-1}$ the exact solution is recovered in one step

Iterative methods use some approximation of the (unknown) inverse, $𝑩\cong {𝑨}^{-1}$. Jacobi iteration uses $𝑩={𝑫}^{-1}$ since

recovering (10). Table 1 shows several common choices for $𝑩$. Two key aspects govern the choice of the inverse approximant:

1. Computational efficiency stated as a requirement that each iteration cost either $𝒪\left(m\right)$ or $𝒪(m\log m)$ operations;

2. Capturing the essential aspects of $𝑨$.

The iteration converges to the solution if $||{𝜹}_k||\to 0$ for increasing $k$. The error at iteration $k+1$ is expressed as

The repeated matrix multiplication indicates that the eigenstructure of the iteration matrix $𝑴=𝑰-𝑩𝑨$ determines iteration convergence. Indeed the above is simply power iteration for $𝑴$ and can be expected to converge as

with $({\mu }_1,{𝒒}_1)$ the eigenpair that corresponds to the largest eigenvalue in absolute value, $𝑴{𝒒}_1={\mu }_1{𝒒}_1$, known as the spectral radius of $𝑴$, denoted by $\rho \left(𝑴\right).$ Clearly, the above iterations will exhibit linear order of convergence when $\rho \left(𝑴\right)<1$. The rate of convergence $s_k$ at iteration $k$ is estimated by the Rayleigh quotient

and is monitored in implementations of iterative methods, and determined by the eigenvalues $\lambda =1-\mu$ of $𝑩𝑨$.

The eigenstructure of $𝑩𝑨$ is difficult for arbitrary matrices $𝑨$, but can be carried out when $𝑨$ has special structure induced by known physical phenomena. The relation between analytical and numerical formulations plays an essential role in convergence analysis. The diffusion equation (5) leads to a symmetric matrix $𝑨$, $𝑨={𝑨}^T$ due to two aspects:

1. the chosen discretization is symmetric using centered finite differences;

2. the operator itself exhibits symmetry.

Insight into the above two aspects is most readily gained from the one-dimensional case $-u_{xx}=f$ with homogeneous Dirichlet boundary conditions $u\left(0\right)=u\left(1\right)=0$. The linear system $𝑨𝒖=𝒇$ obtained from the centered finite difference discretization

has a symmetric tridiagonal system matrix $𝑨=\mathrm{diag}\left(\left[\begin{array}{lll}-1& 2& -1\end{array}\right]\right)$. The ${𝑨}^T=𝑨$ symmetry can be expressed through scalar products in a way that generalizes to differential operators. Recall that a real-valued scalar product $(𝒖,𝒗)$ must satisfy symmetry $(𝒗,𝒖)=(𝒖,𝒗)$. For $𝒖,𝒗\in {ℝ}^m$ the standard inner product $(𝒖,𝒗)={𝒖}^T𝒗$ has this property. Consider the action of the operator $𝑨\in {ℝ}^{m\times m}$ on the two terms. If $(𝑨𝒖,𝒗)=(𝒖,𝑨𝒗)$ the operator $𝑨$ is said to be symmetric. For the inner product

and the two expressions are equal if $𝑨={𝑨}^T$. The same approach extends to the $d=1$ diffusion operator $L=-{\partial }_{xx}$ using the scalar product

Applying integration by parts

For homogeneous boundary conditions $f\left(0\right)=f\left(1\right)=g\left(0\right)=g\left(1\right)=0$, the symmetry condition is satisfied. Note that symmetry involves both the operator anId the boundary conditions of the problem. For $d=2$ the scalar product

is defined on the unit square $\Omega =[0,1]\times [0,1]$, and for homogeneous Dirichlet boundary conditions two applications of Green's formula leads to $({\nabla }^{2}u,v)=(u,{\nabla }^{2}v)$, and the Laplace operator is symmetric.

For $d=1$, $x_j=jh$, $h=1/(m+1)$, the eigenvalues ${\nu }_k$of $𝑨=\mathrm{diag}\left(\left[\begin{array}{ll}1-2& 1\end{array}\right]\right)$ are inferred from those of the $-{\partial }_{xx}$ operator with homogeneous boundary conditions at $x=0$, $x=1$

Positing that an eigenvector $𝒒$ of $𝑨$ is the discretization of the continuum eigenfunction leads to

hypothesis that is verified by the calculation of component $j$ of $𝑨𝒒$

thereby obtaining the eigenvalue

which recovers the analytical eigenvalue in the $h\to 0$ limit

For Jacobi $𝑩={𝑫}^{-1}$, so the eigenvalues of $𝑩𝑨$ are

The eigenvalues of $𝑴=𝑰-𝑩𝑨$ are therefore

Replacing $h=1/(m+1)$ the largest eigenvalue is obtained for $\kappa =1$

For large $m$, ${\mu }_1⪅1$, and slow convergence is expected as verified in the numerical experiment from Fig. 2.

$\circ$

Figure 2. Convergence of Jacobi iteration. Blue: exact solution. Orange, green, red: iterates after 1000, 2000, 3000 iterations.

∴	function ResidualCorrection(A,B,c,u,maxIter) m = length(u0) for k=1:maxIter r = c - A(u) e = B(r) u = u + e end return u end;

∴	function ALaplace1D(u) m = length(u); v = zeros(m,1) v[1] = 2*u[1] - u[2] for k=2:m-1 v[k] = -u[k-1]+2*u[k]-u[k+1] end v[m] = -u[m-1] + 2*u[m] return v end;

∴	function BJacobiLaplace1D(r) beta = 0.5 return beta*r end;

∴	m=100; h=1.0/(m+1); x=(1:m)*h; c = (pi*h)^2*sin.(pi*x);

∴	u0=zeros(m,1); uex=sin.(pi*x);

∴	u1=ResidualCorrection(ALaplace1D,BJacobiLaplace1D,c,u0,1000);

∴	u2=ResidualCorrection(ALaplace1D,BJacobiLaplace1D,c,u0,2000);

∴	u3=ResidualCorrection(ALaplace1D,BJacobiLaplace1D,c,u0,3000);

∴	clf(); plot(x,uex,x,u1,x,u2,x,u3);

∴	grid("on"); xlabel(L"$x$"); ylabel(L"$u(x)$");

∴	title(L"Jacobi iteration convergence for $u_{x x}=\pi^2 sin(\pi x)$");

∴

Jacobi iteration and its variants have limited practical utility by comparison to other iterative procedures due to slow convergence, but the concept of operator spliting has wide applicability. A more consequential example of operator splitting is to consider the advection diffusion equation

which can be interpreted as stating that the time evolution of $q$ is due to the effect of a diffusion operator $𝒜=\alpha \Delta$ and an advection operator $ℬ=-𝒗\cdot \nabla$

Suppose both operators are discretized leading to matrices $𝑨,𝑩$ and the discrete system

with initial condition $𝒒(𝒙,t=0)={𝒒}_0$. Advancing the solution by a time step $\Delta t$ can be written as

and can be separated into two stages

The quantity $\stackrel{\sim }{𝒒}=e^{\Delta t𝑩}{𝒒}^{\left(l\right)}$ captures advection effects and ${𝒒}^{(l+1)}=e^{\Delta t𝑨}\stackrel{\sim }{𝒒}$ is the diffusion correction. Since advection and diffusion are markedly different physical effects describing resolved versus unresolved transport, it can be expected that matrices $𝑨,𝑩$ have different properties that require specific solution procedures.

By contrast, the Jacobi iteration splitting $𝑨=𝑳+𝑫+𝑼$ does not separate physical effects and simply is suggested by the sparsity of $𝑨$ and computational efficiency per iteration. For example, the forward Gauss-Seidel iteration with $𝑩={(𝑫+𝑳)}^{-1}$ leads to

The matrix $𝑫+𝑳$ is (non-strictly) lower-triangular and (14) is easily solved by forward substitution. The implementation is very simple to express while preserving two-dimensional indexing.

In the above implementation a single memory space is used for $𝒖$ with new values taking place of the old, leading to just four floating point additions and one multiplication per grid point. Since the algorithm essentially takes the current average of neighboring values, it is also known as a relaxation method, smoothing out spatial variations in $𝒖$.

Though such simple implementation is desirable, the non-physical splitting $𝑨=𝑳+𝑫+𝑼$ and associated slow convergence usually outweighs ease of coding effort and suggests looking for alternative approaches. The only scenarios where such simply implemented iterations find practical use is in parallel execution and as a preliminary modification of the system prior to use of some other algorithm.