Numerical method for hyperbolic conservation laws via forward backward SDEs

It is well known that for solutions of semi-linear parabolic PDEs, there are equivalent probabilistic interpretations, which yields the so called nonlinear Feymman-Kac formula. By adopting such formula, we consider in this work a novel numerical approach for solutions of hyperbolic conservation laws. Our numerical method consists in efficiently computing the viscosity solutions of conservation laws. However, instead of solving the viscosity problem directly (which is difficult), we find its equivalent probabilistic solution by adopting the Feymman-Kac formula, which relies on solving the equivalent forward backward stochastic differential equations. It is noticed that such framework possesses the following advantages: (i) the viscosity parameter can be chosen sufficiently small (say $10^{-10}$); (ii) the computational procedure on each discretized time level can be \textit{completely parallel}; (iii) the traditional CFL condition is dramatically weakened; (iv) one does not need to handle the transition layers and discertizations of derivatives. Thus, high accuracy viscosity solutions can be efficiently found. Several numerical examples are given to demonstrate the effectiveness of the proposed numerical method.

Introduction
Many problems in physics and engineering are modeled by hyperbolic systems of conservation laws. For example, the shallow water equations of hydrology, the Euler equations for inviscid and compressible flow, and the Magnetohydrodynamic equations of plasma physics, see [7,13]. The general initial value problems for hyperbolic conservation laws yield where U : R d → R m is the vector of unknowns and F j : R m → R m is the flux vector for the jth direction with m being a positive integer. It is well known that solutions of (1.1) develop discontinuities in finite time even when the initial data is smooth. This holds true even for the scalar case and solutions to (1.1) are sought in the weak sense. Furthermore, weak solutions are augmented by additional admissibility criteria or entropy conditions [7] in order to ensure uniqueness.
There has been a large amount of work on the classical numerical methods for solving hyperbolic conservation laws [13,19], starting from the Godunov scheme [14], to the Lax-Wendroff scheme [20] and the upwind scheme [31], to the total variation bounded high-order essentially non-oscillatory (ENO) scheme [15][16][17] and the weighted essentially non-oscillatory method [32], to some recent approaches such us the moving mesh method (cf. [34]) and the (local) discontinuous Galerkin methods (cf. [37,43]), to name a few.
Another classical approach for conservation laws is the viscosity method, in which one introduces a diffusion term into the original equation and obtains a semi-linear parabolic equation which admits a unique smooth solution, precisely, (1.2) Here, ǫ is a small parameter. It is shown that for certain conservation laws [4,6,9] the limit of the smooth solution u ǫ (ǫ → 0) is the viscosity solution of the original problem (1.1). The viscosity method plays an important role in theoretical analysis for hyperbolic conservation laws. However, the direct simulation for (1.2) is tough due to sharp transition layers, and thus, it is hard to compute the viscosity solution with high accuracy. For details about the viscosity method, please refer to [4,6,9,36] and references therein.
In this paper, we will essentially extend the Crank-Nicolson type scheme introduced in [39] for BSDEs to solve the corresponding weakly coupled nonlinear FBSDEs. It is noted that in such framework (compared to an direct solver for (1.2)), one does not need to handle the transition layers and discertizations of the derivatives. Furthermore, such approach possesses the following features: (i) the viscosity parameter can be chosen sufficiently small (say 10 10 ); (ii) the computational procedure on each discretized time level can be completely parallel ; (iii) the traditional CFL condition is dramatically weakened. To demonstrate the effectiveness of our method, several model problems such as the linear advection equation, the Burgers' equation, the Buckley-Leverett equation, and the Euler system are tested.
We summarize here the new features of our discussions: (i) we provide an alternative way to compute the viscosity solution of conservation laws; (ii) we extend the Crank-Nicolson type scheme introduced in [39] for BSDEs to solve the resulting weakly coupled nonlinear FBSDEs. Despite the large amount of the literatures, such problems have not been studied.
The rest of this paper is organized as follows. In Section 2, by introducing the nonlinear Feymman-Kac formula, we describe the equivalent probabilistic solution of the viscosity solution. Such probabilistic solution can be obtained by solving the corresponding FBSDEs. The numerical approach for FBSDEs will be discussed in Section 3. Several numerical tests are given in Section 4 to demonstrate the effectiveness of our numerical method. We finally given some conclusions in Section 5.

Conservation laws, viscosity solution, and the nonlinear Feymman-Kac formula
Without loss of generality, we consider here the following one dimensional hyperbolic system of conservation law where u = (u 1 , ..., u m ) ⊤ is the vector of unknowns. As we stated, the main purpose of this work is to find the following viscosity solution of (2.1) where ǫ is a small positive number. For scalar conservation laws, u ǫ converges to an entropy solution when ǫ goes to zero (see e.g. [4,6,9]). Furthermore, consider the following Hamilton-Jacobi problem with strictly convex hamiltonian H and compactly supported (or periodic) C 2 initial data u ǫ | t=0 = u 0 (x), there holds [21]: Theorem 2.1. Let {u ǫ } ǫ be the family of approximate viscosity solutions of (2.3), and u be the solution of the original problem without the viscosity term. For any fixed T , there exists C(T ) such that the following error estimate holds The viscosity method plays an important role in the theoretical analysis for conservation laws. However, solving equation (2.2) with high accuracy is a difficult work due to the sharp transition layers and the nonlinearity. Therefore, instead of solving (2.2) directly, we consider in this work solving its equivalent probabilistic solution. To do this, we first review some related results of the nonlinear Feymman-Kac formula.

Nonlinear Feymman-Kac formula
The nonlinear Feymman-Kac formula reveals that for solutions of certain semilinear parabolic PDEs, there are equivalent probabilistic interpretations, which yield the representations of the solutions of FBSDEs. To begin, let us first introduce the following coupled FBSDEs Under certain regularity conditions on the data b, σ, f and ϕ, the adapted solution (Y t , Z t ) of (2.4) can be represented in the following form [28,30]: where u(t, x) is the smooth solution of the following parabolic partial differential equation is the solution of the FBSDEs (2.4), then u defined by (2.5) is the viscosity solution of (2.6).
The representations in (2.5) are well known as the nonlinear Feynman-Kac formula.

Probability representation of the viscosity solution
Consider the viscosity equation in the general form of By the Feymman-Kac formula (2.5), the solution u ǫ of the viscosity equation (2.8) can be represented as where (Y τ , Z τ ) is the solution of the following FBSDE . Our aim of this paper is to numerically solve the associated FBSDEs (2.9) to obtain u ǫ (τ, X τ ) = Y τ (Note that we do not need to solve Z τ ). We remark that there has been a few works on using FBSDEs to study PDE problems (cf. [2,11] and references therein). However, the PDE models considered and the corresponding numerical methods derived in these works are totally different with the results presented in this paper.

