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]

The dynamic formulation of optimal transport has attracted growing interest in scientific computing and machine learning, and its computation requires solving a PDE-constrained optimization problem. The classical Eulerian discretization based approaches suffer from the curse of dimensionality, which arises from the approximation of a high-dimensional velocity field. In this work, we propose a deep learning based method for solving the dynamic optimal transport in highdimensional space. Our method contains three main ingredients: a carefully designed representation of the velocity field, the discretization of the PDE constraint along the characteristics, and the computation of a high-dimensional integral by the Monte Carlo method in each time step. Specifically, in the representation of the velocity field, we apply the classical nodal basis function in time and the deep neural networks in the space domain with the H1-norm regularization. This technique promotes the regularity of the velocity field in both time and space such that the discretization along the characteristic remains stable during the training process. Extensive numerical examples have been conducted to test the proposed method. Compared to other solvers of optimal transport, our method could give more accurate results in high-dimensional cases and has very good scalability with respect to dimension. Finally, we extend our method to more complicated cases such as the crowd motion problem. [ABSTRACT FROM AUTHOR]

The computation of f(A)b, the action of a matrix function on a vector, is a task arising in many areas of scientific computing. In many applications, the matrix A is sparse but so large that only a rather small number of Krylov basis vectors can be stored. Here we discuss a new approach to overcome this limitation by randomized sketching combined with an integral representation of f(A)b. Two different approximation methods are introduced, one based on sketched FOM and another based on sketched GMRES. The convergence of the latter method is analyzed for Stieltjes functions of positive real matrices. We also derive a closed-form expression for the sketched FOM approximant and bound its distance to the full FOM approximant. Numerical experiments demonstrate the potential of the presented sketching approaches. [ABSTRACT FROM AUTHOR]

We present the derivation, implementation, and analysis of a multiresolution adaptive grid framework for numerical simulations on octree-based three-dimensional block-structured collocated grids with distributed computational architectures. Our approach provides a consistent handling of nonlifted and lifted interpolating wavelets of arbitrary order demonstrated using second-, fourth-, and sixth-order wavelets, combined with standard finite-difference-based discretization op-erators. We first validate that the wavelet family used provides strict and explicit error control when coarsening the grid, and show that lifting wavelets increase the grid compression rate while conserving discrete moments across levels. Further, we demonstrate that high-order PDE discretization schemes combined with sufficiently high-order wavelets retain the expected convergence order even at resolution jumps. We then simulate the advection of a scalar to analyze convergence for the temporal evolution of a PDE. The results shows that our wavelet-based refinement criterion is successful at controlling the overall error while the coarsening criterion is effective at retaining the relevant information on a compressed grid. Our software exploits a block-structured grid data structure for efficient multilevel operations, combined with a parallelization strategy that relies on a one-sided MPI-RMA communication approach with active post-start-complete-wait synchronization. Using performance tests up to 16,384 cores, we demonstrate that this leads to a highly scalable performance. The associated code is available under a BSD-3 license at https://github.com/vanreeslab/murphy. [ABSTRACT FROM AUTHOR]

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

The article offers information on The Seventeenth Copper Mountain Conference on Iterative Methods (CMCIM), which took place virtually in April 2022.

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

We propose a supervised deep learning algorithm based on low-discrepancy sequences as the training set. By a combination of theoretical arguments and extensive numerical experiments we demonstrate that the proposed algorithm significantly outperforms standard deep learning algorithms that are based on randomly chosen training data for problems in moderately high dimensions. The proposed algorithm provides an efficient method for building inexpensive surrogates for many underlying maps in the context of scientific computing. [ABSTRACT FROM AUTHOR]

In this work, we present an adaptive unfitted finite element scheme that combines the aggregated finite element method with parallel adaptive mesh refinement. We introduce a novel scalable distributed-memory implementation of the resulting scheme on locally adapted Cartesian forest-of-trees meshes. We propose a two-step algorithm to construct the finite element space at hand by means of a discrete extension operator that carefully mixes aggregation constraints of problematic degrees of freedom, which get rid of the small cut cell problem, and standard hanging degree of freedom constraints, which ensure trace continuity on nonconforming meshes. Following this approach, we derive a finite element space that can be expressed as the original one plus well-defined linear constraints. Moreover, it requires minimum parallelization effort, using standard functionality available in existing large-scale finite element codes. Numerical experiments demonstrate its optimal mesh adaptation capability, robustness to cut location, and parallel efficiency, on classical Poisson hp-adaptivity benchmarks. Our work opens the path to functional and geometrical error-driven dynamic mesh adaptation with the aggregated finite element method in large-scale realistic scenarios. Likewise, it can offer guidance for bridging other scalable unfitted methods and parallel adaptive mesh refinement. [ABSTRACT FROM AUTHOR]

