{'fraction_non_alphanumeric': 0.029739175457899166, 'fraction_numerical': 0.011779172869069963, 'mean_word_length': 4.915581929555896, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Distributed Machine Learning for Computational Engineering using MPI

We propose a framework for training neural networks that are coupled with partial differential equations (PDEs) in a parallel computing environment. Unlike most distributed computing frameworks for deep neural networks, our focus is to parallelize both numerical solvers and deep neural networks in forward and adjoint computations. Our parallel computing model views data communication as a node in the computational graph for numerical simulations. The advantage of our model is that data communication and computing are cleanly separated and thus provide better flexibility, modularity, and testability. We demonstrate using various large-scale problems that we can achieve substantial acceleration by using parallel solvers for PDEs in training deep neural networks that are coupled with PDEs.


Introduction
Inverse modeling [1,2], which aims at identifying parameters or functions in a physical model from observations, is one key technique to solving data-driven problems in many applications. With the advent of deep learning techniques [3], we are now able to build neural-network-based datadriven models to solve inverse problems involving more complex physical phenomena. Specifically, by substituting unknown components within a physical system described by partial differential equations (PDEs) with deep neural networks (DNNs) [4,5,6,7,8] , we can leverage the expressive power of DNNs while satisfying the physics to the largest extent.
Gradient-based optimization algorithms [9,10] are usually applied to train DNNs that are coupled with PDEs. A loss function, which measures the discrepancy between true and hypothetical observations, is minimized by iteratively updating DNN weights and biases. The gradients can be computed using automatic differentiation [11,12,13] and adjoint state methods [14]. In our previous work [15,16,5] , we have developed a general framework, ADCME 1 , that expresses both DNNs and numerical PDE solvers (e.g., finite element methods) as computational graphs. Therefore, we can calculate the gradients automatically and with machine accuracy using reversemode automatic differentiation (AD). In ADCME, the automatic differentiation operates on a higher level of abstractions, such as tensor operations and matrix solvers, or even PDE solvers, instead of elementary ones (e.g., arithmetic operations) considered in general purpose AD software packages. This coarser granularity allows us to apply tensor-level performance optimizations and implement algorithms with a more intuitive and mathematical interface for scientific computing applications.
However, as the problem size grows, the memory consumption becomes prohibitive because reverse-mode AD requires saving all intermediate results. For example, a direct AD implementation of a 2D elastic equation double-precision solver with a mesh size 1000×1000 and 2000 steps requires at least 193 gigabytes memory 2 [17], which is impractical for a typical CPU. Parallel computing using MPI [18,19] is the de-facto standard for solving such a large-scale problem on modern distributed memory high performance computing (HPC) architectures. However, integrating a parallel computing capability into a reverse-mode automatic differentiation framework is challenging [20,21,22]. We need to back-propagate gradients through DNNs, numerical solvers, and control flows of parallel communications (Figure 1). These difficulties are magnified by the fact that we need the flexibility for hybrid threaded MPI capabilities (e.g., each MPI processor runs a multithreading task and uses multiple cores on the same CPU) and therefore special efforts are needed to topological design of computational graph to avoid deadlocks.
Currently, distributed computing for reverse-mode AD have been intensively studied by two research communities: parallelizing numerical solvers [23,24,25] in scientific computing, and parallelized training (e.g., DNNs) [26,27] in machine learning. In the scientific computing community, a computational grid is usually split into multiple patches and each MPI processor owns one patch. This allows to distribute degrees of freedoms onto different machines and thus relaxes the memory demands. Examples include libadjoint [28,29] and DAfoam [30]. However, these frameworks are domain specific and are tightly integrated into dolfin [31] and OpenFOAM [32], respectively. This integration makes them difficult to adapt for general purpose data-driven inverse modeling, especially when deep neural networks are present and coupled with PDEs. In the machine learning community, one popular model is data-parallelism distributed computing [33], where each machine has full information of a complete DNN but only a portion of training data [34]. These DNNs share the same parameters, through a synchronization mechanism such as parameter servers ( Figure 2) [35]. In these cases, the computations are mostly independent and there is no communication among different processors within each computational graphs. There are also the so-called modelparallelism distributed computing, in which DNNs themselves are split onto different machines [36]. However, these methods eventually lead to asynchronous training, where the correctness of gradients or the convergence is not guaranteed. Note that the distributed computing models from both two communities are quite different and are designed to meet their own needs. To couple DNNs and PDEs, where both synchronization of DNN parameters and parallel PDE solvers are present, a new distributed computing model that combines their features is desirable.
In this work, we propose a novel framework for developing MPI-based distributed reverse-mode AD algorithms for coupling DNNs and PDEs. We consider treating data communication within both PDE solvers and DNN updates as nodes (operators) in the computational graph ( Figure 3). These operators are non-intrusive to the other operators in the sense that existing implementation of the latter operators also works for parallel execution. The clean separation of data communication and computing offers the flexibility of reusing the existing ADCME framework and having the same implementation for both serial and parallel computing environments.
We use the high level language Julia [37] to specify the dependency of data communication operators and other computation operators. These data communication operators can be either built-in ADCME operators (e.g., mpi bcast, which broadcasts a tensor to all workers), or provided by users in the form of custom operators. All computation and data communication are delegated  to C++ or CUDA kernels at runtime, and the correctness of parallel programs is ensured by the topology of the computational graph. The MPI capability of ADCME relieves us from reasoning about the order of MPI calls, e.g., carefully ordering sends and receives to avoid deadlocks, and parallel communication rules for AD, e.g., forward sends become backward receives. Therefore, ADCME allows us to manage complex numerical simulations and accelerate development cycles. The separation of data communication and computation operators also allows us to test, analyze, and optimize them individually.
The aim of this paper is to describe the mechanism underlying ADCME MPI-based distributed computing capability. We also apply the tool to compute DNN gradients in a coupled system of DNNs and PDEs, such as Poisson's equations and wave equations. We also present some important operators and data structures for scientific computing. For example, we use Hypre [38] as the distributed sparse linear solver backend and use CSR formats (mpi SparseTensor) for storing a sparse matrix on each machine. In gradient back-propagation, the transpose of the original matrix is usually needed and such functionalities are missing in Hypre as of September 2020. The transposition is implemented as a separate operator in ADCME. The numerical examples demonstrate that ADCME can scale up to thousands of cores and solve large-scale scientific machine learning problems.

Machine Learning for Computational Engineering
As an example, we consider the Poisson's equation with a spatially varying diffusivity coefficient Here f (x) is a given source function, and Ω is a domain in R d . Assume that we can observe the state variable u(x) at some locations {x i } i∈I , denoted as {u i } i∈I , where I is the index set. We want to estimate κ(x), whose form is in unspecified, from the observations. This inverse problem can be formulated as an optimization problem min κ∈K i∈I Here u(x i ) are the solutions to Equation (1) for a given κ(x), and K is an appropriate function space for the unknown κ. However, Equation (2) is typically an infinitely dimensional optimization problem if K cannot be parametrized by parameters with a finite dimension. To solve it numerically, we can approximate κ with a deep neural network where θ is the neural network weights and biases. This approach has been shown to be very effective, especially when κ is a high dimensional mapping or lacks smoothness [39], compared to traditional approaches such as piecewise linear functions. Additionally, DNN representation also provides regularization thanks to its spatial correlation (e.g., if x 1 and x 2 are close, κ(x 1 ) and κ(x 2 ) are also close) and appropriate initialization.
Therefore, the optimization problem is reduced to a finite dimensional optimization problem, where the optimization variable is the DNN parameters: The PDE constraint in Equation (3) can be numerically solved using finite difference methods (or finite element methods), which results in a linear system Here u ∈ R n , f ∈ R n , A ∈ R n×n , and n is the degrees of freedoms. The reconstructed numerical solution u(x) depends on θ, although in general, we can not find an explicit expression for u(x) with respect to θ.
To solve the optimization problem Equation (3), we will calculate the gradients ∇ θ J(θ). In the context of reverse-mode automatic differentiation, this calculation requires us to back-propagate gradients through both the numerical solver and the neural network. In our previous work, we proposed a general framework for expressing both numerical simulations and DNNs as computational graphs, and thus unified the implementation of reverse-mode AD for deep neural networks and numerical solvers. This technique enables us to easily extract gradients ∇ θ J(θ) by reverse-mode automatic differentiation.

Parallel Computing
One challenge associated with reverse-mode AD is that all intermediate variables must be saved for gradient back-propagation. This requires additional memory compared to forward computation. Large-scale problems bring another difficulty to memory costs: we may not be able to solve even the forward problem, because it is challenging to save a large-sized discrete solution u on a single CPU.
One solution to ensure performance and scalability of numerical solvers is MPI. When we use MPI to solve Equation (3), the grid is split into multiple patches, and each machine (an MPI processor) owns one patch. We can rearrange the indices so that each MPI processor owns a stripe of rows in the matrix A and a segment of f with continuous indices. We resort to a distributed solver (e.g., from Hypre [38]) to solve the linear system Equation (4). After solving the linear system, each patch owns a segment of u, and a local loss function can be computed. The local loss functions are then aggregated to calculate the global loss function.
However, when we compute the gradients, reverse-mode automatic differentiation with MPI is challenging because reasoning about reverse rules is difficult while avoiding deadlocks, ensuring correctness, and improving performance. There are many researches that focus on adjoint MPI [30,40]. Many of them are tightly coupled with a specific application, and/or are not developed in the context of DNNs and PDEs.
Note that there are also active research in the deep learning community for distributed computing of DNNs. However, the parallel models are usually data-parallelism and do not consider communication within layers. For model parallelism in DNNs [41], usually DNNs are split into several parts and an asynchronous training method is used. These models are not appropriate for the inverse modeling problems that we are considering because communication does not happen within computational graphs. Instead, in our parallel model, the data communication pattern is more complicated: the data communication may occur multiple times within the computational graph, e.g., while performing time-stepping schemes or accumulating the loss function, etc. This difference calls for a new parallel model for both forward computation and gradient back-propagation.

Basic MPI Operators
ADCME provides a subset of the MPI routines to facilitate implementing computational graphs with data communication capabilities. Most of the routines have gradient back-propagation capabilities. For example, the mpi gather function corresponds to MPI Gather, except that mpi gather is differentiable. Internally, this API is implemented by observing that the scatter operation is the reverse of the gather operation in an MPI program. We show an excerpt from the mpi gather implementation in ADCME: Here MPIGather forward implements the forward computation, and broadcast a from the root to all worker nodes buffers out. MPIGather backward is the corresponding gradient back-propagation implementation. Having received grad out from the downstream of the computational graph, this mpi gather operator back-propagates the gradient to grad a by reversing MPI Gather. Table 1 presents parts of the built-in MPI operators in ADCME.
As an example on how to use these primitive operators, consider computing The parameter θ is first broadcasted onto 4 processors using mpi bcast. Processor i computed θ i locally. Then the results are summed on the root processor using mpi sum. An ADCME program is as follows (Figure 4) using ADCME mpi_init() # initialize MPI theta0 = placeholder(1.0) theta = mpi_bcast(theta0) l = theta^mpi_rank() L = mpi_sum(l) # initialize a Session sess = Session(); init(sess) L_value = run(sess, L) # only the value on rank 0 is valid mpi_finalize() # finalize MPI Receive a tensor from a given processor YES mpi send Send a tensor from a given processor YES mpi sendrecv Send a tensor from a given processor to another given processor YES mpi gather Concat tensors from all processors on the root processor YES  Additionally, we can use the following one-liner to extract gradients grad = gradients(L, theta0) This simple line of code hides the details for back-propagating gradients through local operations as well as MPI routines (mpi bcast, mpi sum).

Halo Exchange
In many applications, we use domain decomposition to distribute the workload onto different processors. To reduce the data communication overheads, subdomains usually overlap at the boundaries and exchange boundary data with their adjacent subdomains ( Figure 5). This communication pattern is referred to as halo exchange. Halo exchange allows for communicating only a small portion of data compared to the data used for stencil computation within the subdomain. Halo exchange appears in many applications, and thus is implemented as a standalone operator in ADCME. The implementation of forward computation of halo exchange is straightforward. A direct application of blocking sends and receives will cause deadlocks. Instead, we use nonblocking sends and receives ( Figure 6). For the reverse mode, the communication is also reversed for gradient backpropagation. The idea is that the forward sends correspond to reverse receives and the forward receives correspond to reverse sends. Additionally, the order should also be reversed. Once we work out the one-by-one substitution, the algorithm will work out correctly.

Distributed Sparse Matrices
ADCME uses Hypre as the distributed linear algebra backend for MPI matrix solvers. Each MPI processor owns a stripe of rows with continuous indices. The local matrix is stored in CSR (compressed sparse row) format and ADCME provides an interface mpi SparseTensor to hide the implementation details.
For reverse-mode automatic differentiation, matrix transposition is an operator that are common in gradient back-propagation. For example, assume the forward computation is (x is the input, y is the output, and A is a matrix) y = Ax Here ∂L(y) ∂y is a row vector, and therefore ∂L(y(x))) ∂x T = A T ∂L(y) ∂y T requires a matrix vector multiplication, where the matrix is A T . Given that Equation (6) is ubiquitous in numerical PDE schemes, a distributed implementation of parallel transposition is very important. As mentioned, ADCME distributes a continuous set of rows in a sparse matrix onto different MPI processors. The sub-matrices are stored in the CSR format. To transpose the sparse matrix in a parallel environment, we first split the matrices in each MPI processor into blocks and then use MPI Isend/MPI Irecv to exchange data. Note that this procedure consists of two phases: in the first phase, each block determines the number of nonzero entries to send, and sends the number to its corresponding receiver; in the second phase, each block actually sends all nonzero entries. Since the receiver already knows the number of entries expected to receive, it can prepare a buffer of an appropriate size. Finally, we transpose the matrices in place for each block. Using this method, we obtained a CSR representation of the transposed matrix ( Figure 10). The gradient back-propagation is also implemented accordingly by reversing all the communication steps.