Discretization of the FBSDEs
In this section, we discuss the discretizations of the FBSDEs (2.9). We remark that one can find related information on such numerical approach in e.g. [39][40][41].
We introduce the following equally distributed time partition, We define the density of the grids in D n h by 2. the number of elements in D n h,x is finite and uniformly bounded, that is, there exists a positive integer N e , such that, #D n h,x ≤ N e . We call D n h,x the neighbor grid set in D n h at x.

The reference equations
By the BSDE in (2.9), we have for 0 ≤ n ≤ N − 1 : Then, taking the conditional mathematical expectation E x tn [·] on both sides of (3.2) yields The integrand E x tn [f (s, X s , Y s )] on the right-hand side of (3.3) is a deterministic smooth function of time s under F tn , and we may use some numerical integration methods to approximate the integral in (3.3). In particular, we approximate this integral by the trapezoidal rule, namely, where R n y is the truncation error.
For the FSDE of X t , we introduce the following scheme where R n x is the error introduced by the scheme.
Remark 3.1. The equations of (3.5) and (3.6) are our reference equations of the FBSDEs, and our numerical method will be designed based on them.

The semi-discrete scheme
Based on the reference equations (3.5) and (3.6), let Y n and X n be the numerical approximations of Y tn and X tn respectively, we propose a Crank-Nicolson type numerical scheme to solve Y n : (3.7) From Scheme 1, it is clear that Y n is a function of X n . Note that Scheme 1 is nonlinear for Y n . To solve it, some iteration methods are needed. In this paper, we use the following iteration method: , and set a tolerance parameter δ.
2. For l=0, 1, ..., solve Y l+1 n . Note that we have used Y n+1 (X n ) as the initial gauss for Y 0 n , and it is noted that in such a framework, for a reasonable δ (say δ = 10 −5 ), about 2-3 iterations are enough for convergence.

Remark 3.2.
Although the main purpose of this paper is to provide an alternative way to compute the viscosity solutions of conservation laws, our algorithms Scheme 2 is, to the best of our knowledge, new for solving FBSDEs .

The fully discrete scheme
Note that Scheme 1 is a semi-discrete scheme. For any fixed t n , to solve Y (t n , X tn ), the space partition of X tn is needed, and the conditional expectations should be approximated at each grid point x ∈ D n h . We denote such an approximation by E x,h tn [·|X tn = x] ≈ E tn [·|X tn = x]. Generally speaking, after taking the conditional expectation E tn [·|X tn = x], X n+1 in Scheme 1 does not belong to D n+1 h for x ∈ D n h . Thus, interpolation methods are needed to approximate the value of Y n+1 at X n+1 by the values of Y n+1 on D n+1 h . Here, we will adopt a local interpolation operator I n h,x such that I n h,x g is the interpolation value of the function g at space point x ∈ R q by using the values of g only on D n+1 h,x . Note that any interpolation methods can be used here, however, care should be made if one wants to guarantee the stability and accuracy. Now, we are ready to introduce the following fully discrete scheme: Notice that the iteration method like Scheme 2 is still needed for the fully discrete scheme. Remark 3.3. It is worth to note that in each time level t n , the computations at the spatial grids are independent with each other, and thus, the computational procedure on each discretized time level can be completely parallel. This means the total computational time for each time level is decided by the computational time for one single grid point, provided that enough processors are available, and this is obviously promising for real world simulations.
Remark 3.4. For smooth data b, f and ϕ, it is known that the solution Y t = u(t, X t ) is smooth with respect to t and X t , and furthermore, the conditional expectation E x tn [f (s, X s , Y s )] is a smooth function of s. Thus, the local truncation error R n y has the local estimate which implies that the numerical method can be a high-order numerical scheme.

Gauss-Hermite quadrature rule for E x tn [·]
As mentioned in the last section, numerical approximations for the conditional expectations in Scheme 3 are needed. In this work, we will use the Gauss-Hermite integral to approximate the condition expectations. The Gauss-Hermite quadrature rule is an extension of Gaussian quadrature method for approximating the value of integrals of are the corresponding weights [1]: The truncation error R(g, L) of the Gauss-Hermite quadrature formula (3.10) is where η is a real number in R. The Gauss-Hermite quadrature formula (3.10) is exact for polynomial functions g of degree less than 2L − 1.
For a d-dimensional function g(x), x ∈ R d , the Gauss-Hermite quadrature formula becomes x 2 j , and .
It is well known that, for a standard d-dimensional normal random variable N (0, 1), it holds that Then by (3.12), we deduce where R GH E,L (g) is the truncation error of the Gauss-Hermite quadrature rule for g.
Recall that, in Scheme 3, the conditional expectation E xi is the approximation of E x tn [·], and I n h,x Y n+1 is the interpolation approximation of Y n+1 . By the nonlinear Feynman-Kac formula, Y n+1 has the following explicit representation.  (3.14). Other kinds of conditional expectations can be approximated in a similar way.
Remark 3.5. Note that the Gaussian Hermite approximation for the conditional expectations admits some local properties, namely, only some local informations of Y n+1 (function values around X n ) contribute much for the integral. Then, the scaling idea for the Gaussian Hermite approximation [35] can be used, and in such a way, 6-7 points are enough to guarantee a good approximation.

Numerical experiments
In this section, we provide some numerical examples to demonstrate the effectiveness of our numerical method for solving hyperbolic conservation laws. The numerical examples chosen here will involve all kinds of flow features, that is, shocks, contacts, expansion fans, sonic points, and smooth regions. As we discussed in the last section, we will use an equally spaced grid in space with step ∆x and a uniform time step ∆t. We denote by N x the number of grid points used in the spatial discretizations. In all our computations, 8 points will be used to evaluate conditional expectations.

The advection equation
We first consider the following scalar advection equation where φ(x) is a periodic function with period 2. Notice that the exact solution is In what follows, the computational domain will be fixed to [−1, 1].
To show the convergence rates of our numerical method, we consider here a smooth initial condition φ(x) = − sin(πx). The exact solution is smooth without sonic points.
In Table 1 we list the discrete L 1 -errors and convergence rates (r c ) at t = 2 with different parameters ǫ, ∆t/∆x. The notations Maxvalue and Minvalue represent the maximum and minimum values of the numerical solutions. The results clearly show that our scheme (FBSDEs) is a high-order scheme, and it conserves the Maxvalue principle. Moreover, the convergence rate is maintained even when a large CFL condition is considered (∆t/∆x = 1.6 the second row). Also, the convergence rate is independent of the viscosity parameter ǫ (the twice and third row, ǫ = 10 −6 , 10 −8 ).

The Burgers' equation
We next consider the following scaler nonlinear Burgers' equation The initial function φ(x) is given by (4.3) 9.9991e-1 -9.9991e-1 81 1.5330e-5 9.9999e-1 -9.9999e-1 161 1.9229e-6 1.0000 -1.0000 One can derive the following exact piecewise linear solution for (4.2) and (4.3): Notice that the shock and the expansion fan interact for t > 2/3, which complicates the solution. The exact solution admits at x = −1/3 a sonic expansion fan and at x = 1/3 a steady shock. We   Fig.1 (Left). Also, we provide the numerical results with the parameters ǫ = 10 −10 , ∆t ∆x = 2.0 and 2.8, and N x = 121 in Fig.1 (Right). Both cases show that our scheme behaves very well for such problem with the CFL condition numbers either bigger or smaller than 1.

