Method of lines transpose: an efficient A-stable solver for wave propagation

Building upon recent results obtained in [7,8,9], we describe an efficient second order, A-stable scheme for solving the wave equation, based on the method of lines transpose (MOL$^T$), and the resulting semi-discrete (i.e. continuous in space) boundary value problem. In [7], A-stable schemes of high order were derived, and in [9] a high order, fast $\mathcal{O}(N)$ spatial solver was derived, which is matrix-free and is based on dimensional-splitting. In this work, are interested in building a wave solver, and our main concern is the development of boundary conditions. We demonstrate all desired boundary conditions for a wave solver, including outflow boundary conditions, in 1D and 2D. The scheme works in a logically Cartesian fashion, and the boundary points are embedded into the regular mesh, without incurring stability restrictions, so that boundary conditions are imposed without any reduction in the order of accuracy. We demonstrate how the embedded boundary approach works in the cases of Dirichlet and Neumann boundary conditions. Further, we develop outflow and periodic boundary conditions for the MOL$^T$ formulation. Our solver is designed to couple with particle codes, and so special attention is also paid to the implementation of point sources, and soft sources which can be used to launch waves into waveguides.

Introduction
Numerical solutions of the wave equation have been an area of investigation for many decades. The wave equation is ubiquitous in the physical world, arising in acoustics, electromagnetics, and fluid dynamics. Our main interest is in electromagnetic wave propagation, where traditional finite difference methods such as the finite-difference time-domain (FDTD) method are often used to solve Maxwell's equations. When the classical Yee scheme is used, the Courant-Friedrichs-Lewy (CFL) stability criterion restricts the time step to scale with the smallest cells in the domain. This becomes computationally prohibitive in a variety of interesting problems, such as electromagnetic scattering or waveguide design, where complex geometries need to be embedded in a Cartesian mesh, leading to small spatial cells near the boundaries. Alternatively, problems with multiple and disparate temporal or spatial scales, such as those presented by plasma simulations, require time steps and mesh spacings which are on the order of the shortest temporal and spatial scales. Due to the large number of charged particles in a typical simulation [4], and the high dimensionality of their distribution space, time steps are computationally prohibitive, and the amount of them must be minimized.
As a result, A-stable Maxwell solvers have been developed, which remove the CFL restriction and restore the ability of the user to define a time step which is based on the physical problem, rather than on its spatial discretization. Notable advances include the introduction of the FDTD-ADI [14,15,23,24] algorithm, as well as several semi-implicit time split schemes. However, it remains difficult to preserve both accuracy, and A-stability when non-rectangular domains are required.
Alternatively, Maxwell's equations can be reformulated (e.g., using the scalar and vector potential formulation), so that each field component independently satisfies the second order wave equation. Method of lines (MOL) formulations of the wave equation are well-studied, and often lead to conditionally stable schemes. But when the discretization is first performed in time, following the method of lines transpose (MOL T ), the wave equation is found at discrete time levels by solving a semi-discrete boundary value problem. In [1,2], the MOL T is used to produce an exact integral formulation of the wave equation, where the solution is determined by a convolution over the domain of dependence against a space-time Green's function G d (x, t), in dimensions d = 1, 2, 3. However, more commonly the semi-discrete solution is obtained, in which the modified Helmholtz equation must be solved, and the corresponding semi-discrete Green's function G(x, ∆t) is the Yukawa, or modified Helmholtz kernel. The resulting boundary integrals methods can then be solved at each time step using fast summation algorithms, such as tree-codes [3,12,20,21], or the fast multipole method (FMM) [17,10,13,19,16,20,11]. These methods often scale as O(N ) or O(N log N ), but require substantial storage or precomputing stages.
It is worth noting however that in one spatial dimension, the modified Helmholtz kernel is a simple decaying exponential. This fact has been combined with alternate dimension implicit (ADI) splitting to achieve an A-stable wave solver with computational complexity of O(N log N ) [6,22], and O(N ) [8,7], respectively. The scaling presented in [6,22] is due to the Fourier continuation method, which makes use of the fast Fourier transform to compute a periodic extension of the the convolution integral, where the period is taken sufficiently longer than the domain, so that boundary conditions can be satisfied.
In our approach [7,8,9], the one-dimensional convolution is instead computed strictly over the computational domain, using polynomial interpolation. In [8], the solution is proven to be A-stable, second order accurate, and that the matrix formed by the discrete convolution can be applied in O(N ) operations. However, the discrete convolution matrix was formed using piecewise linear integration, and it was found that a Lax-type correction was required to ensure convergence of the scheme in the semi-discrete limit ∆t → 0, where ∆x is fixed. This issue was addressed in [9], where spatial discretization was extended to orders M ≥ 2, on non-uniform grids while retaining O(N ) complexity. In [7], the scheme was extended to higher orders in time using a novel approach, successive convolution. These higher order schemes were proven to be A-stable.
The purpose of this present work is to develop a robust embedded boundary approach for complex geometry for the MOL T formulation. Because this paper is focused on developing a range of boundary conditions for the wave equation, we limit our attention to the 2nd order accurate solver. In addition to addressing standard closed (i.e., Dirichlet, Neumann and periodic) boundary conditions, we also develop an open, or outflow boundary condition which is suitable for our implicit solver. The method is extended to higher spatial dimensions using the factorization developed in [7], which is similar to but different from the traditional ADI splitting. This means that solutions are computed line-by-line along dimensional sweeps leveraging the O(N ) 1D solver to construct a high dimensional O(N ) implicit method.
The rest of this paper is laid out as follows. In section 2, we use the method of lines transpose (MOL T ) to reformulate the 2D wave equation as an implicit, semi-discrete modified Helmholtz equation. The Helmholtz operator is inverted analytically using dimensional splitting, and we recast the solution as a series of one-dimensional boundary integral equations. We specifically show how to obtain a fully discretized, second order accurate solution by performing spatial quadrature, using a fast O(N ) convolution algorithm.
In section 3, we modify the time-centered scheme of [7] to introduce artificial dissipation. This will be used to stabilize the embedded boundary method in our 2D algorithm, presented in section 4. We demonstrate our solution to be fast, second order accurate, and A-stable, even on non-rectangular domains with numerical results presented in section 5. We conclude the body of the paper with several remarks in section 6.
We also include in Appendix A, a table summarizing the main algorithmic aspects of our wave solver. In Appendix B, we show how point sources, which may represent charged plasma particles, or soft point sources launched into a waveguide, can be implemented.

Integral solution using MOL T
We now develop a dimensionally split algorithm for solving the initial value problem with consistent boundary conditions. We utilize the method of lines transpose (MOL T ) to perform a second order accurate temporal discretization, as was shown in [7]. In particular, the discretization is time-centered, and implicit in the spatial derivatives and β ≤ 2 is chosen to enforce A-stability [7]. Inversion of this operator leads to a boundary integral formulation, where the Green's function is the Yukawa potential, defined for 2D in terms of the modified Bessel function K 0 (r). However we will employ our previously developed [7,8] splitting algorithm, which has O(N ) complexity Here the subscripts denote the spatial component of univariate modified Helmholtz operators, The modified Helmholtz operators are formally inverted using the Green's function. We define convolution with this Green's function by the integral operator so that where the coefficients A and B of the homogeneous solution are determined by applying boundary conditions. A similar definition holds for L −1 y . Formally inverting both operators to the right hand side we find the explicit equation where the multidimensional operator is now . Since the application of L −1 x is done for fixed y, and vice versa for L −1 y , the operator D xy can be constructed in a line-by-line fashion, similar to ADI algorithms. It was also proven in [7] that the scheme (6) is A-stable, for 0 < β ≤ 2. Remark 1. In the continuous case, L x and L y (and their inverses) commute, so D xy = D yx . However in practice the discretize operators will not commute, and some small spatial anisotropy is introduced. This can be controlled by applying both operators D xy and D yx , and then averaging the result.
Remark 2. Similar to a traditional ADI formulation of the wave equation, this factorization produces a fourth order splitting error term, which can be compensated for by adding a term to the right hand side Note our use of the identity (L x − 1)(L y − 1) = ∂ xx ∂ yy /α 4 . Formally inverting both operators to the right hand side, the solution can be rearranged and found as where the convolution operator C is the operator defined in [7] as This form was used in [7] to achieve schemes of higher order through successive convolution, where removing splitting errors is of paramount concern.
2.1. Fast convolution algorithm for the integral solution. In [8], the particular solution (4) was discretized in space using the weighted midpoint and trapezoidal rules, which amounted to replacing u with a piecewise constant and linear approximation, respectively. More recently [9], we have detailed the spatial discretization of I x to arbitrary order, while providing a means for its rapid evaluation. In this work we focus on the second order accurate implementation of this algorithm, and reiterate the relevant details. This particular solution is first decomposed into a left and right oriented integral, split at y = x so that where These quantities can be updated locally, using exponential recursion The recursive updates (11) and (12) are exact (in space), and making δ L and δ R small (typically, δ L = δ R = ∆x) effectively localizes the contribution of the integrals. Based on these observations, we now outline the fast convolution algorithm. In [9], we derived a spatial quadrature of general order M ≥ 2 for irregular grids. Here, we will utilize a second-order accurate method on a uniform grid, which we explicitly develop.
Consider the domain (a, b) discretized by uniform grid points x 1 = a < x 2 < · · · < x N +1 = b of width ∆x = b−a N = x j+1 − x j , j = 1, ..., N . Suppose, given a function f compactly supported in (a, b), we are to evaluate the convolution operator at each grid points, that is compute We proceed by evaluating the local integrals J L j = J L [f ](x j ) and J L j = J L [f ](x j ), as in (11) and (12). The local integrals may be evaluated with quadrature, or, if possible, analytically. A second-order accurate quadrature is given by where defining ν = α∆x and d = e −ν , the quadrature weights are given by We summarize our fast method in Algorithm 1.