Distributed Optimization
The objective function of our problem can be written as a sum of local objective functions where f i is the local objective functions, N is the number of processors. Despite many existing distributed optimization algorithm, in this work we adopt a simple approach: aggregating gradients ∇ θ f i (θ) and updating θ on the root processors. Note this is possible because in scientific computing applications, neural networks are not necessarily very large and therefore the dimension of θ is typically small. This approach allows us to apply sophisticated gradient-based optimization problem easily with existing off-the-shelf optimizers, such as L-BFGS [42], nonlinear conjugate gradient method [43], and ADAM optimizer [44]. Figure 8 shows how we can convert an existing optimizer to an MPI-enabled optimizer. The basic idea is to let the root processor notify worker processors whether to compute the loss function or the gradient. Then the root processor and workers will collaborate on executing the same routines and thus ensuring the correctness of collective MPI calls.

Hybrid Programming
Each MPI processor can communicate data between processes, which do not share memory. Within each process, ADCME also allows for multi-threaded parallelism with a shared-memory model. For example, we can use OpenMP to accelerate matrix vector production. We can also use a threadpool per process 3 to manage more complex and dynamic parallel tasks. However, the hybrid model brings challenges to communicate data using MPI. When we post MPI calls from different threads within the same process, we need to prevent data races and match the corresponding broadcast and collective operators. For example, without any guarantee on the ordering of concurrent MPI calls, we might incorrectly matched a send operator with a gather operator.
In ADCME, we adopt the dependency injection [45] technique: we explicitly serialize the MPI calls by adding ghost dependencies. For example, in the computational graph in Figure 9, originally, Operator 2 and Operator 3 are independent. In a concurrent computing environment, Rank 0 may execute Operator 2 first and then Operator 3, while Rank 1 executes Operator 3 first and then Operator 2. Then there is a mismatch of the MPI call (race condition): Operator 2 in Rank 0 coacts with Operator 3 in Rank 1, and Operator 3 in Rank 0 coacts with Operator 2 in Rank 1. To resolve the data race issue, we can explicitly make Operator 3 depend on Operator 2. In this way, we can ensure that the MPI calls Operator 1, 2, and 3 are executed in order. Note this technique sacrifices some concurrency (Operator 2 and Operator 3 cannot be executed concurrently), but the concurrency of most non-MPI operators is still preserved.

