Search

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]

Surrogate models can reduce computational costs for multivariable functions with an unknown internal structure (black boxes). In a discrete formulation, surrogate modeling is equivalent to restoring a multidimensional array (tensor) from a small part of its elements. The alternating least squares (ALS) algorithm in the tensor train (TT) format is a widely used approach to effectively solve this problem in the case of nonadaptive tensor recovery from a given training set (i.e., tensor completion problem). TT-ALS allows obtaining a low-parametric representation of the tensor, which is free from the curse of dimensionality and can be used for fast computation of the values at arbitrary tensor indices or efficient implementation of algebra operations with the black box (integration, etc.). However, to obtain high accuracy in the presence of restrictions on the size of the train data, a good choice of initial approximation is essential. In this work, we construct the ANOVA representation in the TT-format and use it as an initial approximation for the TT-ALS algorithm. The performed numerical computations for a number of multidimensional model problems, including the parametric partial differential equation, demonstrate a significant advantage of our approach for the commonly used random initial approximation. For all considered model problems we obtained an increase in accuracy by at least an order of magnitude with the same number of requests to the black box. The proposed approach is very general and can be applied in a wide class of real-world surrogate modeling and machine learning problems. [ABSTRACT FROM AUTHOR]

We analyze a new algorithm for the finite element approximation of a family of eigenvalue problems for the curl operator that includes, in particular, the approximation of the helicity of a bounded domain. It exploits a tree-cotree decomposition of the graph relating the degrees of freedom of the Lagrangian finite elements and those of the first family of N\'ed\'elec finite elements to reduce significantly the dimension of the algebraic eigenvalue problem to be solved. The algorithm is well adapted to domains of general topology. Numerical experiments, including a not simply connected domain with a not connected boundary, are presented in order to assess the performance and generality of the method. [ABSTRACT FROM AUTHOR]

The tensor-train (TT) format is a highly compact low-rank representation for highdimensional tensors. TT is particularly useful when representing approximations to the solutions of certain types of parametrized partial differential equations. For many of these problems, computing the solution explicitly would require an infeasible amount of memory and computational time. While the TT format makes these problems tractable, iterative techniques for solving the PDEs must be adapted to perform arithmetic while maintaining the implicit structure. The fundamental operation used to maintain feasible memory and computational time is called rounding, which truncates the internal ranks of a tensor already in TT format. We propose several randomized algorithms for this task that are generalizations of randomized low-rank matrix approximation algorithms and provide significant reduction in computation compared to deterministic TT-rounding algorithms. Randomization is particularly effective in the case of rounding a sum of TT-tensors (where we observe 20\times speedup), which is the bottleneck computation in the adaptation of GMRES to vectors in TT format. We present the randomized algorithms and compare their empirical accuracy and computational time with deterministic alternatives. [ABSTRACT FROM AUTHOR]

We present an extension of the AAA (adaptive Antoulas--Anderson) algorithm for periodic functions, called "AAAtrig."" The algorithm uses the key steps of AAA approximation by (i) representing the approximant in (trigonometric) barycentric form and (ii) selecting the support points greedily. Accordingly, AAAtrig inherits all the favorable characteristics of AAA and is thus extremely flexible and robust, being able to consider quite general sets of sample points in the complex plane. We consider a range of applications with particular emphasis on solving Laplace's equation in periodic domains and compressing periodic conformal maps. These results reproduce the tapered exponential clustering effect observed in other recent studies. The algorithm is implemented in Chebfun. [ABSTRACT FROM AUTHOR]

A selection of algorithms for the rational approximation of matrix-valued functions are discussed, including variants of the interpolatory adaptive Antoulas--Anderson (AAA) method, the rational Krylov fitting (RKFIT) method based on approximate least squares fitting, vector fitting, and a method based on low-rank approximation of a block Loewner matrix. A new method, called the block-AAA algorithm, based on a generalized barycentric formula with matrix-valued weights, is proposed. All algorithms are compared in terms of obtained approximation accuracy and runtime on a set of problems from model order reduction and nonlinear eigenvalue problems, including examples with noisy data. It is found that interpolation-based methods are typically cheaper to run, but they may suffer in the presence of noise for which approximation-based methods perform better. [ABSTRACT FROM AUTHOR]