Artificial dissipation
As will be discussed in section 4.2.2, it is necessary to include some artificial dissipation in the numerical scheme to maintain stability with embedded boundary methods for Neumann boundary conditions. We first present a version of the wave solver based on a backwards difference formula (BDF) time discretization, leading to what we call a diffusive scheme, which is dissipative. This method has a larger truncation error than a centered scheme, does not possess a means to tune the level of dissipation, and also has an implicit source term (at time level n + 1), which is problematic for application in the context of particlein-cell (PIC) methods for the simulation of plasmas. The second-order centered (dispersive) scheme given above is therefore preferable, but are non-dissipative in their original forms. We give a method for adding tunable artificial dissipation terms into the centered scheme, while maintaining A-stability.
3.1. Diffusive wave solver. We substitute the following backward difference formula (BDF) discretization: into the wave equation 1 c 2 u tt − ∇ 2 u = S(x, t). Rearranging, defining α = √ 2/(c∆t) and dividing by α 2 gives the semi-discrete scheme This method is A-stable and dissipative, but does not possess a mechanism for tuning the dissipation, has an inconvenient implicit source term (at time level n + 1), and typically has a larger truncation error compared to the centered scheme.

Artificial dissipation in centered schemes.
3.2.1. Artificial Dissipation in 1D. We give a modified form of the centered version of the implicit wave solver with tunable artificial dissipation that retains the property of unconditional stability. We let denote a small artificial dissipation parameter, and D Ignoring sources, we have the second order scheme with dissipation, We now prove the unconditional stability of this scheme for prescribed values of β. As in [7], we pass to the high-frequency limit.
We obtain the Von Neumann polynomial ρ 2 − (2 − β 2 )ρ + (1 − ). We can check that the roots of this polynomial will be complex if 0 < β ≤ 2 + 2 √ 1 − , and that in this case the roots satisfy which shows both the stability and dissipative nature of the scheme. We note that the maximum allowed value of β is slightly smaller than what is allowed by the corresponding scheme without dissipation. A more detailed analysis shows that the effective damping rate is k 2 k 2 +α 2 2 , meaning that high frequencies are more rapidly damped than low frequencies.
3.2.2. Artificial Dissipation in 2D. Using the notation defined in [7], and again ignoring sources, we have the second order scheme with dissipation, For further details on these operators, see [7]. Numerical experiments indicate that the 2D scheme with artificial dissipation is indeed unconditionally stable with the same maximum value of β as with the 1D schemes.

Boundary conditions
We will now discuss the implementation of boundary conditions. Since our algorithm is dimensionally split, we first develop the boundary conditions in one spatial dimension where the solution (6) reduces to and where we consider the following boundary conditions Periodic Once these boundary conditions have been derived, we use them to build a boundary solver in 2D.

4.1.
Boundary conditions in one dimension. In 1D, this homogeneous solution requires the determination of two coefficients from the imposed boundary conditions and the endpoint values of the particular (integral) solution by solving a 2 × 2 linear system. We now show how to impose several common boundary conditions in 1D. These methods are extended to the 2D case in Section 4.2.
4.1.1. 1D Dirichlet boundary conditions. Let us begin with Dirichlet boundary conditions (25). Evaluating the semi-discrete solution (24) at x = a and b, we find which, after solving for the unknown coefficients can be written as . Homogeneous boundary conditions are recovered upon setting U L (t) = U R (t) = 0. Solving the resulting linear system for the unknowns A n and B n gives 4.1.2. 1D Neumann boundary conditions. For Neumann conditions, first observe that all dependence on x in the integral solution (5) is on the Green's function, which is a simple exponential function. Using this, we obtain the following identities Now, differentiating the semi-discrete solution (24), and applying the Neumann boundary conditions (26) at x = a and b yields which, after solving for the unknown coefficients can be written as Upon solving the linear system we obtain Remark 3. The cases of applying mixed boundary conditions at x = a and b are not considered here, but the details follow from an analogous procedure to that demonstrated above.

4.1.3.
1D Periodic boundary conditions. We impose periodic boundary conditions, by assuming that u n (b) = u n (a), u n x (a) = u n x (b), n ≥ 0. Enforcing this in the semi-discrete solution (24) then yields where we have used the identity (30) applied to derivatives of I. Solving this linear system is accomplished quickly by dividing the second equation by α, and either adding or subtracting it from the first equation, to produce 4.1.4. 1D Outflow boundary conditions. When computing wave phenomena, whether we are interested in finite or infinite domains, it is often the case that we must restrict our attention to some smaller subdomain Ω of the problem, which does not include the physical boundaries. We say that Ω is the computational domain, and that the boundary ∂Ω is the non-physical, or artificial boundary. Under these circumstances, it is necessary to enforce an outflow, or non-reflecting boundary condition, which allows the wave to leave the computational domain, but not incur (non-physical) reflections at the artificial boundary. For this reason, let us consider the free space solution but where we are only interested in evaluating this expression for x ∈ Ω = [a, b]. Then the contributions can be decomposed as where the homogeneous coefficients are which we observe do not depend on x. Since the coefficients A an B are the contributions of the integral to the left and right of [a, b] respectively, they can be thought of as transmission conditions (rather than boundary conditions). We make use of this fact to develop outflow boundary conditions, and it will serve as a key idea in our planed follow-up work on a domain decomposition algorithm and multi-core computing with our implicit wave solver. For the one-dimensional wave equation the exact outflow boundary conditions (28) turn out to be local in space and time. We emphasize that this is only the case in one spatial dimension, but we shall utilize this fact to obtain an outflow boundary integral solution from the integral equation (33). We extend the support of our function to (−∞, ∞), and extend the definition of the outflow boundary conditions to the domains exterior to [a, b] Next, assume the initial conditions have some compact support; for simplicity we will take this support to be Ω 0 = [a, b]. Then after a time t = t n , the domain of dependence of u n (x) is Ω t = [a − ct n , b + ct n ], since the propagation speed is c. Now the free space solution (33) becomes with coefficients At first glance, these coefficients are not at all helpful, as they require computing integrals along spatial domains which not only are outside of the computational domain, but also grow linearly in time. However, we will now make use of the extended boundary conditions to turn these spatial integrals into time integrals, which exist at precisely the endpoints x = a and b respectively. Consider first x > b. By assumption, this region contains only right traveling waves, u(x, t) = u(x − ct), and by tracing backward along a characteristic ray we find Thus, and so B n is equivalently represented by a convolution in time, rather than space. Now, knowing the history of u at x = b is sufficient to impose outflow boundary conditions. Furthermore, we find in analog to equation (12), a temporal recurrence relation due to the exponential where β = αc∆t, by definition (3). Thus, the coefficient B n , which imposes an outflow boundary condition at x = b, can be computed locally in both time and space. To maintain second order accuracy, we fit u with a quadratic interpolant and integrate the expression analytically to arrive at where In this outflow update equation (41), the quantities u n+1 (b) and B n are both unknown. In order to determine these values, we must also evaluate the update equation for u n+1 (24) at We now use these two equations to solve for u n+1 (b), and eliminate it from the outflow update equation (41), so that where Remark 4. While this procedure could be avoided by omitting u n+1 (b) in the interpolation stencil, it turns out to be necessary to obtain convergent outflow boundary conditions.
Likewise, upon considering x < a, we find Solving the resulting linear system produces where w Out Boundary conditions in two dimensions. We now describe our approach for imposing boundary conditions up to second-order accuracy in two dimensions for our MOL T formulation of the wave equation. Boundary conditions must be supplied for the intermediate sweep variable w. Since w = u + O(c 2 ∆t 2 ), for second order accuracy it suffices for w to inherit the boundary condition imposed on the main solution variable u. In the case of rectangular, grid-aligned boundaries, the 1D boundary correction terms can be imposed   Thus, a line y = y k , which is discretized in x, and has endpoints a k and b k , determined as the intersection of the line with the boundary curve. Analogously we define endpoints of the line x = x j as c j and d j , and now have for 1 ≤ k ≤ N y and 1 ≤ j ≤ N x respectively. If the boundary is defined using a level set function C(x, y) = 0, then the points x = a k , b k can be found by solving C(x, y k ) = 0, and following from the analogous approach, the endpoints y = c j , d j corresponding to x = x j are found.
Once the domain has been discretized and the lines have been defined, it remains to compute the discrete form of the scheme (6). Since the one-dimensional convolution algorithm presented in Section 2.1 is formulated for non-uniform grid points, the embedded boundary points do not affect the implementation of the horizontal (45) and vertical (46) sweeps. Thus, the only point of consideration that remains is that of boundary conditions.
First we consider an intermediate variable w (1) and which, according to the implicit scheme (2), can be seen to satisfy Boundary conditions must be applied to w (1) at the boundary points which terminate each line y = y k , defined as (a k , y k ) and (b k , y k ). Since u n+1 , u n , and u n−1 , will be prescribed at these boundary points, we see that the boundary condition will be of the form 1 − ∂ yy α 2 u n+1 − 2u n + u n−1 + β 2 u n (x, y), and the order of the limit and the partial derivatives commute only when partial derivatives of the boundary data can be constructed. But this is the case only when the tangent line at the boundary is in the y-direction, which is precisely why we must restrict our attention to the aforementioned cases for boundary conditions. Upon introducing a second intermediate variable w (2) (x, y), defined by we make the observation that L y [w (2) ] = w (1) . Boundary conditions are now applied to w (2) along the lines x = x j , at the boundary points (x j , c j ) and (x j , d j ), where now no difficulties remain in considering since the right hand side will be fully prescribed. We summarize this procedure in Algorithm 2.

