Search

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]

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]

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

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]

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]

Investigating blood flow in the cardiovascular system is crucial for assessing cardiovascular health. Computational approaches offer some noninvasive alternatives to measure blood flow dynamics. Numerical simulations based on traditional methods such as finite-element and other numerical discretizations have been extensively studied and have yielded excellent results. However, adapting these methods to real-life simulations remains a complex task. In this paper, we propose a method that offers flexibility and can efficiently handle real-life simulations. We suggest utilizing the physics-informed neural network to solve the Navier--Stokes equation in a deformable domain, specifically addressing the simulation of blood flow in elastic vessels. Our approach models blood flow using an incompressible, viscous Navier--Stokes equation in an arbitrary Lagrangian--Eulerian form. The mechanical model for the vessel wall structure is formulated by an equation of Newton's second law of momentum and linear elasticity to the force exerted by the fluid flow. Our method is a mesh-free approach that eliminates the need for discretization and meshing of the computational domain. This makes it highly efficient in solving simulations involving complex geometries. Additionally, with the availability of well-developed open-source machine learning framework packages and parallel modules, our method can easily be accelerated through GPU computing and parallel computing. To evaluate our approach, we conducted experiments on regular cylinder vessels as well as vessels with plaque on their walls. We compared our results to a solution calculated by finite element methods using a dense grid and small time steps, which we considered as the ground truth solution. We report the relative error and the time consumed to solve the problem, highlighting the advantages of our method. [ABSTRACT FROM AUTHOR]

In this paper we consider the filtering of a class of partially observed piecewise deterministic Markov processes. In particular, we assume that an ordinary differential equation (ODE) drives the deterministic element and can only be solved numerically via a time discretization. We develop, based upon the approach in Lemaire, Thieullen, and Thomas [Adv. Appl. Probab., 52 (2020), pp. 138--172], a new particle and multilevel particle filter (MLPF) in order to approximate the filter associated to the discretized ODE. We provide a bound on the mean square error associated to the MLPF which provides guidance on setting the simulation parameters of the algorithm and implies that significant computational gains can be obtained versus using a particle filter. Our theoretical claims are confirmed in several numerical examples. [ABSTRACT FROM AUTHOR]

Tikhonov regularization is a widely used technique in solving inverse problems that can enforce prior properties on the desired solution. In this paper, we propose a Krylov subspace based iterative method for solving linear inverse problems with general-form Tikhonov regularization term xT Mx, where M is a positive semidefinite matrix. An iterative process called the preconditioned Golub--Kahan bidiagonalization (pGKB) is designed, which implicitly utilizes a proper preconditioner to generate a series of solution subspaces with desirable properties encoded by the regularizer xT Mx. Based on the pGKB process, we propose an iterative regularization algorithm via projecting the original problem onto small dimensional solution subspaces. We analyze the regularization properties of this algorithm, including the incorporation of prior properties of the desired solution into the solution subspace and the semiconvergence behavior of the regularized solution. To overcome instabilities caused by semiconvergence, we further propose two pGKB based hybrid regularization algorithms. All the proposed algorithms are tested on both small-scale and large-scale linear inverse problems. Numerical results demonstrate that these iterative algorithms exhibit excellent performance, outperforming other state-of-the-art algorithms in some cases. [ABSTRACT FROM AUTHOR]

Diffusive transport in media with discontinuous properties is a challenging problem that arises in many applications. This paper focuses on one dimensional discontinuous media with generalized permeable boundary conditions at the discontinuity interface. It presents novel analytical expressions from the method of images to simulate diffusive processes, such as mass or thermal transport. The analytical expressions are used to formulate a generalization of the existing Skew Brownian Motion, HYMLA, and Uffink's method, here named as GSBM, GHYMLA, and GUM, respectively, to handle generic interface conditions. The algorithms rely upon the random walk method and are tested by simulating transport in a bimaterial and in a multilayered medium with piecewise constant properties. The results indicate that the GUM algorithm provides the best performance in terms of accuracy and computational cost. The methods proposed can be applied for simulation of a wide range of differential problems. [ABSTRACT FROM AUTHOR]

Small parameters in partial differential equations can give rise to solutions with sharp inner layers that evolve over time. However, the standard model reduction method becomes inefficient when applied to these problems due to the slow decaying Kolmogorov N-width of the solution manifold. To address this issue, a natural approach is to transform the equation in such a way that the transformed solution manifold exhibits a fast decaying Kolmogorov N-width. In this paper, we focus on the Allen--Cahn equation as a model problem. We employ asymptotic analysis to identify slow variables and perform a transformation of the partial differential equations accordingly. Subsequently, we apply the proper orthogonal decomposition method and a QR discrete empirical interpolation method (qDEIM) technique to the transformed equation with the slow variables. Numerical experiments demonstrate that the new model reduction method yields significantly improved results compared to direct model reduction applied to the original equation. Furthermore, this approach can be extended to other equations, such as the convection equation and the Burgers equation. [ABSTRACT FROM AUTHOR]

