ADI splitting schemes for a fourth-order nonlinear partial differential equation from image processing

We present directional operator splitting schemes for the numerical solution of a fourth-order, nonlinear partial differential evolution equation which arises in image processing. This equation constitutes the $H^{-1}$-gradient flow of the total variation and represents a prototype of higher-order equations of similar type which are popular in imaging for denoising, deblurring and inpainting problems. The efficient numerical solution of this equation is very challenging due to the stiffness of most numerical schemes. We show that the combination of directional splitting schemes with implicit time-stepping provides a stable and computationally cheap numerical realisation of the equation.


Introduction
In this paper we propose directional operator splitting methods for the numerical solution of fourth-order nonlinear partial differential equations of the type (1.1) u t = ∇ · (h(u)∇q) q ∈ ∂E(u) in Ω × (0, ∞), where E is the total variation (TV) seminorm (see, for instance, [1]) (1.2) E(u) := |Du|(Ω) = sup p∈C ∞ 0 (Ω;R 2 ), p ≤1 Ω u∇ · p dx. and ∂Eu is the subdifferential of E in u (see [1]). Here Ω ⊂ R 2 is open and bounded with Lipschitz boundary, h : R → R, and u 0 a sufficiently regular initial condition. In our considerations equation (1.1) is endowed with periodic boundary conditions. The elements q of the subdifferential ∂E have the property that, if q ∈ ∂E(u), then (see [52,Proposition 4.1]): The above characterization shows the fourth-differential order and the strong nonlinearity in equation (1.1).
Equations of the form (1.1) typically arise in the study of tension driven flow of thin liquid films (see [41,43]), and have recently found applications in image processing (see, e.g., [44,8,49]). In this paper we limit ourselves to the consideration of one prototype of (1.1), that is we consider the special case where h(u) ≡ 1. The equation we obtain in that case reads (1.4) u t = ∆q, q ∈ ∂E(u), in Ω × (0, ∞), which constitutes a gradient flow of the total variation seminorm in the space H −1 .
In image processing problems -such as image denoising and inpainting -nonlinear partial differential equations of the type (1.4) became popular due to their ability of preserving edges in the process of reconstruction (see [46,11,47] for instance). In the latter works the authors typically deal with second-order partial differential equations, that is L 2 -gradient flows of the total variation functional. Instead, in [38,44,40,33,19] an approach like (1.4) is proposed for image denoising and image decomposition. When applied to image denoising, a given noisy image u 0 is decomposed into its piecewise smooth part u = u(T ), solution of (1.4) at time T , and its oscillatory, noisy part n, i.e., u 0 = u + n. Similarly, in image decomposition the piecewise smooth part u represents the structure/cartoon part of the image, and the oscillatory part n its texture part. The advantage of using an H −1 subgradient of the total variation rather than an L 2 subgradient is that (1.4) shows better separability of the data into oscillatory and piecewise constant components and it reduces unwanted artifacts of the total variation filter such as staircasing, cf. [44,33,36].
In [8,48] equation (1.4) is used for image inpainting, i.e. the problem of filling the missing parts of damaged images using the information from the surrounding areas. Here, the higher-order subgradient is necessary for improving upon the ability of total variationtype inpainting to smoothly connect image structures also over large distances, cf. for instance [14,15,16,12,37,25]. In [8] the inpainted image u is reconstructed from a given image f ∈ L 2 (Ω) which is damaged inside the so-called inpainting domain D ⊂ Ω by evolving it via the flow (1.5) u t = ∆q, q ∈ ∂T V (u), in D, where T V (u) := |Du|(Ω) if u ∞ ≤ 1 a.e. in Ω, +∞ otherwise.
Here, the L ∞ -bound on solutions of (1.5) is a technical assumption needed for the existence analysis in [7] and does not present a restriction when dealing with grayvalue images whose values are always bounded within a fixed interval.
If the function h in (1.1) is the identity function h(u) = u we obtain the following equation (1.6) u t = ∇ · (u∇q), q ∈ ∂E(u), in Ω × (0, ∞), This equation can be formally derived as the L 2 -Wasserstein gradient flow of the total variation E in (1.2) and has been proposed in [7] for density estimation and smoothing. Equation (1.6) has been further investigated in [20], where the authors numerically studied the scale space properties and high-contrasting effects of the equation by exploiting a numerical scheme similar to the ones we consider in this paper. In [6] the authors consider such equation for the problem of denoising, solving an alternative primal-dual formulation of (1.6) with a Newton-type method.
From a computational point of view, finding numerical schemes that solve equations of the form (1.1) is a challenging problem. Dealing with an evolutionary nonlinear fourthorder partial differential equation, we would like to find an efficient and easily applicable method. Note that a naive explicit discretisation in time may restrict the choice of the time-step size ∆t up to an order ∆x 4 , where ∆x denotes the step size of the spatial grid, compare [51]. Moreover, the strong nonlinearity of subgradients of the total variation adds additionally constraints to the stability condition of a discrete time stepping scheme, compare [13,19]. In particular, when approximating the subgradient of the total variation by regularising it -either by a square root ε-regularisation or by a regularisation of its dual formulation -the size of the regularisation parameter encodes the accuracy of this approximation, that is the strength of the nonlinearity in the approximated subgradient. The presence of this nonlinearity together with the fourth differential order of the equation then may result in restrictive stability conditions on numerical time stepping for small values of this regularisation parameter. Having these complications in mind, some numerical methods have been proposed to solve equations like (1.4). Lieu and Vese [33] proposed a numerical method to solve TV-H −1 denoising/decomposition by using the Fourier representation of the H −1 norm on the whole R d , d ≥ 1. This leads to a second-order PDE defined in Fourier space. In [23] and [24] the authors propose an algorithm using a finite element method to solve such equation, while in [2,50] a dual approach similar to the one described in [10] is presented with interesting applications both to denoising and inpainting. In [48] the authors present results of convergence and unconditional stability for a particular numerical splitting method solving equation (1.4), called convexity splitting. Therein, the equation is modeled as the gradient flow in H −1 of the difference of two convex energies. The result is the presence of a linear diffusion term in the numerical scheme which balances the unstable behaviours coming from the nonlinear terms.
In this work we are interested in performing a directional operator splitting for the numerical solution of equation (1.4), i.e. (1.5) for image inpainting. In particular, we consider alternating directional splitting (ADI) methods that reach back to [45] and have been further developed in, e.g, [32,31,4,28,29,55]. We discuss three different splitting methods and their stability properties. Our investigation culminates in the proposal of a fully-implicit, additive multiplicative operator splitting scheme for a slightly modified form of equation (1.4). In our numerical experiments we obtain stable numerical solutions for arbitrary values of the regularization parameter, suggesting unconditional stability of the method. This allows to use fewer iterations (larger time steps) while each iteration is still computationally cheap due to the directional splitting. We want to stress, that in the existing literature so far, most methods have been considered for the solution of second-order linear and nonlinear PDEs. Up to our knowledge the following discussion is the first one that carries out an ADI splitting for a higher-order total variation flow, which is both highly nonlinear and of fourth differential order.
Outline of the paper. After the preliminary Section 2, in Section 3 we will review the general features of the main ADI methods. We will apply them to equation (1.4) in Section 4. In order to come up with new and efficient ADI methods solving numerically (1.4), we exploit at first (see Section 4.1) an auxiliary, linear, fourth-order equation. We next present in Section 4.2 our ideas applied to the nonlinear cases we are interested on. Our strategy is as follows: after linearising the equation in a suitable way, we expand and discretise the appearing differential operators and we build up the ADI schemes where we alternatively apply differential operators acting just along either the x-or the y-direction. According to the different choices of linearisation performed, some mixed derivative terms may appear as well. They need to be controlled properly in order to get stable results. Stability properties appear to be related also with the regularisation of the term |∇u| appearing in both equations. In Sections 4.2.1 and 4.2.2 we will give more details on that, presenting in Section 4.2.4 a modified stable ADI scheme dealing implicitly just with pure derivatives. Section 4.3 presents a quite different approach to solve the so-called primaldual formulation of the TV-H −1 problem. In the concluding Section 5 we present some numerical results obtained by applying the ADI methods above to some problems arising in image processing.

