Search

In this paper we develop and analyze a variational data assimilation method with efficient decoupled iterative numerical algorithms for the Stokes–Darcy equations with the Beavers–Joseph interface condition. By using Tikhonov regularization and formulating the variational data assimilation into an optimization problem, we establish the existence, uniqueness, and stability of the optimal solution. Based on the weak formulation of the Stokes–Darcy equations, the Lagrange multiplier rule is utilized to derive the first order optimality system for both the continuous and discrete variational data assimilation problems, where the discrete data assimilation is based on a finite element discretization in space and the backward Euler scheme in time. By rescaling the optimality system and then analyzing its corresponding bilinear forms, we prove the optimal finite element convergence rate with special attention paid to recovering uncertainties missed in the optimality system. To solve the discrete optimality system efficiently, three decoupled iterative algorithms are proposed to address the computational cost for both well-conditioned and ill-conditioned variational data assimilation problems, respectively. Finally, numerical results are provided to validate the proposed methods. [ABSTRACT FROM AUTHOR]

The Euler elastica (EE) model with surface curvature can generate artifact-free results compared with the traditional total variation regularization model in image processing. However, strong nonlinearity and singularity due to the curvature term in the EE model pose a great challenge for one to design fast and stable algorithms for the EE model. In this paper, we propose a new, fast, hybrid alternating minimization (HALM) algorithm for the EE model based on a bilinear decomposition of the gradient of the underlying image and prove the global convergence of the minimizing sequence generated by the algorithm under mild conditions. The HALM algorithm comprises three subminimization problems and each is either solved in the closed form or approximated by fast solvers, making the new algorithm highly accurate and efficient. We also discuss the extension of the HALM strategy to deal with general curvature-based variational models, especially with a Lipschitz smooth functional of the curvature. A host of numerical experiments are conducted to show that the new algorithm produces good results with much-improved efficiency compared to other state-of-the-art algorithms for the EE model. As one of the benchmarks, we show that the average running time of the HALM algorithm is at most one-quarter of that of the fast operator-splitting-based Deng--Glowinski--Tai algorithm. [ABSTRACT FROM AUTHOR]

The constrained pressure residual (CPR) algorithm is a family of well-known and industry-standard preconditioners for large-scale reservoir simulation. The CPR algorithm is a twostage preconditioner to deal with different blocks stage-by-stage, and is often able to effectively improve the robustness behavior and the convergence speed of linear iterations. Nonetheless, as a linear preconditioner, it is hard for the traditional CPR method to be capable of action at the nonlinear level for effectively solving large sparse nonlinear systems of equations with high nonlinearity. In this paper, we present and study the extension of this linear method to the nonlinearly CPR (NCPR) case for solving fractured reservoir problems, to deal with the difficulty of the slow convergence or stagnation from the nonlinear level. In the proposed nonlinear preconditioning, a subspace nonlinear block system is first built and solved to remove the unbalanced nonlinearities of the pressure and the saturation fields, and the fast convergence can then be restored when a variant of semismooth Newton methods is called after the subspace nonlinear block system is solved. Experiments on two or three dimensional porous media applications are presented to demonstrate the applicability and parallel scalability of the aforementioned NCPR method. We also show that the proposed algorithm is superior to the commonly used nonlinear algorithm in terms of the robustness and the positivity-preserving property. [ABSTRACT FROM AUTHOR]

Stochastic Parareal (SParareal) is a probabilistic variant of the popular parallel-in-time algorithm known as Parareal. Similarly to Parareal, it combines fine- and coarse-grained solutions to an ODE using a predictor-corrector (PC) scheme. The key difference is that carefully chosen random perturbations are added to the PC to try to accelerate the location of a stochastic solution to the ODE. In this paper, we derive superlinear and linear mean-square error bounds for SParareal applied to nonlinear systems of ODEs using different types of perturbations. We illustrate these bounds numerically on a linear system of ODEs and a scalar nonlinear ODE, showing a good match between theory and numerics. [ABSTRACT FROM AUTHOR]

We propose in this paper a new minimization algorithm based on a slightly modified version of the scalar auxialiary variable (SAV) approach coupled with a relaxation step and an adaptive strategy. It enjoys several distinct advantages over popular gradient based methods: (i) it is unconditionally energy diminishing with a modified energy which is intrinsically related to the original energy, and thus no parameter tuning is needed for stability; (ii) it allows the use of large step-sizes which can effectively accelerate the convergence rate. We also present a convergence analysis for some SAV based algorithms, which include our new algorithm without the relaxation step as a special case. We apply our new algorithm to several illustrative and benchmark problems and compare its performance with several popular gradient based methods. The numerical results indicate that the new algorithm is very robust, and its adaptive version usually converges significantly faster than those popular gradient decent based methods. [ABSTRACT FROM AUTHOR]

Policy gradient algorithms have been widely applied to Markov decision processes and reinforcement learning problems in recent years. Regularization with various entropy functions is often used to encourage exploration and improve stability. This paper proposes an approximate Newton method for the policy gradient algorithm with entropy regularization. In the case of Shannon entropy, the resulting algorithm reproduces the natural policy gradient algorithm. For other entropy functions, this method results in brand-new policy gradient algorithms. We prove that all these algorithms enjoy Newton-type quadratic convergence and that the corresponding gradient flow converges globally to the optimal solution. We use synthetic and industrial-scale examples to demonstrate that the proposed approximate Newton method typically converges in single-digit iterations, often orders of magnitude faster than other state-of-the-art algorithms. [ABSTRACT FROM AUTHOR]