We propose a novel algorithm for the approximation of surface quasi-geostrophic (SQG) flows modeled by a nonlinear partial differential equation coupling transport and fractional diffusion phenomena. The time discretization consists of an explicit strong-stability-preserving three-stage Runge--Kutta method while a flux-corrected-transport (FCT) method coupled with Dunford--Taylor representations of fractional operators is advocated for the space discretization. Standard continuous piecewise linear finite elements are employed, and the algorithm does not have restrictions on the mesh structure or on the computational domain. In the inviscid case, we show that the resulting scheme satisfies a discrete maximum principle property under a standard CFL condition and observe, in practice, its second order accuracy in space. The algorithm successfully approximates several benchmarks with sharp transitions and fine structures typical of SQG flows. In addition, theoretical Kolmogorov energy decay rates are observed on a freely decaying atmospheric turbulence simulation. [ABSTRACT FROM AUTHOR]

We design and implement an efficient parallel algorithm for finding a perfect matching in a weighted bipartite graph such that weights on the edges of the matching are large. This problem differs from the maximum weight matching problem, for which scalable approximation algorithms are known. It is primarily motivated by finding good pivots in scalable sparse direct solvers before factorization. Due to the lack of scalable alternatives, distributed solvers use sequential implementations of maximum weight perfect matching algorithms, such as those available in MC64. To overcome this limitation, we propose a fully parallel distributed memory algorithm that first generates a perfect matching and then iteratively improves the weight of the perfect matching by searching for weight-increasing cycles of length 4 in parallel. For most practical problems the weights of the perfect matchings generated by our algorithm are very close to the optimum. An efficient implementation of the algorithm scales up to 256 nodes (17,408 cores) on a Cray XC40 supercomputer and can solve instances that are too large to be handled by a single node using the sequential algorithm. [ABSTRACT FROM AUTHOR]

This paper presents a weakly intrusive strategy for computing a low-rank approximation of the solution of a system of nonlinear parameter-dependent equations. The proposed strategy relies on a Newton-like iterative solver which only requires evaluations of the residual of the parameter-dependent equation and of a preconditioner (such as the differential of the residual) for instances of the parameters independently. The algorithm provides an approximation of the set of solutions associated with a possibly large number of instances of the parameters, with a computational complexity which can be orders of magnitude lower than when using the same Newton-like solver for all instances of the parameters. The reduction of complexity requires efficient strategies for obtaining low-rank approximations of the residual, of the preconditioner, and of the increment at each iteration of the algorithm. For the approximation of the residual and the preconditioner, weakly intrusive variants of the empirical interpolation method are introduced, which require evaluations of entries of the residual and the preconditioner. Then, an approximation of the increment is obtained by using a greedy algorithm for low-rank approximation, and a low-rank approximation of the iterate is finally obtained by using a truncated singular value decomposition. When the preconditioner is the differential of the residual, the proposed algorithm is interpreted as an inexact Newton solver for which a detailed convergence analysis is provided. Numerical examples illustrate the efficiency of the method. [ABSTRACT FROM AUTHOR]

We consider the maximum vertex-weighted matching problem (MVM), in which nonnegative weights are assigned to the vertices of a graph, the weight of a matching is the sum of the weights of the matched vertices, and we are required to compute a matching of maximum weight. We describe an exact algorithm for MVM with O(| V | | E| ) time complexity, and then we design a 2/3-approximation algorithm for MVM on bipartite graphs by restricting the length of augmenting paths to at most three. The latter algorithm has time complexity O(| E| +| V | log | V | ). The approximation algorithm solves two MVM problems on bipartite graphs, each with weights only on one vertex part, and then finds a matching from these two matchings using the Mendelsohn--Dulmage theorem. The approximation ratio of the algorithm is obtained by considering failed vertices, i.e., vertices that the approximation algorithm fails to match but the exact algorithm does. We show that at every step of the algorithm there are two distinct heavier vertices that we can charge each failed vertex to. We have implemented the 2/3-approximation algorithm for MVM and compare it with four other algorithms: an exact maximum edge-weighted matching (MEM) algorithm, the exact MVM algorithm, a 1/2-approximation algorithm for MVM, and a scaling-based (1 - ϵ)-approximation algorithm for MEM. On a test set of nineteen problems with several millions of vertices, we show that the maximum time taken by the exact MEM algorithm is 15 hours, while it is 22 minutes for the exact MVM algorithm, and less than 5 seconds for the 2/3-approximation algorithm. The 2/3-approximation algorithm obtains more than 99.5% of the weight and cardinality of an MVM, whereas the scaling-based approximation algorithms yield lower weights and cardinalities while taking an order of magnitude more time than the former algorithm. We also show that MVM problems should not be first transformed to MEM problems and solved using exact algorithms for the latter, since this transformation can increase runtimes by several orders of magnitude [ABSTRACT FROM AUTHOR]