Preliminaries
In this section we briefly introduce the main notation employed in the paper and the definition of the differential discrete operators we will use in the following to describe the numerical schemes. Dealing with a fourth-order nonlinear PDE, we also need to make precise what we mean exactly by stability of a numerical scheme.

Notation
In this paper we discuss the numerical solution of evolutionary differential equations. We need to distinguish between the exact solution u of the continuous equation and the approximate semi-discrete and fully discrete solutions of the numerical schemes. Let u(x, y, t) the exact solution at point (x, y) ∈ Ω = [a, b] × [c, d] and t ≥ 0. We write u n (x, y) to indicate the approximation in time of the function u(x, y, t n ), where (x, y) ∈ Ω, t n = n∆t, n ≥ 0 and ∆t is the time step size. Further, we approximate Ω by a finite grid {a = x 1 < . . . < x N = b} × {c = y 1 < . . . < y M = d} with equidistant spatial size h = (b − a)/N = (d − c)/M . We then denote by U (t) the approximation of u(x i , y j , t) in the node (x i , y j ) at time t. Finally, the fully discrete approximation of u is denoted by U = (U n ) i,j . When dealing with vectors, we will indicate their components using the superscripts notation: Y = (Y 1 Y 2 ) . Finally, we will indicate by δ * the differential operator acting in the direction * and we will use the notations D ∇ , D div and D ∆ U to indicate the discrete gradient, divergence and laplacian of the approximating solution U , respectively.

Definition of the discrete operators
In the following we define all the discrete operators that we use for approximating the spatial derivatives that appear in our numerical schemes. We use periodic boundary conditions throughout. We also performed the discretisation changing the boundary conditions to homogeneous Neumann boundary conditions, but this did not affect the numerical results significantly.
We approximate the first derivatives u x (x i , y j , ·) by forward differences (δ + x U ) i,j and backward differences (δ − x U ) i,j , where: The first derivatives of u with respect to y are approximated analogously by (δ + y U ) i,j and by (δ − y U ) i,j . The pure second derivatives will be approximated by using the five-point formula, this means the Laplacian operator ∆u = u xx + u yy is approximated by: The mixed derivatives u xy (x i , y j , ·) are approximated by: For the discretisation of the pure and mixed fourth-order derivatives appearing in Section 4.1, we will use a 5 × 5 stencil and approximate by:

Definition of stability
In the following, we introduce the notion of stability of the numerical schemes solving the nonlinear equations (1.4) and (1.5): Definition 2.1. Let u be an element of a suitable function space H defined on Ω × [0, T ], with Ω ⊂ R 2 open and bounded and T > 0. Let F be a real valued function and u t = G(u, D α u) a partial differential equation with all spatial D α with order |α| ≤ 4. A corresponding time stepping method where G n is a suitable approximation of G in u n and u n+1 is unconditionally stable if all the solutions of (2.7) are bounded for all ∆t > 0 and all n such that n∆t ≤ T and stable if, for a given ∆t > 0, the solutions of (2.7) are bounded for all n such that n∆t ≤ T .

ADI splitting schemes
In this section we introduce the numerical method of directional splitting. In particular, we discuss three types of directional splitting: the Peaceman-Rachford scheme, the Douglas-Hundsdorfer scheme, and an additive multiplicative operator splitting scheme.
In everything that follows we consider large systems of generic ordinary differential equations that arise from a semi-discretisation of initial boundary value problems such as for a given function F . The explicit dependence on time of F may not appear, thus considering autonomous systems. Typically, directional splitting methods (see, for instance, [31]) decompose the function F into the sum: where s is the spatial dimension of the problem. The components F j , j = 1, . . . , s, encode the linear action of F along the space direction j = 1, . . . , s, respectively, and F 0 contains contributes coming from mixed directions and/or non stiff nonlinear terms. A generic ADI scheme constitutes a time-stepping method that treats the unidirectional components F j , j ≥ 1 implicitly and the F 0 component, if present, explicitly in time.

The Peaceman-Rachford scheme
The first method we consider in this framework is the second-order Peaceman-Rachford ADI method (see [45] and [31]) where the operator F can be splitted into F = F 1 + F 2 , i.e. no mixed derivative or nonlinear terms are present. Let ∆t > 0 be the time step size of the scheme, then for every n ≥ 1 the numerical solution U n+1 of (3.1) is computed as follows: In (3.3) we can see that forward and backward Euler are applied alternatingly in a symmetrical fashion, thus obtaining second-order accuracy (see [31]). In each half time step, one of the two dimensions is treated implicitly, whereas the other one is explicitly considered. Note that (3.3) can be generalised to problems with nonlinearities f (U ) and by adding 1 2 f (U n ) and 1 2 f (U n+1/2 ) to the right hand side of the first and second equation in (3.3) respectively. The scheme (3.3), however, does not have a natural extension for more than two components F j , so more general ADI methods have been proposed in the literature.

The Douglas-Hundsdorfer scheme
Another ADI method that allows a more general decomposition like the one in (3.2) is the so-called Douglas method (see [30], [31] and [32]). In that, the numerical approximation in each time step is computed by applying at first a forward Euler predictor and then it is stabilised by fractional s steps where just the unidirectional components F j appear, weighted by a parameter θ ∈ [0, 1]. Its size controls the implicit/explicit character of such steps. In other words, the unidirectional operators are applied to the convex combination θU n+1 + (1 − θ)U n , thus considering fully implicit steps for θ = 1, explicit ones for θ = 0 and a Crank-Nicolson type scheme when θ = 1/2. The consistency order in time of the scheme is equal to two whenever F 0 = 0 and θ = 1/2 and it is of order one otherwise. In many applications we need to consider, however, F 0 = 0 (for instance, if we are considering contributions coming from mixed derivative operators) for any given θ. Some extensions of the Douglas scheme have been proposed in order to overcome this drawback. The following scheme proposed in [32] by Hundsdorfer is an extension of the Douglas method where a second stabilising parameter σ > 0 appears: The advantage of this extension is that for any given θ the scheme (3.4) has timeconsistency order two if σ = 1/2 and one otherwise, independently of F 0 . The parameter θ is typically fixed to θ = 1/2. Its choice is discussed in [32]. Larger values of θ give stronger damping of implicit terms, lower values typically better accuracy. Stability properties of this scheme when applied to linear convection-diffusion equations with mixed derivative terms have been investigated in [28,29]. There, the preferable value for θ is θ = 1/2 + √ 3/6. In [55] the authors combine the approach presented above with iterative methods for solving nonlinear systems.

Additive multiplicative operator splitting schemes
As pointed out above, in both schemes (3.3) and (3.4) some explicit terms appear. These may affect stability properties of the methods and, generally, their accuracy. As described in [4], splitting schemes as directional splitting methods belong to the family of multiplicative locally one-dimensional (LOD) schemes. In their general semi-implicit form, when a splitting similar to (3.2) holds, they appear as: where, similarly as before, each operator F i , 1 ≤ i ≤ s, is acting just along the i-th direction, whereas F 0 encodes mixed contributions. As pointed out above, it is generally difficult to deal with such an operator as the matrix (I − ∆tF 0 ) is, generally, not tridiagonal. For this reason, the scheme (3.4) deals explicitly with such a term. Analogously, explicit components appear also when applying (3.3) because of the alternating application of forward and backward Euler. As we will point out later on, these explicit contributions may create stability problems in the methods we are going to present. Following the strategy presented in [4], our attempt is to modify (3.3) such that no explicit contributions appear. The cost of such an operation will be reducing the accuracy of the method to order one, against the second-order achieved with the classical Peaceman-Rachford method (3.3). In order to preserve such accuracy as well as the symmetry of the method, at each time step two calculations are performed: which, written in the same form as (3.3) and (3.4), read as: To get the numerical solution U n+1 we simply average: thus ensuring a symmetric splitting. Due to the nature of such a method, we refer to (3.6)-(3.7) as additive multiplicative operator splitting (AMOS) ADI method. Note, that this scheme is identical to an earlier version of the well-known Strang splitting (see, for instance, [31, (1.12), p.329]). For nonlinear problems, such a scheme is second-order accurate, in contrast to first-order accuracy of the classical Strang splitting scheme (compare [4]). Furthermore, the scheme (3.6)-(3.7) has the advantage of allowing a parallel implementation, as suggested in [34,35,54].