This paper presents an extension of the Adaptive Antoulas-Anderson (AAA) algorithm for rational modeling. Specifically, our new stable multi-input multi-output AAA (smiAAA) algorithm builds rational approximations of multi-input signals with a common set of stable poles. A new methodology is presented for iteratively enforcing the stability constraint on the poles. We demonstrate the strengths of this approach compared to the stability enforcement in the FastAAA algorithm. Results using the smiAAA algorithm are compared with the commonly used Vector Fitting algorithm and the more recently published Rational Krylov fitting (RKFIT) algorithm. The Vector Fitting and RKFIT algorithms both require the user to input the number of poles to use in the approximations. If the final approximation is not accurate enough, the user must restart the Vector Fitting or RKFIT algorithm with a larger number of poles and/or a new starting location for the poles. In contrast, the smiAAA algorithm is designed to allow the user to simply input the desired accuracy of the approximations, and the necessary number of poles is detected automatically. This permits users to produce approximations of a desired accuracy with no knowledge about the underlying order of the system being approximated, preventing the algorithm from ever needing to be rerun. An additional feature for preventing extraneous poles from being returned by the AAA algorithm is also discussed. The cause of these extraneous poles is efficiently detected and removed by our presented methodology. The examples presented demonstrate that the smiAAA algorithm can efficiently produce approximations of similar or better accuracy than the Vector Fitting and RKFIT algorithms while requiring less input from the user. [ABSTRACT FROM AUTHOR]

In this paper, we present the new "asynchronous truncated multigrid-reduction-intime"" (AT-MGRIT) algorithm for introducing time parallelism to the solution of discretized timedependent problems. The new algorithm is based on the multigrid-reduction-in-time (MGRIT) approach, which, in certain settings, is equivalent to another common multilevel parallel-in-time method, Parareal. In contrast to Parareal and MGRIT that both consider a global temporal grid over the entire time interval on the coarsest level, the AT-MGRIT algorithm uses truncated local time grids on the coarsest level, each grid covering certain temporal subintervals. These local grids can be solved completely in an independent way from each other, which reduces the sequential part of the algorithm and, thus, increases parallelism in the method. Here, we study the effect of using truncated local coarse grids on the convergence of the algorithm, both theoretically and numerically, and show, using challenging nonlinear problems, that the new algorithm consistently outperforms classical Parareal/MGRIT in terms of time to solution. [ABSTRACT FROM AUTHOR]

The tensor-train (TT) decomposition expresses a tensor in a data-sparse format used in molecular simulations, high-order correlation functions, and optimization. In this paper, we propose four parallelizable algorithms that compute the TT format from various tensor inputs: (1) Parallel-TTSVD for traditional format, (2) PSTT and its variants for streaming data, (3) Tucker2TT for Tucker format, and (4) TT-fADI for solutions of Sylvester tensor equations. We provide theoretical guarantees of accuracy, parallelization methods, scaling analysis, and numerical results. For example, for a d-dimension tensor in Rn× ... × n, a two-sided sketching algorithm PSTT2 is shown to have a memory complexity of O (n d/2), improving upon O (nd 1) from previous algorithms. [ABSTRACT FROM AUTHOR]

This paper focuses on efficient algorithms for finding the Dantzig selector which was first proposed by Candès and Tao as an effective variable selection technique in the linear regression. This paper first reformulates the Dantzig selector problem as an equivalent convex composite optimization problem and proposes a semismooth Newton augmented Lagrangian (Ssnal) algorithm to solve the equivalent form. This paper also applies a proximal point dual semismooth Newton (PpdSsn) algorithm to solve another equivalent form of the Dantzig selector problem. Comprehensive results on the global convergence and local asymptotic superlinear convergence of the Ssnal and PpdSsn algorithms are characterized under very mild conditions. The computational costs of a semismooth Newton algorithm for solving the subproblems involved in the Ssnal and PpdSsn algorithms can be cheap by fully exploiting the second order sparsity and employing efficient techniques. Numerical experiments on the Dantzig selector problem with synthetic and real data sets demonstrate that the Ssnal and PpdSsn algorithms substantially outperform the state-of-the-art first order algorithms even for the required low accuracy, and the proposed algorithms are able to solve the large-scale problems robustly and efficiently to a relatively high accuracy. [ABSTRACT FROM AUTHOR]

This paper proposes some efficient and accurate adaptive two-grid (ATG) finite element algorithms for linear and nonlinear PDEs. The main idea of these algorithms is to utilize the solutions on the kth-level adaptive meshes to find the solutions on the (k + 1)th-level adaptive meshes which are constructed by performing adaptive element bisections on the k th-level adaptive meshes. These algorithms transform nonsymmetric positive definite (non-SPD) PDEs (resp., nonlinear PDEs) into symmetric positive definite (SPD) PDEs (resp., linear PDEs). The proposed algorithms are both accurate and efficient due to the following advantages: They do not need to solve the nonsymmetric or nonlinear systems; the degrees of freedom are very small, they are easily implemented, and the interpolation errors are very small. Next, this paper constructs residual-type a posteriori error estimators which are shown to be reliable and efficient. The key ingredient in proving the efficiency is to establish an upper bound of the oscillation terms, which may not be higher-order terms due to the low regularity of the numerical solution. Furthermore, the convergence of the algorithms is proved when bisection is used for the mesh refinements. Finally, numerical experiments are provided to verify the accuracy and efficiency of the ATG finite element algorithms compared to regular adaptive finite element algorithms and two-grid finite element algorithms [J. Xu, SIAM J. Numer. Anal., 33 (1996), pp. 1759-1777]. [ABSTRACT FROM AUTHOR]

In this paper, an easy to implement and efficient algorithm is proposed for large-scale sparse least squares problems. It reduces the original large-scale problem into a sequence of small-scale subproblems by utilizing the "sparsity" of the gradient. A well-designed FISTA algorithm is presented to solve the subproblems. With the desired termination criterion, the algorithm is proved to globally converge and converge locally at linear rate. Moreover, a complexity of O(mpln1ε) is given for the presented algorithm. Numerical comparisons between the presented algorithms and a number of state-of-the-art algorithms on solving LASSO problems with large-scale LIBSVM data sets and image deblurring problems demonstrate the robustness and high performance of the presented algorithms. As an example, the presented algorithm only needs 9 seconds to solve a LASSO problem with 19,264,097 samples and 29,890,095 features. [ABSTRACT FROM AUTHOR]