The Buckley-Leverett equation
The Buckley-Leverett equation is where the non-convex flux function f (u) is given by

The Euler System
The Euler system of gas dynamics can be written as with u = (ρ, v, p) ⊤ , where ρ is the density, v is the velocity and p is the pressure. The matrix A is given by with a = γp ρ being the speed of sound. For details of the problem see e.g. [16]. The Euler system can be rewritten into the following equivalent characteristic  Test case 1. A well-known test problem is the shock tube problem proposed by Sod [33]. The initial conditions are given by We simulate the problem up to time t = 1.5 with ∆x = 0.05 and ∆t/∆x = 1/2.4. Test case 2. Another frequently used test problem is the Lax test case proposed by Lax [18], for which the initial states are (ρ L , v L , p L ) = (0.445, 0.698, 3.528), (ρ R , v R , p R ) = (0.5, 0, 0.571).
Compared to the Sod's problem, the contact discontinuity and the shock of such test are stronger. The problem is computed up to t = 1.5 with space step size ∆x = 0.05 and ∆t/∆x = 0.2. The numerical solutions for the density, the velocity and the presser are given in Fig.3, Fig.4 and Fig.5, respectively. The left plots are for the test case 1 while the right plots are for the test case 2. The plots show that our scheme admits no spurious overshoots or oscillations, and captures the expansion fan, the contact discontinuity and the shock very well.
We take ∆x = ∆y = 0.05 and solve the problem up to t = 0.5 with ∆t/∆x = 1/4 and ∆t/∆x = 1/3. The numerical solutions are shown in Fig.6 (a), and we also show the solutions that cut at some x and y, respectively in Fig.6 (b) and Fig.6 (c). It is noted that our scheme solves the problem very well.

Conclusions
In this work, we provide an alternative way to solve the viscosity solution of hyperbolic conservation laws. The method consists in solving the equivalent probability solution (FBSDEs) of the viscosity solution. Numerical approach for the resulting weakly coupled FBSDEs is discussed. It is noticed that in such framework, the viscosity parameter can be chosen sufficiently small (say 10 −10 ), and the computational procedure on each discretized time level can be completely parallel, furthermore, the traditional CFL condition is dramatically weakened. Several numerical tests are given to show the effectiveness of the numerical method. Although the original motivation of this work is to provide an alternative to compute the viscosity solutions of hyperbolic conservation laws, however, we would like to emphasize that our numerical method for solving FBSDEs is, to the best of our knowledge, novel.
We remark that for multi-dimensional problems, we have introduced the tensorized grids in space, which is computational inefficient. In our future studies, we will update this part into the sparse grid setting, which is efficient for high dimensional approaches. Numerical comparisons between our scheme and other methods (e.g., TVD, WENO, etc.) will also be part of our future work.
We address the issue here and want to open up the possibility of designing efficient stochastic approaches (vis FBSDEs) for more general semi-linear parabolic PDEs, such as general convection-dominated diffusion problems, the Schrödinger equations, etc.


extfieldsofstudy: Biology Medicine

text:
CGAL: computing genome assembly likelihoods

Assembly algorithms have been extensively benchmarked using simulated data so that results can be compared to ground truth. However, in de novo assembly, only crude metrics such as contig number and size are typically used to evaluate assembly quality. We present CGAL, a novel likelihood-based approach to assembly assessment in the absence of a ground truth. We show that likelihood is more accurate than other metrics currently used for evaluating assemblies, and describe its application to the optimization and comparison of assembly algorithms. Our methods are implemented in software that is freely available at http://bio.math.berkeley.edu/cgal/.

Background
Genome assembly is the process of merging fragments of a DNA sequence produced by shotgun sequencing in order to reconstruct the original genome. The assembly problem is known to be NP hard for a number of formulations [1][2][3] and is also complicated by the many types of sequencing errors, experimental biases and the volume of data that must be processed. For these reasons, in addition to differences in underlying theory and algorithms, popular assembly methods employ many different heuristics and assemblies produced by existing methods differ substantially from each other [4,5].
Paradoxically, the difficulties of sequence assembly have been compounded by sequencing advances in recent years collectively termed next-generation sequencing technologies. Next-generation sequencing technologies such as 454 pyrosequencing by Applied Sciences [6], Solexa/Illumina sequencing, the SOLiD technology from Applied Biosystems and Helicos single-molecule sequencing [7] produce data of much greater volume at a much lower cost than traditional Sanger sequencing [8]. However, read lengths are considerably shorter and error rates are higher than those in Sanger sequencing. To allow de novo sequencing from short reads from next-generation sequencing machines several assemblers have been developed such as Velvet [9], Euler-sr [10], ABySS [11], Edena [12], SSAKE [13], VCAKE [14], SHARCGS [15], ALLPATHS [16], SOAPdenovo [17], Celera WGA [18], the CLC bio assembler and others [4,5]. A key problem that has arisen is to determine which assembler is 'the best'. In the past this has been done with the help of a number of measures such as N50 scaffold or contig lengths -which is the maximum contig (scaffold) length such that at least half the total length is contained in contigs (scaffolds) of length greater than or equal to that length. Although simulation studies show that simple metrics correlate with assembly quality, the currently used metrics are crude and provide only condensed summaries of the result. They can therefore be very misleading [5,19]. For example, the assembly consisting of simply gluing all reads end-to-end has a very large N50 length, but is obviously a poor assembly. Phillippy et al. presented software called amosvalidate [20] that identifies mis-assembly features and suspicious regions; however, it does not have high specificity and has not been widely adopted. Narzisi et al. used a featureresponse curve [21] to rank assemblies based on features identified by amosvalidate. Studies such as [22][23][24][25] have discussed these issues and produce interesting insights into assembler performance but do not provide an intrinsic direct measure of assembly quality. The recent Assemblathon 1 competition used 10 different metrics [4] in an attempt to reveal more information than just N50 values, but most of the metrics can only be computed when the genome that is being assembled is known, and are therefore not useful in practice on real data.
In this paper we present a computationally efficient approach for computing the likelihood of an assembly, which provides a way to assess assemblies without a ground truth. Intuitively, the likelihood assessment evaluates the uniformity of coverage of the assembly, taking into account errors in the reads, the insert size distribution and the extent of unassembled data. Genome assembly by maximizing likelihood has been proposed previously by Myers [26] and Medvedev and Brudno [1] but their formulations are based on simplified models that do not use important parameters, especially the sequencing error. To demonstrate the power of our approach for assembly quality evaluation, we have implemented our methods in a program called CGAL. We have evaluated assemblies by testing several of them from different programs with varying input parameters in a setting where the desired target genome is known. For each assembly, we compute the likelihood using our tool and then compare our likelihood computation to standard measures such as N50 contig values, sequence similarity with the reference genome as well as values reported by amosvalidate. Although it is beyond the scope of this paper to compare all assemblers and explore all parameters, our results indicate that likelihood is meaningful and useful for evaluating assemblies.

