Search

Recently different formulations of the first-order Lighthill-Whitham-Richards (LWR) model have been identified in different coordinates and state variables. However, relationships between higher-order continuum and car-following traffic flow models are still not well understood. In this study, we first categorize traffic flow models according to their coordinates, state variables, and orders in the three-dimensional representation of traffic flow and propose a unified approach to convert higher-order car-following models into continuum models and vice versa. The conversion method consists of two steps: equivalent transformations between the secondary Eulerian (E-S) formulations and the primary Lagrangian (L-P) formulations, and approximations of L-P derivatives with anisotropic (upwind) finite differences. We use the method to derive continuum models from general second- and third-order car-following models and derive car-following models from second-order continuum models. Furthermore, we demonstrate that corresponding higher-order continuum and car-following models have the same fundamental diagrams, and that the string stability conditions for vehicle-continuous car-following models are the same as the linear stability conditions for the corresponding continuum models. A numerical example verifies the analytical results. In a sense, we establish a weak equivalence between continuum and car-following models, subject to errors introduced by the finite difference approximation. Such an equivalence relation can help us to pick out anisotropic solutions of higher-order models with non-concave fundamental diagrams.

Implicit samplers are algorithms for producing independent, weighted samples from multivariate probability distributions. These are often applied in Bayesian data assimilation algorithms. We use Laplace asymptotic expansions to analyze two implicit samplers in the small noise regime. Our analysis suggests a symmetrization of the algorithms that leads to improved implicit sampling schemes at a relatively small additional cost. Computational experiments confirm the theory and show that symmetrization is effective for small noise sampling problems.© 2016 Wiley Periodicals, Inc.

The single particle energies obtained in a Kohn-Sham density functional theory (DFT) calculation are generally known to be poor approximations to electron excitation energies that are measured in transport, tunneling and spectroscopic experiments such as photo-emission spectroscopy. The correction to these energies can be obtained from the poles of a single particle Green’s function derived from a many-body perturbation theory. From a computational perspective, the accuracy and efficiency of such an approach depends on how a self energy term that properly accounts for dynamic screening of electrons is approximated. The G0W0 approximation is a widely used technique in which the self energy is expressed as the convolution of a noninteracting Green’s function (G0) and a screened Coulomb interaction (W0) in the frequency domain. The computational cost associated with such a convolution is high due to the high complexity of evaluating W0 at multiple frequencies. In this paper, we discuss how the cost of G0W0 calculation can be reduced by constructing a low rank approximation to the frequency dependent part of W0. In particular, we examine the effect of such a low rank approximation on the accuracy of the G0W0 approximation. We also discuss how the numerical convolution of G0 and W0 can be evaluated efficiently and accurately by using a contour deformation technique with an appropriate choice of the contour.

We present several numerical methods and establish their error estimates forthe discretization of the nonlinear Dirac equation in the nonrelativistic limitregime, involving a small dimensionless parameter $0

The single particle energies obtained in a Kohn-Sham density functional theory (DFT) calculation are generally known to be poor approximations to electron excitation energies that are measured in transport, tunneling and spectroscopic experiments such as photo-emission spectroscopy. The correction to these energies can be obtained from the poles of a single particle Green’s function derived from a many-body perturbation theory. From a computational perspective, the accuracy and efficiency of such an approach depends on how a self energy term that properly accounts for dynamic screening of electrons is approximated. The G0W0 approximation is a widely used technique in which the self energy is expressed as the convolution of a noninteracting Green’s function (G0) and a screened Coulomb interaction (W0) in the frequency domain. The computational cost associated with such a convolution is high due to the high complexity of evaluating W0 at multiple frequencies. In this paper, we discuss how the cost of G0W0 calculation can be reduced by constructing a low rank approximation to the frequency dependent part of W0. In particular, we examine the effect of such a low rank approximation on the accuracy of the G0W0 approximation. We also discuss how the numerical convolution of G0 and W0 can be evaluated efficiently and accurately by using a contour deformation technique with an appropriate choice of the contour.

Motivated by classical vortex blob methods for the Euler equations, wedevelop a numerical blob method for the aggregation equation. This provides acounterpoint to existing literature on particle methods. By regularizing thevelocity field with a mollifier or "blob function", the blob method has afaster rate of convergence and allows a wider range of admissible kernels. Infact, we prove arbitrarily high polynomial rates of convergence to classicalsolutions, depending on the choice of mollifier. The blob method conserves massand the corresponding particle system is both energy decreasing for aregularized free energy functional and preserves the Wasserstein gradient flowstructure. We consider numerical examples that validate our predicted rate ofconvergence and illustrate qualitative properties of the method.