The "fast iterative shrinkage-thresholding algorithm,"" a.k.a. FISTA, is one of the most well known first-order optimization scheme in the literature, as it achieves the worst-case O(1/k2) optimal convergence rate for objective function value. However, despite such an optimal theoretical convergence rate, in practice the (local) oscillatory behavior of FISTA often damps its efficiency. Over the past years, various efforts have been made in the literature to improve the practical performance of FISTA, such as monotone FISTA, restarting FISTA, and backtracking strategies. In this paper, we propose a simple yet effective modification to the original FISTA scheme which has two advantages: It allows us to (1) prove the convergence of generated sequence and (2) design a so-called lazy-start strategy which can be up to an order faster than the original scheme. Moreover, we propose novel adaptive and greedy strategies which further improve the practical performance. The advantages of the proposed schemes are tested through problems arising from inverse problems, machine learning, and signal/image processing. [ABSTRACT FROM AUTHOR]

In this paper we consider wave propagation problems in two-dimensional unbounded domains, including dissipative effects, reformulated in terms of space-time boundary integral equations. For their solution, we employ a convolution quadrature (CQ) for the temporal and a Galerkin boundary element method (BEM) for the spatial discretization. It is known that one of the main advantages of the CQ-BEMs is the use of the FFT algorithm to retrieve the discrete time integral operators with an optimal linear complexity in time, up to a logarithmic term. It is also known that a key ingredient for the success of such methods is the efficient and accurate evaluation of all the integrals that define the matrix entries associated to the full space-time discretization. This topic has been successfully addressed when standard Lagrangian basis functions are considered for the space discretization. However, it results that, for such a choice of the basis, the BEM matrices are in general fully populated, a drawback that prevents the application of CQ-BEMs to large-scale problems. In this paper, as a possible remedy to reduce the global complexity of the method, we consider approximant functions of wavelet type. In particular, we propose a numerical procedure that, by taking advantage of the fast wavelet transform, allows us on the one hand to compute the matrix entries associated to the choice of wavelet basis functions by maintaining the accuracy of those associated to the Lagrangian basis ones and, on the other hand, to generate sparse matrices without the need of storing a priori the fully populated ones. Such an approach allows in principle the use of wavelet basis of any type and order, combined with CQ based on any stable ordinary differential equations solver. Several numerical results, showing the accuracy of the solution and the gain in terms of computer memory saving, are presented and discussed. [ABSTRACT FROM AUTHOR]

In this paper, we develop the ball-approximated characteristics (B-char) method, which is an algorithm for efficiently implementing characteristic-based schemes in two and three dimensions. Core to the implementation of numerical schemes is the evaluation of integrals, which in the context of characteristic-based schemes with piecewise constant approximations boils down to computing the intersections between two regions. In the literature, these regions are approximated by polytopes (polygons in two dimensions and polyhedra in three dimensions) and, due to this, the implementation in three dimensions is nontrivial. The main novelty in this paper is the approximation of the regions by balls, whose intersections are much cheaper to compute than those of polytopes. Of course, balls cannot fully tessellate a region, and hence some mass may be lost. We perform some adjustments, and also solve an optimization problem, in order to yield a scheme that is both locally and globally mass conserving. This algorithm can achieve results that are similar to those obtained from an implementation which uses polytopal intersections, with a much cheaper computational cost. [ABSTRACT FROM AUTHOR]

This paper presents a multilevel convergence framework for multigrid-reduction-intime (MGRIT) as a generalization of previous two-grid estimates. The framework provides a priori upper bounds on the convergence of MGRIT V- and F-cycles, with different relaxation schemes, by deriving the respective residual and error propagation operators. The residual and error operators are functions of the time-stepping operator, analyzed directly and bounded in the norm, both numerically and analytically. We present various upper bounds of different computational cost and varying sharpness. These upper bounds are complemented by proposing analytic formulae for the approximate convergence factor of V-cycle algorithms that take the number of fine grid time points, the temporal coarsening factors, and the eigenvalues of the time-stepping operator as parameters. The paper concludes with supporting numerical investigations of parabolic (anisotropic diffusion) and hyperbolic (wave equation) model problems. We assess the sharpness of the bounds and the quality of the approximate convergence factors. Observations from these numerical investigations demonstrate the value of the proposed multilevel convergence framework for estimating MGRIT convergence a priori and for the design of a convergent algorithm. We further highlight that observations in the literature are captured by the theory, including that two-level Parareal and multilevel MGRIT with F-relaxation do not yield scalable algorithms and the benefit of a stronger relaxation scheme. An important observation is that with increasing numbers of levels MGRIT convergence deteriorates for the hyperbolic model problem, while constant convergence factors can be achieved for the diffusion equation. The theory also indicates that L-stable Runge--Kutta schemes are more amendable to multilevel parallel-in-time integration with MGRIT than A-stable Runge--Kutta schemes. [ABSTRACT FROM AUTHOR]

Two asymptotically exact a posteriori error estimates are proposed for eigenvalues by the nonconforming Crouzeix--Raviart and enriched Crouzeix--Raviart elements. The main challenge in the design of such error estimators comes from the nonconformity of the finite element spaces used. Such nonconformity causes two difficulties: the first is the construction of high accuracy gradient recovery algorithms, and the second is a computable high accuracy approximation of a consistency error term. The first difficulty was solved for both nonconforming elements in a previous paper. Two methods are proposed to solve the second difficulty in the present paper. In particular, this solution allows the use of high accuracy gradient recovery techniques. Further, a postprocessing algorithm is designed by utilizing asymptotically exact a posteriori error estimators to construct the weights of a combination of two approximate eigenvalues. This algorithm requires solving only one eigenvalue problem and admits high accuracy eigenvalue approximations both theoretically and numerically. [ABSTRACT FROM AUTHOR]