Results
Our overall approach is simple: we describe a probabilistic generative model for sequencing that captures many aspects of sequencing experiments, and from which we can compute the likelihood of an assembly. This intuitive framework is, however, complicated by one major difficulty, which is the problem we address in this paper: to compute the likelihood of an assembly it is necessary, in principle, to consider the possibility that a read was produced from every single location in the assembly. This results in an intractable computation, which we circumvent by approximating the likelihood via a reduction to a small set of 'likely' sites from which each read originated (using a mapping of the reads to the assembly). This requires an examination of the quality of the approximation, and leads to yet another difficulty, which is how to compute the likelihood for reads that do not map to the assembly at all. These issues are addressed in this paper and their solution is what enables our program for likelihood computation to be efficient and practical.
We begin by describing the statistical model that forms the basis for our likelihood computation. We believe that our model incorporates many aspects of typical sequencing experiments, but it can be easily generalized to accommodate additional parameters if desired.

A generative model for sequencing
Let R = {r 1 , r 2 ,..., r N } be a set of N paired-end (PE) (or mate pair) reads generated from a genome, G (our model can, in principle, be adapted to single-end reads but we do not consider these here). We assume a fragment represented by two paired-end reads r i = (r i1 , r i2 ) is generated according to the following model: • A fragment length l i is selected according to a distribution F. • A site for the 5' end of the fragment s i is selected according to a distribution S.
• The ends of the fragment are read as r i1 and r i2 according to an error model E which comprises mismatches as well as indels.
The generative model is illustrated in Figure 1.

Computing likelihood
Computing the likelihood of an assembly means that the probability of the (observed) set of reads is computed with respect to a proposed assembly using the model described in the previous section. The probability of a sequence of length L generating a paired-end or mate pair read (termed 'read' from now on) r i is where a s ... a s + l -1 is the assembly subsequence starting at s of length l, E denotes all possible ways of obtaining r i from a s ... a s + l -1 . p F (l) is the probability that the fragment length is l, p S (s) is the probability that the 5' end of the fragment is at site s and p E (r i |a s ... a s + l -1 , e) is the probability of obtaining r i from a s ... a s + l -1 with sequencing errors given by e.
Although in theory a read could have been generated from any site (assuming that every base could have been an error), in practice the probability decreases considerably with increasing number of disagreements between the source sequence and the read sequence. We there- fore approximate the probability p(r i ) by mapping the read to the assembly and ignoring mappings with a large number of differences. If M i is the number of such mappings of read r i , the probability is given by where l i, j , s i, j , a i, j and e i, j are the fragment length, start site, assembly subsequence and errors corresponding to the jth mapping of the ith read, respectively. The above equation generalizes to assemblies with more than one contig. Given an assembly A and a set of reads R = {r 1 , r 2 ,..., r N }, the log likelihood is given by In the above equation M i ≥ 1 for all reads r i , and in Methods we explain how we obtain alignments for all reads and how to learn the needed distributions.

Validation with simulated data
To test our implementation, we developed a simulator that generates reads according to given error parameters and fragment lengths distributed according to a Gaussian distribution.
We generated 3 million 35 bp paired-end reads from a strain of Escherichia coli ([NCBI: NC_000913.2]) and an assembly of Grosmannia clavigera ([DDBJ/EMBL/Gen-Bank: ACXQ00000000]) reported in [27]. Table 1 shows the percentage difference in likelihood values computed using true parameters provided to the simulator and using parameters inferred by CGAL.

Performance of assemblers on E. coli reads
We assessed the performance of four assemblers: Velvet, Euler-sr, ABySS and SOAPdenovo on an Escherichia coli dataset ([SRA:SRR 001665] and [SRA:SRR 001666]). We chose E. coli because its assembly is a true 'gold standard' without questions about reliability or accuracy. We assembled the reads using the assemblers mentioned for different hash lengths (k-mer was used for constructing the de Bruijn graph [10]). Likelihood values for assemblies along with the likelihood value for the reference ([NCBI: U00096.2]) are shown in Figure 2.
For this dataset ABySS outperforms the others when likelihood is used as the metric. We also aligned the assemblies to the reference with NUCmer [28] and Figure 3 shows the differences from the reference against the hash lengths. The relations among likelihood, N50 length and similarity are illustrated in Figure 4 and Additional file 1, Figure S1. They suggest that likelihood values are better at capturing sequence similarity than other metrics commonly used for evaluating assemblies, such as the N50 scaffold or contig lengths. We also ran the amosvalidate pipeline to obtain the numbers of mis-assembly of features and suspicious regions ( Figure 5) and plotted the feature response curves (FRCs) [21] of the assemblies (Additional file 1, Figures S4, S5). The FRCs also rank an ABySS assembly as the best one. A similar analysis was performed on a different Escherichia coli dataset downloaded from CLC bio [29]. It consists of approximately 2.6 million 35 bp paired-end Illumina reads (approximately 40 times coverage) along with a reference genome ([NCBI: NC_010473.1]). We noticed that many of the assemblies have a better likelihood than the reference. However, we assembled reads that could not be mapped to the reference and after running BLAST [30] we found another substrain of Escherichia coli strain K-12, MG1655 ([NCBI: NC_000913.2]), which has a better likelihood than all assemblies.  We conjecture that the reads were generated from NC_000913.2. Likelihood values are shown in Figure 6 and relationships among likelihood, similarity and N50 values are illustrated in Additional file 1, Figures S6-S10. The differences between assemblies and the reference are shown on the y-axis where the difference refers to the numbers of bases in the reference not covered by the assembly or differ between the reference and the assembly.   Figure 6 Hash length vs log likelihood for E. coli data from CLC bio. Log likelihoods of assemblies of E. coli reads from CLC bio are shown on the y-axis. Assemblies are generated using different assemblers for varying k-mer length, which is shown on the x-axis.
The yellow dotted line corresponds to the log likelihood of the reference provided and the gray dotted line corresponds to the log likelihood of the strain we believe the reads were generated from.  Figure 7. It also shows likelihood values for assemblies [DDBJ/EMBL/GenBank: ACXQ00000000] and [DDBJ/EMBL/GenBank: ACYC000 00000] reported in [27], which were generated using Sanger and 454 reads as well as Illumina reads. The numbers of mis-assembly features and suspicious regions identified by amosvalidate and the feature response curves are shown in Additional file 1, Figures S14-S15. Figure 8 shows that the assembly with most sequence coverage is produced by ABySS. However, in this case ABySS assemblies are much longer compared to other assemblies and references (Additional file 1, Tables S9-S11). This results in lower likelihoods compared to some assemblies by Velvet and SOAPdenovo. Figure 9 and Additional file 1, Figure S11 show relationships among likelihood, similarity and N50 values. In FRC analysis, genome coverage is estimated using assembly length and so it does not take into account the unassembled sequences and ranks ABySS assemblies above others. It is interesting that assemblies with the best likelihood and sequence similarity are generated for higher values of hash length than are optimal for producing high N50 values.