Algorithm 2 Application of Boundary Conditions in 2 dimensions
(1) Initialize temporary variables w (1) and w (2) , which are the same size u n .
(2) For each horizontal line y = y k , for 1 ≤ k ≤ N y , create the temporary variable w k ](x) = β 2 u n + 1 α 2 S n (x, y k ). The homogeneous coefficients are determined by the fact that which is applied at x = a k , b k .
k (x). The homogeneous coefficients are determined by the fact that which is applied at y = c j , d j .
(5) The dimensional splitting error is corrected by adding

Neumann boundary conditions in two dimensions.
In implementing Neumann boundary conditions for boundary geometries conforming to grid lines, such as a rectangular domain, we can directly impose a two-point boundary correction. One way to extend this method to a general polygonal domain would be to use multiple overset grids, each aligned with a boundary segment, which communicate with the interior grid through interpolation on a ghost cell region, though we do not pursue that approach in this work.
For curved boundaries, an alternative approach is that of an embedded boundary method, which involves determining the Dirichlet values at the endpoints of each x-and y-sweep lines that result in the approximate satisfaction of the Neumann boundary condition (in effect, constructing an approximate Neumann-to-Dirichlet map). We present the implementation of an embedded boundary method for Neumann boundary conditions for the implicit wave solver on a curved boundary geometry. The approach taken here follows the work in [18], which proposes an embedded boundary method for Neumann boundary conditions with a finite difference method for the wave equation. The analysis in that work suggests that, on the continuous level, the modified equations and boundary conditions resulting from typical truncation error terms possess unstable boundary layer solutions, so that the addition of a dissipative term is necessary to achieve a stable method. This is consistent with our experience in the implementation described here, with the embedded boundary method becoming unstable when applied to the non-dissipative dispersive solver, but remaining stable for the diffusive solver, which is dissipative.
In the following, we briefly describe the two-point boundary correction method for a gridaligned rectangular boundary. We then describe the embedded boundary method for a 1D problem. This method requires an iterative procedure, which we show, in the setting of the 1D problem, to be a convergent contraction mapping with a rate of convergence that depends on the CFL number. We describe the implementation of the embedded boundary method in 2D, and finally give numerical results.

Description of the Two-Point Boundary Correction Method.
In a rectangular domain where the boundaries conform to grid lines, it is straightforward to impose the two-point boundary correction terms in a line-by-line fashion, since in this case, the grid lines are not coupled through the normal derivative. As this is a simple extension of the 1D boundary correction algorithm, we do not elaborate further.

Description of the Embedded Boundary Method and Proof of Convergence of the Iterative Solution in 1D.
We consider the situation of a one-dimensional domain {x B < x} with a single boundary point not aligned with the grid points, as displayed in Figure 2. We have grid points x 0 , x 1 , ... with uniform grid spacing x i+1 − x i = ∆x, boundary location x B , and ghost point location x G = x 0 . We define interior points to be any grid points lying within the domain (including the boundary), and exterior points to be any grid points lying outside of the domain. We define a ghost point to be any exterior point for which at least one of the neighboring points x i±1 is an interior point. We neglect the right boundary in the present analysis for simplicity, though it can be extended to the case with both boundaries. We consider applying the diffusive version of the wave solver, having calculated the convolution integral I(x), and now needing to find the value of the coefficient A such that the solution u(x) at the next time step is given by Given the value of the convolution integral and the solution at the ghost point, I G = I(x G ) and u G = u(x G ), respectively, the coefficient may be computed as A = u G − I G . We now describe the procedure for determining the value of the solution at the ghost point, u G , that leads to a solution consistent with homogeneous Neumann boundary conditions to second-order accuracy. We construct a quadratic interpolant using the boundary condition and interior interpolation points x I = x B + ∆s I and x II = x B + 2∆s I , lying in between grid points x m and x m+1 , and x n and x n+1 , respectively. The interpolation distance will be chosen such that ∆x < ∆s I < (3/2)∆x. We define the distances ξ G = x B − x G , ξ I = x I − x B = ∆s I , and ξ II = x II − x B = 2∆s I , and construct a quadratic Hermite-Birkhoff [5] interpolant P (ξ) by imposing the conditions P (0) = 0, P (ξ I ) = u I and P (ξ II ) = u II . We then obtain the following second-order approximation to the ghost point value, given by As the coefficients γ I = , we only need supply second-order accurate approximations to u I and u II to maintain overall second-order accuracy. (The coefficients would be O(1/∆x) in the case of nonhomogeneous Neumann boundary conditions, which would require third-order accurate approximations to u I and u II . For simplicity, we consider only the case of homogeneous Neumann boundary conditions in the present work.) Such approximations may be obtained through linear interpolation, giving , and u j = u(x j ) are the values of the function at the uniform gridpoints for j = m, m + 1, n, n + 1.
x 0 x G x B x 1 x m x m+1 x I x n x n+1 x II Figure 2. Boundary geometry in 1D.
Hence, to determine the ghost point value u G that leads to a solution consistent with homogeneous Neumann boundary conditions, we should solve the following system of equations.
with γ I , γ II , σ I and σ II defined as above, and where I j = I(x j ) is the convolution integral evaluated at uniform grid points for j = m, m + 1, n, n + 1. This system can be solved formally for u G , giving To show that this solution formula is well-defined, we argue that for the relevant values of m, n and ξ G , ξ I , and ξ II . We define d = e −α∆x , and noting that 0 < d < 1, m < n, ξ G < ξ I = ∆s I < ξ II = 2∆s I , γ I > 0, γ II < 0, 0 ≤ σ I ≤ 1, and 0 ≤ σ II ≤ 1, we obtain Now, since ∆x < ∆s I < (3/2)∆x, we can see that it is the case that either m = 1 and n = 2 or 3, or that m = 2 and n = 3. It is then a matter of some simple calculus to check that that the functions f m,n (x) = (4x m − x n+1 )/3 satisfy f m,n (x) < 1 for 0 < x < 1 and the given combinations of m and n. This proves that K < 1 for the relevant values of the parameters, so that the solution for u G is well-defined. We note, however, that K approaches 1 as d approaches 1, that is, as the CFL number becomes large. Thus, we may expect an ill-conditioned system when the CFL number is very large.
To obtain the lower bound on K, we observe In the two-dimensional case, the line-by-line solution method couples the ghost point values, and a general explicit solution formula is impossible to write down. In principle, one may write out and directly solve a linear system to obtain the ghost point values. Instead, we propose an iterative solution method that avoids the formation of a matrix. We now describe this iterative solution method and prove its convergence in the context of the one-dimensional problem described above.
Suppose we have the convolution integral evaluated at the gridpoints, I j , and a k-th iterate for the ghost point value, u k G . Then we may obtain the next iterate by the formulas u k+1 ) where quantities are defined as above. Now, to prove the convergence of the interation, we show it is contractive. Taking the difference of two iterates, we have Hence, the Contraction Mapping Theorem implies that the iteration converges to a unique fixed point (which is the solution given in the formula above). We note again that K approaches 1 as the CFL number becomes large, so that the rate of convergence will become slower for larger CFL numbers.