The need to solve l¹ regularized linear problems can be motivated by various compressive sensing and sparsity related techniques for data analysis and signal or image processing. These problems lead to nonsmooth convex optimization in high dimensions. Theoretical works predict a sharp phase transition for the exact recovery of compressive sensing problems. Our numerical experiments show that state-of-the-art algorithms are not effective enough to observe this phase transition accurately. This paper proposes a simple formalism that enables us to produce an algorithm that computes an l¹ minimizer under the constraints Au = b up to the machine precision. In addition, a numerical comparison with standard algorithms available in the literature is exhibited. The comparison shows that our algorithm compares advantageously with other state-of-the-art methods, both in terms of accuracy and efficiency. With our algorithm, the aforementioned phase transition is observed at high precision. [ABSTRACT FROM AUTHOR]

In this paper, we propose a parallel-in-time algorithm for approximately solving parabolic equations. In particular, we apply the k-step backward differentiation formula and then develop an iterative solver by using the waveform relaxation technique. Each resulting iteration represents a periodic-like system, which could be further solved in parallel by using the diagonalization technique. The convergence of the waveform relaxation iteration is theoretically examined by using the generating function method. The argument could be further applied to the time-fractional subdiffusion equation, whose discretization shares common properties of the standard BDF methods due to the nonlocality of the fractional differential operator. Illustrative numerical results are presented to complement the theoretical analysis. [ABSTRACT FROM AUTHOR]

Least squares form one of the most prominent classes of optimization problems with numerous applications in scientific computing and data fitting. When such formulations aim at modeling complex systems, the optimization process must account for nonlinear dynamics by incorporating constraints. In addition, these systems often incorporate a large number of variables, which increases the difficulty of the problem and motivates the need for efficient algorithms amenable to large-scale implementations. In this paper, we propose and analyze a Levenberg--Marquardt algorithm for nonlinear least squares subject to nonlinear equality constraints. Our algorithm is based on inexact solves of linear least-squares problems that only require Jacobian-vector products. Global convergence is guaranteed by the combination of a composite step approach and a nonmonotone step acceptance rule. We illustrate the performance of our method on several test cases from data assimilation and inverse problems; our algorithm is able to reach the vicinity of a solution from an arbitrary starting point and can outperform the most natural alternatives for these classes of problems. [ABSTRACT FROM AUTHOR]

The paper considers supervised learning problems of labeled data with grouped input features. The groups are nonoverlapped such that the model coefficients corresponding to the input features form disjoint groups. The coefficients have group sparsity structure in the sense that coefficients corresponding to each group shall be simultaneously either zero or nonzero. To make effective use of such group sparsity structure given a priori, we introduce a novel log-composite regularizer, which can be minimized by an iterative algorithm. In particular, our algorithm iteratively solves for a traditional group least absolute shrinkage and selection operator (LASSO) problem that involves summing up the l2 norm of each group until convergent. By updating group weights, our approach enforces a group of smaller coefficients from the previous iterate to be more likely to set to zero compared to the group LASSO. Theoretical results include a minimizing property of the proposed model as well as the convergence of the iterative algorithm to a stationary solution under mild conditions. We conduct extensive experiments on synthetic and real datasets, indicating that our method yields a performance that is superior to that of the state-of-the-art methods in linear regression and binary classification. [ABSTRACT FROM AUTHOR]

When using Laguerre and Hermite spectral methods to numerically solve PDEs in unbounded domains, the number of collocation points assigned inside the region of interest is often insufficient, particularly when the region is expanded or translated in order to safely capture the unknown solution. Simply increasing the number of collocation points cannot ensure a fast convergence to spectral accuracy. In this paper, we propose a scaling technique and a moving technique to adaptively cluster enough collocation points in a region of interest in order to achieve fast spectral convergence. Our scaling algorithm employs an indicator in the frequency domain that both is used to determine when scaling is needed and informs the tuning of a scaling factor to redistribute collocation points in order to adapt to the diffusive behavior of the solution. Our moving technique adopts an exterior-error indicator and moves the collocation points to capture the translation. Both frequency and exterior-error indicators are defined using only the numerical solutions. We apply our methods to a number of different models, including diffusive and moving Fermi-Dirac distributions and nonlinear Dirac solitary waves, and demonstrate recovery of spectral convergence for time-dependent simulations. A performance comparison in solving a linear parabolic problem shows that our frequency scaling algorithm outperforms the existing scaling approaches. We also show our frequency scaling technique is able to track the blowup of average cell sizes in a model for cell proliferation. In addition to the Laguerre and Hermite basis functions with exponential decay at infinity, we also successfully apply the frequency-dependent scaling technique into rational basis functions with algebraic decay at infinity. [ABSTRACT FROM AUTHOR]

Sparsity regularization has been widely adopted in many fields, such as signal and image processing and machine learning. In this paper, we mainly consider nonconvex minimization problems involving three terms, for applications such as sparse signal recovery and low rank matrix recovery. We employ a three-operator splitting proposed by Davis and Yin [Set-Valued Var. Anal., 25 (2017), pp. 829-858] (namely, DYS) to solve the resulting possibly nonconvex problems and develop the convergence theory for this three-operator splitting algorithm in the nonconvex case. We show that when the step size is chosen less than a computable threshold, the whole sequence converges to a stationary point. By defining a new decreasing energy function associated with the DYS method, we establish the global convergence of the whole sequence and a local convergence rate under an additional assumption that F, G, and H are semialgebraic. We also provide sufficient conditions for the boundedness of the generated sequence. Finally, some numerical experiments are conducted to compare the DYS algorithm with some classical efficient algorithms for sparse signal recovery and low rank matrix completion. The numerical results indicate that DYS outperforms the existing methods for these applications. [ABSTRACT FROM AUTHOR]