Benchmarks
The purpose of this section is to present the distributed computing capability of ADCME via MPI. With the MPI operators, ADCME is well suited to parallel applications on clusters with very large numbers of cores. We benchmark individual operators as well as the gradient calculation as a whole. Particularly, we use two metrics for measuring the scaling of the implementation: 1. Weak scaling, i.e., how the solution time varies with the number of processors for a fixed problem size per processor. 2. Strong scaling, i.e., the speedup for a fixed problem size with respect to the number of processors, and is governed by Amdahl's law.
For most operators, ADCME is just a wrapper of existing third-party parallel computing software (e.g., Hypre). However, for gradient back-propagation, some functions may be missing and are implemented anew. For example, in Hypre, distributed sparse matrices split into multiple stripes, where each MPI rank owns a stripe with continuous row indices. In gradient back-propagation, the transpose of the original matrix is usually needed and such functionalities are missing in Hyper. Note that ADCME uses hybrid parallel computing models, i.e., a mixture of multithreading programs and MPI communication; therefore, when we talk about one MPI processor, it may contain multiple CPU cores. The MPI processors may also be distributed in different CPUs/nodes, which are interconnected via network communications.
All the experiments are conducted on a cluster with 133 computational node, and each node has 32 cores in total (Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60GHz, 32 GB RAM). The MPI programs are verified with serial programs.

