I blogged about the Douglas-Rachford method before and in this post I’d like to dig a bit into the history of the method.

As the name suggests, the method has its roots in a paper by Douglas and Rachford and the paper is

Douglas, Jim, Jr., and Henry H. Rachford Jr., “On the numerical solution of heat conduction problems in two and three space variables.” Transactions of the American mathematical Society 82.2 (1956): 421-439.

At first glance, the title does not suggest that the paper may be related to monotone inclusions and if you read the paper you’ll not find any monotone operator mentioned. So let’s start and look at Douglas and Rachford’s paper.

**1. Solving the heat equation numerically **

So let us see, what they were after and how this is related to what is known as Douglas-Rachford splitting method today.

Indeed, Douglas and Rachford wanted to solve the instationary heat equation

with Dirichlet boundary conditions (they also considered three dimensions, but let us skip that here). They considered a rectangular grid and a very simple finite difference approximation of the second derivatives, i.e.

(with modifications at the boundary to accomodate the boundary conditions). To ease notation, we abbreviate the difference quotients as operators (actually, also matrices) that act for a fixed time step

With this notation, our problem is to solve

in time.

Then they give the following iteration:

(plus boundary conditions which I’d like to swipe under the rug here). If we eliminate from the first equation using the second we get

This is a kind of implicit Euler method with an additional small term . From a numerical point of it has one advantage over the implicit Euler method: As equations (1) and (2) show, one does not need to invert in every iteration, but only and . Remember, this was in 1950s, and solving large linear equations was a much bigger problem than it is today. In this specific case of the heat equation, the operators and are in fact tridiagonal, and hence, solving with and can be done by Gaussian elimination without any fill-in in linear time (read Thomas algorithm). This is a huge time saver when compared to solving with which has a fairly large bandwidth (no matter how you reorder).

How do they prove convergence of the method? They don’t since they wanted to solve a parabolic PDE. They were after stability of the scheme, and this can be done by analyzing the eigenvalues of the iteration. Since the matrices and are well understood, they were able to write down the eigenfunctions of the operator associated to iteration (3) explicitly and since the finite difference approximation is well understood, they were able to prove approximation properties. Note that the method can also be seen, as a means to calculate the steady state of the heat equation.

We reformulate the iteration (3) further to see how is actually derived from : We obtain

**2. What about monotone inclusions? **

What has the previous section to do with solving monotone inclusions? A monotone inclusion is

with a *monotone operator*, that is, a multivalued mapping from a Hilbert space to (subsets of) itself such that for all and and it holds that

We are going to restrict ourselves to real Hilbert spaces here. Note that linear operators are monotone if they are positive semi-definite and further note that monotone linear operators need not to be symmetric. A general approach to the solution of monotone inclusions are so-called *splitting methods*. There one splits additively as a sum of two other monotone operators. Then one tries to use the so-called *resolvents* of and , namely

to obtain a numerical method. By the way, the resolvent of a monotone operator always exists and is single valued (to be honest, one needs a regularity assumption here, namely one need *maximal monotone operators*, but we will not deal with this issue here).

The two operators and from the previous section are not monotone, but and are, so the equation is a special case of a montone inclusion. To work with monotone operators we rename

and write the iteration~(4) in terms of monotone operators as

i.e.

Using and we rewrite this in terms of resolvents as

This is not really applicable to a general monotone inclusion since there and may be multi-valued, i.e. the term is not well defined (the iteration may be used as is for splittings where is monotone and single valued, though).

But what to do, when both and and are multivaled? The trick is, to introduce a new variable . Plugging this in throughout leads to

We cancel the outer and use to get

and here we go: This is exactly what is known as Douglas-Rachford method (see the last version of the iteration in my previous post). Note that it is not that converges to a solution, but , so it is convenient to write the iteration in the two variables

The observation, that these splitting method that Douglas and Rachford devised for linear problems has a kind of much wider applicability is due to Lions and Mercier and the paper is

Lions, Pierre-Louis, and Bertrand Mercier. “Splitting algorithms for the sum of two nonlinear operators.” SIAM Journal on Numerical Analysis 16.6 (1979): 964-979.

Other, much older, splitting methods for linear systems, such as the Jacobi method, the Gauss-Seidel method used different properties of the matrices such as the diagonal of the matrix or the upper and lower triangluar parts and as such, do not generalize easily to the case of operators on a Hilbert space.