A comprehensive analysis of the morphology of the solution space for a special type of quadratic quaternion equation is presented. This equation, which arises in a surface construction problem, incorporates linear terms in a quaternion variable and its conjugate with right and left quaternion coefficients, while the quadratic term has a quaternion coefficient placed between the variable and its conjugate. It is proved that, for generic coefficients, the equation has two, one, or no solutions, but in certain special instances the solution set may comprise a circle or a 3-sphere in the quaternion space H. The analysis yields solutions for each case, and intuitive interpretations of them in terms of the four-dimensional geometry of the quaternion space H.

We study the numerical performance of a continuous data assimilation (downscaling) algorithm, based on ideas from feedback control theory, in the context of the two-dimensional incompressible Navier-Stokes equations. Our model problem is to recover an unknown reference solution, asymptotically in time, by using continuous-in-time coarse-mesh nodal-point observational measurements of the velocity field of this reference solution (subsampling), as might be measured by an array of weather vane anemometers. Our calculations show that the required nodal observation density is remarkably less than what is suggested by the analytical study; and is in fact comparable to the number of numerically determining Fourier modes, which was reported in an earlier computational study by the authors. Thus, this method is computationally efficient and performs far better than the analytical estimates suggest.

We study the Riemannian optimization methods on the embedded manifold of low rank matrices for the problem of matrix completion, which is about recovering a low rank matrix from its partial entries. Assume $m$ entries of an $n\times n$ rank $r$ matrix are sampled independently and uniformly with replacement. We first prove that with high probability the Riemannian gradient descent and conjugate gradient descent algorithms initialized by one step hard thresholding are guaranteed to converge linearly to the measured matrix provided \begin{align*} m\geq C_\kappa n^{1.5}r\log^{1.5}(n), \end{align*} where $C_\kappa$ is a numerical constant depending on the condition number of the underlying matrix. The sampling complexity has been further improved to \begin{align*} m\geq C_\kappa nr^2\log^{2}(n) \end{align*} via the resampled Riemannian gradient descent initialization. The analysis of the new initialization procedure relies on an asymmetric restricted isometry property of the sampling operator and the curvature of the low rank matrix manifold. Numerical simulation shows that the algorithms are able to recover a low rank matrix from nearly the minimum number of measurements.

We prove a finiteness principle for interpolation of data by nonnegative Cm functions. Our result raises the hope that one can start to understand constrained interpolation problems in which e.g. the interpolating function F is required to be nonnegative.

The first order condition of the constrained minimization problem leads to asaddle point problem. A multigrid method using a multiplicative Schwarzsmoother for saddle point problems can thus be interpreted as a successivesubspace optimization method based on a multilevel decomposition of theconstraint space. Convergence theory is developed for successive subspaceoptimization methods based on two assumptions on the space decomposition:stable decomposition and strengthened Cauchy-Schwarz inequality, andsuccessfully applied to the saddle point systems arising from mixed finiteelement methods for Poisson and Stokes equations. Uniform convergence isobtained without the full regularity assumption of the underlying partialdifferential equations. As a byproduct, a V-cycle multigrid method fornon-conforming finite elements is developed and proved to be uniform convergentwith even one smoothing step.

Given a set of mixtures, blind source separation attempts to retrieve the source signals without or with very little information of the mixing process. We present a geometric approach for blind separation of nonnegative linear mixtures termed facet component analysis. The approach is based on facet identification of the underlying cone structure of the data. Earlier works focus on recovering the cone by locating its vertices (vertex component analysis) based on a mutual sparsity condition which requires each source signal to possess a stand-alone peak in its spectrum. We formulate alternative conditions so that enough data points fall on the facets of a cone instead of accumulating around the vertices. To find a regime of unique solvability, we make use of both geometric and density properties of the data points and develop an efficient facet identification method by combining data classification and linear regression. For noisy data, total variation technique may be employed. We show computational results on nuclear magnetic resonance spectroscopic data to substantiate our method.