Applications to higher-order PDEs
As we are dealing with the particular case of regular domains in R 2 , we will consider in the following s = 2. Thus, the F 1 component (F 2 , respectively) will contain operators acting just along the x-direction (y-direction, respectively). When appearing, the term F 0 will deal, instead, with the mixed xy-direction. Our aim is to adapt the ADI schemes

A linear example: the biharmonic equation
We start by considering the auxiliary linear, fourth-order parabolic biharmonic equation: Such equation appears in many applied mathematical models such as the Cahn-Hilliard models describing phase transitions and separation in binary mixtures (see, for instance, [22] and [42]). Some work on ADI schemes applied to the biharmonic equation already exists in literature. It dates back to [18] where the authors consider the equation to model vibrational modes for thin plates and it has been analysed also in [55,28] where linear stability results are proved as well. We are looking for the solution U of the following semi-discretised version of (4.1): According to the rules of ADI schemes, we decompose the function F into the sum where the components F i , i = 0, 1, 2, are defined by and the differential operators have been discretised as discussed in Section 2. With this choice, we can find for every n the approximating solution U n+1 of (4.1) using the Hundsdorfer ADI scheme (3.4).
Simplifying the problem by splitting the fourth-order equation into a mathematically equivalent autonomous system of two partial differential equations of order two produces the system: Then, the Hundsdorfer scheme applied to (4.2) can be equivalently written as a coupled ADI scheme for approximate solutions (U n , V n ) of (4.3). For positive parameters θ, σ this gives: where the functions F, F 1 , F 2 and G, G 1 , G 2 are given by: To verify that (4.4) gives the same solution as the Hundsdorfer scheme applied to (4.2) we consider as example the first implicit step in (4.4) providing the approximations (Y 1 1 , Y 2 1 ). For the approximation Y 1 1 of U n+1 we have: Using now the expression of Y 2 1 given by the implicit step relating to the equation for V we have: which, compared to the respective step performed applying the Hundsdorfer scheme (3.4) to equation (4.2), gives exactly the same result. We perform the same technique for the other implicit steps. Moreover, as the reader may note, in both the explicit steps of the scheme above we swapped the order of application of the method for consistency issues. Namely, we first find consistent approximations for V n+1 using them to get consistent approximations of U n+1 . In the application of both these methods we get stable solutions both for ∆t = C(∆x) 3 and also for ∆t = C(∆x) 2 , improving largely upon the condition ∆t = C(∆x) 4 that a naive explicit discretisation would give. Scheme (4.4) is indeed unconditionally stable, as proved earlier in [28] and confirmed in our numerical experiments in Section 5. However, note that considering bigger time steps such as ∆t = C(∆x) the numerical accuracy suffers not producing sensible solutions.

Applications to the primal formulation of TV-H −1 equation
We now aim to derive an ADI method solving the TV-H −1 equation (1.4). Expanding directly the differential operators appearing in the equation generates an intractable number of nonlinear terms of various differential orders. This makes a direct application of the ADI scheme to (1.4) impractical. Therefore, following the ideas presented in Section 4.1, we reduce the original fourth-order equation to an autonomous system of two secondorder equations. In the following we consider the primal formulation of (1.4), in contrast to the primal-dual formulation presented in Section 4.3. To do so, we first regularise the subgradient of the total variation in (1.4) characterised by (1.3). We use the standard square root ε-regularisation that replaces |∇u| by |∇u| ε := |∇u| 2 + ε, 1 ε > 0, see for instance [9]. This results in the following regularised version of (1.4) (4.6) u t = −∆∇ · ∇u |∇u| ε and, equivalently, In the following we present two different linearisations of the problem above. Heuristically, such a choice is important from two different points of view, intrinsically related to each other. The former is the accuracy of the scheme we are considering: rough linearisations (i.e. linearisations which consider most of the nonlinear terms explicitly evaluating them in one or more given approximations of the solution in previous time steps) are likely to present poor accuracy as well as stability issues. This is a general consideration in the numerical solution of every partial differential equation and it has to be taken into account and balanced with the choice of linearisation which might be more accurate and precise, but which could present, on the other hand, difficulties in its implementation and application. The latter point of view is, in some sense, peculiar to our choice of performing a directional splitting scheme. As pointed out above, our purpose is splitting our partial differential operator into the sum of components which are considered both explicitly (see F 0 above) and implicitly (see F 1 and F 2 ), as in (3.2). As we are going to present in the following, the choice of the linearisation affects such a splitting as the F 0 component and the linearised quantities multiplying the differential operators acting in x and y may change accordingly. For instance, the F 0 component might not appear changing the choice of the ADI scheme we want to use.
In the following we proceed by presenting two ADI schemes of the form (3.3) and (3.4) for two different linearisations of (4.6). For a given initial condition (U 0 , V 0 ), our problem consists in finding an approximation (U n+1 , V n+1 ) of the solution (u(t n+1 ), v(t n+1 )) to (4.6) for every n ≥ 0. In the following we always linearise around the solution at the previous time step (U n , V n ).