Transposition
This example benchmarks transposition of a distributed sparse matrix. The transposition operation is very important for gradient back-propagation. We consider solving the Poisson's equation Here f (x) ≡ 1, and κ(x) is approximated by a deep neural network where θ is the neural network weights and biases. The equation is discretized using finite difference method on a uniform grid and the discretization leads to a linear system Here u is the solution vector and f is the source vector. Note that A is a sparse matrix and its entries depend on θ.
The domain is split into multiple patches and each patch correspond to a strip of rows in A. The results are shown in Figure 10. In the plots, we show the strong scaling for a fixed matrix size of 25,000,000×25,000,000 as well as the weak scaling, where each MPI processor owns a patch of the whole mesh and has a size 300 × 300, i.e., 300 2 = 90000 rows. The run time for a fixed problem size is effectively reduced as we increase the number of processors, and no stagnation is observed for the current setting. For the weak scaling results, we see only 2∼3 times increase in runtime when we scale up to hundreds of processors.

Poisson's Equation
This example demonstrates the capability of ADCME to use Hypre for solving linear systems in parallel. We show that the overhead of ADCME is quite small compared to the linear solver time. We use the same problem setting as Section 4.1. We formulate the loss function using the numerical solution Here I consists of 80% randomly picked degrees of freedom for each patch. Figure 11 shows the computational model for solving the Poisson's equation. In the strong scaling experiments, we consider a fixed problem size 1800 × 1800 (mesh size, which implies the matrix size is around 32 million × 32 million). In the weak scaling experiments, each MPI processor owns a 300 × 300 block. For example, a problem with 3600 processors has the problem size 90000 × 3600 ≈ 0.3 billion. We first consider the weak scaling. Figure 12 shows the runtime for forward computation as well as gradient back-propagation. There are two important observations: 1. By using more cores per processor, the runtime is reduced significantly. For example, the runtime for the backward is reduced to around 10 seconds from 30 seconds by switching from 1 core to 4 cores per processor. 2. The runtime for the backward is typically less than twice the forward computation. Although the backward requires solve two linear systems (one of them is in the forward computation), the AMG (algebraic multigrid) linear solver in the back-propagation may converge faster, and therefore costs less than the forward.
Additionally, we show the overhead in Figure 13, which is defined as the difference between total runtime and Hypre linear solver time, of both the forward and backward calculation. We see that the overhead is quite small compared to the total time, especially when the problem size is large. This indicates that the ADCME MPI implementation is very effective. Now we consider the strong scaling. In this case, we fixed the whole problem size and split the mesh onto different MPI processors. Figure 14 shows the runtime for the forward computation and the gradient back-propagation. We can reduce the runtime by more than 20 times for the expensive gradient back-propagation by utilizing more than 100 MPI processors. Figure 15 shows the speedup and efficiency. We can see that the 4 cores have smaller runtime compared to 1 core. Finally, we apply the distributed L-BFGS optimizer to train the deep neural network. The distributed optimizer is constructed from a serial distributed optimizer using the approach described in Section 3.4. Figure 16 illustrate the results at 400-th iteration with 4 and 9 MPI processors respectively. Each processor owns a 300 × 300 patch. We can see that the estimated DNN-based κ(x) is very similar to the exact κ(x). Figure 17 also shows the corresponding loss function profiles.