In this paper, we present a fast divergence-free spectral algorithm for the curl-curl problem. Divergence-free bases in two and three dimensions are constructed by using the generalized Jacobi polynomials. An accurate spectral method with exact preservation of the divergence-free constraint pointwisely is then proposed, and its corresponding error estimate is established. We then present a highly efficient solution algorithm based on a combination of the matrix-free preconditioned Krylov subspace iterative method and a fully diagonalizable auxiliary problem, which is derived from the spectral discretizations of generalized eigenvalue problems of Laplace and biharmonic operators. We rigorously prove that the dimensions of the invariant subspace of the preconditioned linear system resulting from the divergence-free spectral method with respect to the dominant eigenvalue 1 are (N - 3)² and 2(N - 3)³ for two- and three-dimensional problems with (N - 1)² and 2(N - 1)³ unknowns, respectively. Thus, the proposed method usually takes only several iterations to converge, and, astonishingly, as the problem size (polynomial order) increases, the number of iterations will decrease, even for highly indefinite system and oscillatory solutions. As a result, the computational cost of the solution algorithm is only a small multiple of N³ and N4 floating number operations for two- and three-dimensional problems, respectively. Plenty of numerical examples for solving the curlcurl problem with both constant and variable coefficients in two and three dimensions are presented to demonstrate the accuracy and efficiency of the proposed method. [ABSTRACT FROM AUTHOR]

In this paper, we propose a deep learning framework for solving high-dimensional partial integro-differential equations (PIDEs) based on the temporal difference learning. We introduce a set of Lévy processes and construct a corresponding reinforcement learning model. To simulate the entire process, we use deep neural networks to represent the solutions and nonlocal terms of the equations. Subsequently, we train the networks using the temporal difference error, the termination condition, and properties of the nonlocal terms as the loss function. The relative error of the method reaches O (10- 3) in 100-dimensional experiments and O (10- 4) in one-dimensional pure jump problems. Additionally, our method demonstrates the advantages of low computational cost and robustness, making it well-suited for addressing problems with different forms and intensities of jumps. [ABSTRACT FROM AUTHOR]

The multigrid V-cycle method is a popular method for solving systems of linear equations. It computes an approximate solution by using smoothing on fine levels and solving a system of linear equations on the coarsest level. Solving on the coarsest level depends on the size and difficulty of the problem. If the size permits, it is typical to use a direct method based on LU or Cholesky decomposition. In settings with large coarsest-level problems, approximate solvers such as iterative Krylov subspace methods, or direct methods based on low-rank approximation, are often used. The accuracy of the coarsest-level solver is typically determined based on the experience of the users with the concrete problems and methods. In this paper, we present an approach to analyzing the effects of approximate coarsest-level solves on the convergence of the V-cycle method for symmetric positive definite problems. Using these results, we derive coarsest-level stopping criterion through which we may control the difference between the approximation computed by a V-cycle method with approximate coarsest-level solver and the approximation which would be computed if the coarsestlevel problems were solved exactly. The coarsest-level stopping criterion may thus be set up such that the V-cycle method converges to a chosen finest-level accuracy in (nearly) the same number of V-cycle iterations as the V-cycle method with exact coarsest-level solver. We also utilize the theoretical results to discuss how the convergence of the V-cycle method may be affected by the choice of a tolerance in a coarsest-level stopping criterion based on the relative residual norm. [ABSTRACT FROM AUTHOR]

The integration of graph neural networks (GNNs) and neural ordinary and partial differential equations has been extensively studied in recent years. GNN architectures powered by neural differential equations allow us to reason about their behavior, and develop GNNs with desired properties such as controlled smoothing or energy conservation. In this paper we take inspiration from Turing instabilities in a reaction diffusion (RD) system of partial differential equations, and propose a novel family of GNNs based on neural RD systems, called RDGNN. We show that our RDGNN is powerful for the modeling of various data types, from homophilic, to heterophilic, and spatiotemporal datasets. We discuss the theoretical properties of our RDGNN, its implementation, and show that it improves or offers competitive performance to state-of-the-art methods. [ABSTRACT FROM AUTHOR]