Description of the Method in 2D.
We now describe the implementation of the embedded Neumann boundary condition in the 2D case. We consider the situation displayed in Figure 3, in which we need to determine the value of our unknown u G at the ghost point location (x G , y G ). In the 2D case, we define a ghost point to be any exterior point (x i , y j ) for which at least one of the neighboring points (x i±1 , y j ) or (x i , y j±1 ) is an interior point. Similarly to the 1D case, we will construct a quadratic Hermite-Birkhoff boundary interpolant P (ξ) along the direction normal to the boundary, which intersects the boundary curve Γ at location (x B , y B ), and supply the interior interpolation point values u I and u II , at points (x I , y I ) and (x II , y II ), respectively, by further interpolation from interior grid points. These points are selected along the normal, in analogy to the 1D case, such that ξ I = |(x I , y I ) − (x B , y B )| = ∆s I and ξ II = |(x II , y II ) − (x B , y B )| = 2∆s I , where we will typically take ∆s I = √ 2∆x. We construct a quadratic Hermite-Birkhoff interpolant P (ξ) by imposing the conditions P (0) = 0, P (ξ I ) = u I and P (ξ II ) = u II . Defining further the distance from the boundary to the ghost point ξ G = |(x G , y G ) − (x B , y B )|, we obtain, as in the 1D case, the following second-order approximation to the ghost point value, given by where the coefficients are defined as γ I = In the 2D case, we find approximations to u I and u II through bilinear interpolation. This is in contrast to [18], who find the intersection of the normal with grid lines, then interpolate along the grid lines. We have also implemented a second-order accurate version of this approach and compared to the bilinear interpolation scheme proposed here. We have found that the two schemes behave similarly, however the bilinear interpolation scheme is slightly more accurate and simpler to code, not requiring to handle separate cases of intersection with horizontal and vertical grid lines. The bilinear interpolation scheme is standard, but we give it here for completeness. If the interpolation point u I lies in a cell with corners (x i , y j ), (x i+1 , y j ), (x i+1 , y j+1 ) and (x i , y j+1 ), then we have the following approximation for u I : . With this interpolation scheme established, we now outline the algorithm for the 2D dimensionallysplit wave solver.
Γ (x G , y G ) (x I , y I ) (x II , y II ) (x B , y B ) Figure 3. Boundary geometry in 2D.
The above interpolation procedure applies regardless of the variety of the wave solver that it is used with, provided that the wave solver has sufficient dissipation to maintain stability. We now describe the rest of the embedded boundary algorithm in the context of the diffusive wave solver, though it it may be similarly applied to the the dispersive scheme with artificial dissipation described above. In analogy to the iteration presented in the 1D case, the 2D iterative algorithm proceeds by using the interpolation scheme to provide values at the ghost points, which in turn provide new values for the boundary correction coefficients, which are then used to update the values at the interior grid points, comprising one full iteration. It should be noted that not all interior grid points need be updated in the iteration, only those near the boundary that lie within the boundary interpolation stencils.
Using values from previous time steps, the initial guess for the interior grid points in the boundary interpolation stencils may be given by extrapolation in time, as u n+1,0 = 2u n −u n−1 (linear extrapolation) or u n+1,0 = 3u n − 3u n−1 + u n−2 (quadratic extrapolation). Either extrapolated initial guess provides a modest reduction in the number of iterations required versus a zero initial guess, with only a slight further reduction in the number of iterations going from linear to quadratic extrapolation. An effective stopping criterion for iteration is |u n+1,l+1 − u n+1,l | ∞ < , where is some chosen tolerance, which may be chosen to be quite small, as the iteration is a fixed point interation. In the numerical example, we choose a tolerance of 10 −15 , and we achieve convergence in less than 40 iterations at a CFL number of 2.
In applying the diffusive version of the wave solver, we assume we have previous time steps u n , u n−1 and u n−2 . We have to solve the modified Helmholtz equation with homogeneous Neumann boundary conditions, c∆t . We apply dimensional splitting to find so we define w = 1 − 1 α 2 ∂ yy u, and noting that w = u + O (c∆t) 2 so that ∂w ∂n = ∂u ∂n + O (c∆t) 2 , we obtain the following approximate system We now suppose our domain is embedded in a uniform Cartesian grid, with horizontal grid lines corresponding to y = y k , 1 ≤ k ≤ N y and vertical grid lines corresponding to x = x j , 1 ≤ j ≤ N x . The embedded boundary algorithm will be applied when calculating the intermediate variable w n+1 in horizontal line sweeps as well as the solution variable u n+1 in vertical line sweeps. The iterations for these two variables are separate; first, the iterative procedure is applied to w to convergence, and then this value of w is used to compute u, and the iterative procedure is applied to u to convergence. However, in each iteration, the grid lines are coupled through the interpolation scheme, so that all grid lines must be iterated together. The overall iterative algorithm is described in 3, with details specified for the iteration on w. The iteration on u is very similar, and so we omit the details.

Algorithm 3 Application of Neumann Boundary Conditions in 2 dimensions
(1) (Initialization of ghost points) Perform the interpolation scheme described above to obtain the values of u n , u n−1 , and u n−2 at the ghost points, which are the endpoints of the horizontal and vertical grid lines.
(2) (Evaluation of particular solution) For each horizontal line y = y k , for 1 ≤ k ≤ N y , with ghost (end) points x = a k and b k find the particular solution w n+1 p,k for the intermediate variable w n+1 k (x) by evaluating the discrete convolution operator for each grid point x j in the horizontal line, including the ghost points. to calculate the particular solution for u n+1 , then apply the bounday correction interation.

Double Circle Cavity.
In this example, we solve the wave equation with homogeneous Dirichlet boundary conditions on a 2D domain Ω which is, as in Figure 4, the union of two overlapping disks, with centers P 1 = (−γ, 0) and P 2 = (γ, 0), respectively, and each with radius R: where | (x, y) | = x 2 + y 2 is the usual Euclidean vector norm, and γ < R.
This geometry is of interest due to, for example, its similarity to that of the radio frequency (RF) cavities used in the design of linear particle accelerators, and presents numerical difficulties due to the curvature of, and presence of corners in, the boundary. Our method avoids the staircase approximation used in typical finite difference methods to handle curved boundaries, which reduces accuracy to first order and may introduce spurious numerical diffraction. As initial conditions, we choose |(x,y)−P 2 | 0.8γ 2 | (x, y) − P 2 | < 0.8γ 0 otherwise and u t (x, y, 0) = 0 for (x, y) ∈ Ω. Selected snapshots of the evolution are given in Figure 5, and the results of a refinement study are given in Table 1. The discrete L 2 error was computed against a well-refined numerical reference solution (∆x = 2.1875 × 10 −4 ); the error displayed in the table is the maximum over time steps with t ∈ [0.28, 0.29]. For this example, R = 0.3, γ = 0.2, c = 1, and the CFL is 2. In both cases, the solution converges to second order. An overlay of the two are shown in Figure 6, demonstrating the close agreement.  Figure 6. Two separate numerical constructions of a Bessel mode are superimposed, demonstrating that the solution on the quarter circle using Neumann boundary conditions is equivalent to that of the full circle.