When solving partial differential equations (PDEs) with random fields as coefficients, the efficient sampling of random field realizations can be challenging. In this paper we focus on the fast sampling of Gaussian fields using quasi-random points in a finite element and multilevel quasi Monte Carlo (MLQMC) setting. Our method uses the stochastic PDE (SPDE) approach of Lindgren et al. combined with a new fast algorithm for white noise sampling which is tailored to (ML)QMC. We express white noise as a wavelet series expansion that we divide into two parts. The first part is sampled using quasi-random points and contains a finite number of terms in order of decaying importance to ensure good quasi Monte Carlo (QMC) convergence. The second part is a correction term which is sampled using standard pseudo-random numbers. We show how the sampling of both terms can be performed in linear time and memory complexity in the number of mesh cells via a supermesh construction, yielding an overall linear cost. Furthermore, our technique can be used to enforce the MLQMC coupling even in the case of nonnested mesh hierarchies. We demonstrate the efficacy of our method with numerical experiments. [ABSTRACT FROM AUTHOR]

In many multiphysics applications, an ultimate quantity of interest can be written as a linear functional of the solution to the discretized governing nonlinear partial differential equations and finding a sufficiently accurate pointwise solution may be regarded as a step toward that end. In this paper, we derive a posteriori approximate error bounds for linear functionals corresponding to quantities of interest using two kinds of nonlinear preconditioning techniques. Nonlinear preconditioning, such as the inexact Newton with backtracking and nonlinear elimination algorithm and the multiplicative Schwarz preconditioned inexact Newton algorithm, may be effective in improving global convergence for Newton's method. It may prevent stagnation of the nonlinear residual norm and reduce the number of solutions of large ill-conditioned linear systems involving a global Jacobian required at each nonlinear iteration. We illustrate the effectiveness of the new bounds using canonical nonlinear PDE models: a flame sheet model and a nonlinear coupled lid-driven cavity problem. [ABSTRACT FROM AUTHOR]

In ptychography experiments, redundant scanning is usually required to guarantee the stable recovery, such that a huge number of frames are generated, and thus it poses a great demand on parallel computing to solve this large-scale inverse problem. In this paper, we propose the overlapping domain decomposition methods to solve the nonconvex optimization problem in ptychographic imaging. They decouple the problem defined on the whole domain into subproblems only defined on the subdomains with synchronizing information in the overlapping regions of these subdomains, thus leading to highly parallel algorithms with good load balance. More specifically, for the nonblind recovery (with known probe in advance), by enforcing the continuity of the overlapping regions for the image (sample), the nonlinear optimization model is established based on a novel smooth-truncated amplitude-Gaussian metric (ST-AGM). Such a metric allows for fast calculation of the proximal mapping with closed form, and meanwhile provides the possibility for the convergence guarantee of the first-order nonconvex optimization algorithm due to its Lipschitz smoothness. Then the alternating direction method of multipliers is utilized to generate an efficient overlapping domain decomposition based ptychography algorithm (OD2P) for the two-subdomain domain decomposition (DD), where all subproblems can be computed with closed-form solutions. Due to the Lipschitz continuity for the gradient of the objective function with ST-AGM, the convergence of the proposed OD2P is derived under mild conditions. Moreover, it is extended to more general cases including multiple-subdomain DD and blind recovery. Numerical experiments are further conducted to show the performance of proposed algorithms, demonstrating good convergence speed and robustness to the noise. Especially, we report the virtual wall-clock time of proposed algorithm up to 10 processors, which shows potential for upcoming massively parallel computations. [ABSTRACT FROM AUTHOR]

In this paper, we consider a nonlinear and nonlocal parabolic model for multispecies ionic fluids and introduce a semi-implicit finite volume scheme, which is second-order accurate in space, is first-order in time, and satisfies the positivity preserving and mass conservation properties. The scheme also satisfies an energy dissipation relation at the semidiscrete level. In addition, our scheme involves a fast algorithm on the convolution terms with singular but integrable kernels on the uniform grid, which otherwise impedes the accuracy and efficiency of the whole scheme. Error estimates on the fast convolution algorithm are shown next. Numerous numerical tests are provided to demonstrate the properties, such as unconditional stability, the order of convergence, energy dissipation, and the complexity of the fast convolution algorithm. Furthermore, extensive numerical experiments are carried out to explore the modeling effects in specific examples, such as the steric repulsion, concentration of ions at the boundary, and blowup phenomenon of Keller--Segel equations. [ABSTRACT FROM AUTHOR]

Large-scale optimization problems arising from the discretization of problems involving PDEs sometimes admit solutions that can be well approximated by low-rank matrices. In this paper, we will exploit this low-rank approximation property by solving the optimization problem directly over the set of low-rank matrices. In particular, we introduce a new multilevel algorithm, where the optimization variable is constrained to the Riemannian manifold of fixed-rank matrices. In contrast to most other multilevel algorithms, where the rank is chosen adaptively on each level in order to control the perturbation due to the low-rank truncation, we can keep the ranks (and thus the computational complexity) fixed throughout the iterations. Furthermore, classical implementations of line searches based on Wolfe conditions enable computing a solution where the numerical accuracy is limited to about the square root of the machine epsilon. Here, we propose an extension to Riemannian manifolds of the line search of Hager and Zhang, which uses approximate Wolfe conditions that enable computing a solution on the order of the machine epsilon. Numerical experiments demonstrate the computational efficiency of the proposed framework. [ABSTRACT FROM AUTHOR]

Topology optimization problems often support multiple local minima due to a lack of convexity. Typically, gradient-based techniques combined with continuation in model parameters are used to promote convergence to more optimal solutions; however, these methods can fail even in the simplest cases. In this paper, we present an algorithm to perform a systematic exploratory search for the solutions of the optimization problem via second order methods without a good initial guess. The algorithm combines the techniques of deflation, barrier methods, and primal-dual active set solvers in a novel way. We demonstrate this approach on several numerical examples, observe mesh independence in certain cases and show that multiple distinct local minima can be recovered. [ABSTRACT FROM AUTHOR]