GAGE results
We computed likelihoods for the assemblies generated in the GAGE project [5]. In Additional file 1, Tables S12-S14 show likelihoods of Library 1 and the number     Figure 9 Log likelihood vs N50 scaffold length for G. clavigera. Log likelihoods are shown on the x-axis and N50 scaffold lengths are shown on the y-axis. Each circle corresponds to an assembly generated using an assembler for some hash length and the sizes of the circles correspond to similarity with reference. The R 2 values are: (i) log likelihood vs similarity: 0.4545793, (ii) log likelihood vs N50 scaffold length: 0.002397233, (iii) N50 scaffold length vs similarity: 0.006084032.
of reads mapped to assemblies by Bowtie 2 [31]. We found that the likelihood values of Library 1 are dominated by coverage and contiguity does not affect these values greatly. However, contiguity has more effect on the likelihoods of Library 2, which has a longer insert size (Additional file 1, Tables S12-S14), as might be expected. The total likelihood along with coverage and N50 values are shown in Tables 2, 3, 4, 5. For human chromosome 14, we computed Library 2 likelihoods for assemblies with the best three likelihoods of Library 1.
The likelihood values of Library 2 for bumble bee assemblies were not computed as only a small fraction of the reads could be mapped to the assemblies.

Assemblathon 1 results
We also analyzed the assemblies submitted for Assemblathon 1 [4]. The likelihoods for a library of an insert size of mean 200 bp for all assemblies are given in Additional file 1, Table S15. Figure 10 shows the relationship between likelihood and coverage. We took the entries with the highest likelihood for the top ten participants and computed the likelihoods for the libraries of insert sizes of means 3,000 bp and 10,000 bp. Table 6 shows the total likelihoods of the top ten participants along with their Assemblathon 1 rankings.

E. coli
We found that for both E. coli datasets, the assemblies with the best likelihoods were constructed by ABySS. They also have most similarity with the references (assuming [NCBI: NC 000913.2] is the reference for the CLC bio dataset). The R 2 values (Figures 4, 5 and Additional file 1, Figure S1) reveal that the likelihoods reflect sequence similarity better than contiguity statistics such as N50 values as well as the numbers of mis-assembly features and suspicious regions identified by amosvalidate. The analysis of the two different E. coli datasets also reveal that for assemblers like Velvet and SOAPdenovo higher likelihood values are achieved for different values of the k-mer length used to construct the de Bruijn graph during assembly.

G. clavigera
For the G. clavigera dataset, one of the Velvet assemblies has the highest likelihood. Although ABySS assemblies have more coverage, they have lower likelihood because of the much longer total length. Despite this we see from the R 2 values that likelihood values reflect sequence similarity better than the N50 values ( Figure 9 and Additional file 1, Figures S11, S14) and the numbers

GAGE
For the GAGE Staphylococcus aureus dataset, we find that the assembly generated using Velvet has the best likelihood but likelihoods of a few other assemblies are close. For Rhodobacter sphaeroides, the ALLPATHS-LG assembly has the best likelihood, which is also the assembly with the highest coverage and N50 scaffold length. The CABOG assembly of human chromosome 14 is the one with the best likelihood. The CABOG assembly also has the highest coverage and N50 contig length among the assemblies. In all three cases, we find that the reference sequences have the highest likelihoods and the highest number of reads mapped to them by Bowtie 2.
For the bumblebee data, the assembly using CABOG has the best likelihood of the three (the likelihood of the SOAPdenovo assembly could not be computed as reads could not be mapped to it using Bowtie 2).
Assemblathon 1 Figure 10 reveals that for the Assemblathon 1 dataset, the likelihoods for a small fragment library correlate well with coverages. Overall, we find that participants with the ten highest likelihoods were ranked within the top eleven by the Assemblathon 1 organizers but there are differences between the two rankings. The entry with the highest likelihood is from the Beijing Genomics Institute (BGI), which was ranked two in the original paper. The differences in rankings are primarily due to the emphasis on contiguity made by the Assemblathon 1 organizers while our likelihood model implicitly places more importance on coverage. This brings up the issue that better contiguity statistics can be achieved by not reporting hard-to-assemble regions and these values may be misleading if they are not used in conjunction with an indicator of coverage.

Applications
Currently, assembly evaluation projects rely mostly on simulated data or data from genomes that have been sequenced previously [4,5]. Having a tool that can assess   the quality of assembly without the need for a reference will allow researchers who work with real data from genomes that have not been sequenced before to assess the performance of different assemblers on their data, and to optimize the parameters in the programs they are using. The analysis of two different datasets from E. coli reveals that the performance of some assemblers varies significantly depending on the k-mer chosen for constructing the de Bruijn graph. Moreover, the 'optimal' value depends on read length and sequence coverage. Likelihood values can therefore guide selection of parameter values.
The concept of maximum likelihood genome assembly was introduced by Medvedev and Brudno [1] but they do not consider sequencing errors or paired-end reads. A likelihood model taking into account these may be the next step towards genome assemblers for real data that try to maximize likelihood.

Conclusions
In this paper we presented a tool for computing the likelihood of an assembly. The result can be used as a metric for evaluating and comparing assemblies. In the past this has been done using many different criteria including N50 lengths, total sequence length and number of contigs. The likelihood model incorporates these directly or indirectly in addition to other important factors such as genome coverage and assembly accuracy and combines them into a single metric for evaluation.
We have also used our tool to assess the performance of assemblers using different datasets. Our results indicate that likelihood reflects sequence similarity, which is missed by other metrics commonly used and will be a valuable tool for evaluating assemblies generated by different assemblers and for different values of the input parameters.

Mapping reads
The first step in computing the likelihood is mapping reads to the assembly. A number of tools are available for this such as Bowtie [31,32], MAQ [33], BWA [34] and BFAST [35]. Our present implementation can use either BFAST or Bowtie 2 for mapping reads as they support mapping with indels and report multiple alignments in a way that gives all the required information without accessing the assembly sequence. But any tool that reports multiple alignments of reads and allows for insertions and deletions can be used with some minor modifications.
However, existing tools do not usually map all reads, and for the likelihood computation it is necessary to assign probabilities to reads that are unmapped. We found that mapping tools were unable to map a large fraction of reads in our experiments. One option is to assign probabilities to these reads, assuming that they could have been generated from any site, using the number and types of errors not handled by the mapping tool. But it is then often the case that unmapped reads are deemed more probable than mapped ones, which we believe is anomalous. Furthermore, in our analyses we determined that the resulting probabilities were inaccurate (results not shown). Therefore, we chose to directly align the reads not mapped by BFAST or Bowtie 2 using an adaptation of the Smith-Waterman algorithm. We adapted the striped implementation of the Smith-Waterman algorithm by Farrar [36]. This step is time consuming, so we align only a random subset of reads with the number specified by the user and approximate probabilities using these.

