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.

In this paper, we introduce the ensemble Kaiman--Stein gradient descent (EnKSGD) class of algorithms. The EnKSGD class of algorithms builds on the ensemble Kalman filter (EnKF) line of work, applying techniques from sequential data assimilation to unconstrained optimization and parameter estimation problems. An essential idea is to exploit the EnKF as a black box (i.e., derivative-free, zeroth order) optimization tool if iterated to convergence. In this paper, we return to the foundations of the EnKF as a sequential data assimilation technique, including its continuous-time and mean-field limits, with the goal of developing faster optimization algorithms suited to noisy black box optimization and inverse problems. The resulting EnKSGD class of algorithms can be designed to both maintain the desirable property of affine-invariance and employ the well-known backtracking line search. Furthermore, EnKSGD algorithms are designed to not necessitate the subspace restriction property and to avoid the variance collapse property of previous iterated EnKF approaches to optimization, as both these properties can be undesirable in an optimization context. EnKSGD also generalizes beyond the L² loss and is thus applicable to a wider class of problems than the standard EnKF. Numerical experiments with empirical risk minimization type problems, including both linear and nonlinear least squares problems, as well as maximum likelihood estimation, demonstrate the faster empirical convergence of EnKSGD relative to alternative EnKF approaches to optimization. [ABSTRACT FROM AUTHOR]

In this paper, we address the numerical solution of the quadratic optimal transport problem in its dynamical form, the so-called Benamou-Brenier formulation. When solved using interior point methods, the main computational bottleneck is the solution of large saddle point linear systems arising from the associated Newton-Raphson scheme. The main purpose of this paper is to design efficient preconditioners to solve these linear systems via iterative methods. Among the proposed preconditioners, we introduce one based on the partial commutation of the operators that compose the dual Schur complement of these saddle point linear systems, which we refer to as the BB-preconditioner. A series of numerical tests show that the BB-preconditioner is the most efficient among those presented, despite a performance deterioration in the last steps of the interior point method. It is in fact the only one having a CPU time that scales only slightly worse than linearly with respect to the number of unknowns used to discretize the problem. [ABSTRACT FROM AUTHOR]

For a nonlinear dynamical system that depends on parameters, this paper introduces a novel tensorial reduced-order model (TROM). The reduced model is projection-based, and for systems with no parameters involved, it resembles proper orthogonal decomposition (POD) combined with the discrete empirical interpolation method (DEIM). For parametric systems, TROM employs low-rank tensor approximations in place of truncated SVD, a key dimension-reduction technique in POD with DEIM. Three popular low-rank tensor compression formats are considered for this purpose: canonical polyadic, Tucker, and tensor train. The use of multilinear algebra tools allows the incorporation of information about the parameter dependence of the system into the reduced model and leads to a POD-DEIM type ROM that (i) is parameter-specific (localized) and predicts the system dynamics for out-of-training set (unseen) parameter values, (ii) mitigates the adverse effects of high parameter space dimension, (iii) has online computational costs that depend only on tensor compression ranks but not on the full-order model size, and (iv) achieves lower reduced space dimensions compared to the conventional POD-DEIM ROM. This paper explains the method, analyzes its prediction power, and assesses its performance for two specific parameter-dependent nonlinear dynamical systems. [ABSTRACT FROM AUTHOR]

The article introduces a new section structure for the SIAM Journal on Scientific Computing (SISC), including a focus on Scientific Machine Learning as a major development in the field of scientific computing. Topics include the development of accurate, efficient, and scalable machine learning algorithms for solving mathematical modeling problems, emphasizing the importance of reproducibility in the review process.

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]

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]

The classic exact simulation scheme for the Ornstein--Uhlenbeck driven stochastic volatility model is designed for the single volatility factor case. Extension to the multifactor case results in a cumbersome procedure requiring multiple numerical inversions of Laplace transforms and subsequent random sampling through numerical methods, resulting in it being perceptively slow to run. Moreover, for each volatility factor, the error is controlled by two parameters, ensuring difficult control of the bias. In this paper, we propose a new exact simulation scheme for the multifactor Ornstein--Uhlenbeck driven stochastic volatility model that is easier to implement, faster to run, and allows for an improved control of the error, which, in contrast to the existing method, is controlled by only one parameter, regardless of the number of volatility factors. Numerical results show that the proposed approach is three times faster than the original approach when one volatility factor is considered and 11 times faster when four volatility factors are considered, while still being theoretically exact. [ABSTRACT FROM AUTHOR]