Acoustic Wave Equation
In this example, we consider the acoustic wave equation with perfectly matched layer (PML) [46,17]. The governing equation for the acoustic equation is where u is the displacement, f is the source term, and c is the spatially varying acoustic velocity. Figure 18 shows the snapshots for the acoustic wave propagation and Figure 19 shows the true velocity model we have used in the acoustic and elastic wave equation examples.
In the inverse problem, only the wavefield u on the surface is observable, and we want to use this information to estimate c. The problem is usually ill-posed, so regularization techniques are usually used to constrain c. One approach is to represent c by a deep neural network c(x, y) = N θ (x, y)  where θ is the neural network weights and biases. The loss function is formulated by the square loss for the wavefield on the surface. The computational model is shown in Figure 21.
To implement an MPI version of the acoustic wave equation propagator, we use mpi halo exchange, which is implemented using MPI and performs the halo exchange mentioned in the last example for both wavefields and axilliary fields. This function communicates the boundary information for each block of the mesh. The following plot shows the computational graph for the numerical simulation of the acoustic wave equation Figure 21 shows the strong scaling and weak scaling of our implementation. Each processor consists of 32 processors, which are used at the discretion of ADCME's backend, i.e., TensorFlow. The strong scaling result is obtained by using a 1000 × 1000 grid and 100 times steps. For the weak scaling result, each MPI processor owns a 100 × 100 grid, and the total number of steps is 2000. It is remarkable that even though we increase the number of processors from 1 to 100, the total time only increases 2 times in the weak scaling.
We also show the speedup and efficiency for the strong scaling case in Figure 22. We can achieve more than 20 times acceleration by using 100 processors (3200 cores in total) and the trend is not slowing down at this scale.