Randomness extractors and error correcting codes are fundamental objects in computer science. Recently, there have been several natural generalizations of these objects, in the context and study of tamper-resilient cryptography. These are seeded nonmalleable extractors, introduced by Dodis and Wichs (STOC 2009); seedless nonmalleable extractors, introduced by Cheraghchi and Guruswami (TCC 2014); and nonmalleable codes, introduced by Dziembowski, Pietrzak, and Wichs (J. ACM, 2018). Besides being interesting on their own, they also have important applications in cryptography, e.g., privacy amplification with an active adversary, explicit nonmalleable codes, etc., and often have unexpected connections to their nontampered analogues. However, the known constructions are far behind their nontampered counterparts. Indeed, the best known seeded nonmalleable extractor requires min-entropy rate at least 0.49 [X. Li, in Proceedings of the 53rd Annual IEEE Symposium on Foundations of Computer Science, 2012, pp. 688--697], while explicit construction of nonmalleable two-source extractors was not known even if both sources have full min-entropy and was left as an open problem in [M. Cheraghchi and V. Guruswami, J. Cryptology, 30 (2017), pp. 191--241]. In this paper we make progress towards solving the above problems and other related generalizations. Our contributions are as follows: (i) We construct an explicit seeded nonmalleable extractor for min-entropy k \geq log2 n. This dramatically improves all previous results and gives a simpler two-round privacy amplification protocol with optimal entropy loss, matching the best known result in [X. Li, in Theory of Cryptography (TCC 2015), Springer, 2015, pp. 502--531]. In fact, we construct more general seeded nonmalleable extractors (that can handle multiple adversaries) which were used in the recent construction of explicit two-source extractors for polylogarithmic min-entropy [E. Chattopadhyay and D. Zuckerman, Ann. of Math. (2), 189 (2019), pp. 653--705]. (ii) We construct the first explicit nonmalleable two-source extractor for min-entropy k \geq n n\Omega (1), with output size n\Omega (1) and error 2 n\Omega (1), thus resolving the open question in [M. Cheraghchi and V. Guruswami, J. Cryptology, 30 (2017), pp. 191--241]. (iii) We motivate and initiate the study of two natural generalizations of seedless nonmalleable extractors and nonmalleable codes, where the sources or the codeword may be tampered many times. For this, we construct the first explicit nonmalleable two-source extractor with tampering degree t up to n\Omega (1). By using the connection in [M. Cheraghchi and V. Guruswami, J. Cryptology, 30 (2017), pp. 191--241] and providing efficient sampling algorithms, we obtain the first explicit nonmalleable codes with tampering degree t up to n\Omega (1). We call these stronger notions one-many and many-many nonmalleable codes. This provides a stronger information theoretic analogue of a primitive known as continuous nonmalleable codes. Our basic technique used in all of our constructions can be seen as inspired, in part, by the techniques previously used to construct cryptographic nonmalleable commitments. [ABSTRACT FROM AUTHOR]

Computing units that carry out a fused multiply-add (FMA) operation with matrix arguments, referred to as tensor units by some vendors, have great potential for use in scientific computing. However, these units are inherently mixed precision, and existing rounding error analyses do not support them. We consider a mixed precision block FMA that generalizes both the usual scalar FMA and existing tensor units. We describe how to exploit such a block FMA in the numerical linear algebra kernels of matrix multiplication and LU factorization and give detailed rounding error analyses of both kernels. An important application is to GMRES-based iterative refinement with block FMAs, about which our analysis provides new insight. Our framework is applicable to the tensor core units in the NVIDIA Volta and Turing GPUs. For these we compare matrix multiplication and LU factorization with TC16 and TC32 forms of FMA, which differ in the precision used for the output of the tensor cores. Our experiments on an NVDIA V100 GPU confirm the predictions of the analysis that the TC32 variant is much more accurate than the TC16 one, and they show that the accuracy boost is obtained with almost no performance loss. [ABSTRACT FROM AUTHOR]

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