We studied the least-squares ReLU neural network (LSNN) method for solving a linear advection-reaction equation with discontinuous solution in [Z. Cai et al., J. Comput. Phys., 443 (2021), 110514]. The method is based on a least-squares formulation and uses a new class of approximating functions: ReLU neural network (NN) functions. A critical and additional component of the LSNN method, differing from other NN-based methods, is the introduction of a properly designed and physics preserved discrete differential operator. In this paper, we study the LSNN method for problems with discontinuity interfaces. First, we show that ReLU NN functions with depth [log2(d+1)]+1 can approximate any d-dimensional step function on a discontinuity interface generated by a vector field as streamlines with any prescribed accuracy. By decomposing the solution into continuous and discontinuous parts, we prove theoretically that the discretization error of the LSNN method using ReLU NN functions with depth [log2(d+1)]+1 is mainly determined by the continuous part of the solution provided that the solution jump is constant. Numerical results for both two- and three-dimensional test problems with various discontinuity interfaces show that the LSNN method with enough layers is accurate and does not exhibit the common Gibbs phenomena along discontinuity interfaces. [ABSTRACT FROM AUTHOR]

The spectrum of a kernel matrix significantly depends on the parameter values of the kernel function used to define the kernel matrix. This makes it challenging to design a preconditioner for a regularized kernel matrix that is robust across different parameter values. This paper proposes the adaptive factorized Nyström (AFN) preconditioner. The preconditioner is designed for the case where the rank k of the Nyström approximation is large, i.e., for kernel function parameters that lead to kernel matrices with eigenvalues that decay slowly. AFN deliberately chooses a well-conditioned submatrix to solve with and corrects a Nyström approximation with a factorized sparse approximate matrix inverse. This makes AFN efficient for kernel matrices with large numerical ranks. AFN also adaptively chooses the size of this submatrix to balance accuracy and cost. [ABSTRACT FROM AUTHOR]

Quantum dynamics, typically expressed in the form of a time-dependent Schrödinger equation with a Hermitian Hamiltonian, is a natural application for quantum computing. However, when simulating quantum dynamics that involves the emission of electrons, it is necessary to use artificial boundary conditions (ABCs) to confine the computation within a fixed domain. The introduction of ABCs alters the Hamiltonian structure of the dynamics, and existing quantum algorithms cannot be directly applied since the evolution is no longer unitary. The current paper utilizes a recently introduced Schrödingerization method that converts non-Hermitian dynamics into a Schrödinger form for the artificial boundary problems [S. Jin, N. Liu, and Y. Yu, Quantum Simulation of Partial Differential Equations via Schrödingerisation, preprint, arXiv:2212.13969, 2022], [S. Jin, N. Liu, and Y. Yu, Phys. Rev. A, 108 (2023), 032603]. We implement this method for three types of ABCs, including the complex absorbing potential technique, perfectly matched layer methods, and Dirichlet-to-Neumann approach. We analyze the query complexity of these algorithms and perform numerical experiments to demonstrate the validity of this approach. This helps to bridge the gap between available quantum algorithms and computational models for quantum dynamics in unbounded domains. [ABSTRACT FROM AUTHOR]

The present study is an extension of the work done by Peddle, Haut, and Wingate [SIAM J. Sci. Comput., 41 (2019), pp. A3476--A3497] and Haut and Wingate [SIAM J. Sci. Comput., 36 (2014), pp. A693--A713], where a two-level Parareal method with mapping and averaging is examined. The method proposed in this paper is a multilevel Parareal method with arbitrarily many levels, which is not restricted to the two-level case. We give an asymptotic error estimate which reduces to the two-level estimate for the case when only two levels are considered. Introducing more than two levels has important consequences for the averaging procedure, as we choose separate averaging windows for each of the different levels, which is an additional new feature of the present study. The different averaging windows make the proposed method especially appropriate for nonlinear multiscale problems, because we can introduce a level for each intrinsic scale of the problem and adapt the averaging procedure such that we reproduce the behavior of the model on the particular scale resolved by the level. The method is applied to nonlinear differential equations. The nonlinearities can generate a range of frequencies in the problem. The computational cost of the new method is investigated and studied on several examples. [ABSTRACT FROM AUTHOR]

This paper describes a class of shape optimization problems for optical metamaterials comprised of periodic microscale inclusions composed of a dielectric, low-dimensional material suspended in a nonmagnetic bulk dielectric. The shape optimization approach is based on a homogenization theory for time-harmonic Maxwell's equations that describes effective material parameters for the propagation of electromagnetic waves through the metamaterial. The control parameter of the optimization is a deformation field representing the deviation of the microscale geometry from a reference configuration of the cell problem. This allows for describing the homogenized effective permittivity tensor as a function of the deformation field. We show that the underlying deformed cell problem is well-posed and regular. This, in turn, proves that the shape optimization problem is well-posed. In addition, a numerical scheme is formulated that utilizes an adjoint formulation with either gradient descent or BFGS as optimization algorithms. The developed algorithm is tested numerically on a number of prototypical shape optimization problems with a prescribed effective permittivity tensor as the target. [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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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

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

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]

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]

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]

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]

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

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]

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]

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]