Elastic Wave Equation
In the last example, we consider the elastic wave equation [47] ρ ∂v i ∂t where v is the velocity, σ is the stress tensor, ρ is the density, and λ and µ are the Lamé constants. Similar to the acoustic equation, we use the PML boundary conditions and have observations on the surface. However, the inversion parameters are now spatially varying ρ, λ and µ. Figure 23 shows snapshots of the wave propagation.
As an example, we approximate λ by a deep neural network λ(x, y) = N θ (x, y) and the other two parameters are kept fixed. We use the same geometry settings as the acoustic wave equation case. Note that elastic wave equation has more state variables as well as auxilliary fields, and thus is more memory demanding. The huge memory cost calls for a distributed framework, especially for large-scale problems.    Additionally, we use fourth-order finite difference scheme for discretizing Equation 2. This scheme requires us to exchange two layers on the boundaries for each block in the mesh. This data communication is implemented using MPI, i.e., mpi halo exchange2. Figure 24 shows the strong and weak scaling. Again, we see that the weak scaling of the implementation is quite effective because the runtime increases mildly even if we increase the number of processors from 1 to 100.  Figure 25 shows the speedup and efficiency for the strong scaling. We can achieve more than 20 times speedup when using 100 processors.

Conclusion
In this paper we proposed a framework for distributed optimization in scientific machine learning for solving inverse problems. The key is to express numerical simulation using a computational graph. The data communication operations are treated as nodes in the comptuational graph. This view makes the implementation of inverse modeling algorithms quite flexible (modular and testable) and conceptually simple. We also proposed a method to convert existing gradient-based optimization algorithms to MPI-enabled optimizers. These ideas are implemented in the ADCME library, which provides a set of MPI primitives, which can back-propagate gradients. We can also develop customized data communication operators to tailor to specific applications for better performance. To demonstrate the effectiveness of the proposed algorithm, we have trained a neural network that is coupled with either a Poisson's equation or a wave equation. The PDE is solved numerically and in a distributed way. The results show that our method provides a promising approach towards scalable inverse modeling.