This paper proposes an adaptive space-time algorithm based on domain decomposition for the large-scale simulation of a recently developed thermodynamically consistent reservoir problem. In the approach, the bound constraints are represented by means of a minimum-type complementarity function to enforce the positivity of the reservoir model, and a space-time mixed finite element method is applied for the parallel-in-time monolithic discretization. In particular, we propose a time-adaptive strategy using the improved backward differencing formula of second order, to take full advantage of the high degree of space-time parallelism. Moreover, the complicated dynamics with higher nonlinearity of space-time discretization require some innovative nonlinear and linear solution strategies. Therefore, we present a class of modified semismooth Newton algorithms to enhance the convergence rate of nonlinear iterations. Multilevel space-time restricted additive Schwarz algorithms, whose subdomains cover both space and time variables, are also studied for domain decomposition-based preconditioning. Numerical experiments demonstrate the robustness and parallel scalability of the proposed adaptive space-time algorithm on a supercomputer with tens of thousands of processor cores. [ABSTRACT FROM AUTHOR]

In this paper we present a framework for the construction and implementation of general virtual element spaces based on projections built from constrained least squares problems. Building on the triples used for finite element spaces, we introduce the concept of a virtual element method (VEM) tuple which encodes the necessary building blocks to construct these projections. Using this approach, a wide range of virtual element spaces can be defined. We discuss Hk-conforming spaces for k = 1, 2 as well as divergence and curl free spaces. This general framework has the advantage of being easily integrated into any existing finite element package, and we demonstrate this within the open source software package Dune. [ABSTRACT FROM AUTHOR]

In this paper, we propose a cubic-regularized Riemannian optimization method (RDRSOM), which partially exploits the second-order information and achieves the iteration complexity of O(1/ε³/²). In order to reduce the periteration computational cost, we further propose a practical version of RDRSOM which is an extension of the well-known Barzilai--Borwein method, which enjoys the worst-case iteration complexity of O(1/ε²). Moreover, under more stringent conditions, RDRSOM achieves the iteration complexity of O(1/ε1.5). We apply our method to solve a nonlinear formulation of the wireless sensor network localization problem whose feasible set is a Riemannian manifold that has not been considered in the literature before. Numerical experiments are conducted to verify the high efficiency of our algorithm compared to state-of-the-art Riemannian optimization methods and other nonlinear solvers. [ABSTRACT FROM AUTHOR]

In this paper, we introduce a new algorithm for rare event estimation based on adaptive importance sampling. We consider a smoothed version of the optimal importance sampling density, which is approximated by an ensemble of interacting particles. The particle dynamics is governed by a McKean--Vlasov stochastic differential equation, which was introduced and analyzed in [Carrillo et al., Stud. Appl. Math., 148 (2022), pp. 1069-1140] for consensus-based sampling and optimization of posterior distributions arising in the context of Bayesian inverse problems. We develop automatic updates for the internal parameters of our algorithm. This includes a novel time step size controller for the exponential Euler method, which discretizes the particle dynamics. The behavior of all parameter updates depends on easy to interpret accuracy criteria specified by the user. We show in numerical experiments that our method is competitive to state-of-the-art adaptive importance sampling algorithms for rare event estimation, namely a sequential importance sampling method and the ensemble Kalman filter for rare event estimation. [ABSTRACT FROM AUTHOR]