Bessel Mode with Neumann Boundary Conditions.
Here we present a numerical example of the embedded boundary method for homogeneous Neumann boundary conditions given in Section 4.2.2. We apply the method to a circular domain, for which analytical solutions exist. We consider a radially-symmetric Bessel mode with homogeneous Neumann boundary conditions, with an analytic solution given by where J 0 is the Bessel function of the first kind of order 0, r = x 2 + y 2 , R is the radius of the domain, and Z 0 ≈ 3.8317 is the smallest nonzero root of J 0 (so that ∂u ∂n (R, t) = ∂u ∂r (R, t) = Z 0 R J 0 (Z 0 ) cos Z 0 ct R = 0). In this example, we take radius R = π/2 and wave speed c = 1. An example of the embedded boundary grid used is given in Figure 7. We perform a refinement study with a fixed CFL number of 2, with the results in Figure 8 indicating the expected second-order convergence. We set the iteration tolerance to 10 −15 , and we see convergence of the boundary correction iteration in fewer than 40 iterations. We note some oscillation of the L ∞ error about the line giving second-order accuracy, which we believe to be due to the grid points moving with respect to the boundary through the refinement, causing some variation in the maximum error. 5.4. Periodic Slit Diffraction Grating. In this example, we apply our method to model an infinite, periodic diffraction grating under an incident plane wave. Diffraction gratings are periodic structures used in optics to separate different wavelengths of light, much like a prism. The high resolution that can be achieved with diffraction gratings makes them useful in spectroscopy, for example, in the determination of atomic and molecular spectra. Our numerical experiment, depicted in Figure 9, demonstrates the use of our method with multiple boundary conditions and nontrivial geometry in a single simulation to capture complex wave phenomena.
An idealized slit diffraction grating consists of a reflecting screen of vanishing thickness, with open slits of aperture width a, spaced distance d apart, measured from the end of one slit to the beginning of the next (that is, the periodicity of the grating is d). We impose an incident plane wave of the form u inc (x, y, t) = cos (ωt + ky), where k = 2π/a and ω = k/c, where c is the wave speed. Periodic BCs at x = ±d/2 (determining the periodicity of the grating), and homogeneous Dirichlet BCs are imposed at the screen. We also test outflow boundary conditions in multiple dimensions, which are imposed at y = ±L y /2.   In Figure 10, we observe the time evolution of the incident plane wave passing through the aperture, and the resulting interference patterns as the diffracted wave propagates across the periodic boundaries. The outflow boundary conditions allow the waves to propagate outside the domain. While a rigorous analysis of the efficacy of our outflow BCs is the subject of future work, the results look quite reasonable, as no spurious reflections are seen at the artificial boundaries. 5.5. Point Sources and Outflow Boundary Conditions. As a final example, we illustrate some of the more interesting features of our solver. We launch two point sources (from points not on the mesh), and employ Dirichlet, periodic and outflow boundary conditions along the edges of the domain. As can be seen in Figure 11, the point sources propagate perfectly, despite not being placed at grid points.

Conclusion
In this paper we have presented a fast, A-stable, second order scheme for solving the wave equation. Using the method of lines transpose (MOL T ), the solution can be interpreted in the semi-discrete sense as a boundary integral solution, posed as a convolution against an exponential kernel. We have exploited this fact to develop a matrix-free, O(N ) fast spatial convolution algorithm, capable of embedding boundary points in a regular Cartesian mesh, without affecting the accuracy or stability.
In addition to demonstrating second order accuracy for wave propagation in a variety of non-Cartesian geometries, with time steps in excess of the usual CFL restriction, we have also developed a novel method for implementing outflow boundary conditions, as well as methods for launching waves, using soft sources, from points located at arbitrary locations (e.g., not located at a mesh point) inside the domain, which is of interest in particle-wave simulations such as those required in studying plasma. Several topics warrant further investigation. We are developing a domain decomposition approach to multi-core computing with our implicit wave solver based on the properties of the exponential kernel, where the subdomains require only pointwise data communication from adjoining edges. Future work will investigate the implementation of higher-order accurate boundary conditions for complex boundary geometries and outflow boundary conditions, and further, apply these methods to Maxwell's equations both in electromagnetic scattering and plasma physics problems.
Appendix A. Summary Table For Second-Order Wave Solver
We are predominantly interested in the case where S(x, t) consists of a large number of time dependent point sources. However, it is often the case that in electromagnetics problems, a soft source is prescribed to excite waves of a prescribed frequency, or range of frequencies, within the domain. A soft source is so named because, although incident fields are generated at a prescribed fixed spatial location, no scattered fields are generated.
The implementation of a soft source σ(t) at x = x s is accomplished by prescribing the source condition u(x s , t) = σ(t).
(48) However, it can be shown that if we set and insert it into the wave equation, then the soft source condition (48) is satisfied, and the solutions are equivalent. Thus, a soft source is nothing more than a point source, whose time-varying field is integrated by the wave equation. Upon convolving this source term with the Green's function according to (4), we find where the definition of α = β/(c∆t) has been utilized.
Remark 5. It is often the case that taking the analytical derivative σ (t n ) is to be avoided, for various reasons. In this case, any finite difference approximation which is of the desired order of accuracy can be substituted.
Likewise for general point sources, the corresponding form of the source term is Therefore, it suffices to consider delta functions both for the implementation of soft sources, as well as including time dependent point sources.


Fully Connected Deep Structured Networks

Convolutional neural networks with many layers have recently been shown to achieve excellent results on many high-level tasks such as image classification, object detection and more recently also semantic segmentation. Particularly for semantic segmentation, a two-stage procedure is often employed. Hereby, convolutional networks are trained to provide good local pixel-wise features for the second step being traditionally a more global graphical model. In this work we unify this two-stage process into a single joint training algorithm. We demonstrate our method on the semantic image segmentation task and show encouraging results on the challenging PASCAL VOC 2012 dataset.

Introduction
In the past few years, Convolutional Neural Networks (CNNs) have revolutionized computer vision. They have been shown to achieve state-of-the-art performance in a variety of vision problems, including image classification [19,31], object detection [11], human pose estimation [32], stereo [36], and caption generation [15,24,35,8,14,10]. This is mainly due to their high representational power achieved by learning complex, non-linear dependencies.
It is only very recently that convolutional nets have proven also very effective for semantic segmentation [12,30,21,41,3]. This is perhaps due to the fact that to achieve invariance, pooling operations are performed, often reducing the dimensionality of the prediction. A Markov random field (MRF) is then used as a refinement step in order to obtain segmentations that respect well segment boundaries. The seminal work of [17] showed that inference in fully connected MRFs is possible if the smoothness potentials are Gaussian. Impressive performance was demonstrated in semantic segmentation with hand craft features. Later, [3] extended the unary potentials to incorporate convolutional network features. However, these current approaches train the segmentation models in a piece-wise fashion, fixing the unary weights during learning of the parameters of the pairwise terms which enforce smoothness.
In this paper we present an algorithm that is able to train jointly the parameters of the convolutional network defining the unary potentials as well as the smoothness terms taking into account the dependencies between the random variables. We demonstrate the effectiveness of our approach using the dataset of the PASCAL VOC 2012 challenge [9].

Background
We begin by describing how to learn probabilistic deep networks which take into account correlations between multiple output variables y = (y 1 , . . . , y N ) that are of interest to us. Moreover, a valid configuration y ∈ Y = N i=1 Y i is assumed to lie in the product space of the discrete variable domains Y i = {1, . . . |Y i |}.
For a given data sample x ∈ X , and a parameter vector w ∈ R A , the score F of a configuration y ∈ Y is generally modeled by the mapping F : X × Y × R A → R.  The prediction task amounts to finding the configuration which maximizes the score F (x,ŷ; w). Note that the best scoring configuration y * is equivalently given as the maximizer of the probability distribution since the exponential function is a monotone increasing function and the normalization constant is independent of the configurationŷ ∈ Y, i.e., it is constant indeed.
The learning task is concerned with finding a parameter vector which maximizes the likelihood of a given training set D = {(x, y)}. The training set consists of input-output pairs (x, y) which are assumed to be independent and identically distributed. Note that maximizing the likelihood is equivalent to maximizing the cross entropy between the modeled distribution p(ŷ | x, w) and a target distribution which places all its mass on the groundtruth configuration y. Throughout this work we make no further assumptions about the dependence of the scoring function F (x,ŷ; w) on the parameter vector w, i.e., F (x,ŷ; w) is generally neither convex nor smooth.
For problems where the output-space size |Y| = N i=1 |Y i | is in the thousands, we can exactly solve the inference task given in Eq. (1) by searching over all possible output space configurationsŷ ∈ Y. In such a setting, those different configurations are typically referred to as different classes. Similarly, we normalize the distribution p(ŷ | x, w) by summing up the exponentiated score exp F (x,ŷ; w) over all possibilitiesŷ ∈ Y. This is often referred to as a soft-max computation. Non-convexity and non-smoothness of the learning objective w.r.t. the parameters w is answered with stochastic gradient ascent. For efficiency, the gradient is often computed on a small subset of the training data, i.e., a mini-batch.
We summarize the resulting training algorithm in Fig. 1. On a high level it consists of four steps which are iterated until a stopping criterion is met: (i) the forward pass to compute the scoring function F (x,ŷ; w) for all output space configurationsŷ ∈ Y. (ii) normalizing the scoring function via a soft-max computation to obtain the probability distribution p(ŷ | x, w). (iii) computation and back-propagation of the gradient of the loss function, i.e., often the log-likelihood or equivalently the cross-entropy. (iv) an update of the parameters.
However, solving the inference task given in Eq. (1) or the learning problem stated in Eq. (2) is computationally challenging if we consider more complex output spaces Y, e.g., those arising from tasks like image tagging. The situation is even more severe if we target image segmentation where the exponential number of possible output space configurations prevents even storage of F (x,ŷ; w) ∀ŷ ∈ Y. Note that this is required in the first line of the algorithm summarized in Fig. 1.
Given an exponential amount of possible configurations |Y| = N i=1 |Y i |, how do we represent the scoring function F (x,ŷ; w) efficiently? Assuming we have an efficient representation, how can we effectively normalize the probability p(ŷ | x, w)? One possible answer to those questions was given by Chen et al. [4], who discussed extending log-linear models, i.e., those with a scoring function of the form F (x,ŷ; w) = w φ(x,ŷ), to the more general setting, i.e., an arbitrary dependence of the scoring function F (x,ŷ; w) on the parameter vector w.
In short, [4] assumed the global scoring function F (x,ŷ; w) to decompose into a sum of local scoring functions f r , each depending on a small subset r ⊆ {1, . . . , N } of variablesŷ r = (ŷ i ) i∈r . All restrictions r required to compute the global function via are subsumed in the set R. If the size of each and every local restriction set r ∈ R is small, To compute the gradient of the log-likelihood cost function, we require a properly normalized distribution p(ŷ | x, w), or more specifically its marginals b (x,y),r (ŷ r ) for each restriction r ∈ R. To this end, message passing type algorithms were employed by [4]. Such an approach is exact if the distribution p(ŷ | x, w) is of low tree-width. Otherwise computational complexity is prohibitively large and approximations like loopy belief propagation [26], convex belief propagation [39] or treereweighted message passing [37] are alternatives that were successfully applied.
The resulting iterative method of [4] is summarized in Fig. 2. In a first step the forward pass computes all outputs of every local scoring function. Afterwards (approximate) marginals are obtained in a second step, and utilized to compute the derivative of the (approximated) maximum likelihood cost function w.r.t. the parameters w. The following backward pass computes the gradient of the parameters by repeatedly applying the chain-rule according to the definition of the scoring function F (x,ŷ; w). The gradient is then utilized during the final parameter update.
Not only does the approach presented by [4] fail if the decomposition assumed in Eq. (3) is not available. But it is also computationally challenging to obtain the required marginals if too many local functions are required. I.e., computation is slow if the number of restrictions |R| is large, e.g., when working with densely connected image segmentation models where every pixel is possibly correlated to every other pixel in the image.