This paper presents a recently developed particle simulation code package PIFE-PIC, which is a novel three-dimensional (3-D) parallel immersed finite element (IFE) particle-in-cell (PIC) simulation model for particle simulations of plasma-material interactions. This framework is based on the recently developed nonhomogeneous electrostatic IFE-PIC algorithm, which is designed to handle complex plasma-material interface conditions associated with irregular geometries using a Cartesian mesh-based PIC. Three-dimensional domain decomposition is utilized for both the electrostatic field solver with IFE and the particle operations in PIC to distribute the computation among multiple processors. A simulation of the orbital motion-limited (OML) sheath of a dielectric sphere immersed in a stationary plasma is carried out to validate parallel IFE-PIC and profile the parallel performance of the code package. Furthermore, a large-scale simulation of plasma charging at a lunar crater containing 2 million PIC cells (10 million FE/IFE cells) and about 1 billion particles, running for 20,000 PIC steps in about 154 wall-clock hours, is presented to demonstrate the high-performance computing capability of PIFE-PIC. [ABSTRACT FROM AUTHOR]

In this paper, we address the numerical solution of the optimal transport problem on undirected weighted graphs, taking the shortest path distance as transport cost. The optimal solution is obtained from the long-time limit of the gradient descent dynamics. Among different time stepping procedures for the discretization of this dynamics, a backward Euler time stepping scheme combined with the inexact Newton-Raphson method results in a robust and accurate approach for the solution of the optimal transport problem on graphs. It is found experimentally that the algorithm requires solving between O(1) and O(m0:36) linear systems involving weighted Laplacian matrices, where m is the number of edges. These linear systems are solved via algebraic multigrid methods, resulting in an efficient solver for the optimal transport problem on graphs. [ABSTRACT FROM AUTHOR]

This paper proposes a parallel domain decomposition Rayleigh--Ritz projection scheme to compute a selected number of eigenvalues (and, optionally, associated eigenvectors) of large and sparse symmetric pencils. The projection subspace associated with interface variables is built by computing a few of the eigenvectors and associated leading derivatives of a zeroth-order approximation of the nonlinear matrix-valued interface operator. On the other hand, the projection subspace associated with interior variables is built independently in each subdomain by exploiting local eigenmodes and matrix resolvent approximations. The sought eigenpairs are then approximated by a Rayleigh--Ritz projection onto the subspace formed by the union of these two subspaces. Several theoretical and practical details are discussed, and upper bounds of the approximation errors are provided. Our numerical experiments demonstrate the efficiency of the proposed technique on sequential/distributed memory architectures as well as its competitiveness against schemes such as shift-and-invert Lanczos and automated multilevel substructuring combined with p-way vertex-based partitionings. [ABSTRACT FROM AUTHOR]

Randomized iterative algorithms have recently been proposed to solve large-scale linear systems. In this paper, we present a simple randomized extended average block Kaczmarz algorithm that exponentially converges in the mean square to the unique minimum norm least squares solution of a given linear system of equations. The proposed algorithm is pseudoinversefree and therefore different from the projection-based randomized double block Kaczmarz algorithm of Needell, Zhao, and Zouzias [Linear Algebra Appl., 484 (2015), pp. 322--343]. We emphasize that our method works for all types of linear systems (consistent or inconsistent, overdetermined or underdetermined, full-rank or rank-deficient). Moreover, our approach can be implemented for parallel computation, yielding remarkable improvements in computational time. Numerical examples are given to show the efficiency of the new algorithm. [ABSTRACT FROM AUTHOR]

A framework, PhaseLift, was recently proposed to solve the phase retrieval problem. In this framework, the problem is solved by optimizing a cost function over the set of complex Hermitian positive semidefinite matrices. This approach to phase retrieval motivates a more general consideration of optimizing cost functions on semidefinite Hermitian matrices where the desired minimizers are known to have low rank. This paper considers an approach based on an alternative cost function defined on a union of appropriate manifolds. It is related to the original cost function in a manner that preserves the ability to find a global minimizer and is significantly more efficient computationally. A rank-based optimality condition for stationary points is given and optimization algorithms based on state-of-the-art Riemannian optimization and dynamically reducing rank are proposed. Empirical evaluations are performed using the PhaseLift problem. The new approach is shown to be an effective method of phase retrieval with computational efficiency increased substantially compared to a state-of-the-art algorithm, the Wirtinger flow algorithm. A preliminary version of this paper can be found in [W. Huang, K. A. Gallivan, and X. Zhang, Procedia Comput. Sci., 80 (2016), pp. 1125-1134]. [ABSTRACT FROM AUTHOR]

In this paper, we present a novel numerical method that redistributes unevenly given points on an evolving closed curve to satisfy equi-arclength(-like) condition. Without substantial difficulty, it is also capable of remeshing or employing adaptive mesh refinement. The key idea is to find the discrete inverse of the arclength(-like) function in the framework of the Fourier spectral method to obtain overall spectral accuracy. Both equi-arclength and curvature-dependent redistributions are extensively studied, and their spectral accuracy is verified by application to smoothly perturbed points on various curves. We further confirm that our method converges even for the points being perturbed nonsmoothly and randomly. To leverage the robustness of our method, a remeshing technique is applied in which the accuracy is not affected. Application to a periodic planar curve without any modification of our algorithm is also discussed. Then, to show the practical applicability, an evolving curve with large deformation is studied by coupling with point redistribution and remeshing in various flows such as mean curvature flow, Willmore flow, and Stokes flow. [ABSTRACT FROM AUTHOR]

We consider the problem of optimizing an unknown function over a multidimensional continuous domain, where function evaluation is noisy and expensive. We assume that a globally accurate model of the function is not available, but there exist some parametric models that can well approximate the function in local regions. In this paper, we propose an algorithm in the optimal learning framework that learns the shape of the function and finds the optimal design with a limited number of measurements. We construct belief functions using a radial basis function--based local approximation technique, and use the knowledge gradient policy to decide where to measure, aiming at maximizing the value of information from each measurement. Experiments on both synthetic test problems and a real materials science application show the strong performance of our algorithm. [ABSTRACT FROM AUTHOR]