The first linearisation
Indicating byŨ the value of the solution U n in the previous time step, our first choice of linearisation is the following (compare with [20]): where (U, V ) is the semi-discrete approximation to a solution of (4.6). Here, the linearisation of V constitutes a semi-implicit approximation of the second-order nonlinear diffusion, evaluating all the first-order derivatives of the expansion in the previous time step. As before, we use the following notation for the system (4.7): for suitable matrices A, B, C and D in R N M ×N M . We split F and G into the sum of three different terms: F 0 and G 0 containing the mixed derivative term and F 1 , G 1 and F 2 , G 2 containing the derivatives with respect to x and y only, respectively. This produces the splitting: with: Due to the presence of a mixed derivative operator, a simple Peaceman-Rachford scheme cannot be applied. The Hundsdorfer scheme (3.4) is needed. The ADI scheme we use has the same form as for the biharmonic equation, i.e. we employ (4.4) with F, F i and G, G i given by (4.9), (4.10).

The second linearisation
Another possibility is to linearise the system (4.6) in the following way: where againŨ = U n and the spatial quantities are discretised as above and the discrete operators δ up x and δ up y are defined below. We observe that with this choice no mixed derivative operator acting on U appears. Mixed terms are encoded and considered in the previous time step. On the other hand, we get first derivative operators and not just second-order ones as in (4.7). Writing again the system (4.11) as in (4.8), we now split F and G in the following way: We note that the splitting (4.12) is a two-components splitting as the one provided for the Peaceman-Rachford ADI scheme (3.3) applied to the system (4.11) and this gives: (4.14) For the discretisation of the first derivative operators in the scheme -present in the equation for V in (4.11) -we use the standard numerical technique of upwinding, i.e. the sign of the coefficients in front of the first derivatives terms affects in which direction the finite differences are computed. More precisely, we use: and 1 S is the indicator function for the set S. A numerical discussion pointing out the differences of the ADI schemes resulting from the two linearisations (4.7) and (4.11) follows in Section 5.

Discussion of stability restrictions for the Hundsdorfer scheme
As we are going to illustrate numerically in Section 5, a stable application of the ADI schemes to equation (1.4) depends on the choice of the regularising parameter ε. In order to use reasonably large time steps ∆t, this parameter has to be taken sufficiently large to get stable results for the numerical solution of (1.4). For the following stability consideration we use the terminology introduced in Definition 2.1 where we consider the solution continuous in space and discrete in time.
Fully explicit numerical schemes solving TV gradient flows turn out to show restrictive stability conditions related to the strength of the nonlinearity in the TV subgradient, cf. [13,19]. On the other hand, fully implicit schemes solving (4.11) without any operator splitting are unconditionally stable. In particular, we have the following stability theorem: Theorem 4.1. Let u 0 be a sufficiently regular initial condition and u n the solution of (4. 16) u n+1 = u n − ∆t∆∇ · ∇u n+1 |∇u n | ε .
Then, the following stability estimate holds Proof. Multiplying equation (4.16) by ∆ −1 (u n+1 − u n ) (where ∆ −1 is the inverse of the negative laplacian with zero boundary conditions) and integrating over Ω we get: where ·, · denotes the L 2 inner product. We can rewrite the left hand side of the equation above using the properties of ∆ −1 and applying the divergence theorem, thus finding: We can now apply the result provided in [23] and summing over all t n = n∆t up to T = N ∆t, finding the following stability estimate: , which gives (4.17). In particular, estimate (4.17) does not depend on the size of ε.
These considerations about explicit and implicit schemes solving directly (i.e. without any splitting) our problem (1.4) serve as a motivation for the following estimates. We focus on the Hundsdorfer scheme (4.4) applied to the TV-H −1 equation with the choice (4.9), (4.10). In each iteration the numerical solution is computed from equations consisting of a combination of explicit and implicit quantities. In particular, the explicit quantities might affect the stability properties of the scheme. To motivate this, we focus in the following just on the first three stages of the scheme (4.4) applied to the TV-H −1 equation with the choice (4.9), (4.10) and θ = 1/2. Considering the first three stages of (4.4) only can be justified by the fact that the subsequent three stages of the scheme have a similar structure and are not expected to change the stability properties drastically.
Combining the three steps of the scheme (4.4) with (4.9), (4.10) we find the following expression: where ∂ x denotes the continuous derivative with respect to x and the positive quantities C 1 (u n ) and C 2 (u n ) come from the linearisation described in Section 4.2.1. They read: We observe that a mixed, eighth-order operator appears, both on the left and on the right hand side of (4.18). In the following stability discussion we neglect these high-order terms which only represent second-order in time contributions. Now, we multiply equation (4.18) by u n+1 and integrate over the domain Ω. By applying integration by parts twice with respect to the x and the y variables to the second and the third terms of the left hand side of the equation, respectively, we get: Here K(ε) is a suitable constant that depends on the regularising parameter ε only. A similar strategy is applied to the analogous terms on the right hand side, where we have also used Young's inequality with weights δ 1 and δ 2 . We obtain: The first term on the left hand side of (4.18) can be dealt with by using Young's inequality again. It remains to consider the first, nonlinear, term on the right hand side of (4.18). By applying the divergence theorem, integration by parts and Young's inequality with weight δ 3 , we get for this term the following estimate: The second term on the right hand side of the inequality above can be merged with the corresponding ones in (4.19), choosing δ 1 small enough. Denoting by C 3 (u n ) = ∂xun∂yun |∇un| 3 ε , for the curvature term in (4.20) we observe that: where we used Cauchy's inequality and upper bounds on C 1 , C 2 and C 3 . Defining K 1 (ε), K 2 (ε) and K 3 (ε) as and using an estimate proved in [17] for the mixed derivative term, we get the following bound Collecting the previous estimates, choosing δ 1 , δ 2 and δ 3 small enough and applying once more upper bounds on C 1 and C 2 we get the following stability estimate for scaled constantsK 1 ,K 2 andK 3 which tend to infinity as the regularising parameter ε 0. For this limit the estimate then blows up, thus indicating possible unstable behaviour when choosing ε small. This will be discussed in more detail in Section 5.