We present a numerical method for polynomial approximation of multivariate functions. The method utilizes Gauss quadrature in tensor product form, which is known to be inefficient in high dimensions. Here we demonstrate that by using a new randomized algorithm and taking advantage of the tensor structure of the grids, a highly efficient algorithm can be constructed. The new method does not require prior knowledge/storage of the entire data set at all the tensor grid points, whose total number of points is excessively large in high dimensions. Instead the method utilizes one data point at a time and iteratively conducts the approximation. This feature allows the use of the method irrespective of the size of the data set. We establish the rate of convergence of this iterative algorithm and show that its operational counts can be lower than the standard methods, when applicable, such as least squares. Numerical examples in up to hundreds of dimensions are presented to verify the theoretical analysis and demonstrate the effectiveness of the method. [ABSTRACT FROM AUTHOR]

In this paper we develop a two-grid convergence theory for the parallel-in-time scheme known as multigrid reduction in time (MGRIT), as it is implemented in the open-source package [XBraid: Parallel Multigrid in Time, http://llnl.gov/casc/xbraid]. MGRIT is a scalable and multilevel approach to parallel-in-time simulations that nonintrusively uses existing time-stepping schemes, and in a specific two-level setting it is equivalent to the widely known parareal algorithm. The goal of this paper is twofold. First, we present a two-level MGRIT convergence analysis for linear problems where the spatial discretization matrix can be diagonalized, and then apply this analysis to our two basic model problems, the heat equation and the advection equation. One important assumption is that the coarse and fine time-grid propagators can be diagaonalized by the same set of eigenvectors, which is often the case when the same spatial discretization operator is used on the coarse and fine time grids. In many cases, the MGRIT algorithm is guaranteed to converge, and we demonstrate numerically that the theoretically predicted convergence rates are sharp in practice for our model problems. Second, we explore how the convergence of MGRIT compares to the stability of the chosen time-stepping scheme. In particular, we demonstrate that a stable time-stepping scheme does not necessarily imply convergence of MGRIT, although MGRIT with FCF-relaxation always converges for the diffusion dominated problems considered here. [ABSTRACT FROM AUTHOR]

We propose a projection-based model order reduction method for the solution of parameter-dependent dynamical systems. The proposed method relies on the construction of timedependent reduced spaces generated from evaluations of the solution of the full-order model at some selected parameter values. The approximation obtained by Galerkin projection is the solution of a reduced dynamical system with a modified ux that takes into account the time dependency of the reduced spaces. An a posteriori error estimate is derived, and a greedy algorithm using this error estimate is proposed for the adaptive selection of parameter values. The resulting method can be interpreted as a dynamical low-rank approximation method with a subspace point of view and a uniform control of the error over the parameter set. [ABSTRACT FROM AUTHOR]

A crude but commonly used technique for compressing ordered scientific data consists of simply retaining every sth datum (with a value of s = 10 generally the default) and discarding the remainder. Should the value of a discarded datum be required afterwards, an approximation is generated by linear interpolation of the two nearest retained values. Despite the widespread use of this and similar techniques, there is little by way of theoretical analysis of their expected performance. First, we quantify the accuracy achieved by linear interpolation when approximating values discarded by decimation, obtaining both deterministic bounds in terms of appropriate smoothness measures of the data and probabilistic bounds in terms of statistics of the data. Second, we investigate the efficiency of the lossless compression scheme consisting of decimation coupled with encoding of the interpolation errors. In particular, we bound the expected compression ratio in terms of the appropriate measures of the data. Finally, we provide numerical illustrations of the practical performance of the algorithm on some real datasets. [ABSTRACT FROM AUTHOR]

Functions with jumps and kinks typically arising from parameter dependent or stochastic hyperbolic PDEs are notoriously difficult to approximate. If the jump location in physical space is parameter dependent or random, standard approximation techniques like reduced basis methods, PODs, polynomial chaos, etc., are known to yield poor convergence rates. In order to improve these rates, we propose a new approximation scheme. It relies on snapshots for the reconstruction of parameter dependent functions so that it is efficiently applicable in a PDE context. However, the physical domain of each individual snapshot is transformed before its use in the recon- struction, which allows realignng the moving discontinuities and yields high convergence rates. The transforms are automatically computed by minimizing a training error. In order to show feasibility of this approach, it is tested by numerical experiments with one and two physical dimensions and one parameter dimension. [ABSTRACT FROM AUTHOR]

We propose and compare goal-oriented projection-based model order reduction meth- ods for the estimation of vector-valued functionals of the solution of parameter-dependent equations. The first projection method is a generalization of the classical primal-dual method to the case of vector-valued variables of interest. We highlight the role played by three reduced spaces: the ap- proximation space and the test space associated to the primal variable, and the approximation space associated with the dual variable. Then we propose a Petrov-Galerkin projection method based on a saddle point problem involving an approximation space for the primal variable and an approximation space for an auxiliary variable. A goal-oriented choice of the latter space, defined as the sum of two spaces, allows us to improve the approximation of the variable of interest compared to a primal-dual method using the same reduced spaces. Then, for both approaches, we derive computable error estimates for the approximations of the variable of interest and we propose greedy algorithms for the goal-oriented construction of reduced spaces. The performance of the algorithms is illustrated on numerical examples and compared to standard (non-goal-oriented) algorithms. [ABSTRACT FROM AUTHOR]

Proper generalized decompositions (PGDs) are a family of methods for efficiently solving high-dimensional PDEs, which seek to find a low-rank approximation to the solution of the PDE a priori. Convergence of PGD algorithms can only be proven for problems which are continuous, symmetric, and strongly coercive. In the particular case of problems which are only weakly coercive we have the additional issue that weak coercivity estimates are not guaranteed to be inherited by the low-rank PGD approximation. This can cause stability issues when employing a Galerkin PGD approximation of weakly coercive problems. In this paper we propose the use of PGD algorithms based on least-squares formulations which always lead to symmetric and strongly coercive problems and hence provide stable and provably convergent algorithms. Taking the Stokes problem as a prototypical example of a weakly coercive problem, we develop and compare rigorous leastsquares PGD algorithms based on continuous least-squares estimates for two different reformulations of the problem. We show that these least-squares PGDs provide a much stabler algorithm than an equivalent Galerkin PGDs, and we provide proofs of convergence of the algorithms. [ABSTRACT FROM AUTHOR]

We construct and analyze acceleration techniques for adaptive Monte Carlo simulations for general multivariate probability laws when the sample average approximation is employed for optimal parameter search. Our goal is to accelerate the adaptive Monte Carlo estimation by lead- ing the parameter search line based on the sample average approximation to reach a nearly optimal realm at its small sample size stage. First, we introduce an auxiliary parameter into the parameter search line and update it wisely to aim for small-sample convergence. All three lines of the algorithm (the Monte Carlo averaging, the importance sampling parameter search, and the auxiliary parame- ter updating) run concurrently in a fully automated manner with a common set of random vectors without ad hoc tuning of the simulation system. We next propose and examine various criteria for the auxiliary parameter updating, all under which the asymptotic normality of the estimator of the desired mean and of the importance sampling parameter hold true. To illustrate the applicability and effectiveness of the proposed methods, we provide numerical results throughout the paper with particular attention toward the very early stages of the algorithm. In particular, direct use of the importance sampling parameter as the auxiliary parameter has the great potential to accelerate the parameter search line, and thus the Monte Carlo averaging as a consequence, by orders of magnitude without further restriction of the parameter domain or additional computing effort. [ABSTRACT FROM AUTHOR]

We describe a half-approximation algorithm, b-Suitor, for computing a b-Matching of maximum weight in a graph with weights on the edges. b-Matching is a generalization of the well-known Matching problem in graphs, where the objective is to choose a subset of M edges in the graph such that at most a specified number b(v) of edges in M are incident on each vertex v. Subject to this restriction we maximize the sum of the weights of the edges in M. We prove that the b-Suitor algorithm computes the same b-Matching as the one obtained by the Greedy algorithm for the problem. We implement the algorithm on serial and shared-memory parallel processors and compare its performance against a collection of approximation algorithms that have been proposed earlier. Our results show that the b-Suitor algorithm outperforms the Greedy and locally dominant edge algorithms by one to two orders of magnitude on a serial processor. The b-Suitor algorithm has a high degree of concurrency, and it scales well up to 240 threads on a shared-memory multiprocessor. The b-Suitor algorithm outperforms the locally dominant edge algorithm by a factor of 14 on 16 cores of an Intel Xeon multiprocessor. [ABSTRACT FROM AUTHOR]

Under some weak conditions, the first-passage time of the Brownian motion to a continuous curved boundary is an almost surely finite stopping time. Its probability density function (pdf) is explicitly known only in few particular cases. Several mathematical studies proposed to approximate the pdf in a quite general framework or even to simulate this hitting time using a discrete time approximation of the Brownian motion. The authors study a new algorithm which permits one to simulate the first-passage time using an iterating procedure. The convergence rate presented in this paper suggests that the method is very efficient. [ABSTRACT FROM AUTHOR]

This paper compares the performance of penalty and Lagrange multiplier approaches for the necessary unit-length constraint in the computation of liquid crystal equilibrium configurations. Building on previous work in [SIAM J. Sci. Comput., 37 (2015), pp. S157-S176; SIAM J. Numer. Anal., 53 (2015), pp. 2226-2254], the penalty method is derived and well-posedness of the linearizations within the nonlinear iteration is discussed. In addition, the paper considers the effects of tailored trust-region methods in the context of finite-element discretizations and nested iteration for both formulations. Such methods are aimed at increasing the efficiency and robustness of each algorithm's nonlinear iterations. Three representative elastic equilibrium problems are considered to examine each method's performance. The first two configurations have analytical expressions for their exact solutions and, therefore, convergence to the true solution is considered. The third problem considers complicated boundary conditions, relevant in ongoing research, simulating surface nano-patterning. Finally, a novel multigrid scheme is introduced and tested for electrically and flexoelectrically coupled models to establish scalability for highly complicated applications. The Lagrange multiplier method is found to outperform the penalty method in a number of measures, the developed trust regions are shown to improve robustness, and nested iteration proves highly effective at reducing computational costs. [ABSTRACT FROM AUTHOR]

We address an important research area in stochastic multiscale modeling, namely, the propagation of uncertainty across heterogeneous domains characterized by partially correlated processes with vastly different correlation lengths. This class of problems arises very often when computing stochastic PDEs and particle models with stochastic/stochastic domain interaction but also with stochastic/deterministic coupling. The domains may be fully embedded, adjacent, or partially overlapping. The fundamental open question we address is the construction of proper transmission boundary conditions that preserve global statistical properties of the solution across different subdomains. Often, the codes that model different parts of the domains are black box and hence a domain decomposition technique is required. No rigorous theory or even effective empirical algorithms have yet been developed for this purpose, although interfaces defined in terms of functionals of random fields (e.g., multipoint cumulants) can overcome the computationally prohibitive problem of preserving sample-path continuity across domains. The key idea of the different methods we propose relies on combining local reduced-order representations of random fields with multilevel domain decomposition. Specifically, we propose two new algorithms: The first one enforces the continuity of the conditional mean and variance of the solution across adjacent subdomains by using Schwarz iterations. The second algorithm is based on PDE-constrained multiobjective optimization, and it allows us to set more general interface conditions. The effectiveness of these new algorithms is demonstrated in numerical examples involving elliptic problems with random diffusion coefficients, stochastically advected scalar fields, and nonlinear advection-reaction problems with random reaction rates. [ABSTRACT FROM AUTHOR]

Reduced order models, in particular the reduced basis method, rely on empirically built and problem dependent basis functions that are constructed during an off-line stage. In the online stage, the precomputed problem-dependent solution space, that is spanned by the basis functions, can then be used in order to reduce the size of the computational problem. For complex problems, the number of basis functions required to guarantee a certain error tolerance can become too large in order to benefit computationally from the model reduction. To overcome this, the present work introduces a framework where local approximation spaces (in parameter space) are used to define the reduced order approximation in order to have explicit control over the on-line cost. This approach also adapts the local approximation spaces to local anisotropic behavior in the parameter space. We present the algorithm and numerous numerical tests. [ABSTRACT FROM AUTHOR]

New algorithms are proposed for the Tucker approximation of a 3-tensor accessed only through a tensor-by-vector-by-vector multiplication subroutine. In the matrix case, the Krylov methods are methods of choice to approximate the dominant column and row subspaces of a sparse or structured matrix given through a matrix-by-vector operation. Using the Wedderburn rank reduction formula, we propose a matrix approximation algorithm that computes the Krylov subspaces and can be generalized to 3-tensors. The numerical experiments show that on quantum chemistry data the proposed tensor methods outperform the minimal Krylov recursion of Savas and Eldén. [ABSTRACT FROM AUTHOR]