In this paper, we propose a new efficient preconditioner for iteratively solving the large-scale indefinite saddle-point sparse linear system, which arises from discretizing the optimality system in optimal control problems of wave equations with a one-shot second-order finite difference scheme in both space and time. The proposed preconditioner can be implemented in a parallel- in-time (PinT) manner via a carefully designed unitary diagonalization decomposition. Such an explicit unitary diagonalization is rarely seen in the literature. We also analyze the eigenvalue bounds of the preconditioned system, which are shown to be highly clustered around one. Moreover, a simple splitting algorithm that alternates between a linear complementarity problem (LCP) and a quasi-Newton iteration is discussed for handling the case with control constraints. Within the quasi- Newton iteration, our proposed PinT preconditioner can be directly used in preconditioning the Jacobian system of the same structure. Both 1D and 2D numerical examples are given to illustrate the promising convergence performance of our proposed PinT preconditioner in comparison with a recently proposed matching Schur complement (MSC) preconditioner. [ABSTRACT FROM AUTHOR]

In this paper, we develop a novel spectral method for solving the time-dependent Kohn{Sham equations in the semiclassical regime. Our strategy is to use a single-step predictor- corrector algorithm, where the propagator is derived using a Fourier integral operator commonly known as the frozen Gaussian approximation (FGA) ansatz. In the case of laser potentials, we derived a simplified FGA to avoid the high dimensional integration. Numerical examples are provided to verify applicability, as well as efficiency and accuracy of the scheme. [ABSTRACT FROM AUTHOR]

The Strang splitting method, formally of order two, can suffer from order reduction when applied to semilinear parabolic problems with inhomogeneous boundary conditions. The recent work [L. Einkemmer and A. Ostermann, SIAM J. Sci. Comput., 37, 2015; SIAM J. Sci. Comput., 38, 2016] introduces a modification of the method to avoid the reduction of order based on the nonlinearity. In this paper we introduce a new correction constructed directly from the flow of the nonlinearity and which requires no evaluation of the source term or its derivatives. The goal is twofold. One, this new modification requires only one evaluation of the diffusion flow and one evaluation of the source term flow at each step of the algorithm and it reduces the computational effort to construct the correction. Second, numerical experiments suggest it is well suited in the case where the nonlinearity is stiff. We provide a convergence analysis of the method for a smooth nonlinearity and perform numerical experiments to illustrate the performances of the new approach. [ABSTRACT FROM AUTHOR]

Abstract. The discrete empirical interpolation method (DEIM) is a popular technique for nonlinear model reduction, and it has two main ingredients: an interpolating basis that is computed from a collection of snapshots of the solution, and a set of indices which determine the nonlinear components to be simulated. The computation of these two ingredients dominates the overall cost of the DEIM algorithm. To specifically address these two issues, we present randomized versions of the DEIM algorithm. There are three main contributions of this paper. First, we use randomized range finding algorithms to efficiently find an approximate DEIM basis. Second, we develop randomized subset selection tools, based on leverage scores, to efficiently select the nonlinear components. Third, we develop several theoretical results that quantify the accuracy of the randomization on the DEIM approximation. We also present numerical experiments that demonstrate the benefits of the proposed algorithms. [ABSTRACT FROM AUTHOR]

This paper introduces a fast and numerically stable algorithm for the solution of fourth-order linear boundary value problems on an interval. This type of equation arises in a variety of settings in physics and signal processing. Our method reformulates the equation as a collection of second-kind integral equations defined on local subdomains. Each such equation can be stably discretized and solved. The boundary values of these local solutions are matched by solving a banded linear system. The method of deferred corrections is then used to increase the accuracy of the scheme. Deferred corrections require applying the integral operator to a function on the entire domain, for which we provide an algorithm with linear cost. We illustrate the performance of our method on several numerical examples. [ABSTRACT FROM AUTHOR]

Accurate interpolation and approximation techniques for functions with discontinuities are key tools in many applications, such as medical imaging. In this paper, we study a radial basis function type of method for scattered data interpolation that incorporates discontinuities via a variable scaling function. For the construction of the discontinuous basis of kernel functions, information on the edges of the interpolated function is necessary. We characterize the native space spanned by these kernel functions and study error bounds in terms of the fill distance of the node set. To extract the location of the discontinuities, we use a segmentation method based on a classi- fication algorithm from machine learning. The results of the conducted numerical experiments are in line with the theoretically derived convergence rates in case that the discontinuities are a priori known. Further, an application to interpolation in magnetic particle imaging shows that the presented method is very promising in order to obtain edge-preserving image reconstructions in which ringing artifacts are reduced. [ABSTRACT FROM AUTHOR]

In many inverse problems, a nonnegativity constraint is natural. Moreover, in some cases, we expect the vector of unknown parameters to have zero components. When a Bayesian approach is taken, this motivates a desire for prior probability density (and hence posterior probability density) functions that have positive mass at the boundary of the set {x∊RN|x&#88040}. Unfortunately, it is difficult to define a prior with this property that yields computationally tractable inference for large-scale inverse problems. In this paper, we use nonnegativity constrained optimization to define such prior and posterior density functions when the measurement error is either Gaussian or Poisson distributed. The numerical optimization methods we use are highly efficient, and hence our approach is computationally tractable even in large-scale cases. We embed our nonnegativity constrained optimization approach within a hierarchical framework, obtaining Gibbs samplers for both Gaussian and Poisson distributed measurement cases. Finally, we test the resulting Markov chain Monte Carlo methods on examples from both image deblurring and positron emission tomography. [ABSTRACT FROM AUTHOR]

In the context of big data, high-order tensor decompositions have to face a new challenge in terms of storage and computational costs. The tensor train (TT) decomposition provides a very useful graph-based model reduction, whose storage cost grows linearly with the tensor order D. The computation of the TT-core tensors and TT-ranks can be done in a stable sequential (i.e., noniterative) way thanks to the popular TT-SVD algorithm. In this paper, we exploit the ideas developed for the hierarchical/tree Tucker decomposition in the context of the TT decomposition. Specifically, a new efficient estimation scheme, called TT-HSVD (Tensor-Train Hierarchical SVD), is proposed as a solution to compute the TT decomposition of a high-order tensor. The new algorithm simultaneously delivers the TT-core tensors and their TT-ranks in a hierarchical way. It is a stable (i.e., noniterative) and computationally more efficient algorithm than TT-SVD, which is very important when dealing with large-scale data. The TT-HSVD algorithm uses a new reshaping strategy and a tailored partial SVD, which allows us to deal with smaller matrices compared to those of the TT-SVD. In addition, TT-HSVD is well suited for a parallel processing architecture. An algebraic analysis of the two algorithms is carried out, showing that TT-SVD and TT-HSVD compute the same TT-ranks and TT-core tensors up to specific bases. Simulation results for different tensor orders and dimensions corroborate the effectiveness of the proposed algorithm. [ABSTRACT FROM AUTHOR]

Many PDEs involving fractional Laplacian are naturally set in unbounded domains with underlying solutions decaying slowly and subject to certain power law. Their numerical solutions are underexplored. This paper aims at developing accurate spectral methods using rational basis (or modified mapped Gegenbauer functions) for such models in unbounded domains. The main building block of the spectral algorithms is the explicit representations for the Fourier transform and fractional Laplacian of the rational basis, derived from some useful integral identities related to modified Bessel functions. With these at our disposal, we can construct rational spectral-Galerkin and direct collocation schemes by precomputing the associated fractional differentiation matrices. We obtain optimal error estimates of rational spectral approximation in the fractional Sobolev spaces and analyze the optimal convergence of the proposed Galerkin scheme. We also provide ample numerical results to show that the rational method outperforms the Hermite function approach. [ABSTRACT FROM AUTHOR]

Orthogonal maps are the solutions of the mathematical model of paper-folding, also called the origami problem. They consist of a system of first-order fully nonlinear equations involving the gradient of the solution. The Dirichlet problem for orthogonal maps is considered here. A variational approach is advocated for the numerical approximation of the maps. The introduction of a suitable objective function allows us to enforce the uniqueness of the solution. A strategy based on a splitting algorithm for the corresponding flow problem is presented and leads to decoupling the time-dependent problem into a sequence of local nonlinear problems and a global linear variational problem at each time step. Numerical experiments validate the accuracy and the efficiency of the method for various domains and meshes. [ABSTRACT FROM AUTHOR]

In this paper, a fast multipole method (FMM) is proposed to compute long-range interactions of wave sources embedded in 3-dimensional (3-D) layered media. The layered media Green's function for the Helmholtz equation, which satisfies the transmission conditions at material interfaces, is decomposed into a free space component and four types of reaction field components arising from wave reflections and transmissions through the layered media. The proposed algorithm is a combination of the classic FMM for the free space component and FMMs specifically designed for the four types of reaction components, made possible by new multipole expansions (MEs) and local expansions (LEs) as well as the multipole-to-local translation (M2L) operators for the reaction field components. Moreover, an equivalent polarization source can be defined for each reaction component based on the convergence analysis of its ME. The FMMs for the reaction components, implemented with the target particles and equivalent polarization sources, are found to be much more efficient than the classic FMM for the free space component due to the fact that the equivalent polarization sources and the target particles are always separated by a material interface. As a result, the FMM algorithm developed for layered media has a similar computational cost as that for the free space. Numerical results validate the fast convergence of the MEs and the O(N) complexity of the FMM for interactions of low-frequency wave sources in 3-D layered media. [ABSTRACT FROM AUTHOR]

In this paper, a novel pre- and post-processing algorithm is presented that can double the convergence rate of finite element approximations for linear wave problems. In particular, it is shown that a q-step pre- and post-processing algorithm can improve the convergence rate of the finite element approximation from order p+1 to order p+1+q in the L2-norm and from order p to order p + q in the energy norm, in both cases up to a maximum of order 2p, with p the polynomial degree of the finite element space. The q-step pre- and post-processing algorithms only need to be applied once and require solving at most q linear systems of equations. The biggest advantage of the proposed method compared to other post-processing methods is that it does not suffer from convergence rate loss when using unstructured meshes. Other advantages are that this new pre- and post-processing method is straightforward to implement, incorporates boundary conditions naturally, and does not lose accuracy near boundaries or strong inhomogeneities in the domain. Numerical examples illustrate the improved accuracy and higher convergence rates when using this method. In particular, they confirm that 2p-order convergence rates in the energy norm are obtained, even when using unstructured meshes or when solving problems involving heterogeneous materials and curved boundaries. [ABSTRACT FROM AUTHOR]

Gaussian kernels can be an eficient and accurate tool for multivariate interpolation. For smooth functions, high accuracies are often achieved near the at limit where the interpolation matrix becomes increasingly ill-conditioned. Stable evaluation algorithms for isotropic Gaussians (Gaussian radial basis functions) have been proposed based on a Chebyshev expansion of the Gaus- sians by Fornberg, Larsson, and Flyer and based on a Mercer expansion with Hermite polynomials by Fasshauer and McCourt. In this paper, we propose a new stabilization algorithm for the multi- variate interpolation with isotropic or anisotropic Gaussians for an arbitrary number of dimensions derived from the generating function of the Hermite polynomials. We also derive and analyze a new analytic cutoff criterion for the generating function expansion that allows us to automatically adjust the number of stabilizing basis functions. [ABSTRACT FROM AUTHOR]

In this paper, we demonstrate that many of the computational tools for univariate orthogonal polynomials have analogues for a family of bivariate orthogonal polynomials on the triangle, including Clenshaw's algorithm and sparse differentiation operators. This allows us to derive a practical spectral method for solving linear partial differential equations on triangles with sparse discretizations. We can thereby rapidly solve partial differential equations using polynomials with degrees in the thousands, resulting in sparse discretizations with as many as several million degrees of freedom. [ABSTRACT FROM AUTHOR]