Approach
Densely connected models were previously considered by [17,33,34,18] and shown to yield impressive results for the image segmentation task. Learning the parameters of densely connected models was considered by Krähenbühl and Koltun [18] in the context of the log-linear setting. Following [4] we aim at extending those fully connected log-linear models to the more general setting of an arbitrary function F (x,ŷ; w), e.g., a deep convolutional neural network. Note that a similar approach has been recently discussed by [41] in independent work.
Let us consider within this section how to efficiently combine deep structured prediction [4] with densely connected probabilistic models [17,33,34,18]. Before getting into the details we note that the presented approach trades computational complexity of the general method of [4] with a restriction on the pairwise functions f ij (i.e., r = {i, j}). Concretely, the local functions f ij are assumed to be mixtures of kernels in a feature space as detailed below. For simplicity we assume that local functions of order higher than two are not required to represent our global scoring function F (x,ŷ; w). Generalizations have however been presented, e.g., by Vineet et al. [34].

Inference
We begin our discussion by considering the inference task. To obtain a computationally efficient prediction algorithm we use a mean field approximation of the model distribution p(ŷ | x, w) for every sample (x, y). More formally, we assume our approximation to factor according to Given some parameters w, we employ a forward pass to obtain our local function representations f r (x,ŷ r ; w). Next we compute the single variable marginals q (x,y),i (ŷ i ) by minimizing the Kullback-Leibler (KL) divergence w.r.t. to the assumed factorization of the mean field distribution q (x,y) (ŷ), i.e., Hereby q ∈ ∆ requires q to be a valid probability distribution. Due to non-convexity, only convergence to a stationary point of the KL divergence cost function is guaranteed for sequential block-coordinate updates [38,16]. More precisely, iterating until convergence through the variables i ∈ {1, . . . , N } using the closed form update which assumes all marginals but q (x,y),i to be fixed, retrieves a stationary point for the cost function of the program given in Eq. (4). The set of variables neighboring i is denoted N (i).
In the case of densely connected variables, the computational bottleneck arises from the second summand which involves j∈N (i) |Y j | additions. The sum ranges over |N (i)| = N − 1 terms for densely connected structured models. Hence the complexity of an update for a single marginal is of O(N ), and updating all N marginals therefore requires O(N 2 ) operations as also discussed by Krähenbühl and Koltun [18].
Importantly, Krähenbühl and Koltun [17] observed that a high dimensional Gaussian filter can be applied to concurrently update all marginals in O(N ). This is achievable when constraining ourselves to pairwise functions being mixtures of M kernels in the feature space as mentioned before.
Formally, we require where µ (m) is a label compatibility function, k (m) is a kernel function, andf i (x) are features of variable i depending on the data x.
However, to ensure convergence to a stationary point of the KL divergence cost function for this parallel update, further restrictions on the form of the pairwise functions f ij apply. Formally, if the label compatibility functions µ (m) are negative semi-definite ∀m, and the kernels k (m) are positive definite ∀m, the KL divergence is readily given as the difference between a concave and a convex term [18]. Hence the concave-convex procedure (CCCP) [40] is directly applicable. We therefore proceed iteratively by first linearizing the concave term at the current location and second minimizing the resulting linearized but convex program.
As detailed by Krähenbühl and Koltun [18], and as discussed above, finding the linearization is equivalently solved via filtering in time linear in N . Solving the convex program in its original form requires solving a non-linear system of equations independently for each marginal q (x,y),i (ŷ i ), e.g., via Newton's method. A further approximation to the cross-entropy term of the KL-divergence relates the efficient filtering based mean field update of the marginals q (x,y),i (ŷ i ) to the corresponding cost function for which a stationary point is found.

Learning
Having observed that mean-field inference can be efficiently addressed with Gaussian filtering, given restrictions on the pairwise functions f ij , we now turn our attention to the learning task. As mentioned before we aim at finding a parameter vector w that maximizes the likelihood objective function. Since the exact likelihood is computationally expensive, we use the log-likelihood based on

Algorithm: Learning Fully Connected Deep Structured Models
Repeat until stopping criteria 1. Forward pass to compute f r (x,ŷ r ; w) ∀r ∈ R, y r ∈ Y r 2. Computation of marginals q t (x,y),i (ŷ i ) via filtering for t ∈ {1, . . . , T } 3. Backtracking through the marginals q t (x,y),i (ŷ i ) from t = T − 1 down to t = 1 4. Backward pass through definition of function via chain rule 5. Parameter update the mean-field marginals. Hence our surrogate loss function L (x,y) for a sample (x, y) with corresponding annotated ground truth labeling y is given by To perform a parameter update step we need the gradient of the surrogate loss function w.r.t. the parameters, i.e., The gradient of the surrogate loss function L (x,y) w.r.t. the marginals is easily obtained from Eq. (6). It is given by where the Iverson bracket ŷ i = y i equals one ifŷ i = y i , and returns zero otherwise.
To perform a gradient step during learning, we additionally require the derivatives of the marginals w.r.t. the parameters, i.e., ∂q (x,y),i (ŷi) ∂w .
More carefully investigating the mean-field update given in Eq. (5) reveals a recursive definition.
More concretely, the derivative ∂q t (x,y),i (ŷi) ∂w of the marginal q t (x,y),i (ŷ i ) after t iterations depends on the results from earlier iterations. Hence, we obtain the desired result by successively back-tracking through the mean-field iterations from the last iteration back to the first. This direct computation is however computationally expensive. Fortunately, back-substitution into the loss gradient yields an algorithm which requires a total of T back-tracking steps, independent of the number of parameters. We refer the interested reader to [18] for additional details regarding the computation of the gradient But contrasting [18], we no longer assume the unaries to be given by a logistic regression model. Contrasting [3], we don't assume the unaries to be fixed during CRF parameter updates. Generalizing the gradient of the marginals w.r.t. parameters to arbitrary unaries is straightforward since the gradients are directly given by the marginals. Combined with the gradient of the log-likelihood loss function w.r.t. the marginals, given in Eq. (8), we obtain ∂L (x,y) ∂w as the difference between the ground-truth and the predicted marginals. This result is then used for back-propagation through any functional structure which provides the unary scoring functions f i , e.g., convolutional neural networks.
Derivatives w.r.t. to label compatibility and kernel shape parameters are readily given in [18]. The resulting algorithm is summarized in Fig. 3. In short, we first obtain again our functional representation via a forward pass through any functional network. Subsequently we compute our meanfield marginals via filtering. Afterwards we obtain the gradient of the loss function via an efficient back-tracking. In the next step the gradient of the parameters is computed by back-propagating the gradient of the loss-function using the chain-rule dictated by the definition of the scoring function. In a final step we update the parameters.

