Search

The article focuses on research reproducibility in scientific computing journals, Journal on Scientific Computing (SISC), Journal of Scientific Computing (JSC), and Journal of Computational Physics (JCP), by requiring authors to publicly share their data and credibility in scientific research.

The weighted essentially nonoscillatory (WENO) schemes are popular high-order numerical methods for hyperbolic conservation laws. When dealing with hyperbolic systems, WENO schemes are usually used in cooperation with the local characteristic decomposition, as the componentwise WENO reconstruction/interpolation procedure often produces oscillatory approximations near shocks. In this paper, we investigate local characteristic decomposition-free WENO schemes for a special class of hyperbolic systems endowed with a coordinate system of Riemann invariants. We apply the WENO procedure to the coordinate system of Riemann invariants instead of the local characteristic fields to save the expensive computational cost on local characteristic decomposition but meanwhile maintain the essentially nonoscillatory performance. Due to the nonlinear algebraic relation between the Riemann invariants and conserved variables, it is difficult to obtain the cell averages of Riemann invariants directly from those of conserved variables and vice versa; thus, we do not use the finite volume WENO schemes in this work. The same difficulty is also faced in the traditional Shu-Osher lemma [C.-W. Shu and S. Osher, J. Comput. Phys., 83 (1989), pp. 32-78]-based finite difference schemes, as the computation of fluxes is based on reconstruction as well. Therefore, we adopt the alternative formulation of the finite difference WENO scheme [Y. Jiang, C.-W. Shu, and M. Zhang, SIAM J. Sci. Comput., 35 (2013), pp. A1137-A1160, C.-W. Shu and S. Osher, J. Comput. Phys., 77 (1988), pp. 439-471] in this paper, which is based on interpolation for nodal values. The efficiency and good performance of our method are demonstrated by extensive numerical tests which indicate that the coordinate system of Riemann invariants is a good alternative of local characteristic fields for the WENO procedure. [ABSTRACT FROM AUTHOR]

This paper presents a fast high-order method for the solution of two-dimensional problems of scattering by penetrable inhomogeneous media, with application to high-frequency configurations containing (possibly) discontinuous refractivities. The method relies on a hybrid direct/iterative combination of (1) a differential volumetric formulation (which is based on the use of appropriate Chebyshev differentiation matrices enacting the Laplace operator) and (2) a secondkind boundary integral formulation (which, once again, utilizes Chebyshev discretization, but, in this case, in the boundary integral context). The approach enjoys low dispersion and high-order accuracy for smooth refractivities, as well as second-order accuracy (while maintaining low dispersion) in the discontinuous refractivity case. The solution approach proceeds by application of impedance-to-impedance (ItI) maps to couple the volumetric and boundary discretizations. The volumetric linear algebra solutions are obtained by means of a multifrontal solver, and the coupling with the boundary integral formulation is achieved via an application of the iterative linear algebra solver GMRES. In particular, the existence and uniqueness theory presented in the present paper provides an affirmative answer to an open question concerning the existence of a uniquely solvable second-kind ItI-based formulation for the overall scattering problem under consideration. Relying on a modestly demanding scatterer-dependent precomputation stage (requiring in practice a computing cost of the order of O(Nα ) operations, with α ≈ 1.07, for an N-point discretization and for the relevant Chebyshev accuracy orders q used), together with fast (O(N)-cost) single-core runs for each incident field considered, the proposed algorithm can effectively solve scattering problems for large and complex objects possibly containing discontinuities and strong refractivity contrasts. [ABSTRACT FROM AUTHOR]

This paper tackles efficient methods for Bayesian inverse problems with priors based on Whittle–Matérn Gaussian random fields. The Whittle–Matérn prior is characterized by a mean function and a covariance operator that is taken as a negative power of an elliptic differential operator. This approach is flexible in that it can incorporate a wide range of prior information including nonstationary effects, but it is currently computationally advantageous only for integer values of the exponent. In this paper, we derive an efficient method for handling all admissible noninteger values of the exponent. The method first discretizes the covariance operator using finite elements and quadrature, and uses preconditioned Krylov subspace solvers for shifted linear systems to efficiently apply the resulting covariance matrix to a vector. This approach can be used for generating samples from the distribution in two different ways: by solving a stochastic partial differential equation, and by using a truncated Karhunen–Loève expansion. We show how to incorporate this prior representation into the infinite-dimensional Bayesian formulation, and show how to efficiently compute the maximum a posteriori estimate, and approximate the posterior variance. Although the focus of this paper is on Bayesian inverse problems, the techniques developed here are applicable to solving systems with fractional Laplacians and Gaussian random fields. Numerical experiments demonstrate the performance and scalability of the solvers and their applicability to model and real-data inverse problems in tomography and a time-dependent heat equation. [ABSTRACT FROM AUTHOR]

This paper is focused on the approximation of the Euler equations of compressible fluid dynamics on a staggered mesh. With this aim, the flow parameters are described by the velocity, the density, and the internal energy. The thermodynamic quantities are described on the elements of the mesh, and thus the approximation is only in L2, while the kinematic quantities are globally continuous. The method is general in the sense that the thermodynamic and kinetic parameters are described by an arbitrary degree of polynomials. In practice, the difference between the degrees of the kinematic parameters and the thermodynamic ones is set to 1. The integration in time is done using the forward Euler method but can be extended straightforwardly to higher-order methods. In order to guarantee that the limit solution will be a weak solution of the problem, we introduce a general correction method in the spirit of the Lagrangian staggered method described in [R. Abgrall and S. Tokareva, SIAM J. Sci. Comput., 39 (2017), pp. A2345--A2364; R. Abgrall, K. Lipnikov, N. Morgan, and S. Tokareva, SIAM J. Sci. Comput., 2 (2020), pp. A343--A370; V. A. Dobrev, T. V. Kolev, and R. N. Rieben, SIAM J. Sci. Comput., 34 (2012), pp. B606--B641], and we prove a Lax--Wendroff theorem. The proof is valid for multidimensional versions of the scheme, even though most of the numerical illustrations in this work, on classical benchmark problems, are one-dimensional because we have easy access to the exact solution for comparison. We conclude by explaining that the method is general and can be used in different settings, for example, finite volume or discontinuous Galerkin method, not just the specific one presented in this paper. [ABSTRACT FROM AUTHOR]