Learning distributions
To compute the likelihood from mapped reads, we need to learn the distribution of fragment lengths, their distribution across the genome and error characteristics. Since they differ with library preparation methods and sequencing instruments, we chose to learn these from sequencing data generated in the experiment. We do this by mapping reads to the assembly and using reads that map uniquely. However, this can be easily extended to take into account all reads by using the expectation-maximization (EM) algorithm at the expense of more iterations. We explain each distribution in more detail below.

Fragment length distribution
The distribution of fragment lengths depends on the method used for size selection and may not be approximated well by common distributions [37]. So, we use the empirical distribution.

Distribution of fragments along genome
In our implementation, we assume that fragments are distributed uniformly across the genome. We leave incorporating sequencing bias as future work.

Error model
In the error model used at present, we have made the assumption that sequencing errors are independent of one another. We learn an error rate for each position in the read since error rates are known to be different across positions in reads [38]. We also learn separate error rates for each type of base and substitution type. Although errors are known to depend on sequence context [38], we have ignored them for the sake of simplicity.
To account for varying indel rates across positions in reads, we learn an insertion rate and a deletion rate for each position in the read. Since short indels are more likely than longer ones, we also count the number of insertions and deletions by length.

Implementation
As mentioned earlier, we use BFAST or Bowtie 2 to map reads to assemblies. The parameters are set so that they report all alignments of a read found.
The remaining code for computing likelihood is written in C++ and it consists of three parts: convert: This converts the output generated by BFAST or Bowtie 2 to an internal format. It also separates reads with no end or one end mapped and reads with ends mapped to different scaffolds if needed. Separating this module also allow us to support other mapping tools by writing a conversion routine.
align: To align the reads not mapped by the mapping tool, we adapted the striped implementation of the Smith-Waterman algorithm by Farrar [36]. As this step is time consuming, we align a random subset of reads with the number determined by the user. This step is multithreaded to speed up the process.
CGAL: This learns the fragment length distribution and parameters for the error model using uniquely mapped reads and then uses these to compute the likelihood value.

Assembling genomes
To assemble reads, we varied the k-mer length used to construct the de Bruijn graph to obtain different assemblies for each assembly tool. For other parameters, the default values or values suggested in manuals were used.

Data analysis
Likelihoods were computed by running CGAL with the default parameters and aligning between 300 and 1,000 randomly chosen reads not mapped by the mapping tool used. The running time for CGAL was approximately 1/3 of the time taken to map reads using Bowtie 2.
To compute the difference between an assembly and the reference, we aligned the assembly to the reference using NUCmer [28]. The difference refers to the number of bases in the reference that are either not covered by the assembly or differ in the reference and assembly. Contigs were generated by splitting scaffolds at sites with 25 or more N's (character representing any base).


Joint Space Neural Probabilistic Language Model for Statistical Machine Translation

A neural probabilistic language model (NPLM) provides an idea to achieve the better perplexity than n-gram language model and their smoothed language models. This paper investigates application area in bilingual NLP, specifically Statistical Machine Translation (SMT). We focus on the perspectives that NPLM has potential to open the possibility to complement potentially `huge' monolingual resources into the `resource-constraint' bilingual resources. We introduce an ngram-HMM language model as NPLM using the non-parametric Bayesian construction. In order to facilitate the application to various tasks, we propose the joint space model of ngram-HMM language model. We show an experiment of system combination in the area of SMT. One discovery was that our treatment of noise improved the results 0.20 BLEU points if NPLM is trained in relatively small corpus, in our case 500,000 sentence pairs, which is often the case due to the long training time of NPLM.

Introduction
A neural probabilistic language model (NPLM) [3,4] and the distributed representations [25] provide an idea to achieve the better perplexity than n-gram language model [47] and their smoothed language models [26,9,48]. Recently, the latter one, i.e. smoothed language model, has had a lot of developments in the line of nonparametric Bayesian methods such as hierarchical Pitman-Yor language model (HPYLM) [48] and Sequence Memoizer (SM) [51,20], including an application to SMT [36,37,38]. A NPLM considers the representation of data in order to make the probability distribution of word sequences more compact where we focus on the similar semantical and syntactical roles of words. For example, when we have two sentences "The cat is walking in the bedroom" and "A dog was running in a room", these sentences can be more compactly stored than the n-gram language model if we focus on the similarity between (the, a), (bedroom, room), (is, was), and (running, walking). Thus, a NPLM provides the semantical and syntactical roles of words as a language model. A NPLM of [3] implemented this using the multi-layer neural network and yielded 20% to 35% better perplexity than the language model with the modified Kneser-Ney methods [9].
There are several successful applications of NPLM [41,11,42,10,12,14,43]. First, one category of applications include POS tagging, NER tagging, and parsing [12,7]. This category uses the features provided by a NPLM in the limited window size. It is often the case that there is no such long range effects that the decision cannot be made beyond the limited windows which requires to look carefully the elements in a long distance. Second, the other category of applications include Semantic Role Labeling (SRL) task [12,14]. This category uses the features within a sentence. A typical element is the predicate in a SRL task which requires the information which sometimes in a long distance but within a sentence. Both of these approaches do not require to obtain the best tag sequence, but these tags are independent. Third, the final category includes MERT process [42] and possibly many others where most of them remain undeveloped. The objective of this learning in this category is not to search the best tag for a word but the best sequence for a sentence. Hence, we need to apply the sequential learning approach. Although most of the applications described in [11,10,12,14] are monolingual tasks, the application of this approach to a bilingual task introduces really astonishing aspects, which we can call "creative words" [50], automatically into the traditional resource constrained SMT components. For example, the training corpus of word aligner is often strictly restricted to the given parallel corpus. However, a NPLM allows this training with huge monolingual corpus. Although most of this line has not been even tested mostly due to the problem of computational complexity of training NPLM, [43] applied this to MERT process which reranks the n-best lists using NPLM. This paper aims at different task, a task of system combination [1,29,49,15,13,35]. This category of tasks employs the sequential method such as Maximum A Posteriori (MAP) inference (Viterbi decoding) [27,44,33] on Conditional Random Fields (CRFs) / Markov Random Fields (MRFs).
Although this paper discusses an ngram-HMM language model which we introduce as one model of NPLM where we borrow many of the mechanism from infinite HMM [19] and hierarchical Pitman-Yor LM [48], one main contribution would be to show one new application area of NPLM in SMT. Although several applications of NPLM have been presented, there have been no application to the task of system combination as far as we know.
The remainder of this paper is organized as follows. Section 2 describes ngram-HMM language model while Section 3 introduces a joint space model of ngram-HMM language model. In Section 4, our intrinsic experimental results are presented, while in Section 5 our extrinsic experimental results are presented. We conclude in Section 5.