Experiments
We evaluate our approach summarized in Fig. 3 on the dataset of the Pascal VOC 2012 challenge [9]. The task is semantic image segmentation of 21 object classes (including background). The original dataset contains 1464 training, 1449 validation and 1456 test images. In addition to this data we make use of the annotations provided by Hariharan et al. [13], resulting in a total of 10582 training instances. The reported performance is measured using the intersection-over-union metric. Note that we conduct our tests on the 1449 validation set images which were neither used during training nor for fine-tuning.

Model
Our model setup follows [3], i.e., we employ the 16 layer DeepNet model [31]. Just like [3] we first convert the fully connected layers into convolutions as first discussed in [12,30]. This is useful since we are not interested in a single variable output prediction, but rather aim at learning probability masks. To obtain a larger probability mask we skip downsampling during the last two max-pooling operations. To take into account the skipped downsampling during subsequent convolutions we employ the 'à trous (with hole) algorithm' [23]. It takes care of the fact that data is stored in an interleaved way, i.e., in our case convolutions sub-sample the input data by a factor of two or four respectively. To adapt to the 21 object classes we also replace the top layer of the DeepNet model to yield 21 classes for each pixel.
Similar to [3] we assume the input size of our network to be of dimension 306 × 306 which results in a 40 × 40 sized spatial output of the DeepNet which is in our case an intermediate result however.
Contrasting [3], we jointly optimize for both unary and CRF parameters using the algorithm presented in Fig. 3. To this end, given images downsampled to a size of 306 × 306, our algorithm first performs a forward pass through the convolutional DeepNet to obtain the 40 × 40 × 21 sized class probability maps in an intermediate stage. These intermediate class probability maps are directly up-sampled to the original image dimension using a bi-linear interpolation layer. This yields the actual output of our augmented DeepNet network defining the scoring function F (x,ŷ, w). Note that the number N of variablesŷ = (ŷ 1 , . . . ,ŷ N ) ∈ Y is therefore equal to the number of pixels of the original image.
For the second step of our algorithm we perform 5 iterations of mean field updates to compute the marginals q (x,y),i (ŷ i ) of the fully connected CRF. Those are then compared to the original groundtruth image segmentations, using as our loss function the sum of cross-entropy terms, i.e., the log-likelihood loss, as specified in Eq. (6). In the third step we back-track through the marginals to obtain a gradient of the loss function. Afterwards we back-propagate the derivatives w.r.t. the unary term through both the bi-linear interpolation and the 16-layer convolutional network. The shape and compatibility parameters of the CRF, detailed below, are updated directly.  Table 1: Performance of our approach for individual classes. In the last two columns of the lower panel we compare our mean to the recently presented baseline by Chen et al. [3].
It was shown independently by many authors [31,4], that successively increasing the number of parameters during training typically yields better performance due to better initialization of larger models. We therefore train our model in two stages. First, we assume no pairwise connections to be present, i.e., we fine-tune the weights obtained from the DeepNet ImageNet model [31,29] to the Pascal dataset [9]. Standard parameter settings for a momentum of 0.9, a weight decay of 0.0005 and learning rates of 0.01 and 0.001 for the top and all other layers are employed respectively. Due to the 12GB memory restrictions on the Tesla K40 GPU we use a mini-batch size of 20 images.
In a second stage we jointly train the convolutional network parameters as well as the compatibility and shape parameters of the dense CRF arising from the pairwise functions Hereby, we employ the Potts potential µ(y i , y j ) = y i = y j and the Gaussian kernels given by As indicated in Eq. (9), we use M = 2 kernels, both with diagonal covariance matrix Σ m . One containing as featuresf i (x) the two-dimensional pixel positions, the other one containing as features the two dimensional pixel positions as well as the three color channels. Hence we obtain a total of nine parameters, i.e., two compatibility parameters w 1 and w 2 and 2+5 = 7 kernel shape parameters for the diagonal covariance matrices Σ m .

Results
As mentioned before, all our results were computed on the validation set of the Pascal VOC dataset. This part of the data was neither used for training nor for fine-tuning.
Unary performance: We first investigate the performance of the first training stage of the proposed approach, i.e., fine-tuning of the 16 layer DeepNet parameters on the Pascal VOC data. The validation set accuracy is plotted over the number of iterations in Fig. 4 (a). We observe the performance to peak at around 4000 iterations with a mean intersection over union measure of 61.476%. The result reported by [3] for this experiment is 59.80%, i.e., we outperform their unary model by 1.5%.
Joint training: Next we illustrate the performance of the second step, i.e., joint training of both convolutional network parameters and CRF compatibility and shape parameters. In Fig. 4 (b) we indicate the best obtained unary performance from the first step and visualize the validation and training set performance over the number of iterations. We observe the results to peak quickly after around 20 iterations and remain largely stable thereafter.
Details: In Tab. 1 we provide the training and test set accuracies for the 21 individual classes. We observe the 'bike' and 'chair' class to be particularly difficult. For both categories the validation set performance is roughly half of the training set accuracy.
Comparison to baseline: As provided in Tab. 1, the peak validation set performance of our approach is 64.060%, which slightly outperforms the separate training result of 63.74% reported by Chen et al. [3]. Visual results: We illustrate visual results of our approach in Fig. 5. Our method successfully segments the object if the images are clearly apparent. Noisy images and objects with many variations pose challenges to the presented approach as visualized in Fig. 6. Also, we observe our learnt parameters to generally over-smooth results while being noisy on the boundaries.

Discussion
We presented a first method that jointly trains convolutional neural networks and fully connected conditional random fields for semantic image segmentation. To this end we generalize [3] to joint training. Note that a method along those lines has also been recently made publicly available in independent work [41]. Whereas the latter combines dense conditional random fields [17] with the fully convolutional networks presented by Long et al. [21], we employ and modify the 16 layer DeepNet architecture presented in work by Simonyan and Zisserman [31].
Ideas along the lines of joint training were discussed within machine learning and computer vision as early as the 90's in work done by Bridle [2] and Bottou [1]. More recently [5,27,22,6,28,25] incorporate non-linearities into unary potentials but generally assume exact inference to be tractable. Even more recently, Li and Zemel [20] investigate training with hinge-loss objectives using nonlinear unaries, but the pairwise potentials remain fixed, i.e., no joint training. Domke [7] decomposes the learning objective into logistic regressors which will be computationally expensive in our setting. Tompson et al. [32] propose joint training for pose estimation based on a heuristic approximation which ignores the normalization constant of the model distribution. Joint training of conditional random fields and deep networks was also discussed recently by [4] for graphical models in general. Techniques based on convex and non-convex approximations were described for obtaining marginals in the general non-linear setting.

Conclusion
We discussed a method for semantic image segmentation that jointly trains convolutional neural networks and conditional random fields. Our approach combines techniques from deep convolutional neural networks with variational mean-field approximations from the graphical model literature. We obtain good results on the challenging Pascal VOC 2012 dataset.
In the future we plan to train our method on larger datasets. Additionally we want to investigate training with weakly labeled data.


Hadoop, MapReduce and HDFS: A Developers Perspective☆

Abstract The applications running on Hadoop clusters are increasing day by day. This is due to the fact that organizations have found a simple and efficient model that works well in distributed environment. The model is built to work efficiently on thousands of machines and massive data sets using commodity hardware. HDFS and MapReduce is a scalable and fault-tolerant model that hides all the complexities for Big Data analytics. Since Hadoop is becoming increasingly popular, understanding technical details becomes essential. This fact inspired us to explore Hadoop and its components in-depth. The process of analysing, examining and processing huge amount of unstructured data to extract required information has been a challenge. In this paper we discuss Hadoop and its components in detail which comprise of MapReduce and Hadoop Distributed File System (HDFS). MapReduce engine uses JobTracker and TaskTracker that handle monitoring and execution of job. HDFS a distributed file-system which comprise of NameNode, DataNode and Secondary NameNode for efficient handling of distributed storage purpose. The details provided can be used for developing large scale distributed applications that can exploit computational power of multiple nodes for data and compute intensive applications.

Introduction
Nowadays data is being generated with a very high rate from different fields like business, scientific data, e-mails, blogs, etc. To analyse and process this huge amount of data and to extract meaningful information for users there is a need of deploying data intensive application and storage clusters 1 . Requirements for this type of application are fault tolerance; parallel processing, data-distribution, load balancing, scalability and highly availability. To deal with such type of problems Google introduced the MapReduce programming model 2 . Apache Hadoop is an open source implementation of MapReduce system 3 .