A stable AMOS ADI scheme
In order to counteract the dependence of the stability properties of the ADI schemes (4.4) and (4.14) on the size of ε, we now consider as an alternative the AMOS operator splitting scheme (3.6)-(3.7) for solving (4.6). Due to the fully implicit character of the scheme, the hope is that stability properties improve.
To see this, both ADI methods (4.4) and (4.14) can be represented in a vectorial, multiplicative form similar to (3.5). For (4.4), the operator F 0 appearing in (3.5) is taken explicitly in time. Thus, when writing the correspondent numerical scheme in a multiplicative form, we have additional explicit terms on the right hand side. Writing (4.14) in the form (3.5) we again obtain additional explicit terms appearing due to the forward Euler steps in (4.14). This, together with the stability estimates from the previous section, indicates possible stability restrictions for these schemes that will be confirmed by the numerical results in the following Section 5.
In the following we present an AMOS scheme for solving a slightly modified version of (1.4). In system form, this new equation reads: This equation is more anisotropic than the original (1.4) in the sense that the nonlinear diffusion in x and y directions are considered separately. Only the diffusion weighting involves the whole image gradient taking both x and y variations into account. In this way, linearising v 1 and v 2 by considering the diffusion weighting 1/|∇ũ| ε for a givenũ, results in an equation with only pure x and y derivatives. Considering such an equation reduces the explicit components appearing in the scheme, allowing a fully implicit treatment of the operators, though still exploiting the advantage of directional splitting by solving along the two directions x and y separately. Applying the AMOS scheme (3.6)-(3.7) to the linearisation of (4.22), we obtain: where the operators are defined exactly as in (4.13) and upwinding is used for the first derivatives as in (4.15). We recall that the alternating application of the scheme first in the y − x direction and subsequently in the x − y direction allows to achieve order two of accuracy, as explained in [4]. As we will see in Section 5, the scheme (4.23) has better stability properties where the choice of the time step does not seem to depend on the size of ε, thus suggesting -at least empirically -unconditional stability.

Primal-dual formulation of TV-H −1 equation with penalty term
An interesting alternative to the linearisations of the primal formulation of equation (1.4) is motivated by earlier work on numerically characterising elements in the subdifferential of the total variation seminorm by primal-dual iterations, see [5,39,6] for instance. In what follows we briefly outline such a strategy combined with ADI splitting for solving equation (1.4), that is (4.24) u t = ∆q in Ω × (0, ∞), q ∈ ∂|Du|(Ω).
Here, q ∈ ∂|Du|(Ω) by definition of the subdifferential means which is typically known as the primal-dual formulation of the problem (4.24). The constraint on p appearing in (4.27) can be relaxed, for instance, by a penalty method. That is, we remove the constraint from the minimisation in (4.27) and instead add a term that penalises the functional if p ∞ > 1. A typical example for such a penalty term F is F (s) = 1 2 max{s, 0} 2 2 .
With these considerations we reformulate (4.27) as where the parameter 1 ε > 0 is the weight of our penalisation. We can then find the optimality conditions for solutions p and u of (4.28) which, merged with equation (4.24), allow us to consider the following, approximate formulation of (1.4) In the system above we have indicated by H the derivative of the penalty term F (|p| − 1), i.e.
which we linearise via its first-order Taylor approximation, that is Here, H denotes the Jacobian of H. In order to guarantee the invertibility of the now linear operator that defines the system (4.29) we add an additional damping term in p as suggested, for instance, in [39]. Collecting everything, we propose the following numerical scheme for solving (4.29), which consists of two nested iterations, the inner damped Newton iteration with index k and an outer time stepping with index n for the evolution in time of U . Here, τ k is a sequence of parameters controlling the damping of the Newton iteration in every time step. The authors in [39,6] suggest to start with a large value τ 0 and then decrease it in the inner iterations to ensure efficient convergence.
System (4.30) could now be discretised in space and either solved directly using Newton's iterations (compare [6] where the authors used this strategy to solve equation (1.6)) or fitted into a numerical ADI scheme. A detailed presentation and analysis of the resulting scheme is beyond the scope of the present paper and a matter of future research.