Ngram-HMM Language Model
Generative model Figure 1 depicted an example of ngram-HMM language model, i.e. 4-gram-HMM language model in this case, in blue (in the center). We consider a Hidden Markov Model (HMM) [40,21,2] of size K which emits n-gram word sequence w i , . . . , w i−K+1 where h i , . . . , h i−K+1 denote corresponding hidden states. The arcs from w i−3 to w i , · · · , w i−1 to w i show the back-off relations appeared in language model smoothing, such as Kneser-Ney smoothing [26], Good-Turing smoothing [24], and hierarchical Pitman-Yor LM smoothing [48]. In the left side in Figure 1, we place one Dirichlet Process prior DP(α, H), with concentration parameter α and base measure H, for the transition probabilities going out from each hidden state. This construction is borrowed from the infinite HMM [2,19]. The observation likelihood for the hidden word h t are parameterized as in w t |h t ∼ F (φ ht ) since the hidden variables of HMM is limited in its representation power where φ ht denotes output parameters. This is since the observations can be regarded as being generated from a dynamic mixture model [19] as in (1), the Dirichlet priors on the rows have a shared parameter.
In the right side in Figure 1, we place Pitman-Yor prior PY, which has advantage in its power-law behavior as our target is NLP, as in (2): where α is a concentration parameter, θ is a strength parameter, and G i is a base measure. This construction is borrowed from hierarchical Pitman-Yor language model [48].
Inference We compute the expected value of the posterior distribution of the hidden variables with a beam search [19]. This blocked Gibbs sampler alternate samples the parameters (transition matrix, output parameters), the state sequence, hyper-parameters, and the parameters related to language model smoothing. As is mentioned in [19], this sampler has characteristic in that it adaptively truncates the state space and run dynamic programming as in (3): where u t is only valid if this is smaller than the transition probabilities of the hidden word sequence h 1 , . . . , h K . Note that we use an auxiliary variable u i which samples for each word in the sequence from the distribution u i ∼ Uniform(0, π (hi−1,hi) ). The implementation of the beam sampler consists of preprocessing the transition matrix π and sorting its elements in descending order.
Initialization First, we obtain the parameters for hierarchical Pitman-Yor process-based language model [48,23], which can be obtained using a block Gibbs sampling [32].
Second, in order to obtain a better initialization value h for the above inference, we perform the following EM algorithm instead of giving the distribution of h randomly. This EM algorithm incorporates the above mentioned truncation [19]. In the E-step, we compute the expected value of the posterior distribution of the hidden variables. For every position h i , we send a forward message α(h i−n+1:i−1 ) in a single path from the start to the end of the chain (which is the standard forward recursion in HMM; Hence we use α). Here we normalize the sum of α considering the truncated variables u i−n+1:i−1 .
Then, for every position h j , we send a message β(h i−n+2:i , h j ) in multiple paths from the start to the end of the chain as in (5), This step aims at obtaining the expected value of the posterior distribution (Similar construction to use expectation can be seen in factored HMM [22]). In the M-step, using this expected value of the posterior distribution obtained in the E-step to evaluate the expectation of the logarithm of the complete-data likelihood.

Joint Space Model
In this paper, we mechanically introduce a joint space model. Other than the ngram-HMM language model obtained in the previous section, we will often encounter the situation where we have another hidden variables h 1 which is irrelevant to h 0 which is depicted in Figure 2. Suppose that we have the ngram-HMM language model yielded the hidden variables suggesting semantic and syntactical role of words. Adding to this, we may have another hidden variables suggesting, say, a genre ID. This genre ID can be considered as the second context which is often not closely related to the first context. This also has an advantage in this mechanical construction that the resulted language model often has the perplexity smaller than the original ngram-HMM language model. Note that we do not intend to learn this model jointly using the universal criteria, but we just concatenate the labels by different tasks on the same sequence. By this formulation, we intend to facilitate the use of this language model. It is noted that those two contexts may not be derived in a single learning algorithm. For example, language model with the sentence context may be derived in the same way with that with the word context. In the above example, a hidden semantics over sentence is not a sequential object. Hence, this can be only considering all the sentence are independent. Then, we can obtain this using, say, LDA.

Extrinsic Evaluation: Task of System Combination
We applied ngram-HMM language model to the task of system combination. For given multiple Machine Translation (MT) outputs, this task essentially combines the best fragments among given MT outputs to recreate a new MT output. The standard procedure consists of three steps: Minimum Bayes Risk decoding, monolingual word alignment, and monotonic consensus decoding. Although these procedures themselves will need explanations in order to understand the following, we keep the main text in minimum, moving some explanations (but not sufficient) in appendices. Note that although this experiment was done using the ngram-HMM language model, any NPLM may be sufficient for this purpose. In this sense, we use the term NPLM instead of ngram-HMM language model.