Hadoop is the most popular and open source implementation of MapReduce programming model. Apache
Hadoop is a software framework for reliable, scalable, parallel and distributed computing. In lieu of relying on costly hardware and expensive systems for processing and storing data, Apache Hadoop empowers parallel processing on Big Data on commodity hardware 4 . HDFS and MapReduce programming standard grant itself to these Big Data intensive analytic jobs because of its scale-out architecture and its ability to process data in parallel manner in multi-node clusters. MapReduce, and its existing open-source project called Hadoop, enables parallel processing of huge amount of data and automatic partition of data, data distribution, fault tolerance and load balancing management which finally yields reliable and scalable computing. The rapid and high growth rate of data first posed challenges on big companies like Facebook, Google, Amazon and Yahoo. These companies need to execute and carry out terabytes and petabytes of data on daily basis to deduce demands and queries of their users. Existing and traditional tools and application becomes deficient to process large amount of data. Hadoop explained and answered the problem of handling and processing such terabytes and petabytes of data 5 . Many enterprises, industries and universities works on parallel and distributed computing but for neophytes of Big Data and its processing and analysis is not an easy task and it requires much effort to manage due to the lack of experience and insufficient money to invest in super computers and servers. For such neophytes it is important to utilize limited resources in effective and efficient manner.
In this paper the section II provide the explanation and relation of cloud computing with hadoop and why and how big organization exploiting hadoop along with cloud computing. Section III comprised of Hadoop Architecture which explains all the daemons of Hadoop including JobTracker, TaskTracker, NameNode, DataNode and Secondary NameNode in detail followed by conclusion.

Cloud computing with hadoop
Enterprises like Google, Facebook other Internet colossal organizations process user demands, queries and workloads. Millions and billions of queries and requests are served every hour on Internet. Therefore, migrating towards cloud computing to take advantage of reduced costs and improved performance is becoming need of an hour 6 . The key is to select the applications and technology that is best suited for organization's mission, infrastructure, and long-term service. Apache Hadoop offers a broad selection of effective cloud deployable solutions to assist organizations in achieving its goal and in successfully migrating to the cloud.
These organizations depend on Linux servers and cloud computing for getting the adequate performance in terms of scalability and availability 7 . The pliability of Linux merged with smooth and consistent scalability of cloud environment makes it capable of supplying the ideal framework for analysing and processing Big Data, while eliminating the need for costly hardware and software. Hadoop is seen as a preferred choice in open source cloud computing community for providing an efficient platform for Big Data processing. It has become obvious that Apache Hadoop in cloud computing is now an interesting topic because cloud computing is regarded as the next milestone of the IT industry 8 .

Hadoop Architecture
Apache Hadoop comprise of five different daemons and each of these daemons run its own JVM-1) NameNode, 2) DataNode, 3) Secondary NameNode, 4) JobTracker and 5)TaskTracker 9 . Daemons which stores data and metadata i.e., NameNode and DataNode, comes under HDFS layer and JobTracker and TaskTracker, as shown in Fig. 1, which keeps track and actually executes the job, comes under MapReduce layer. Hadoop cluster comprised of multiple slave nodes and a master node. The master node runs the master daemons for each layer i.e., NameNode for the HDFS storage layer, and JobTracker for the MapReduce processing layer. Rest of the machines will run the "slave" daemons: DataNode for the HDFS layer and TaskTracker for MapReduce layer. Master node can also play a role of slave. Thus, in addition to master daemons, master node can run the slave daemons as well. Fundamentally, the daemons running on master node take responsibility for coordinating and managing the slave daemons on all nodes which carryout work for data storage and processing 10,11 .
Hadoop core architecture comprise of -

HDFS (Hadoop Distributed File System)
-HDFS is a self-healing, distributed file system that provides reliable, scalable and fault tolerant data storage on commodity hardware 12,13 It works closely with MapReduce by distributing storage and computation across large clusters by combining storage resources that can scale depending upon requests and queries while remaining inexpensive and in budget 14,15 . HDFS accepts data in any format like text, images, videos, etc regardless of architecture and automatically optimizes for high bandwidth streaming.
The foremost advantage of HDFS is fault tolerance. By provisioning fast data transfer between the nodes and enabling Hadoop to continue to provide service even in event of node failures decreases the risk of catastrophic failure 16 . HDFS is also capable of providing scale-out storage solution for Hadoop. a) NameNode: A single NameNode daemon runs on master node. NameNode stores and manages metadata about the file system in a file named fsimage. This metadata is cached in main memory to provide faster access to the clients on read/write requests. The NameNode controls the slave DataNode daemons to execute the I/O tasks.
The NameNode manages and controls how files are broken down into blocks, identifies which slave node should store these blocks, along with the overall condition and fitness of the distributed file system. The tasks carried out by the NameNode are memory and I/O intensive. b) DataNodes: Hadoop cluster comprise of DataNode daemon that runs on every slave node. DataNodes are primary storage elements of HDFS that store data blocks and service read/write requests on files stored on HDFS. These are controlled by NameNode. Blocks stored in DataNodes are replicated as per the configuration to provide reliability and high availability. These replicated blocks are distributed across the cluster to provide rapid computation. c) Secondary NameNode: Secondary NameNode is not a backup for the NameNode. The Secondary NameNode's job is to periodically read the file system, log the changes and apply them into the fsimage file. This helps in updating NameNode to start up faster next time as shown in Fig. 3.  A MapReduce job generally breaks and divides the input data into chunks which are first processed by "Map phase" in parallel and then by "Reduce phase" 19,20,21 . Hadoop framework sorts out the output of the Map phase which are then given as an input to Reduce phase to initiate parallel reduce tasks as shown in Fig. 4. These input and output files are stored in file system. By default, the MapReduce framework gets input datasets from HDFS 22,23 . It is not necessary that Map and Reduce tasks proceed in a sequential manner i.e., reduce tasks can begin as soon as any of the Map task completes its assigned work. It is also not necessary that all Map tasks completes before any reduce tasks starts working. MapReduce works on key-value pairs. Conceptually, a MapReduce task takes input data set as key-value pair and gives output in the form of key-value pair only by processing input data sets through MapReduce phases. Output generated by the Map phase is called intermediate results which are given as an input to reduce phase as shown in Fig. 5.

Map
Reduce Fig 5: MapReduce key-value pairs Similar to HDFS, MapReduce also exploits master/slave architecture in which JobTracker daemon runs on master node and TaskTracker daemon runs on each salve node as shown in Fig. 6. MapReduce processing layer comprised of two daemonsa) JobTracker: The JobTracker service runs on master node and monitors MapReduce tasks executed by TaskTracker on salve nodes. User in interaction with the Master node submits the job to the JobTracker. JobTracker then asks NameNode for the actual location of data in HDFS to be processed. JobTracker locates TaskTracker on slave nodes and submits the jobs to TaskTracker on slave nodes. The TaskTracker sends heartbeat message back to JobTracker periodically to ensure that TaskTracker on a particular slave node is alive and executing its allotted job. If the heartbeat message is not received within a stipulated time then TaskTracker on a particular salve node is considered to be non-functional and the assigned job is scheduled on another TaskTracker. JobTracker and TaskTracker are known as the MapReduce engine. JobTracker is a single point-of-failure for Hadoop MapReduce service, if it goes down all executing jobs will be stopped. b) TaskTracker: A TaskTracker daemon runs on slave nodes of a cluster. It accepts jobs from JobTracker and executes MapReduce operations. Each TaskTracker has a finite number of "task slots" based on the ability of a node. The heartbeat protocol permits JobTracker to know how many "task slots" are available in TaskTracker on a slave node. It is the task of JobTracker to allocate appropriate jobs to the particular TaskTracker depending upon how many free task slots are available.
TaskTracker governs the execution of every MapReduce operation on each slave node. Although, there is only one TaskTracker per slave node, each TaskTracker can start multiple JVM's to execute MapReduce operation in parallel. TaskTracker from every slave node sends heartbeat message to JobTracker, master node, periodically to ensure JobTracker that TaskTracker is still alive.

Conclusion
Hadoop MapReduce programming paradigm and HDFS are increasingly being used for processing large and unstructured data sets. Hadoop enables interacting with the MapReduce programming model while hiding the complexity of deploying, configuring and running the software components in the public or private cloud. Hadoop enables users to create cluster of commodity servers. MapReduce has been modelled as an independent platform-asa-service layer suitable for different requirement by cloud providers. It also enables users to understand the data processing and analysing.