Numerical results
In this section we discuss the numerical performance of the ADI methods proposed in this paper. To do so, we present numerical experiments for a Gaussian, an oscillatory function made up of sine and cosine functions and grayscale images as initial conditions. The paper is furnished with numerical examples for image inpainting.
We start presenting the numerical results obtained applying the Hundsdorfer scheme (3.4) to compute the numerical solution of the biharmonic equation (4.2) and to the equivalent system (4.3). Next, we report on numerical experiments for the TV-H −1 equation (1.4) using the Hundsdorfer scheme (4.4) with (4.9)-(4.10). The choice of the time step size ∆t of this scheme is constrained by very strong stability restrictions related to the size of the regularising parameter ε. A possible dependence of this type for stability was already discussed in Section 4. Our numerical tests for the application of the Peaceman-Rachford scheme (4.14) to the TV-H −1 equation show the same stability behaviour and are therefore not included in the paper. Finally, we show the numerical results for the solution of the modified TV-H −1 system (4.22) solved with the AMOS scheme (4.23). This scheme shows stable behaviour independent of the size of ∆t and ε.

The biharmonic equation: numerical results
The linear system that arises from the application of the Hundsdorfer scheme (4.4) with (4.5) to the biharmonic equation (4.3) is solved by the Schur complement method.
We consider a grid of 100 × 100 grid points for the discretisation of the spatial domain Ω being the unit square. We analyse the example of the evolution of the biharmonic equation having as initial condition U 0 the Gaussian density U 0 ij = exp (−((x i − 1/2) 2 +(y j + 1/2) 2 )/γ 2 ) where the variance γ 2 is equal to 100. Figure 1 shows two iterates of the scheme with θ = σ = 1/2 for a time-step size ∆t = C(∆x) 2 , where C now and for the rest of the discussion is equal to 0.1. Even for such a big choice of ∆t the result is stable and the convergence of the solution to the steady state is quick. Considering another initial condition U 0 and taking, for instance, some very oscillatory function does not effect the performance of the method. This confirms the unconditional stability of this scheme when applied to a linear fourth-order equation such as the biharmonic equation (4.2), compare [28]. In the following section we will see that the unconditional stability of the Hundsdorfer scheme breaks down when applied to the nonlinear fourth-order diffusion equation TV-H −1 (1.4).

TV-H −1 equation: numerical results
In what follows we provide numerical discussion for applying the directional splitting schemes introduced in Section 4 for the numerical solution of the regularised TV-H −1 equation (4.6). In particular we consider Hundsdorfer ADI scheme (4.4) with (4.9), (4.10) and the AMOS scheme (4.23). Once again we use the Schur complement technique to solve the linear systems that arise in the numerical solution of these schemes. The expected edge-preserving behaviour of the TV-H −1 equation (1.4) which is due to the subgradient of the TV functional is closely related to the size of the regularising parameter ε used in (4.6). This parameter, in fact, becomes a "measure" of how close the nonlinear diffusion is to the linear, biharmonic one. Choosing ε too big the smoothing behaviour of (4.6) is similar to the one of the biharmonic equation (4.2). In this case the stability properties of the Hundsdorfer scheme applied to the nonlinear equation are close to the ones discussed for the biharmonic equation in the previous section. On the other hand, small values of ε keep the regularised version of the subgradient of the total variation close to its exact characterisation and hence solution show edge-preserving features. However, stability issues that disturb and worsen the performance of the methods appear. The explicit treatment of some of the terms in the Hundsdorfer and the Peaceman-Rachford schemes (namely, the mixed derivative term in (4.4) and the half-direction forward Euler steps in (4.14)) seem to influence the stability of the schemes in a negative way. In particular, the time step sizes ∆t have to be decreased significantly with small values of ε. Moreover, the choice of the initial condition also influences the stability properties. For instance, for smooth Gaussian initial conditions with large support the time steps can be chosen larger than for an oscillatory initial condition, see Figures 2-3. These issues are resolved by the AMOS scheme (4.23) which did not show dependence of ∆t on the size of ε nor on the type of initial condition in order to get stable results.
In the following we present the numerical results obtained considering the linearisation of the system (4.6) given by (4.7) and solved by the Hundsdorfer ADI method (4.4). We consider as initial conditions both the Gaussian density U 0 ij = exp (−((x i − 1/2) 2 +(y j + 1/2) 2 )/γ 2 ) with γ 2 = 100 and the oscillatory function U 0 ij = sin(8πx i ) + cos(8πy j ). Figure 2 shows iterates of the Hundsdorfer scheme with θ = σ = 1/2, ε = 5 and ∆t = C(∆x) 3 applied to the Gaussian datum. The value ε = 5 is the smallest possible value that can be used in order to get stable solutions of the Hundsdorfer scheme with ∆t = C(∆x) 3 . Figure 3 shows the evolution of the process for the oscillatory datum with the same choice of ε as before. In this case to get stable results smaller time steps are needed, thus showing stability dependence on the initial condition.  In Figures 4-5 we describe this stability dependence in more detail. For the two different choices of the initial conditions considered above, the smallest values ε needed for stability as the size of ∆t increases are plotted. To compare the different setups of the Hundsdorfer scheme with respect to θ, these tests were performed for different values of the stabilising parameter θ = 0, 1/2, 1/2 + √ 3/6, 1 thus resulting into four graphs per test. Note that different values of θ (compare [31,32,28,29] for such choices) result into different weighting of implicit and explicit terms, that means considering a fully explicit method for θ = 0 and implicit contributions in the unidirectional steps of (4.4) in the other cases (see Section 3.2). For each of these graphs, their epigraph corresponds to the region of stability of the method (according to Definition 2.1). To obtain such graphs we have considered different values of the constant C and of the regularising parameter ε for time step sizes of the order (∆x) k , k = 2, 3, 4. For the choice θ = 0 stable solutions are only obtained by very restrictive choices of ε. This situation improves for θ > 0. In particular, for values of θ close to 1 the stability constraint on the size of ε is reduced. However, in all cases, these plots show a clear dependence of the stability of the Hundsdorfer scheme on the strength of the nonlinearity (encoded by the size of ε). This creates numerical difficulties in the attempt of increasing the time step size to get an efficient numerical scheme for solving the approximated problem (4.6) for sufficiently small values of ε.   Figure  4 shows dependence on the initial condition for admissible values of ε providing stability.
We do not present here the numerics related to the application of the Peaceman-Rachford method (4.14) to the TV-H −1 equation (4.22) as the stability issues resemble the ones described above. In order to overcome such problems, we present in the following the results related to the application of the ADI AMOS scheme (4.23) solving the slightly modified system (4.22).
Due to the implicit character of the scheme (4.23), improved stability properties are expected while keeping the advantages of the directional splitting strategy, that is each can be solved very efficiently. In the Figures 7-8 we show the evolution of the TV-H −1 equation (4.22) for the Gaussian initial condition as above for ∆t = C(∆x) 3 and ∆t = C(∆x) 2 with fixed ε = 0.001. We observe that the time discretisation provides stable results even for large ∆t. However, as clearly visible in Figure 8, choosing ∆t too large badly affects time accuracy. The convergence rate of the AMOS scheme (4.23) for the choice of ∆t = C(∆x) 3 and regularising parameter ε = 0.001 when applied to the Gaussian initial condition is presented in Figure 9. We observe an exponential-type decay for the ∞ norm of the difference between the iterative solution U n of (4.23) and the steady state U ∞ which has been computed numerically beforehand by iterating the scheme till the quantity U n+1 − U n ∞ /(M N ) ≤ 10 −10 . For comparison, the exponential function (7·10 −6 )e −0.05x , has been plotted. Figure 10 shows the decay of total variation energy (1.2) towards the energy of the steady state U ∞ against the number of iterations. This shows that although the modified TV-H −1 equation (4.22) is not as the same as the gradient flow of the total variation in the space H −1 , the total variation energy is still decreasing in every iteration.   Motivated by our original purposes of applying the ADI schemes in the imaging framework, we present in Figure 11 the scale space properties of system (4.22) solved with the AMOS scheme (4.23) for a 120 × 120 image. The time step size is ∆t = C(∆x) 3    Finally, we present some numerical results obtained by using the AMOS scheme (4.23) to solve the H −1 -inpainting model (1.5) with the modified TV-H −1 equation given by (4.22). The constraint u = f outside of the inpainting domain is approximately enforced by adding the fidelity term λ·1 Ω\D (f −u) to the equation, where λ measures how close the reconstructed image is to the original one. In Figure 12 we show the result for inpainting a 150 × 150 image of a cross: the initial condition is inpainted in 1000 iterations using a time step size of∆t = C(∆x) 3 and the regularising parameter ε is chosen to be ε = 0.001, thus allowing the preservation of edges additionally to fulfilling the connectivity principle that is a consequence of the fourth-order of the method. Figure 12: Solution of inpainting problem obtained with 1000 iterations of the ADI AMOS scheme (4.23), ∆t = C(∆x) 3 , ε = 0.001.