Features in Joint Space
The first feature of NPLM is the semantically and syntactically similar words of roles, which can be derived from the original NPLM. We introduce the second feature in this paragraph, which is a genre ID.
The motivation to use this feature comes from the study of domain adaptation for SMT where it becomes popular to consider the effect of genre in testset. This paper uses Latent Dirichlet Allocation (LDA) [5,46,6,45,33] to obtain the genre ID via (unsupervised) document classification since our interest here is on the genre of sentences in testset. And then, we place these labels on a joint space.
LDA represents topics as multinomial distributions over the W unique word-types in the corpus and represents documents as a mixture of topics. Let C be the number of unique labels in the corpus. Each label c is represented by a W -dimensional multinomial distribution φ c over the vocabulary. For document d, we observe both the words in the document w (d) as well as the document labels c (d) . Given the distribution over topics θ d , the generation of words in the document is captured by the following generative model.The parameters α and β relate to the corpus level, the variables θ d belong to the document level, and finally the variables z dn and w dn correspond to the word level, which are sampled once for each word in each document.
Using topic modeling in the second step, we propose the overall algorithm to obtain genre IDs for testset as in (5). 1. Fix the number of clusters C, we explore values from small to big where the optimal value will be searched on tuning set. Modified Process in System Combination Given a joint space of NPLM, we need to specify in which process of the task of system combination among three processes use this NPLM. We only discuss here the standard system combination using confusion-network. This strategy takes the following three steps (Very brief explanation of these three is available in Appendix): • Minimum Bayes Risk decoding [28] (with Minimum Error Rate Training (MERT) process [34] • Monolingual word alignment • (Monotone) consensus decoding (with MERT process) Similar to the task of n-best reranking in MERT process [43], we consider the reranking of nbest lists in the third step of above, i.e. (monotone) consensus decoding (with MERT process). We do not discuss the other two processes in this paper.
On one hand, we intend to use the first feature of NPLM, i.e. the semantically and syntactically similar role of words, for paraphrases. The n-best reranking in MERT process [43] alternate the probability suggested by word sense disambiguation task using the feature of NPLM, while we intend to add a sentence which replaces the words using NPLM. On the other hand, we intend to use the second feature of NPLM, i.e. the genre ID, to split a single system combination system into multiple system combination systems based on the genre ID clusters. In this perspective, the role of these two feature can be seen as independent. We conducted four kinds of settings below.
(A) -First Feature: N-Best Reranking in Monotonic Consensus Decoding without Noise -NPLM plain In the first setting for the experiments, we used the first feature without considering noise. The original aim of NPLM is to capture the semantically and syntactically similar words in a way that a latent word depends on the context. We will be able to get variety of words if we condition on the fixed context, which would form paraphrases in theory.
We introduce our algorithm via a word sense disambiguation (WSD) task which selects the right disambiguated sense for the word in question. This task is necessary due to the fact that a text is natively ambiguous accommodating with several different meanings. The task of WSD [14] can be written as in (6): where k ranges over all possible features, f (feature k i ) is an indicator function whose value is 1 if the feature exists, and 0 otherwise, g(synset i , k) is a parameter for a given synset and feature, θ is a collection of all these parameters in g(synset i , k), and Z is a normalization constant. Note that we use the term "synset" as an analogy of the WordNet [30]: this is equivalent to "sense" or "meaning". Note also that NPLM will be included as one of the features in this equation. If features include sufficient statistics, a task of WSD will succeed. Otherwise, it will fail. We do reranking of the outcome of this WSD task.
On the one hand, the paraphrases obtained in this way have attractive aspects that can be called "a creative word" [50]. This is since the traditional resource that can be used when building a translation model by SMT are constrained on parallel corpus. However, NPLM can be trained on huge monolingual corpus. On the other hand, unfortunately in practice, the notorious training time of NPLM only allows us to use fairly small monolingual corpus although many papers made an effort to reduce it [31]. Due to this, we cannot ignore the fact that NPLM trained not on a huge corpus may be affected by noise. Conversely, we have no guarantee that such noise will be reduced if we train NPLM on a huge corpus. It is quite likely that NPLM has a lot of noise for small corpora. Hence, this paper also needs to provide the way to overcome difficulties of noisy data. In order to avoid this difficulty, we limit the paraphrase only when it includes itself in high probability.
(B)-First Feature: N-Best Reranking in Monotonic Consensus Decoding with Noise -NPLM dep In the second setting for our experiment, we used the first feature considering noise. Although we modified a suggested paraphrase without any intervention in the above algorithm, it is also possible to examine whether such suggestion should be adopted or not. If we add paraphrases and the resulted sentence has a higher score in terms of the modified dependency score [39] (See Figure 3), this means that the addition of paraphrases is a good choice. If the resulted score decreases, we do not need to add them. One difficulty in this approach is that we do not have a reference which allows us to score it in the usual manner. For this reason, we adopt the naive way to deploy the above and we deploy this with pseudo references. (This formulation is equivalent that we decode these inputs by MBR decoding.) First, if we add paraphrases and the resulted sentence does not have a very bad score, we add these paraphrases since these paraphrase are not very bad (naive way). Second, we do scoring between the sentence in question with all the other candidates (pseudo references) and calculate an average of them. Thus, our second algorithm is to select a paraphrase which may not achieve a very bad score in terms of the modified dependency score using NPLM.
(C) -Second Feature: Genre ID -DA (Domain Adaptation) In the third setting of our experiment, we used only the second feature. As is mentioned in the explanation about this feature, we intend to splits a single module of system combination into multiple modules of system combi-  Figure 3: By the modified dependency score [39], the score of these two sentences, "John resigned yesterday" and "Yesterday John resigned", are the same. Figure shows c-structure and f-structure of two sentences using Lexical Functional Grammar (LFG) [8].
nation according to the genre ID. Hence, we will use the module of system combination tuned for the specific genre ID, 1 .
(D) -First and Second Feature -COMBINED In the fourth setting we used both features. In this setting, (1) we used modules of system combination which are tuned for the specific genre ID, and (2) we prepared NPLM whose context can be switched based on the specific genre of the sentence in test set. The latter was straightforward since these two features are stored in joint space in our case.
Experimental Results ML4HMT-2012 provides four translation outputs (s1 to s4) which are MT outputs by two RBMT systems, APERTIUM and LUCY, PB-SMT (MOSES) and HPB-SMT (MOSES), respectively. The tuning data consists of 20,000 sentence pairs, while the test data consists of 3,003 sentence pairs.
Our experimental setting is as follows. We use our system combination module [16,17,35], which has its own language modeling tool, MERT process, and MBR decoding. We use the BLEU metric as loss function in MBR decoding. We use TERP 2 as alignment metrics in monolingual word alignment. We trained NPLM using 500,000 sentence pairs from English side of EN-ES corpus of EUROPARL 3 .
The results show that the first setting of NPLM-based paraphrased augmentation, that is NPLM plain, achieved 25.61 BLEU points, which lost 0.39 BLEU points absolute over the standard system combination. The second setting, NPLM dep, achieved slightly better results of 25.81 BLEU points, which lost 0.19 BLEU points absolute over the standard system combination. Note that the baseline achieved 26.00 BLEU points, the best single system in terms of BLEU was s4 which achieved 25.31 BLEU points, and the best single system in terms of METEOR was s2 which achieved 0.5853. The third setting achieved 26.33 BLEU points, which was the best among our four settings. The fourth setting achieved 25.95, which is again lost 0.05 BLEU points over the standard system combination.
Other than our four settings where these settings differ which features to use, we run several different settings of system combination in order to understand the performance of four settings. Standard system combination using BLEU loss function (line 5 in Table 2), standard system combination using TER loss function (line 6), system combination whose backbone is unanamously taken from the RBMT outputs (MT input s2 in this case; line 11), and system combination whose backbone is selected by the modified dependency score (which has three variations in the figure; modDep preci-sion, recall and Fscore; line 12, 13 and 14). One interesting characteristics is that the s2 backbone (line 11) achieved the best score among all of these variations. Then, the score of the modified dependency measure-selected backbone follows. From these runs, we cannot say that the runs related to NPLM, i.e. (A), (B) and (D), were not particularly successful. The possible reason for this was that our interface with NPLM was only limited to paraphrases, which was not very successfuly chosen by reranking.  Table 2: This table shows single best performance, the performance of the standard system combination (BLEU and TER loss functions), the performance of four settings in this paper ((A),. . .,(D)), the performance of s2 backboned system combination, and the performance of the selection of sentences by modified dependency score (precision, recall, and F-score each).

Conclusion and Perspectives
This paper proposes a non-parametric Bayesian way to interpret NPLM, which we call ngram-HMM language model. Then, we add a small extension to this by concatenating other context in the same model, which we call a joint space ngram-HMM language model. The main issues investigated in this paper was an application of NPLM in bilingual NLP, specifically Statistical Machine Translation (SMT). We focused on the perspectives that NPLM has potential to open the possibility to complement potentially 'huge' monolingual resources into the 'resource-constraint' bilingual resources. We compared our proposed algorithms and others. One discovery was that when we use a fairly small NPLM, noise reduction may be one way to improve the quality. In our case, the noise reduced version obtained 0.2 BLEU points better.
Further work would be to apply this NPLM in various other tasks in SMT: word alignment, hierarchical phrase-based decoding, and semantic incorporated MT systems in order to discover the merit of 'depth' of architecture in Machine Learning.