In this paper, we propose a conservative low rank tensor method to approximate nonlinear Vlasov solutions. The low rank approach is based on our earlier work [W. Guo and J.-M. Qiu, A Low Rank Tensor Representation of Linear Transport and Nonlinear Vlasov Solutions and Their Associated Flow Maps, preprint, https://arxiv.org/abs/2106.08834, 2021]. It takes advantage of the fact that the differential operators in the Vlasov equation are tensor friendly, based on which we propose to dynamically and adaptively build up low rank solution basis by adding new basis functions from discretization of the differential equation, and removing basis from a singular value decomposition (SVD)-type truncation procedure. For the discretization, we adopt a high order finite difference spatial discretization together with a second order strong stability preserving multistep time discretization. While the SVD truncation will remove the redundancy in representing the high dimensional Vlasov solution, it will destroy the conservation properties of the associated full conservative scheme. In this paper, we develop a conservative truncation procedure with conservation of mass, momentum, and kinetic energy densities. The conservative truncation is achieved by an orthogonal projection onto a subspace spanned by 1, v, and v2 in the velocity space associated with a weighted inner product. Then the algorithm performs a weighted SVD truncation of the remainder, which involves a scaling, followed by the standard SVD truncation and rescaling back. The algorithm is further developed in high dimensions with hierarchical Tucker tensor decomposition of high dimensional Vlasov solutions, overcoming the curse of dimensionality. An extensive set of nonlinear Vlasov examples are performed to show the effectiveness and conservation property of proposed conservative low rank approach. Comparison is performed against the nonconservative low rank tensor approach on conservation history of mass, momentum, and energy. [ABSTRACT FROM AUTHOR]

Low Reynolds number fluid flows are governed by the Stokes equations. In two dimensions, Stokes flows can be described by two analytic functions, known as Goursat functions. Brubeck and Trefethen [SIAM J. Sci. Comput., 44 (2022), pp. A1205-A1226] recently introduced a lightning Stokes solver that uses rational functions to approximate the Goursat functions in polyg- onal domains. In this paper, we present the LARS algorithm (lightning-AAA rational Stokes) for computing two-dimensional (2D) Stokes flows in domains with smooth boundaries and multiply con- nected domains using lightning and AAA rational approximation [Y. Nakatsukasa, O. Sète, and L. N. Trefethen, SIAM J. Sci. Comput., 40 (2018), pp. A1494-A1522]. After validating our solver against known analytical solutions, we solve a variety of 2D Stokes flow problems with physical and engineering applications. Using these examples, we show rational approximation can now be used to compute 2D Stokes flows in general domains. The computations take less than a second and give solutions with at least 6-digit accuracy. [ABSTRACT FROM AUTHOR]

In this paper, we demonstrate the efficiency of using semi-Lagrangian discontinuous Galerkin methods to solve the drift-kinetic equation using graphic processing units (GPUs). In this setting we propose a second order splitting scheme and a two-dimensional semi-Lagrangian scheme in the poloidal plane. The resulting method is able to conserve mass up to machine precision, allows us to take large time steps due to the absence of a CFL condition, and provides local data dependency which is essential to obtain good performance on state-of-the-art high-performance computing systems. We report simulations of a drift-kinetic ion temperature gradient instability and show that our implementation achieves a performance of up to 600 GB/s on an A100 GPU. [ABSTRACT FROM AUTHOR]

In this paper, we propose a direct probing method to solve the inverse problem of the Eikonal equation. This problem involves the determination of the inhomogeneous wave-speed distribution from first-arrival time data at measurement surfaces corresponding to distributed point sources. The viscosity solution of the point-source Eikonal equation represents the least traveltime of wave fields from the source to the point at the high-frequency limit. We show that this inverse problem is highly ill-posed. To address this issue, we develop a direct probing method that incorporates solution analysis of the Eikonal equation and several aspects of the velocity models. Specifically, we use the filtered back-projection method to reconstruct the inhomogeneous wave-speed distribution when it has a small variation from the homogeneous medium. For the scenarios involving highcontrast media, we assume a background medium and develop an adjoint-based back-projection method to identify the variations of the medium from the assumed background. [ABSTRACT FROM AUTHOR]

In this paper, we consider an approach to improve the performance of exponential Runge–Kutta integrators and Lawson schemes in cases where the solution of a related, but usually much simpler, problem can be computed efficiently. While for implicit methods such an approach is common (e.g., by using preconditioners), for exponential integrators this has proven more challenging. Here we propose to extract a constant coefficient differential operator from the semilinear advection-diffusion-reaction equation for which, in many situations, efficient methods are known to compute the required matrix functions. Both a linear stability analysis and extensive numerical experiments show that the resulting schemes can be unconditionally stable. In fact, we find that exponential integrators of Runge–Kutta type and Lawson schemes can have better stability properties than similarly constructed implicit-explicit schemes. We also derive two new Lawson-type integrators that further improve on these stability properties. The overall effectiveness of the approach is highlighted by a number of performance comparisons on examples in two and three space dimensions. [ABSTRACT FROM AUTHOR]

Hamilton--Jacobi partial differential equations (HJ PDEs) have deep connections with a wide range of fields, including optimal control, differential games, and imaging sciences. By considering the time variable to be a higher dimensional quantity, HJ PDEs can be extended to the multitime case. In this paper, we establish a novel theoretical connection between specific optimization problems arising in machine learning and the multitime Hopf formula, which corresponds to a representation of the solution to certain multitime HJ PDEs. Through this connection, we increase the interpretability of the training process of certain machine learning applications by showing that when we solve these learning problems, we also solve a multitime HJ PDE and, by extension, its corresponding optimal control problem. As a first exploration of this connection, we develop the relation between the regularized linear regression problem and the linear quadratic regulator (LQR). We then leverage our theoretical connection to adapt standard LQR solvers (namely, those based on the Riccati ordinary differential equations) to design new training approaches for machine learning. Finally, we provide some numerical examples that demonstrate the versatility and possible computational advantages of our Riccati-based approach in the context of continual learning, post-training calibration, transfer learning, and sparse dynamics identification. [ABSTRACT FROM AUTHOR]

In this paper, we propose two methods for multivariate Hermite interpolation of manifold-valued functions. On the one hand, we approach the problem via computing suitable weighted Riemannian barycenters. To satisfy the conditions for Hermite interpolation, the sampled derivative information is converted into a condition on the derivatives of the associated weight functions. It turns out that this requires the solution of linear systems of equations, but no vector transport is necessary. This approach treats all given sample data points equally and is intrinsic in the sense that it does not depend on local coordinates or embeddings. As an alternative, we consider Hermite interpolation in a tangent space. This is a straightforward approach, where one designated point, for example, one of the sample points or (one of) their center(s) of mass, is chosen to act as the base point to which the tangent space is attached. The remaining sampled locations and sampled derivatives are mapped to said tangent space. This requires a vector transport between different tangent spaces. The actual interpolation is then conducted via classical vector space operations. The interpolant depends on the selected base point. The validity and performance of both approaches is illustrated by means of numerical examples. [ABSTRACT FROM AUTHOR]

In this paper, we develop a new type of Runge–Kutta (RK) discontinuous Galerkin (DG) method for solving hyperbolic conservation laws. Compared with the original RKDG method, the new method features improved compactness and allows simple boundary treatment. The key idea is to hybridize two different spatial operators in an explicit RK scheme, utilizing local projected derivatives for inner RK stages and the usual DG spatial discretization for the final stage only. Limiters are applied only at the final stage for the control of spurious oscillations. We also explore the connections between our method and Lax–Wendroff DG schemes and ADER-DG schemes. Numerical examples are given to confirm that the new RKDG method is as accurate as the original RKDG method, while being more compact, for problems including two-dimensional Euler equations for compressible gas dynamics. [ABSTRACT FROM AUTHOR]

In this paper, we present the two-dimensional unstructured grids extension of the a posteriori local subcell correction of discontinuous Galerkin (DG) schemes introduced in [F. Vilar, J. Comput. Phys., 387 (2018), pp. 245–279]. The technique is based on the reformulation of the DG scheme as a finite-volume (FV)-like method through the definition of some specific numerical fluxes referred to as reconstructed fluxes. A high-order DG numerical scheme combined with this new local subcell correction technique ensures positivity preservation of the solution, along with a low oscillatory and sharp shocks representation. The main idea of this correction procedure is to retain as much as possible of the high accuracy and the very precise subcell resolution of DG schemes, while ensuring the robustness and stability of the numerical method, as well as preserving physical admissibility of the solution. Consequently, an a posteriori correction will only be applied locally at the subcell scale where it is needed, but still ensuring the scheme conservativity. Practically, at each time step, we compute a DG candidate solution and check if this solution is admissible (for instance positive, non-oscillating, …). If it is the case, we go further in time. Otherwise, we return to the previous time step and correct locally, at the subcell scale, the numerical solution. To this end, each cell is subdivided into subcells. Then, if the solution is locally detected as bad, we substitute the DG reconstructed flux on the subcell boundaries by a robust first-order numerical flux. For a subcell detected as admissible, we keep the high-order DG reconstructed flux which allows us to retain the very highly accurate resolution and conservation of the DG scheme. As a consequence, only the solution inside troubled subcells and its first neighbors will have to be recomputed; elsewhere, the solution remains unchanged. Another technique blending in a convex combination fashion DG reconstructed fluxes and first-order FV fluxes for admissible subcells in the vicinity of troubled areas will also be presented and prove to improve results in comparison to the original algorithm introduced in [F. Vilar, J. Comput. Phys., 387 (2018), pp. 245–279]. Numerical results on various type of problems and test cases will be presented to assess the very good performance of the designed correction algorithm. [ABSTRACT FROM AUTHOR]

Probabilistic domain decomposition (PDD) is an alternative paradigm for solving boundary value problems (BVPs) in parallel with excellent scalability properties, thanks to its reliance on stochastic representations of the BVP. However, there are cases when the latter is less numerically convenient, or unknown. Semilinear elliptic BVPs and the Helmholtz equation are prominent examples of either class. In this paper, we overcome this issue by designing suitable iterative schemes for either problem. These schemes not only retain the desirable properties of PDD but also are optimally suited for pathwise variance reduction, resulting in a systematic, nearly cost-free reduction of the statistical error through the iterations. Numerical tests carried out on the supercomputer Marconi100 are presented. [ABSTRACT FROM AUTHOR]

The infamous numerical sign problem poses a fundamental obstacle to particle-based stochastic Wigner simulations in high-dimensional phase space. Although the existing particle annihilation (PA) via uniform mesh significantly alleviates the sign problem when dimensionality D ≤ 4, the mesh size grows dramatically when D ≥ 6 due to the curse of dimensionality and consequently makes the annihilation very inefficient. In this paper, we propose an adaptive PA algorithm, termed sequential-clustering particle annihilation via discrepancy estimation (SPADE), to overcome the sign problem. SPADE follows a divide-and-conquer strategy: adaptive clustering of particles via controlling their number-theoretic discrepancies and independent random matching in each cluster. The target is to alleviate the oversampling problem induced by the overpartitioning of phase space and to capture the nonclassicality of the Wigner function simultaneously. Combining SPADE with the variance reduction technique based on the stationary phase approximation, we attempt to simulate the proton-electron couplings in six- and 12-dimensional phase space. A thorough performance benchmark of SPADE is provided with the reference solutions in six-dimensional phase space produced by a characteristic-spectral-mixed scheme under a 73³ × 80³ uniform grid, which fully explores the limit of grid-based deterministic Wigner solvers. [ABSTRACT FROM AUTHOR]

Most numerical approximations of frequency-domain wave propagation problems suffer from the so-called dispersion error, which is the fact that plane waves at the discrete level oscillate at a frequency different from the continuous one. In this paper, we introduce a new technique to reduce the dispersion error in general finite difference (FD) schemes for frequency-domain wave propagation using the Helmholtz equation as a guiding example. Our method is based on the introduction of a shifted wavenumber in the FD stencil which we use to reduce the numerical dispersion for large enough numbers of grid points per wavelength (or for small enough meshsize), and thus we call the method asymptotic dispersion correction. The advantage of this technique is that the asymptotically optimal shift can be determined in closed form by computing the extrema of a function over a compact set. For one-dimensional Helmholtz equations, we prove that the standard 3-point stencil with shifted wavenumber does not have any dispersion error, and that the so-called pollution effect is completely suppressed. For higher dimensional Helmholtz problems, we give easy-to-use closed form formulas for the asymptotically optimal shift associated to the second-order 5-point scheme and a sixth-order 9-point scheme in two dimensions, and the 7-point scheme in three dimensions that yield substantially less dispersion error than their standard (unshifted) version. We illustrate this also with numerical experiments. [ABSTRACT FROM AUTHOR]

Entropy inequalities are crucial to the well-posedness of hyperbolic conservation laws, which help to select the physically meaningful one from among the infinite many weak solutions. Recently, several high order discontinuous Galerkin (DG) methods satisfying entropy inequalities were proposed; see [T. Chen and C.-W. Shu, J. Comput. Phys., 345 (2017), pp. 427-461; J. Chan, J. Comput. Phys., 362 (2018), pp. 346-374; T. Chen and C.-W. Shu, CSIAM Trans. Appl. Math., 1 (2020), pp. 1-52] and the references therein. However, high order numerical methods typically generate spurious oscillations in the presence of shock discontinuities. In this paper, we construct a high order entropy stable essentially oscillation-free DG (OFDG) method for hyperbolic conservation laws. With some suitable modification on the high order damping term introduced in [J. Lu, Y. Liu, and C.-W. Shu, SIAM J. Numer. Anal., 59 (2021), pp. 1299-1324; Y. Liu, J. Lu, and C.-W. Shu, SIAM J. Sci. Comput., 44 (2022), pp. A230-A259], we are able to construct an OFDG scheme with dissipative entropy. It is challenging to make the damping term compatible with the current entropy stable DG framework, that is, the damping term should be dissipative for any given entropy function without compromising high order accuracy. The key ingredient is to utilize the convexity of the entropy function and the orthogonality of the projection. Then the proposed method maintains the same properties of conservation, error estimates, and entropy dissipation as the original entropy stable DG method. Extensive numerical experiments are presented to validate the theoretical findings and the capability of controlling spurious oscillations of the proposed algorithm. [ABSTRACT FROM AUTHOR]

Hierarchical matrices provide a powerful representation for significantly reducing the computational complexity associated with dense kernel matrices. For example, the fast multipole method (FMM) and its variants are highly efficient when the kernel function is related to fundamental solutions of classical elliptic PDEs. For general kernel functions, interpolation-based methods are widely used for the efficient construction of hierarchical matrices. In this paper, we present a fast hierarchical data reduction (HiDR) procedure with O(n) complexity for the memoryefficient construction of hierarchical matrices with nested bases where n is the number of data points. HiDR aims to reduce the given data in a hierarchical way so as to obtain O(1) representations for all nearfield and farfield interactions. Based on HiDR, a linear complexity H² matrix construction algorithm is proposed. The use of data-driven methods enables better efficiency than other general-purpose methods and flexible computation without accessing the kernel function. Experiments demonstrate significantly improved memory efficiency of the proposed data-driven method compared to interpolation-based methods over a wide range of kernels. For the Coulomb kernel, the proposed general-purpose algorithm offers competitive performance compared to FMM and its variants, such as PVFMM. The data-driven approach not only works for general kernels but also leads to much smaller precomputation costs compared to PVFMM. [ABSTRACT FROM AUTHOR]

This paper concerns numerical solution of the diffusion equation with strong anisotropy on meshes not aligned with the anisotropic vector field. In order to resolve the numerical pollution for simulations on a non-anisotropy-aligned mesh and reduce the associated high computational cost we propose an effective preconditioner, extending our previous work [D. Green et al., Comput. Phys. Commun., 9 (2022), 108333]. Similar to the anisotropy-aligned mesh case, we apply the auxiliary space preconditioning framework to design a preconditioner where a continuous finite element space is used as the auxiliary space for the discontinuous finite element space. The key component is an effective line smoother that can mitigate the high-frequency errors perpendicular to the magnetic field. We design a graph-based approach to find such a line smoother that is approximately perpendicular to the vector fields when the mesh does not align with the anisotropy. Numerical experiments for several benchmark problems are presented, demonstrating the effectiveness and robustness of the proposed preconditioner when applied to Krylov iterative methods. [ABSTRACT FROM AUTHOR]

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 GMRES algorithm of Saad and Schultz [SIAM J. Sci. Stat. Comput., 7 (1986), pp. 856-869] is an iterative method for approximately solving linear systems Ax = b, with initial guess x0 and residual r0 = b Ax0. The algorithm employs the Arnoldi process to generate the Krylov basis vectors (the columns of Vκ). It is well known that this process can be viewed as a QR factorization of the matrix Bκ = [ r0, AVκ] at each iteration. Despite an O (ε)κ (Bκ) loss of orthogonality, for unit roundoff ε and condition number κ, the modified Gram-Schmidt formulation was shown to be backward stable in the seminal paper by Paige et al. [SIAM J. Matrix Anal. Appl., 28 (2006), pp. 264-284]. We present an iterated Gauss-Seidel formulation of the GMRES algorithm (IGS-GMRES) based on the ideas of Ruhe [Linear Algebra Appl., 52 (1983), pp. 591-601] and Świrydowicz et al. [Numer. Linear Algebra Appl., 28 (2020), pp. 1-20]. IGS-GMRES maintains orthogonality to the level O (ε)κ (Bκ) or O (ε), depending on the choice of one or two iterations; for two Gauss-Seidel iterations, the computed Krylov basis vectors remain orthogonal to working accuracy and the smallest singular value of Vκ remains close to one. The resulting GMRES method is thus backward stable. We show that IGS-GMRES can be implemented with only a single synchronization point per iteration, making it relevant to large-scale parallel computing environments. We also demonstrate that, unlike MGS-GMRES, in IGS-GMRES the relative Arnoldi residual corresponding to the computed approximate solution no longer stagnates above machine precision even for highly nonnormal systems. [ABSTRACT FROM AUTHOR]

We address the problem of approximating an unknown function from its discrete samples given at arbitrarily scattered sites. This problem is essential in numerical sciences, where modern applications also highlight the need for a solution to the case of functions with manifold values. In this paper, we introduce and analyze a combination of kernel-based quasi-interpolation and multiscale approximations for both scalar- and manifold-valued functions. While quasi-interpolation provides a powerful tool for approximation problems if the data is defined on infinite grids, the situation is more complicated when it comes to scattered data. Here, higher-order quasi-interpolation schemes either require derivative information or become numerically unstable. Hence, this paper principally studies the improvement achieved by combining quasi-interpolation with a multiscale technique. The main contributions of this paper are as follows. First, we introduce the multiscale quasi-interpolation technique for scalar-valued functions. Second, we show how this technique can be carried over using moving least-squares operators to the manifold-valued setting. Third, we give a mathematical proof that converging quasi-interpolation will also lead to converging multiscale quasi-interpolation. Fourth, we provide ample numerical evidence that multiscale quasi-interpolation has superior convergence to quasi-interpolation. In addition, we will provide examples showing that the multiscale quasi-interpolation approach offers a powerful tool for many data analysis tasks, such as denoising and anomaly detection. It is especially attractive for cases of massive data points and high dimensionality. [ABSTRACT FROM AUTHOR]

This paper presents the use of element-based algebraic multigrid (AMGe) hierarchies, implemented in the Parallel Element Agglomeration Algebraic Multigrid Upscaling and Solvers (ParELAG) library, to produce multilevel preconditioners and solvers for H(curl) and H(div) formulations. ParELAG constructs hierarchies of compatible nested spaces, forming an exact de Rham sequence on each level. This allows the application of hybrid smoothers on all levels and the Auxiliary-Space Maxwell Solver or the Auxiliary-Space Divergence Solver on the coarsest levels, obtaining complete multigrid cycles. Numerical results are presented, showing the parallel performance of the proposed methods. As a part of the exposition, this paper demonstrates some of the capabilities of ParELAG and outlines some of the components and procedures within the library. [ABSTRACT FROM AUTHOR]

This paper studies large-scale optimization problems on Riemannian manifolds whose objective function is a finite sum of negative log-probability losses. Such problems arise in various machine learning and signal processing applications. By introducing the notion of Fisher information matrix in the manifold setting, we propose a novel Riemannian natural gradient method, which can be viewed as a natural extension of the natural gradient method from the Euclidean setting to the manifold setting. We establish the almost-sure global convergence of our proposed method under standard assumptions. Moreover, we show that if the loss function satisfies certain convexity and smoothness conditions and the input-output map satisfies a Riemannian Jacobian stability condition, then our proposed method enjoys a local linear--or, under the Lipschitz continuity of the Riemannian Jacobian of the input-output map, even quadratic--rate of convergence. We then prove that the Riemannian Jacobian stability condition will be satisfied by a two-layer fully connected neural network with batch normalization with high probability, provided that the width of the network is sufficiently large. This demonstrates the practical relevance of our convergence rate result. Numerical experiments on applications arising from machine learning demonstrate the advantages of the proposed method over state-of-the-art ones. [ABSTRACT FROM AUTHOR]

This paper proposes and analyzes a novel efficient high-order finite volume method for the ideal magnetohydrodynamics (MHD). As a distinctive feature, the method simultaneously preserves two critical physical constraints: a discretely divergence-free (DDF) constraint on the magnetic field and the positivity-preserving (PP) property, which ensures the positivity of density, pressure, and internal energy. To enforce the DDF condition in each cell, we design a new discrete projection approach that projects the reconstructed point values at the cell interface into a DDF space, without using any approximation polynomials. This projection method is highly efficient, easy to implement, and particularly suitable for the high-order finite volume methods that return only the point values (no explicit approximation polynomials) in the reconstruction. Moreover, we also develop a new finite volume framework for constructing provably PP schemes for the ideal MHD system. The framework comprises the discrete projection technique, a suitable approximation to the Godunov--Powell source terms, and a simple PP limiter. We provide rigorous analysis of the PP property of the proposed finite volume method, demonstrating that the DDF condition and the proper approximation to the source terms eliminate the impact of magnetic divergence terms on the PP property. The analysis is challenging due to the internal energy function's nonlinearity and the intricate relationship between the DDF and PP properties. To address these challenges, we adopt the recently developed geometric quasilinearization approach [K. Wu and C.-W. Shu, SIAM Rev., 65 (2023), pp. 1031--1073], which transforms a nonlinear constraint into a family of linear constraints. Finally, we validate the effectiveness of the proposed method through several benchmark and demanding numerical examples. The results demonstrate that the proposed method is robust, accurate, and highly effective, confirming the significance of the proposed DDF projection and PP techniques. [ABSTRACT FROM AUTHOR]

Randomized matrix algorithms have become workhorse tools in scientific computing and machine learning. To use these algorithms safely in applications, they should be coupled with posterior error estimates to assess the quality of the output. To meet this need, this paper proposes two diagnostics: a leave-one-out error estimator for randomized low-rank approximations and a jackknife resampling method to estimate the variance of the output of a randomized matrix computation. Both of these diagnostics are rapid to compute for randomized low-rank approximation algorithms such as the randomized SVD and randomized Nystr\"om approximation, and they provide useful information that can be used to assess the quality of the computed output and guide algorithmic parameter choices. [ABSTRACT FROM AUTHOR]

A unified discrete gradient structure of the second order nonuniform integral averaged approximations for the Caputo fractional derivative and the Riemann--Liouville fractional integral is established in this paper. The required constraint of the step-size ratio is weaker than that found in the literature. With the proposed discrete gradient structure, the energy stability of the variable step Crank--Nicolson type numerical schemes is derived immediately, which is essential to the longtime simulations of the time fractional gradient flows and the nonlinear integro-differential models. The discrete energy dissipation laws fit seamlessly into their classical counterparts as the fractional indexes tend to one. In particular, we provide a framework for the stability analysis of variable step numerical schemes based on the scalar auxiliary variable type approaches. The time fractional Swift--Hohenberg model and the time fractional sine-Gordon model are taken as two examples to elucidate the theoretical results at great length. Extensive numerical experiments using the adaptive time-stepping strategy are provided to verify the theoretical results in the time multiscale simulations. [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]

In this paper, we present a second-order accurate and linear numerical scheme for the phase field crystal equation and prove its convergence in the discrete L\infty sense. The key ingredient of the error analysis is to justify the boundedness of the numerical solution, so that the nonlinear term, treated explicitly in the scheme, can be bounded appropriately. Benefiting from the existence of the sixth-order dissipation term in the model, we first estimate the discrete H2 norm of the numerical error. The error estimate in the supremum norm is then obtained by the Sobolev embedding, so that the uniform bound of the numerical solution is available. We also show the mass conservation and the energy stability in the discrete setting. The proposed scheme is linear with constant coefficients so that it can be solved efficiently via some fast algorithms. Numerical experiments are conducted to verify the theoretical results, and long-time simulations in two- and three-dimensional spaces demonstrate the satisfactory and high effectiveness of the proposed numerical scheme. [ABSTRACT FROM AUTHOR]

Balanced singular perturbation approximation (SPA) is a model order reduction method for linear time-invariant systems that guarantees asymptotic stability and for which there exists an a priori error bound. In that respect, it is similar to balanced truncation (BT). However, the reduced models obtained by SPA generally introduce better approximation in the lower frequency range and near steady-states, whereas BT is better suited for the higher frequency range. Even so, independently of the frequency range of interest, BT and its variants are more often applied in practice, since there exist more efficient algorithmic realizations thereof. In this paper, we aim at closing this practically relevant gap for SPA. We propose two novel and efficient algorithms that are adapted for different settings. First, we derive a low-rank implementation of SPA that is applicable in the large-scale setting. Second, a data-driven reinterpretation of the method is proposed that only requires input-output data and thus is realization-free. A main tool for our derivations is the reciprocal transformation, which induces a distinct view on implementing the method. While the reciprocal transformation and the characterization of SPA are not new, their significance for the practical algorithmic realization has been overlooked in the literature. Our proposed algorithms have well-established counterparts for BT and, as such, a comparable computational complexity. The numerical performance of the two novel implementations is tested for several numerical benchmarks, and comparisons to their counterparts for BT as well as to existing implementations of SPA are made. [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]

Deep learning (DL), in particular deep neural networks, by default is purely datadriven and in general does not require physics. This is the strength of DL but also one of its key limitations when applied to science and engineering problems in which underlying physical properties--such as stability, conservation, and positivity--and accuracy are required. DL methods in their original forms are often not capable of respecting the underlying mathematical models or achieving desired accuracy even in big-data regimes. On the other hand, many data-driven science and engineering problems, such as inverse problems, typically have limited experimental or observational data, and DL would overfit the data in this case. Leveraging information encoded in the underlying mathematical models, we argue, not only compensates for missing information in low data regimes but also provides opportunities to equip DL methods with the underlying physics, hence promoting better generalization. This paper develops a model-constrained DL approach and its variant TNet--a Tikhonov neural network--which are capable of learning not only information hidden in the training data but also in the underlying mathematical models to solve inverse problems governed by partial differential equations in low data regimes. We provide the constructions and some theoretical results for the proposed approaches for both linear and nonlinear inverse problems. Since TNet is designed to learn inverse solutions with Tikhonov regularization, it is interpretable: in fact it recovers Tikhonov solutions for linear cases while potentially approximating Tikhonov solutions for nonlinear inverse problems. We also prove that data randomization can enhance not only the smoothness of the networks but also their generalizations. Comprehensive numerical results confirm the theoretical findings and show that with even as little as 1 training data sample for one-dimensional (1D) deconvolution, 5 for an inverse 2D heat conductivity problem, 100 for inverse initial conditions for a time-dependent 2D Burgers's equation, and 50 for inverse initial conditions for 2D Navier--Stokes equations, TNet solutions can be as accurate as Tikhonov solutions while being several orders of magnitude faster. This is possible owing to the model-constrained term, replications, and randomization. [ABSTRACT FROM AUTHOR]

The accurate and efficient evaluation of Newtonian potentials over general twodimensional domains is important for the numerical solution of Poisson's equation and volume integral equations. In this paper, we present a simple and efficient high-order algorithm for computing the Newtonian potential over a planar domain discretized by an unstructured mesh. The algorithm is based on the use of Green's third identity for transforming the Newtonian potential into a collection of layer potentials over the boundaries of the mesh elements, which can be easily evaluated by the Helsing--Ojala method. One important component of our algorithm is the use of high-order (up to order 20) bivariate polynomial interpolation in the monomial basis, for which we provide extensive justification. The performance of our algorithm is illustrated through several numerical experiments. [ABSTRACT FROM AUTHOR]

This paper proposes a sketching strategy with spherical designs to equip the classical spherical basis function (SBF) approach for massive spherical data fitting. We conduct theoretical analysis and numerical verifications to demonstrate the feasibility of the proposed sketching strategy. From the theoretical side, we prove that sketching based on spherical designs can reduce the computational burden of the SBF approach without sacrificing its approximation capability. In particular, we provide upper and lower bounds for the proposed sketching strategy to fit noisy data on spheres. From the experimental side, we numerically illustrate the feasibility of the sketching strategy by showing its comparable fitting performance with the SBF approach. These interesting findings show that the proposed sketching strategy is effective in fitting massive and noisy data on spheres. [ABSTRACT FROM AUTHOR]

Parametric optimal control problems governed by partial differential equations (PDEs) are widely found in scientific and engineering applications. Traditional grid-based numerical methods for such problems generally require repeated solutions of PDEs with different parameter settings, which is computationally prohibitive, especially for problems with high-dimensional parameter spaces. Although recently proposed neural network methods make it possible to obtain the optimal solutions simultaneously for different parameters, challenges still remain when dealing with problems with complex constraints. In this paper, we propose the AONN, an adjoint-oriented neural network method, to overcome the limitations of existing approaches in solving parametric optimal control problems. In an AONN, the neural networks are served as parametric surrogate models for the control, adjoint, and state functions to get the optimal solutions all at once. In order to reduce the training difficulty and handle complex constraints, we introduce an iterative training framework inspired by the classical direct-adjoint looping method so that penalty terms arising from the Karush--Kuhn--Tucker system can be avoided. Once the training is done, parameter-specific optimal solutions can be quickly computed through the forward propagation of the neural networks, which may be further used for analyzing the parametric properties of the optimal solutions. The validity and efficiency of the AONN is demonstrated through a series of numerical experiments with problems involving various types of parameters. [ABSTRACT FROM AUTHOR]

Two-phase flow with moving contact lines (MCLs) is an unsolved problem in fluid dynamics. It is challenging to solve the problem numerically due to its intrinsic multiscale property that the microscopic slip must be taken into account in a macroscopic model. It is even more difficult when the solid substrate has microscopic inhomogeneity or roughness. In this paper, we propose a novel unified numerical framework for two-phase flows with MCLs. The framework cover some typical sharp-interface models for MCLs and can deal with the contact angle hysteresis (CAH) naturally. We prove that all the models, including the nonlinear Cox model and a CAH model, are thermodynamically consistent in the sense that an energy dissipation relation is satisfied. We further derive a new variational formula which leads to a stable and consistent numerical method independent of the choice of the slip length and the contact line frictions. This enables us to efficiently solve the macroscopic models for MCLs without resolving very small scale flow field in the vicinity of the contact line. We prove the well-posedness of the fully decoupled scheme which is based on a stabilized extended finite element discretization and a level-set representation for the free interface. Numerical examples are given to show the efficiency of the numerical framework. [ABSTRACT FROM AUTHOR]

In this paper, based on the physics-informed neural networks (PINNs) framework, a meshfree method using the deep neural network approach is developed for solving two kinds of two-phase interface problems governed by different dynamic partial differential equations on either side of the stationary interface with the jump and high-contrast coefficients. The first type of two-phase interface problem is the fluid-fluid (two-phase flow) interface problem modeled by Navier--Stokes equations with high-contrast physical parameters across the interface. The second one is the fluid-structure interaction problem modeled by Navier--Stokes equations on one side of the interface and the structural equation on the other side, where the fluid and the structure interact with each other via the kinematic and dynamic interface conditions across the interface. Following the PINNs framework, the DNN/meshfree method is respectively developed for two kinds of two-phase interface problems by approximating the solutions using different DNN's structures in different subdomains and reformulating the interface problems as least-squares minimization problems based on a spacetime sampling-point set (as the training dataset). Mathematically, the approximation error analyses are carried out for both interface problems, revealing an intrinsic strategy for efficiently sampling points to improve the accuracy. In addition, compared with traditional discretization approaches (e.g., finite element/volume/difference methods), the proposed DNN/meshfree method and its error analysis technique can be smoothly extended to many other dynamic interface problems with stationary interfaces. Numerical experiments illustrate the accuracy of the proposed method for the presented two-phase interface problems and validate theoretical results to some extent through two numerical examples. [ABSTRACT FROM AUTHOR]

Convection-diffusion equations arise in a variety of applications such as particle transport, electromagnetics, and magnetohydrodynamics. Simulation of the convection-dominated regime for these problems, even with high-fidelity techniques, is particularly challenging due to the presence of sharp boundary layers and shocks causing jumps and discontinuities in the solution, and numerical issues such as loss of the maximum principle in the discretization. These complications cause instabilities, admitting large oscillations in the numerical solution when using traditional methods. Drawing connections to the simplex-averaged finite-element method [S. Wu and J. Xu, SIAM J. Numer. Anal., 58 (2020), pp. 884--906], this paper develops a mimetic finite-difference (MFD) discretization using exponentially averaged coefficients to overcome instability of the numerical solution as the diffusion coefficient approaches zero. The finite-element framework allows for transparent analysis of the MFD, such as proving well-posedness and deriving error estimates. Numerical tests are presented confirming the stability of the method and verifying the error estimates. [ABSTRACT FROM AUTHOR]

It is challenging to numerically solve oscillatory Hamiltonian systems due to the stiffness of the problems and the requirement of highly stable and energy-preserving schemes. The previously constructed numerical schemes are generally fully implicit and thus result in a considerable computational cost in long-time integrations. In this paper, a family of explicit and energy-conserving schemes are presented for solving oscillatory Hamiltonian systems. These schemes are developed by using the idea of the construction of exponential Rosenbrock-type (ER) methods and relaxation techniques. The novel relaxation methods can be high-order accurate and have better long-time numerical behavior than the corresponding ER methods. Several numerical experiments on typical models are given to demonstrate the efficiency of the proposed methods. [ABSTRACT FROM AUTHOR]

In this paper we generalize the polynomial time integration framework to additively partitioned initial value problems. The framework we present is general and enables the construction of many new families of additive integrators with arbitrary order-of-accuracy and varying degree of implicitness. In this first work, we focus on a new class of implicit-explicit polynomial block methods that are based on fully implicit Runge--Kutta methods with Radau nodes and that possess high stage order. We show that the new fully implicit-explicit (FIMEX) integrators have improved stability compared to existing IMEX Runge--Kutta methods, while also being more computationally efficient due to recent developments in preconditioning techniques for solving the associated systems of nonlinear equations. For PDEs on periodic domains where the implicit component is trivial to invert, we will show how parallelization of the right-hand side evaluations can be exploited to obtain significant speedup compared to existing serial IMEX Runge--Kutta methods. For parallel (in space) finite element discretizations, the new methods can achieve orders of magnitude better accuracy than existing IMEX Runge--Kutta methods and/or achieve a given accuracy several times times faster in terms of computational runtime. [ABSTRACT FROM AUTHOR]

Gradient-enhanced Kriging (GE-Kriging) is a well-established surrogate modeling technique for approximating expensive computational models. However, it tends to get impractical for high-dimensional problems due to the size of the inherent correlation matrix and the associated high-dimensional hyperparameter tuning problem. To address these issues, a new method, called sliced GE-Kriging (SGE-Kriging), is developed in this paper for reducing both the size of the correlation matrix and the number of hyperparameters. We first split the training sample set into multiple slices, and invoke Bayes' theorem to approximate the full likelihood function via a sliced likelihood function, in which multiple small correlation matrices are utilized to describe the correlation of the sample set rather than one large one. Then, we replace the original high-dimensional hyperparameter tuning problem with a low-dimensional counterpart by learning the relationship between the hyperparameters and the derivative-based global sensitivity indices. The performance of SGE-kriging is finally validated by means of numerical experiments with several benchmarks and a high-dimensional aerodynamic modeling problem. The results show that the SGE-Kriging model features an accuracy and robustness that is comparable to the standard one but comes at much less training costs. The benefits are most evident for high-dimensional problems with tens of variables. [ABSTRACT FROM AUTHOR]

In this paper we present an algebraic dimension-oblivious two-level domain decomposition solver for discretizations of elliptic partial differential equations. The proposed parallel domain decomposition solver is based on a space-filling curve partitioning approach that is applicable to any discretization, i.e., it directly operates on the assembled matrix equations. Moreover, it allows for the effective use of arbitrary processor numbers independent of the dimension of the underlying partial differential equation while maintaining its optimal convergence behavior. This is the core property required to attain a combination technique solver with extreme scalability which can utilize exascale parallel systems efficiently. Moreover, this approach provides a basis for the development of a faulttolerant solver for the numerical treatment of high-dimensional Problems. To achieve the required data redundancy we are therefore concerned with large overlaps of our domain decomposition, which we construct in any dimension via space-filling curves. In this paper, we propose a space-filling curve based domain decompositoin solver and present its convergence properties and scaling behavior. The results of our numerical experiments clearly show that our approach provides optimal convergence and scaling behavior in an arbitrary dimension utilizing arbitrary processor numbers. [ABSTRACT FROM AUTHOR]

In this paper, uniformly unconditionally stable first and second order finite difference schemes are developed for kinetic transport equations in the diffusive scaling. We first derive an approximate evolution equation for the macroscopic density, from the formal solution of the distribution function, which is then discretized by following characteristics for the transport part with a backward finite difference semi-Lagrangian approach, while the diffusive part is discretized implicitly. After the macroscopic density is available, the distribution function can be efficiently solved even with a fully implicit time discretization since all discrete velocities are decoupled, resulting in a low-dimensional linear system from spatial discretizations at each discrete velocity. Both first and second order discretizations in space and in time are considered. The resulting schemes can be shown to be asymptotic preserving (AP) in the diffusive limit. Uniformly unconditional stabilities are verified from a Fourier analysis based on eigenvalues of corresponding amplification matrices. Numerical experiments, including high-dimensional problems, have demonstrated the corresponding orders of accuracy both in space and in time, uniform stability, AP property, and good performance of our proposed approach. [ABSTRACT FROM AUTHOR]

This paper is concerned with a blood flow problem coupled with a slow plaque growth at the artery wall. In the model, the micro (fast) system is the Navier-Stokes equation with a periodically applied force, and the macro (slow) system is a fractional reaction equation, which is used to describe the plaque growth with memory effect. We construct an auxiliary temporal periodic problem and an effective time-average equation to approximate the original problem and analyze the approximation error of the corresponding linearized PDE (Stokes) system, where the simple front-tracking technique is used to update the slow moving boundary. An effective multiscale method is then designed based on the approximate problem and the front-tracking framework. We also present a temporal finite difference scheme with a spatial continuous finite element method and analyze its temporal discrete error. Furthermore, a fast iterative procedure is designed to find the initial value of the temporal periodic problem, and its convergence is analyzed as well. Our designed front-tracking framework and the iterative procedure for solving the temporal periodic problem make it easy to implement the multiscale method on existing PDE solving software. The numerical method is implemented by a combination of the finite element platform COMSOL Multiphysics and the mainstream software MATLAB, which significantly reduce the programming effort and easily handle the fluid-structure interaction, especially moving boundaries with more complex geometries. We present two numerical examples of ODEs and a two-dimensional Navier-Stokes system to demonstrate the effectiveness of the multiscale method. Finally, we have a numerical experiment on the plaque growth problem and discuss the physical implication of the fractional order parameter. [ABSTRACT FROM AUTHOR]

In this paper we analyze a zeroth-order proximal stochastic gradient method suitable for the minimization of weakly convex stochastic optimization problems. We consider nonsmooth and nonlinear stochastic composite problems, for which (sub)gradient information might be unavailable. The proposed algorithm utilizes the well-known Gaussian smoothing technique, which yields unbiased zeroth-order gradient estimators of a related partially smooth surrogate problem (in which one of the two nonsmooth terms in the original problem's objective is replaced by a smooth approximation). This allows us to employ a standard proximal stochastic gradient scheme for the approximate solution of the surrogate problem, which is determined by a single smoothing parameter, and without the utilization of first-order information. We provide state-of-the-art convergence rates for the proposed zeroth-order method using minimal assumptions. The proposed scheme is numerically compared against alternative zeroth-order methods as well as a stochastic subgradient scheme on a standard phase retrieval problem. Further, we showcase the usefulness and effectiveness of our method in the unique setting of automated hyperparameter tuning. In particular, we focus on automatically tuning the parameters of optimization algorithms by minimizing a novel heuristic model. The proposed approach is tested on a proximal alternating direction method of multipliers for the solution of L1/L2-regularized PDE-constrained optimal control problems, with evident empirical success. [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]

In this paper we construct a quadrature scheme to numerically solve the nonlocal diffusion equation (Aα + bI)u = f with Aα the αth power of the regularly accretive operator A. Rigorous error analysis is carried out and sharp error bounds (up to some negligible constants) are obtained. The error estimates include a wide range of cases in which the regularity index and spectral angle of A, the smoothness of f, and the size of b and α are all involved. The quadrature scheme is exponentially convergent with respect to the step size and is root-exponentially convergent with respect to the number of solves. Some numerical tests are presented in the last section to verify the sharpness of our estimates. Furthermore, both the scheme and the error bounds can be utilized directly to solve and analyze time-dependent problems. [ABSTRACT FROM AUTHOR]