Conclusion
In this paper we have studied the applicability of directional splitting techniques to fourthorder nonlinear diffusion equations. In particular we have considered the TV-H −1 equation which is the H −1 gradient flow of the total variation that emerges from imaging applications such image denoising and inpainting.
The numerical challenges when solving this equation are both its fourth-differential order and the strong nonlinearity due to the total variation subgradient in the equation. In former work on numerical schemes for this equation these issues have been either addressed by using explicit time stepping schemes with tiny step sizes, semi-implicit schemes with heavy damping, e.g. [48], or by fully implicit schemes, e.g. [19], which are computationally expensive to solve. Having this in mind directional splitting seems to be a promising compromise. The computational cost of each of its iterations is very small due to the triangular form of the system matrices, while the stability conditions are improved when compared to an explicit in time discretisation.
We have proposed three different ADI methods applied to two linearisations of the TV-H −1 equation and presented numerical results for Gaussian-type initial conditions, image smoothing and image inpainting. Our investigation of stability properties of the schemes results into the following observation: any explicit terms, due to the type of splitting used, result into heavy stability conditions on the size of the time steps. For fourth-order equations these restrictions usually turn out to make a choice of ∆t as small as (∆x) 4 necessary, thus preventing an efficient numerical solution of the equation. In particular, even if each iteration of the numerical scheme is very cheap -due to its "onedimensional" character -a large number of them have to be computed in order to see any significant change of the solution in time. The AMOS scheme provides a stable and cheap numerical scheme for solving a slightly modified version of the TV-H −1 equation, which turns this fourth-order nonlinear equation numerically tractable and hence makes it more attractive for imaging applications.
Exchanges Award IE110314 for the project High-order Compressed Sensing for Medical Imaging, the EPSRC first grant Nr. EP/J009539/1 Sparse & Higher-order Image Restoration, and the EPSRC / Isaac Newton Trust Small Grant on Non-smooth geometric reconstruction for high resolution MRI imaging of fluid transport in bed reactors. Further, this publication is based on work supported by Award No. KUK-I1-007-43 , made by King Abdullah University of Science and Technology (KAUST).