We show that there are CNF formulas which can be refuted in resolution in both small space and small width, but for which any small-width proof must have space exceeding by far the linear worst-case upper bound. This significantly strengthens the space-width trade-offs in [E. Ben-Sasson, SIAM J. Comput., 38 (2009), pp. 2511{2525], and provides one more example of tradeo ffs in the \supercritical" regime above worst case recently identified by [A.A. Razborov, J. ACM, 63 (2016), 16]. We obtain our results by using Razborov's new hardness condensation technique and combining it with the space lower bounds in [E. Ben-Sasson and J. Nordstrom, Short proofs may be spacious: An optimal separation of space and length in resolution, in Proceedings of the 49th Annual IEEE Symposium on Foundations of Computer Science (FOCS '08), 2008, pp. 709{718]. [ABSTRACT FROM AUTHOR]

Set functions with convenient properties (such as submodularity) appear in application areas of current interest, such as algorithmic game theory, and allow for improved optimization algorithms. It is natural to ask (e.g., in the context of data driven optimization) how robust such properties are, and whether small deviations from them can be tolerated. We consider two such questions in the important special case of linear set functions. One question that we address is whether any set function that approximately satisfies the modularity equation (linear functions satisfy the modularity equation exactly) is close to a linear function. The answer to this is positive (in a precise formal sense) as shown by Kalton and Roberts [Trans. Amer. Math. Soc., 278 (1983), pp. 803--816] (and further improved by Bondarenko, Prymak, and Radchenko [J. Math. Anal. Appl., 402 (2013), pp. 234--241]). We revisit their proof idea that is based on expander graphs and provide significantly stronger upper bounds by combining it with new techniques. Furthermore, we provide improved lower bounds for this problem. Another question that we address is that of how to learn a linear function h that is close to an approximately linear function f, while querying the value of f on only a small number of sets. We present a deterministic algorithm that makes only linearly many (in the number of items) nonadaptive queries, and thus improve upon a previous algorithm of Chierichetti, Das, Dasgupta, and Kumar [Proceedings of the 56th Symposium on Foundations of Computer Science, 2015, pp. 1143--1162] that is randomized and makes more than a quadratic number of queries. Our learning algorithm is based on the Hadamard transform. [ABSTRACT FROM AUTHOR]

Multilevel partitioning methods that are inspired by principles of multiscaling are the most powerful practical hypergraph partitioning solvers. Hypergraph partitioning has many applications in disciplines ranging from scientific computing to data science. In this paper we introduce the concept of algebraic distance on hypergraphs and demonstrate its use as an algorithmic component in the coarsening stage of multilevel hypergraph partitioning solvers. The algebraic distance is a vertex distance measure that extends hyperedge weights for capturing the local connectivity of vertices which is critical for hypergraph coarsening schemes. The practical effectiveness of the proposed measure and corresponding coarsening scheme is demonstrated through extensive computational experiments on a diverse set of problems. Finally, we propose a benchmark of hypergraph partitioning problems to compare the quality of other solvers. [ABSTRACT FROM AUTHOR]

In this work, we present a modification of the traditional integral deferred correction (IDC) approach that significantly reduces the storage requirements of the methods. These methods, which we call low-storage IDC or LS-IDC methods, require storing only one copy of each stage vector throughout the iteration process, whereas traditional IDC methods require two copies of each vector. We prove that LS-IDC methods converge with the same formal order of accuracy as traditional IDC methods as the timestep size approaches zero for the case of linear, constant-coefficient systems. A collection of numerical tests are used to evaluate the stability properties of the low-storage methods, and a nonlinear ODE and a linear transport equation are used to compare the accuracy and storage requirements of LS-IDC integrators with other fully implicit schemes. The results demonstrate that LS-IDC methods have similar accuracy but significantly reduced memory requirements compared to other fully implicit methods. [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.

There are three common parallel sparse matrix-vector multiply algorithms: 1D row-parallel, 1D column-parallel, and 2D row-column-parallel. The 1D parallel algorithms offer the advantage of having only one communication phase. On the other hand, the 2D parallel algorithm is more scalable, but it suffers from two communication phases. Here, we introduce a novel concept of heterogeneous messages where a heterogeneous message may contain both input-vector entries and partially computed output-vector entries. This concept not only leads to a decreased number of messages but also enables fusing the input- and output-communication phases into a single phase. These findings are exploited to propose a 1.5D parallel sparse matrix-vector multiply algorithm which is called local row-column-parallel. This proposed algorithm requires a constrained fine-grain partitioning in which each fine-grain task is assigned to the processor that contains either its input- vector entry, its output-vector entry, or both. We propose two methods to carry out the constrained fine-grain partitioning. We conduct our experiments on a large set of test matrices to evaluate the partitioning qualities and partitioning times of these proposed 1.5D methods. [ABSTRACT FROM AUTHOR]

Computation of the trace of a matrix function plays an important role in many scientific computing applications, including applications in machine learning, computational physics (e.g., lattice quantum chromodynamics), network analysis, and computational biology (e.g., protein folding), just to name a few application areas. We propose a linear-time randomized algorithm for approximating the trace of matrix functions of large symmetric matrices. Our algorithm is based on coupling function approximation using Chebyshev interpolation with stochastic trace estimators (Hutchinson's method), and as such requires only implicit access to the matrix, in the form of a function that maps a vector to the product of the matrix and the vector. We provide rigorous approximation error in terms of the extremal eigenvalue of the input matrix, and the Bernstein ellipse that corresponds to the function at hand. Based on our general scheme, we provide algorithms with provable guarantees for important matrix computations, including log-determinant, trace of matrix inverse, Estrada index, Schatten p-norm, and testing positive definiteness. We experimentally evaluate our algorithm and demonstrate its effectiveness on matrices with tens of millions dimensions. [ABSTRACT FROM AUTHOR]

JuMP is an open-source modeling language that allows users to express a wide range of optimization problems (linear, mixed-integer, quadratic, conic-quadratic, semidefinite, and nonlinear) in a high-level, algebraic syntax. JuMP takes advantage of advanced features of the Julia programming language to offer unique functionality while achieving performance on par with commercial modeling tools for standard tasks. In this work we will provide benchmarks, present the novel aspects of the implementation, and discuss how JuMP can be extended to new problem classes and composed with state-of-the-art tools for visualization and interactivity. [ABSTRACT FROM AUTHOR]

Bridging cultures that have often been distant, Julia combines expertise from the diverse fields of computer science and computational science to create a new approach to numerical computing. Julia is designed to be easy and fast and questions notions generally held to be "laws of nature" by practitioners of numerical computing: 1. High-level dynamic programs have to be slow. 2. One must prototype in one language and then rewrite in another language for speed or deployment. 3. There are parts of a system appropriate for the programmer, and other parts that are best left untouched as they have been built by the experts. We introduce the Julia programming language and its design—a dance between specialization and abstraction. Specialization allows for custom treatment. Multiple dispatch, a technique from computer science, picks the right algorithm for the right circumstance. Abstraction, which is what good computation is really about, recognizes what remains the same after differences are stripped away. Abstractions in mathematics are captured as code through another technique from computer science, generic programming. Julia shows that one can achieve machine performance without sacrificing human convenience. [ABSTRACT FROM AUTHOR]

Reduced basis methods are projection-based model order reduction techniques for reducing the computational complexity of solving parametrized partial differential equation problems. In this work we discuss the design of pyMOR, a freely available software library of model order reduction algorithms, in particular reduced basis methods, implemented with the Python programming language. As its main design feature, all reduction algorithms in pyMOR are implemented generically via operations on well-defined vector array, operator, and discretization interface classes. This allows for an easy integration with existing open-source high-performance partial differential equation solvers without adding any model reduction specific code to these solvers. Besides an in-depth discussion of pyMOR's design philosophy and architecture, we present several benchmark results and numerical examples showing the feasibility of our approach. [ABSTRACT FROM AUTHOR]

MADNESS (multiresolution adaptive numerical environment for scientific simulation) is a high-level software environment for solving integral and differential equations in many dimensions that uses adaptive and fast harmonic analysis methods with guaranteed precision that are based on multiresolution analysis and separated representations. Underpinning the numerical capabilities is a powerful petascale parallel programming environment that aims to increase both programmer productivity and code scalability. This paper describes the features and capabilities of MADNESS and briefly discusses some current applications in chemistry and several areas of physics. [ABSTRACT FROM AUTHOR]

CUDA, OpenCL, and OpenMP are popular programming models for the multicore architectures of CPUs and many-core architectures of GPUs or Xeon Phis. At the same time, computational scientists face the question of which programming model to use to obtain their scientific results. We present the linear algebra library ViennaCL, which is built on top of all three programming models, thus enabling computational scientists to interface to a single library, yet obtain high performance for all three hardware types. Since the respective compute back end can be selected at runtime, one can seamlessly switch between different hardware types without the need for errorprone and time-consuming recompilation steps. We present new benchmark results for sparse linear algebra operations in ViennaCL, complementing results for the dense linear algebra operations in ViennaCL reported in earlier work. Comparisons with vendor libraries show that ViennaCL provides better overall performance for sparse matrix-vector and sparse matrix-matrix products. Additional benchmark results for pipelined iterative solvers with kernel fusion and preconditioners identify the respective sweet spots for CPUs, Xeon Phis, and GPUs. [ABSTRACT FROM AUTHOR]

With the growing need of big-data processing in diverse application domains, MapReduce (e.g., Hadoop) has become one of the standard computing paradigms for large-scale computing on a cluster system. Despite its popularity, the current MapReduce framework suffers from in exibility and inefficiency inherent to its programming model and system architecture. In or- der to address these problems, we propose Vispark, a novel extension of Spark for GPU-accelerated MapReduce processing on array-based scientific computing and image processing tasks. Vispark provides an easy-to-use, Python-like high-level language syntax and a novel data abstraction for MapReduce programming on a GPU cluster system. Vispark introduces a programming abstraction for accessing neighbor data in the mapper function, which greatly simplifies many image processing tasks using MapReduce by reducing memory footprints and bypassing the reduce stage. Vispark provides socket-based halo communication that synchronizes between data partitions transparently from the users, which is necessary for many scientific computing problems in distributed systems. Vispark also provides domain-specific functions and language supports specifically designed for high- performance computing and image processing applications. We demonstrate the performance of our prototype system on several visual computing tasks, such as image processing, volume rendering, K-means clustering, and heat transfer simulation. [ABSTRACT FROM AUTHOR]

We present PuLP (partitioning using label propagation), a parallel and memory- efficient graph partitioning method specifically designed to partition low-diameter networks with skewed degree distributions on shared-memory multicore platforms. Graph partitioning is an im- portant problem in scientific computing because it impacts the execution time and energy efficiency of computations on distributed-memory platforms. Partitioning determines the in-memory layout of a graph, which affects locality, intertask load balance, communication time, and overall memory utilization. A novel feature of our PuLP method is that it optimizes for multiple objective metrics simultaneously, while satisfying multiple partitioning constraints. Using our method, we are able to partition a web crawl with billions of edges on a single compute server in under a minute. For a collection of test graphs, we show that PuLP uses up to 7:8× less memory than state-of-the-art partitioners and is 5:0× faster, on average, than alternate approaches (with 16-way parallelism). We also achieve better partitioning quality results for the multiobjective scenario. [ABSTRACT FROM AUTHOR]

The use of composable abstractions allows the application of new and established algorithms to a wide range of problems, while automatically inheriting the benefits of well-known performance optimizations. This work highlights the composition of the PETSc DMPlex domain topology abstraction with the Firedrake automated finite element system to create a PDE solving environment that combines expressiveness, flexibility, and high performance. We describe how Firedrake utilizes DMPlex to provide the indirection maps required for finite element assembly, while supporting various mesh input formats and runtime domain decomposition. In particular, we describe how DMPlex and its accompanying data structures allow the generic creation of user-defined discretizations, while utilizing data layout optimizations that improve cache coherency and ensure overlapped communication during assembly computation. [ABSTRACT FROM AUTHOR]

Given the square matrices A,B,D,E and the matrix C of conforming dimensions, we consider the linear matrix equation AXE + DXB = C in the unknown matrix X. Our aim is to provide an overview of the major algorithmic developments that have taken place over the past few decades in the numerical solution of this and related problems, which are producing reliable numerical tools in the formulation and solution of advanced mathematical models in engineering and scientific computing. [ABSTRACT FROM AUTHOR]

Multi-time-scale stochastic dynamical systems are ubiquitous in science and engineering, and the reduction of such systems and their models to only their slow components is often essential for scientific computation and further analysis. Rather than being available in the form of an explicit analytical model, often such systems can only be observed as a data set which embodies dynamics on several time scales. We focus on applying and adapting data-mining and manifold learning techniques to detect the slow components in a class of such multiscale data. Traditional data-mining methods are based on metrics (and thus, geometries) which are not informed of the multiscale nature of the underlying system dynamics; such methods cannot successfully recover the slow variables. Here, we present an approach which utilizes both the local geometry and the local noise dynamics within the data set through a metric which is both insensitive to the fast variables and more general than simple statistical averaging. Our analysis of the approach provides conditions for successfully recovering the underlying slow variables, as well as an empirical protocol guiding the selection of the method parameters. Interestingly, the recovered underlying variables are gauge invariant--they are insensitive to the measuring instrument/observation function. [ABSTRACT FROM AUTHOR]

The forest-of-octrees approach to parallel adaptive mesh refinement and coarsening has recently been demonstrated in the context of a number of large-scale PDE-based applications. Efficient reference software has been made freely available to the public both in the form of the standalone p4est library and more indirectly by the general-purpose finite element library deal.II, which has been equipped with a p4est backend. Although linear octrees, which store only leaf octants, have an underlying tree structure by definition, it is not fully exploited in previously published mesh-related algorithms. This is because tree branches are not explicitly stored and because the topological relationships in meshes, such as the adjacency between cells, introduce dependencies that do not respect the octree hierarchy. In this work, we combine hierarchical and topological relationships between octants to design efficient recursive algorithms that operate on distributed forests of octrees. We present three important algorithms with recursive implementations. The first is a parallel search for leaves matching any of a set of multiple search criteria, such as leaves that contain points or intersect polytopes. The second is a ghost layer construction algorithm that handles arbitrarily refined octrees that are not covered by previous algorithms, which require a 2:1 condition between neighboring leaves. The third is a universal mesh topology iterator. This iterator visits every cell in a partition, as well as every interface (face, edge, and corner) between these cells. The iterator calculates the local topological information for every interface that it visits, taking into account the nonconforming interfaces that increase the complexity of describing the local topology. To demonstrate the utility of the topology iterator, we use it to compute the numbering and encoding of higher-order C0 nodal basis functions used for finite elements. We analyze the complexity of the new recursive algorithms theoretically and assess their performance, both in terms of single-processor efficiency and in terms of parallel scalability, demonstrating good weak and strong scaling up to 458,000 cores of the JUQUEEN supercomputer. [ABSTRACT FROM AUTHOR]

Nonnegative matrix factorization (NMF) is powerful in data analysis by recovering hidden features from redundant data. Recently, various extensions of NMF have been proposed to further enhance the classification results. In this paper, we present a novel proximal alternating nonnegative least squares (PANLS) framework for NMF. We show the globally convergent property that any accumulation point of PANLS is a stationary point of NMF. A new active set method based on the framework is proposed to solve each subproblem, which reduces to unconstrained optimization after finite steps and converges to the unique minimizer of the subproblem. The PANLS framework and the associated new active set method can be easily extended to variants of NMF. Extensive computational experiments using both synthetic data and real data show that our proposed method is attractive in terms of computational speed, the solution quality, and the flexibility to extensions of NMF. [ABSTRACT FROM AUTHOR]

We present a new method for solving parametrized linear systems. Under certain assumptions on the parametrization, solutions to the linear systems for all parameters can be accurately approximated by linear combinations of solutions to linear systems for a small set of fixed parameters. Combined with either direct solvers or preconditioned iterative solvers for each linear system with a fixed parameter, the method is particularly suitable for situations when solutions to a large number of distinct parameters or a large number of right-hand sides are required. The method is also simple to parallelize. We demonstrate the applicability of the method to the calculation of the response functions in electronic structure theory. We demonstrate the numerical performance of the method using a benzene molecule and a DNA molecule. [ABSTRACT FROM AUTHOR]

We study the cross-entropy method for diffusions. One of the results is a versatile cross-entropy algorithm that can be used to design efficient importance sampling strategies for rare events or to solve optimal control problems. The approach is based on the minimization of a suitable cross-entropy functional, with a parametric family of exponentially tilted probability distributions. We illustrate the new algorithm with several numerical examples and discuss algorithmic issues and possible extensions of the method. [ABSTRACT FROM AUTHOR]

In parametric equations--stochastic equations are a special case--one may want to approximate the solution such that it is easy to evaluate its dependence on the parameters. Interpolation in the parameters is an obvious possibility--in this context often labeled as a collocation method. In the frequent situation where one has a "solver" for a given fixed parameter value, this may be used "nonintrusively" as a black-box component to compute the solution at all the interpolation points independently of each other. By extension, all other methods, and especially simple Galerkin methods, which produce some kind of coupled system, are often classed as "intrusive." We show how, for such "plain vanilla" Galerkin formulations, one may solve the coupled system in a nonintrusive way, and even the simplest form of block-solver has comparable efficiency. This opens at least two avenues for possible speed-up: first, to benefit from the coupling in the iteration by using more sophisticated block-solvers and, second, the possibility of nonintrusive successive rank-one updates as in the proper generalized decomposition (PGD). [ABSTRACT FROM AUTHOR]

The foliage of a plant performs vital functions. As such, leaf models are required to be developed for modeling the plant architecture from a set of scattered data captured using a scanning device. The leaf model can be used for purely visual purposes or as part of a further model, such as a fluid movement model or biological process. For these reasons, an accurate mathematical representation of the surface and boundary is required. This paper compares three approaches for fitting a continuously differentiable surface through a set of scanned data points from a leaf surface, with a technique already used for reconstructing leaf surfaces. The techniques which will be considered are discrete smoothing D²-splines [R. Arcangéli, M. C. Lopez de Silanes, and J. J. Torrens, Multidimensional Minimising Splines, Springer, New York, 2004], the thin plate spline finite element smoother [S. Roberts, M. Hegland, and I. Altas, SIAM J. Numer. Anal., 1 (2003), pp. 208-234], and the radial basis function Clough-Tocher method [M. N. Oqielat, I. W. Turner, and J. A. Belward, Appl. Math. Model., 33 (2009), pp. 2582-2595]. Numerical results show that discrete smoothing D²-splines produce reconstructed leaf surfaces which better represent the original physical leaf. [ABSTRACT FROM AUTHOR]

The aim of this work is to design monotonicity-preserving stabilized finite element techniques for transport problems as a blend of linear and nonlinear (shock-capturing) stabilization. As linear stabilization, we consider and analyze a novel symmetric projection stabilization technique based on a local Scott--Zhang projector. Next, we design a weighting of the aforementioned linear stabilization such that, when combined with a finite element discretization enjoying a discrete maximum principle (usually attained via nonlinear stabilization), it does not spoil these monotonicity properties. Then, we propose novel nonlinear stabilization schemes in the form of an artificial viscosity method where the amount of viscosity is proportional to gradient jumps at either finite element boundaries or nodes. For the nodal scheme, we prove a discrete maximum principle for time-dependent multidimensional transport problems. Numerical experiments support the numerical analysis and we show that the resulting methods provide excellent results. In particular, we observe that the proposed nonlinear stabilization techniques do an excellent job eliminating oscillations around shocks. [ABSTRACT FROM AUTHOR]

We study the convergence of multigrid schemes for the Helmholtz equation, focusing in particular on the choice of the coarse scale operators. Let Gc denote the number of points per wavelength at the coarse level. If the coarse scale solutions are to approximate the true solutions, then the oscillatory nature of the solutions implies the requirement Gc > 2. However, in examples the requirement is more like Gc ≳ 10, in a trade-off involving also the amount of damping present and the number of multigrid iterations. We conjecture that this is caused by the difference in phase speeds between the coarse and fine scale operators. Standard 5-point finite differences in two dimensions are our first example. A new coarse scale 9-point operator is constructed to match the fine scale phase speeds. We then compare phase speeds and multigrid performance of standard schemes with a scheme using the new operator. The required Gc is reduced from about 10 to about 3.5, with less damping present so that waves propagate over > 100 wavelengths in the new scheme. Next, we consider extensions of the method to more general cases. In three dimensions, comparable results are obtained with standard 7-point differences and optimized 27-point coarse grid operators, leading to an order of magnitude reduction in the number of unknowns for the coarsest scale linear system. Finally, we show how to include perfectly matched layers at the boundary, using a regular grid finite element method. Matching coarse scale operators can easily be constructed for other discretizations. The method is therefore potentially useful for a large class of discretized high-frequency Helmholtz equations. [ABSTRACT FROM AUTHOR]

We show that if G is a 4-critical graph embedded in a fixed surface Σ so that every contractible cycle has length at least 5, then G can be expressed as G = G'∪G1∪2∪ ... ∪Gk, where |V(G')| and k are bounded by a constant (depending linearly on the genus of ∑) and G1, ..., Gk are graphs (of unbounded size) whose structure we describe exactly. The proof is computer assisted--we use a computer to enumerate all plane 4-critical graphs of girth 5 with a precolored cycle of length at most 16 that are used in the basic case of the inductive proof of the statement. [ABSTRACT FROM AUTHOR]

We investigate the problem of symmetrically permuting a square sparse matrix into a block diagonal form with overlap. This permutation problem arises in the parallelization of an explicit formulation of the multiplicative Schwarz preconditioner and a more recent block overlapping banded linear solver as well as its application to general sparse linear systems. In order to formulate this permutation problem as a graph theoretical problem, we define a constrained version of the multiway graph partitioning by vertex separator (GPVS) problem, which is referred to as the ordered GPVS (oGPVS) problem. However, existing graph partitioning tools are unable to solve the oGPVS problem. So, we also show how the recursive bipartitioning framework can be utilized for solving the oGPVS problem. For this purpose, we propose a left-to-right bipartitioning approach together with a novel vertex fixation scheme so that existing 2-way GPVS tools that support fixed vertices can be effectively and efficiently utilized in the recursive bipartitioning framework. Experimental results on a wide range of matrices confirm the validity of the proposed approach. [ABSTRACT FROM AUTHOR]

Measuring the connection strength between a pair of vertices ill a graph is one of the most important concerns in many graph applications. Simple measures such as edge weights may not be sufficient for capturing the effects associated with short paths of lengths greater than one. In this paper, we consider an iterative process that smooths an associated value for nearby vertices, and we present a measure of the local connection strength (called the algebraic distance; see [D. Ron, I. Safro, and A. Brandt, Multiscale Model. Simul., 9 (2011), pp. 407-423]) based on this process. The proposed measure is attractive in that the process is simple, linear, and easily parallelized. An analysis of the convergence property of the process reveals that the local neighborhoods play an important role in determining the connectivity between vertices. We demonstrate the practical effectiveness of the proposed measure through several combinatorial optimization problems on graphs and hypergraphs. [ABSTRACT FROM AUTHOR]

A typical first step of a direct solver for the linear system Mx = b is reordering of the symmetric matrix M to improve execution time and space requirements of the solution process. In this work, we propose a novel nested-dissection-based ordering approach that utilizes hypergraph partitioning. Our approach is based on the formulation of graph partitioning by vertex separator (GPVS) problem as a hypergraph partitioning problem. This new formulation is immune to deficiency of GPVS in a multilevel framework and hence enables better orderings. In matrix terms, our method relies on the existence of a structural factorization of the input M matrix in the form of M = AAT (or M = AD²AT). We show that the partitioning of the row-net hypergraph representation of the rectangular matrix A induces a GPVS of the standard graph representation of matrix M. In the absence of such factorization, we also propose simple, yet effective structural factorization techniques that are based on finding an edge clique cover of the standard graph representation of matrix M, and hence applicable to any arbitrary symmetric matrix M. Our experimental evaluation has shown that the proposed method achieves better ordering in comparison to state-of-the-art graph-based ordering tools even for symmetric matrices where structural M = AATM factorization is not provided as an input. For matrices coming from linear programming problems, our method enables even faster and better orderings. [ABSTRACT FROM AUTHOR]

We present scalable algorithms for parallel adaptive mesh refinement and coarsening (AMR), partitioning, and 2:1 balancing on computational domains composed of multiple connected two-dimensional quadtrees or three-dimensional octrees, referred to as a forest of octrees. By distributing the union of octants from all octrees in parallel, we combine the high scalability proven previously for adaptive single-octree algorithms with the geometric flexibility that can be achieved by arbitrarily connected hexahedral macromeshes. in which each macroelement is the root of an adapted octree. A key concept of our approach is an encoding scheme of the interoctree connectivity that permits arbitrary relative orientations between octrees. Based on this encoding we develop interoctree transformations of octants. These form the basis for high-level parallel octree algorithms, which are designed to interact with an application code such as a numerical solver for partial differential equations. We have implemented and tested these algorithms in the p4est software library. We demonstrate the parallel scalability of p4est on its own and in combination with two geophysics codes. Using p4est we generate and adapt multioctree meshes with up to 5.13 x 1011 octants on as many as 220,320 CPU cores and execute the 2:1 balance algorithm in less than 10 seconds per million octants per process. [ABSTRACT FROM AUTHOR]

Data flow acyclic directed graphs (digraphs) can be applied to accurately describe the data dependency for a wide range of grid-based scientific computing applications ranging from numerical algebra to realistic applications of radiation or neutron transport. The parallel computing of these applications is equivalent to the parallel execution of digraphs. This paper presents a framework of scalable heuristic algorithms for the parallel execution of digraphs. This framework consists of three components: the heuristic partitioning method of a digraph, the parallel sweeping algorithm for a partitioned digraph, and the heuristic strategy for vertex scheduling and vertex packing. Evaluation rules of heuristic algorithms are presented for better theoretical understanding and performance optimization. Parallel benchmarks for the multigroup neutron or radiation Sn transport using processors from 100 to 2048 on two massively parallel machines show that these heuristic algorithms scale well. [ABSTRACT FROM AUTHOR]

It has been shown that combinatorial optimization of matrix-vector multiplication can lead to faster evaluation of finite element stiffness matrices. Based on a graph model characterizing relationships between rows, an efficient set of operations can be generated to perform matrix-vector multiplication for this problem. We improve the graph model by extending the set of binary row relationships and solve this combinatorial optimization problem optimally for the binary row relationships implemented, yielding significantly improved results over previously published graph models. We also extend the representation by using hypergraphs to model more complicated row relationships, expressing a three-row relationship with a three-vertex hyperedge, for example. Our initial greedy algorithm for this hypergraph model has yielded significantly better results than the graph model for many matrices. [ABSTRACT FROM AUTHOR]

A novel type of microfluidic biochip employs acoustic surface waves as a pumping mechanism for fluids. Piezoelectric substrate materials in the chip efficiently convert an electromagnetic input signal into an elastic wave confined to the surface layer of the substrate, hence the surface acoustic wave. The interaction of these waves with an adjoining fluid volume leads to streaming patterns in the fluid or of the motion of the fluid as a whole. In this work we study this acoustic streaming effect. We present a two-scale mathematical model, some theoretical results on the well-posedness of the model problem, as well as numerical simulations based on the finite element method. The model will include the case of free capillary fluid boundaries. [ABSTRACT FROM AUTHOR]

A letter to the editor is presented in response to the article "Full Waveform Inversion and the Truncated Newton Method" by L. Métivier and colleagues in the previous issue.