Developing large-scale distributed methods that are robust to the presence of adversarial or corrupted workers is an important part of making such methods practical for real-world problems. In this paper, we propose an iterative approach that is adversary-tolerant for convex optimization problems. By leveraging simple statistics, our method ensures convergence and is capable of adapting to adversarial distributions. Additionally, the efficiency of the proposed methods for solving convex problems is shown in simulations with the presence of adversaries. Through simulations, we demonstrate the efficiency of our approach in the presence of adversaries and its ability to identify adversarial workers with high accuracy and tolerate varying levels of adversary rates. [ABSTRACT FROM AUTHOR]

Computing equilibrium shapes of crystals (ESCs) is a challenging problem in materials science that involves minimizing an orientation-dependent (i.e., anisotropic) surface energy functional subject to a prescribed mass constraint. The highly nonlinear and singular anisotropic terms in the problem make it very challenging from both analytical and numerical aspects. Especially when the strength of anisotropy is very strong (i.e., strongly anisotropic cases), the ESCs will form some singular, sharp corners even if the surface energy function is smooth. Traditional numerical approaches, such as the H-1 gradient flow, are unable to produce true sharp corners due to the necessary addition of a high-order regularization term that penalizes sharp corners and rounds them off. In this paper, we propose a new numerical method based on the Davis--Yin splitting (DYS) optimization algorithm to predict the ESCs instead of using gradient flow approaches. We discretize the infinite-dimensional phase-field energy functional in the absence of regularization terms and transform it into a finite-dimensional constraint minimization problem. The resulting optimization problem is solved using the DYS method, which automatically guarantees the mass-conservation and bound-preserving properties. We also prove the global convergence of the proposed algorithm. These desired properties are numerically observed. In particular, the proposed method can produce real sharp corners with satisfactory accuracy. Finally, we present numerous numerical results to demonstrate that the ESCs can be well simulated under different types of anisotropic surface energies, which also confirms the effectiveness and efficiency of the proposed method. [ABSTRACT FROM AUTHOR]

We present a novel solver technique for the anisotropic heat flux equation, aimed at the high level of anisotropy seen in magnetic confinement fusion plasmas. Such problems pose two major challenges: (i) discretization accuracy and (ii) efficient implicit linear solvers. We simultaneously address each of these challenges by constructing a new finite element discretization with excellent accuracy properties, tailored to a novel solver approach based on algebraic multigrid (AMG) methods designed for advective operators. We pose the problem in a mixed formulation, introducing the directional temperature gradient as an auxiliary variable. The temperature and auxiliary fields are discretized in a scalar discontinuous Galerkin space with upwinding principles used for discretizations of advection. We demonstrate the proposed discretization's superior accuracy over other discretizations of anisotropic heat flux, achieving error 1000 x smaller for anisotropy ratio of 109, for closed field lines. The block matrix system is reordered and solved in an approach where the two advection operators are inverted using AMG solvers based on approximate ideal restriction, which is particularly efficient for upwind discontinuous Galerkin discretizations of advection. To ensure that the advection operators are nonsingular, in this paper we restrict ourselves to considering open (acyclic) magnetic field lines for the linear solvers. We demonstrate fast convergence of the proposed iterative solver in highly anisotropic regimes where other diffusion-based AMG methods fail. [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]

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]

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 rounding (SR) offers an alternative to the deterministic IEEE-754 floating-point rounding modes. In some applications such as PDEs, ODEs, and neural networks, SR empirically improves the numerical behavior and convergence to accurate solutions while the theoretical background remains partial. Recent works by Ipsen, Zhou, Higham, and Mary have computed SR probabilistic error bounds for basic linear algebra kernels. For example, the inner product SR probabilistic bound of the forward error is proportional to √nu instead of nu for the default rounding mode. To compute the bounds, these works show that the errors accumulated in computation form a martingale. This paper proposes an alternative framework to characterize SR errors based on the computation of the variance. We pinpoint common error patterns in numerical algorithms and propose a lemma that bounds their variance. For each probability and through the Bienaymé-Chebyshev inequality, this bound leads to a better probabilistic error bound in several situations. Our method has the advantage of providing a tight probabilistic bound for all algorithms fitting our model. We show how the method can be applied to give SR error bounds for the inner product and Horner polynomial evaluation. [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]

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 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 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]

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 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 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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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 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 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]

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]

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]