Sparse matrix-matrix multiplication (or SpGEMM) is a key primitive for many high-performance graph algorithms as well as for some linear solvers, such as algebraic multigrid. The scaling of existing parallel implementations of SpGEMM is heavily bound by communication. Even though 3D (or 2.5D) algorithms have been proposed and theoretically analyzed in the at MPI model on Erd}os{Rffenyi matrices, those algorithms had not been implemented in practice and their complexities had not been analyzed for the general case. In this work, we present the first implementation of the 3D SpGEMM formulation that exploits multiple (intranode and internode) levels of parallelism, achieving significant speedups over the state-of-the-art publicly available codes at all levels of concurrencies. We extensively evaluate our implementation and identify bottlenecks that should be subject to further research.

A new method for generating variance reduction parameters for strongly anisotropic, deep-penetration radiation shielding studies is presented. This method generates an alternate form of the adjoint scalar flux quantity, φ†Ω which is used by both CADIS and FW-CADIS to generate variance reduction parameters for local and global response functions, respectively. The new method, called CADIS-Ω, was implemented in the Denovo/ADVANTG software suite, and results are presented for a concrete labyrinth test problem. Results indicate that the flux generated by CADIS-Ω incorporates localized angular anisotropics in the flux effectively. CADIS-Ω outperformed CADIS in the test problem while obtaining correct results. This initial work indicates that CADIS-Ω may be highly useful for shielding problems with strong angular anisotropics. A future test plan to fully characterize the new method is proposed, which should reveal more about the types of realistic problems for which the CADIS-Ω will be suited.

We compare two approaches to the predictive modeling of dynamical systems from partial observations at discrete times. The first is continuous in time, where one uses data to infer a model in the form of stochastic differential equations, which are then discretized for numerical solution. The second is discrete in time, where one directly infers a discrete-time model in the form of a nonlinear autoregression moving average model. The comparison is performed in a special case where the observations are known to have been obtained from a hypoelliptic stochastic differential equation. We show that the discrete-time approach has better predictive skills, especially when the data are relatively sparse in time. We discuss open questions as well as the broader significance of the results.

We develop an immersed boundary (IB) method for modeling flows around fixed or moving rigid bodies that is suitable for a broad range of Reynolds numbers, including steady Stokes flow. The spatio-temporal discretization of the fluid equations is based on a standard staggered-grid approach. Fluid-body interaction is handled using Peskin's IB method; however, unlike existing IB approaches to such problems, we do not rely on penalty or fractional-step formulations. Instead, we use an unsplit scheme that ensures the no-slip constraint is enforced exactly in terms of the Lagrangian velocity field evaluated at the IB markers. Fractionalstep approaches, by contrast, can impose such constraints only approximately, which can lead to penetration of the flow into the body, and are inconsistent for steady Stokes flow. Imposing no-slip constraints exactly requires the solution of a large linear system that includes the fluid velocity and pressure as well as Lagrange multiplier forces that impose the motion of the body. The principal contribution of this paper is that it develops an efficient preconditioner for this exactly constrained IB formulation which is based on an analytical approximation to the Schur complement. This approach is enabled by the near translational and rotational invariance of Peskin's IB method. We demonstrate that only a few cycles of a geometric multigrid method for the fluid equations are required in each application of the preconditioner, and we demonstrate robust convergence of the overall Krylov solver despite the approximations made in the preconditioner. We empirically observe that to control the condition number of the coupled linear system while also keeping the rigid structure impermeable to fluid, we need to place the immersed boundary markers at a distance of about two grid spacings, which is significantly larger from what has been recommended in the literature for elastic bodies. We demonstrate the advantage of our monolithic solver over split solvers by computing the steady state flow through a two-dimensional nozzle at several Reynolds numbers. We apply the method to a number of benchmark problems at zero and finite Reynolds numbers, and we demonstrate first-order convergence of the method to several analytical solutions and benchmark computations.

The Bethe-Salpeter eigenvalue problem is a dense structured eigenvalue problem arising from discretized Bethe-Salpeter equation in the context of computing exciton energies and states. A computational challenge is that at least half of the eigenvalues and the associated eigenvectors are desired in practice. We establish the equivalence between Bethe-Salpeter eigenvalue problems and real Hamiltonian eigenvalue problems. Based on theoretical analysis, structure preserving algorithms for a class of Bethe-Salpeter eigenvalue problems are proposed. We also show that for this class of problems all eigenvalues obtained from the Tamm-Dancoff approximation are overestimated. In order to solve large scale problems of practical interest, we discuss parallel implementations of our algorithms targeting distributed memory systems. Several numerical examples are presented to demonstrate the efficiency and accuracy of our algorithms.

A comprehensive analysis of the morphology of the solution space for a special type of quadratic quaternion equation is presented. This equation, which arises in a surface construction problem, incorporates linear terms in a quaternion variable and its conjugate with right and left quaternion coefficients, while the quadratic term has a quaternion coefficient placed between the variable and its conjugate. It is proved that, for generic coefficients, the equation has two, one, or no solutions, but in certain special instances the solution set may comprise a circle or a 3-sphere in the quaternion space H. The analysis yields solutions for each case, and intuitive interpretations of them in terms of the four-dimensional geometry of the quaternion space H.

The Bethe-Salpeter eigenvalue problem is a dense structured eigenvalue problem arising from discretized Bethe-Salpeter equation in the context of computing exciton energies and states. A computational challenge is that at least half of the eigenvalues and the associated eigenvectors are desired in practice. We establish the equivalence between Bethe-Salpeter eigenvalue problems and real Hamiltonian eigenvalue problems. Based on theoretical analysis, structure preserving algorithms for a class of Bethe-Salpeter eigenvalue problems are proposed. We also show that for this class of problems all eigenvalues obtained from the Tamm-Dancoff approximation are overestimated. In order to solve large scale problems of practical interest, we discuss parallel implementations of our algorithms targeting distributed memory systems. Several numerical examples are presented to demonstrate the efficiency and accuracy of our algorithms.

We describe a method for calculating the roots of special functions satisfying second order linear ordinary differential equations. It exploits the recent observation that the solutions of a large class of such equations can be represented via nonoscillatory phase functions, even in the high-frequency regime. Our algorithm achieves near machine precision accuracy and the time required to compute one root of a solution is independent of the frequency of oscillations of that solution. Moreover, despite its great generality, our approach is competitive with specialized, state-of-the-art methods for the construction of Gaussian quadrature rules of large orders when it used in such a capacity. The performance of the scheme is illustrated with several numerical experiments and a Fortran implementation of our algorithm is available at the author's website.

Motivated by classical vortex blob methods for the Euler equations, wedevelop a numerical blob method for the aggregation equation. This provides acounterpoint to existing literature on particle methods. By regularizing thevelocity field with a mollifier or "blob function", the blob method has afaster rate of convergence and allows a wider range of admissible kernels. Infact, we prove arbitrarily high polynomial rates of convergence to classicalsolutions, depending on the choice of mollifier. The blob method conserves massand the corresponding particle system is both energy decreasing for aregularized free energy functional and preserves the Wasserstein gradient flowstructure. We consider numerical examples that validate our predicted rate ofconvergence and illustrate qualitative properties of the method.

We establish theoretical recovery guarantees of a family of Riemannian optimization algorithms for low rank matrix recovery, which is about recovering an $m\times n$ rank $r$ matrix from $p < mn$ number of linear measurements. The algorithms are first interpreted as iterative hard thresholding algorithms with subspace projections. Based on this connection, we show that provided the restricted isometry constant $R_{3r}$ of the sensing operator is less than $C_\kappa /\sqrt{r}$, the Riemannian gradient descent algorithm and a restarted variant of the Riemannian conjugate gradient algorithm are guaranteed to converge linearly to the underlying rank $r$ matrix if they are initialized by one step hard thresholding. Empirical evaluation shows that the algorithms are able to recover a low rank matrix from nearly the minimum number of measurements necessary.

Alternating projection (AP) of various forms, including the Parallel AP (PAP), Real-constrained AP (RAP) and the Serial AP (SAP), are proposed to solve phase retrieval with at most two coded diffraction patterns. The proofs of geometric convergence are given with sharp bounds on the rates of convergence in terms of a spectral gap condition. To compensate for the local nature of convergence, the null initialization is proposed for initial guess and proved to produce asymptotically accurate initialization for the case of Gaussian random measurement. Numerical experiments show that the null initialization produces more accurate initial guess than the spectral initialization and that AP converges faster to the true object than other iterative schemes for non-convex optimization such as the Wirtinger Flow. In numerical experiments, AP with the null initialization converges globally to the true object.

Many applications have benefited remarkably from low-dimensional models inthe recent decade. The fact that many signals, though high dimensional, areintrinsically low dimensional has given the possibility to recover them stablyfrom a relatively small number of their measurements. For example, incompressed sensing with the standard (synthesis) sparsity prior and in matrixcompletion, the number of measurements needed is proportional (up to alogarithmic factor) to the signal's manifold dimension. Recently, a new natural low-dimensional signal model has been proposed: thecosparse analysis prior. In the noiseless case, it is possible to recoversignals from this model, using a combinatorial search, from a number ofmeasurements proportional to the signal's manifold dimension. However, if weask for stability to noise or an efficient (polynomial complexity) solver, allthe existing results demand a number of measurements which is far removed fromthe manifold dimension, sometimes far greater. Thus, it is natural to askwhether this gap is a deficiency of the theory and the solvers, or if thereexists a real barrier in recovering the cosparse signals by relying only ontheir manifold dimension. Is there an algorithm which, in the presence ofnoise, can accurately recover a cosparse signal from a number of measurementsproportional to the manifold dimension? In this work, we prove that there is nosuch algorithm. Further, we show through numerical simulations that even in thenoiseless case convex relaxations fail when the number of measurements iscomparable to the manifold dimension. This gives a practical counter-example tothe growing literature on compressed acquisition of signals based on manifolddimension.

Douglas-Rachford (DR) algorithm is analyzed for Fourier phase retrieval with a single random phase mask. Local, geometric convergence to a unique fixed point is proved with numerical demonstration of global convergence.

Compared with the traditional spherical harmonics, the spherical needlets are a new generation of spherical wavelets that possess several attractive properties. Their double localization in both spatial and frequency domains empowers them to easily and sparsely represent functions with small spatial scale features. This paper is divided into two parts. First, it reviews the spherical harmonics and discusses their limitations in representing functions with small spatial scale features. To overcome the limitations, it introduces the spherical needlets and their attractive properties. In the second part of the paper, a Matlab package for the spherical needlets is presented. The properties of the spherical needlets are demonstrated by several examples using the package.

Many physical systems are described by nonlinear differential equations that are too complicated to solve in full. A natural way to proceed is to divide the variables into those that are of direct interest and those that are not, formulate solvable approximate equations for the variables of greater interest, and use data and statistical methods to account for the impact of the other variables. In the present paper we consider time-dependent problems and introduce a fully discrete solution method, which simplifies both the analysis of the data and the numerical algorithms. The resulting time series are identified by a NARMAX (nonlinear autoregression moving average with exogenous input) representation familiar from engineering practice. The connections with the Mori-Zwanzig formalism of statistical physics are discussed, as well as an application to the Lorenz 96 system.

Classical iterative methods for tomographic reconstruction include the class of Algebraic Reconstruction Techniques (ART). Convergence of these stationary linear iterative methods is however notably slow. In this paper we propose the use of Krylov solvers for tomographic linear inversion problems. These advanced iterative methods feature fast convergence at the expense of a higher computational cost per iteration, causing them to be generally uncompetitive without the inclusion of a suitable preconditioner. Combining elements from standard multigrid (MG) solvers and the theory of wavelets, a novel wavelet-based multi-level (WMG) preconditioner is introduced, which is shown to significantly speed-up Krylov convergence. The performance of the WMG-preconditioned Krylov method is analyzed through a spectral analysis, and the approach is compared to existing methods like the classical Simultaneous Iterative Reconstruction Technique (SIRT) and unpreconditioned Krylov methods on a 2D tomographic benchmark problem. Numerical experiments are promising, showing the method to be competitive with the classical Algebraic Reconstruction Techniques in terms of convergence speed and overall performance (CPU time) as well as precision of the reconstruction.

Classical iterative methods for tomographic reconstruction include the class of Algebraic Reconstruction Techniques (ART). Convergence of these stationary linear iterative methods is however notably slow. In this paper we propose the use of Krylov solvers for tomographic linear inversion problems. These advanced iterative methods feature fast convergence at the expense of a higher computational cost per iteration, causing them to be generally uncompetitive without the inclusion of a suitable preconditioner. Combining elements from standard multigrid (MG) solvers and the theory of wavelets, a novel wavelet-based multi-level (WMG) preconditioner is introduced, which is shown to significantly speed-up Krylov convergence. The performance of the WMG-preconditioned Krylov method is analyzed through a spectral analysis, and the approach is compared to existing methods like the classical Simultaneous Iterative Reconstruction Technique (SIRT) and unpreconditioned Krylov methods on a 2D tomographic benchmark problem. Numerical experiments are promising, showing the method to be competitive with the classical Algebraic Reconstruction Techniques in terms of convergence speed and overall performance (CPU time) as well as precision of the reconstruction.

We propose a new blow-up criterion for the 3D Euler equations ofincompressible fluid flows, based on the 3D Euler-Voigt inviscidregularization. This criterion is similar in character to a criterion proposedin a previous work by the authors, but it is stronger, and better adapted forcomputational tests. The 3D Euler-Voigt equations enjoy global well-posedness,and moreover are more tractable to simulate than the 3D Euler equations. Amajor advantage of these new criteria is that one only needs to simulate the 3DEuler-Voigt, and not the 3D Euler equations, to test the blow-up criteria, forthe 3D Euler equations, computationally.

We present a new method to compare the shapes of genus-zero surfaces. We introduce a measure of mutual stretching, the symmetric distortion energy, and establish the existence of a conformal diffeomorphism between any two genus-zero surfaces that minimizes this energy. We then prove that the energies of the minimizing diffeomorphisms give a metric on the space of genus-zero Riemannian surfaces. This metric and the corresponding optimal diffeomorphisms are shown to have properties that are highly desirable for applications.

The Immersed Boundary method is a simple, efficient, and robust numerical scheme for solving PDE in general domains, yet it only achieves first-order spatial accuracy near embedded boundaries. In this paper, we introduce a new high-order numerical method which we call the Immersed Boundary Smooth Extension (IBSE) method. The IBSE method achieves high-order accuracy by smoothly extending the unknown solution of the PDE from a given smooth domain to a larger computational domain, enabling the use of simple Cartesian-grid discretizations (e.g. Fourier spectral methods). The method preserves much of the flexibility and robustness of the original IB method. In particular, it requires minimal geometric information to describe the boundary and relies only on convolution with regularized delta-functions to communicate information between the computational grid and the boundary. We present a fast algorithm for solving elliptic equations, which forms the basis for simple, high-order implicit-time methods for parabolic PDE and implicit-explicit methods for related nonlinear PDE. We apply the IBSE method to solve the Poisson, heat, Burgers', and Fitzhugh-Nagumo equations, and demonstrate fourth-order pointwise convergence for Dirichlet problems and third-order pointwise convergence for Neumann problems.

We propose a modification of the weak Galerkin methods and show itsequivalence to a new version of virtual element methods. We also show theoriginal weak Galerkin method is equivalent to the non-conforming virtualelement method. As a consequence, ideas and techniques used for one method canbe transferred to another. The key of the connection is the degree of freedoms.

Polynomial chaos expansions are used to reduce the computational cost in the Bayesian solutions of inverse problems by creating a surrogate posterior that can be evaluated inexpensively. We show, by analysis and example, that when the data contain significant information beyond what is assumed in the prior, the surrogate posterior can be very different from the posterior, and the resulting estimates become inaccurate. One can improve the accuracy by adaptively increasing the order of the polynomial chaos, but the cost may increase too fast for this to be cost effective compared to Monte Carlo sampling without a surrogate posterior.

Polynomial chaos expansions are used to reduce the computational cost in the Bayesian solutions of inverse problems by creating a surrogate posterior that can be evaluated inexpensively. We show, by analysis and example, that when the data contain significant information beyond what is assumed in the prior, the surrogate posterior can be very different from the posterior, and the resulting estimates become inaccurate. One can improve the accuracy by adaptively increasing the order of the polynomial chaos, but the cost may increase too fast for this to be cost effective compared to Monte Carlo sampling without a surrogate posterior.

Implicit sampling is a weighted sampling method that is used in data assimilation to sequentially update state estimates of a stochastic model based on noisy and incomplete data. Here we apply implicit sampling to sample the posterior probability density of parameter estimation problems. The posterior probability combines prior information about the parameter with information from a numerical model, e.g., a partial differential equation (PDE), and noisy data. The result of our computations are parameters that lead to simulations that are compatible with the data. We demonstrate the usefulness of our implicit sampling algorithm with an example from subsurface flow. For an efficient implementation, we make use of multiple grids, BFGS optimization coupled to adjoint equations, and Karhunen-Loève expansions for dimensional reduction. Several difficulties of Markov chain Monte Carlo methods, e.g., estimation of burn-in times or correlations among the samples, are avoided because the implicit samples are independent.

We survey results that have been obtained for self-adjoint operators, and especially Schrödinger operators, associated with mathematical models of quasicrystals. After presenting general results that hold in arbitrary dimensions, we focus our attention on the one-dimensional case, and in particular on several key examples. The most prominent of these is the Fibonacci Hamiltonian, for which much is known by now and to which an entire section is devoted here. Other examples that are discussed in detail are given by the more general class of Schrödinger operators with Sturmian potentials. We put some emphasis on the methods that have been introduced quite recently in the study of these operators, many of them coming from hyperbolic dynamics. We conclude with a multitude of numerical calculations that illustrate the validity of the known rigorous results and suggest conjectures for further exploration.

We survey results that have been obtained for self-adjoint operators, and especially Schrödinger operators, associated with mathematical models of quasicrystals. After presenting general results that hold in arbitrary dimensions, we focus our attention on the one-dimensional case, and in particular on several key examples. The most prominent of these is the Fibonacci Hamiltonian, for which much is known by now and to which an entire section is devoted here. Other examples that are discussed in detail are given by the more general class of Schrödinger operators with Sturmian potentials. We put some emphasis on the methods that have been introduced quite recently in the study of these operators, many of them coming from hyperbolic dynamics. We conclude with a multitude of numerical calculations that illustrate the validity of the known rigorous results and suggest conjectures for further exploration.

Three complementary methods have been implemented in the code Denovo that accelerate neutral particle transport calculations with methods that use leadership-class computers fully and effectively: a multigroup block (MG) Krylov solver, a Rayleigh quotient iteration (RQI) eigenvalue solver, and a multigrid in energy preconditioner. The multigroup Krylov solver converges more quickly than Gauss Seidel and enables energy decomposition such that Denovo can scale to hundreds of thousands of cores. The new multigrid in energy preconditioner reduces iteration count for many problem types and takes advantage of the new energy decomposition such that it can scale efficiently. These two tools are useful on their own, but together they enable the RQI eigenvalue solver to work. Each individual method has been described before, but this is the first time they have been demonstrated to work together effectively. RQI should converge in fewer iterations than power iteration (PI) for large and challenging problems. RQI creates shifted systems that would not be tractable without the MG Krylov solver. It also creates ill-conditioned matrices that cannot converge without the multigrid in energy preconditioner. Using these methods together, RQI converged in fewer iterations and in less time than all PI calculations for a full pressurized water reactor core. It also scaled reasonably well out to 275,968 cores.

We present analysis for the reduction of an inertial description offluid-structure interactions subject to thermal fluctuations. We show how theviscous coupling between the immersed structures and the fluid can besimplified in the regime where this coupling becomes increasingly strong. Manydescriptions in fluid mechanics and in the formulation of computational methodsaccount for fluid-structure interactions through viscous drag terms to transfermomentum from the fluid to immersed structures. In the inertial regime, thiscoupling often introduces a prohibitively small time-scale into the temporaldynamics of the fluid-structure system. This is further exacerbated in thepresence of thermal fluctuations. We discuss here a systematic reductiontechnique for the full inertial equations to obtain a simplified descriptionwhere this coupling term is eliminated. This approach also accounts for theeffective stochastic equations for the fluid-structure dynamics. The analysisis based on use of the Infinitesmal Generator of the SPDEs and a singularperturbation analysis of the Backward Kolomogorov PDEs. We also discuss thephysical motivations and interpretation of the obtained reduced description ofthe fluid-structure system. Working paper currently under revision. Pleasereport any comments or issues to tabak.gil@gmail.com.

The maximum entropy principle is a powerful tool for solving underdeterminedinverse problems. This paper considers the problem of discretizing a continuousdistribution, which arises in various applied fields. We obtain theapproximating distribution by minimizing the Kullback-Leibler information(relative entropy) of the unknown discrete distribution relative to an initialdiscretization based on a quadrature formula subject to some momentconstraints. We study the theoretical error bound and the convergence of thisapproximation method as the number of discrete points increases. We prove that(i) the theoretical error bound of the approximate expectation of any boundedcontinuous function has at most the same order as the quadrature formula westart with, and (ii) the approximate discrete distribution weakly converges tothe given continuous distribution. Moreover, we present some numerical examplesthat show the advantage of the method and apply to numerically solving anoptimal portfolio problem.

The kinetic Monte Carlo (kMC) method is used in many scientific fields in applications involving rare-event transitions. Due to its discrete stochastic nature, efforts to parallelize kMC approaches often produce unbalanced time evolutions requiring complex implementations to ensure correct statistics. In the context of parallel kMC, the sequential update technique has shown promise by generating high quality distributions with high relative efficiencies for short-range systems. In this work, we provide an extension of the sequential update method in a parallel context that rigorously obeys detailed balance, which guarantees exact equilibrium statistics for all parallelization settings. Our approach also preserves nonequilibrium dynamics with minimal error for many parallelization settings, and can be used to achieve highly precise sampling. © 2014 Elsevier B.V. All rights reserved.

We describe an algorithm for the numerical solution of second order linear differential equations in the highly-oscillatory regime. It is founded on the recent observation that the solutions of equations of this type can be accurately represented using nonoscillatory phase functions. Unlike standard solvers for ordinary differential equations, the running time of our algorithm is independent of the frequency of oscillation of the solutions. We illustrate the performance of the method with several numerical experiments.

We observe that solutions of a large class of highly oscillatory second order linear ordinary differential equations can be approximated using nonoscillatory phase functions. In addition, we describe numerical experiments which illustrate important implications of this fact. For example, that many special functions of great interest --- such as the Bessel functions $J_\nu$ and $Y_\nu$ --- can be evaluated accurately using a number of operations which is $O(1)$ in the order $\nu$. The present paper is devoted to the development of an analytical apparatus. Numerical aspects of this work will be reported at a later date.

We introduce a version of the asymptotic expansions for Bessel functions $J_\nu(z)$, $Y_\nu(z)$ that is valid whenever $|z| > \nu$ (which is deep in the Fresnel regime), as opposed to the standard expansions that are applicable only in the Fraunhofer regime (i.e. when $|z| > \nu^2$). As expected, in the Fraunhofer regime our asymptotics reduce to the classical ones. The approach is based on the observation that Bessel's equation admits a non-oscillatory phase function, and uses classical formulas to obtain an asymptotic expansion for this function; this in turn leads to both an analytical tool and a numerical scheme for the efficient evaluation of $J_\nu(z)$, $Y_\nu(z)$, as well as various related quantities. The effectiveness of the technique is demonstrated via several numerical examples. We also observe that the procedure admits far-reaching generalizations to wide classes of second order differential equations, to be reported at a later date.

This work considers a computationally and statistically efficient parameter estimation method for a wide class of latent variable models-including Gaussian mixture models, hidden Markov models, and latent Dirichlet allocation-which exploits a certain tensor structure in their low-order observable moments (typically, of second- and third-order). Specifically, parameter estimation is reduced to the problem of extracting a certain (orthogonal) decomposition of a symmetric tensor derived from the moments; this decomposition can be viewed as a natural generalization of the singular value decomposition for matrices. Although tensor decompositions are generally intractable to compute, the decomposition of these specially structured tensors can be efficiently obtained by a variety of approaches, including power iterations and maximization approaches (similar to the case of matrices). A detailed analysis of a robust tensor power method is provided, establishing an analogue of Wedin's perturbation theorem for the singular vectors of matrices. This implies a robust and computationally tractable estimation approach for several popular latent variable models.

Rapid and reliable detection and identification of unknown chemical substances are critical to homeland security. It is challenging to identify chemical components from a wide range of explosives. There are two key steps involved. One is a non-destructive and informative spectroscopic technique for data acquisition. The other is an associated library of reference features along with a computational method for feature matching and meaningful detection within or beyond the library. In this paper, we develop a new iterative method to identify unknown substances from mixture samples of Raman spectroscopy. In the first step, a constrained least squares method decomposes the data into a sum of linear combination of the known components and a non-negative residual. In the second step, a sparse and convex blind source separation method extracts components geometrically from the residuals. Verification based on the library templates or expert knowledge helps to confirm these components. If necessary, the confirmed meaningful components are fed back into step one to refine the residual and then step two extracts possibly more hidden components. The two steps may be iterated until no more components can be identified. We illustrate the proposed method in processing a set of the so called swept wavelength optical resonant Raman spectroscopy experimental data by a satisfactory blind extraction of a priori unknown chemical explosives from mixture samples. We also test the method on nuclear magnetic resonance (NMR) spectra for chemical compounds identification. © 2013 Published by Elsevier B.V.

In this paper, we provide local and global convergence guarantees forrecovering CP (Candecomp/Parafac) tensor decomposition. The main step of theproposed algorithm is a simple alternating rank-$1$ update which is thealternating version of the tensor power iteration adapted for asymmetrictensors. Local convergence guarantees are established for third order tensorsof rank $k$ in $d$ dimensions, when $k=o \bigl( d^{1.5} \bigr)$ and the tensorcomponents are incoherent. Thus, we can recover overcomplete tensordecomposition. We also strengthen the results to global convergence guaranteesunder stricter rank condition $k \le \beta d$ (for arbitrary constant $\beta >1$) through a simple initialization procedure where the algorithm isinitialized by top singular vectors of random tensor slices. Furthermore, theapproximate local convergence guarantees for $p$-th order tensors are alsoprovided under rank condition $k=o \bigl( d^{p/2} \bigr)$. The guarantees alsoinclude tight perturbation analysis given noisy tensor.

Differential optical absorption spectroscopy (DOAS) is a powerful tool for detecting and quantifying trace gases in atmospheric chemistry [22]. DOAS spectra consist of a linear combination of complex multi-peak multi-scale structures. Most DOAS analysis routines in use today are based on least squares techniques, for example, the approach developed in the 1970s [18, 19, 20, 21] uses polynomial fits to remove a slowly varying background (broad spectral structures in the data), and known reference spectra to retrieve the identity and concentrations of reference gases [23]. An open problem [] is that fitting residuals for complex atmospheric mixtures often still exhibit structure that indicates the presence of unknown absorbers. In this work, we develop a novel three step semi-blind source separation method. The first step uses a multi-resolution analysis called empirical mode decomposition (EMD) to remove the slow-varying and fast-varying components in the DOAS spectral data matrix X. This has the advantage of avoiding user bias in fitting the slow varying signal. The second step decomposes the preprocessed data X in the first step into a linear combination of the reference spectra plus a remainder, or X = A S + R, where columns of matrix A are known reference spectra, and the matrix S contains the unknown non-negative coefficients that are proportional to concentration. The second step is realized by a convex minimization problem S = argminnorm(X - AS), where the norm is a hybrid I1fl2 norm (Huber estimator) that helps to maintain the non-negativity of S. Non-negative coefficients are necessary in order for the derived proportional concentrations to make physical sense. The third step performs a blind independent component analysis of the remainder matrix R to extract remnant gas components. This step demonstrates the ability of the new fitting method to extract orthogonal components without the use of reference spectra. We illustrate utility of the proposed method in processing a set of DOAS experimental data by a satisfactory blind extraction of an a-priori unknown trace gas (ozone) from the remainder matrix. Numerical results also show that the method can identify trace gases from the residuals. ©2014 American Institute of Mathematical Sciences.

This work considers a computationally and statistically efficient parameter estimation method for a wide class of latent variable models-including Gaussian mixture models, hidden Markov models, and latent Dirichlet allocation-which exploits a certain tensor structure in their low-order observable moments (typically, of second- and third-order). Specifically, parameter estimation is reduced to the problem of extracting a certain (orthogonal) decomposition of a symmetric tensor derived from the moments; this decomposition can be viewed as a natural generalization of the singular value decomposition for matrices. Although tensor decompositions are generally intractable to compute, the decomposition of these specially structured tensors can be efficiently obtained by a variety of approaches, including power iterations and maximization approaches (similar to the case of matrices). A detailed analysis of a robust tensor power method is provided, establishing an analogue of Wedin's perturbation theorem for the singular vectors of matrices. This implies a robust and computationally tractable estimation approach for several popular latent variable models.