Accurate and Efficient Expression Evaluation and Linear Algebra

We survey and unify recent results on the existence of accurate algorithms for evaluating multivariate polynomials, and more generally for accurate numerical linear algebra with structured matrices. By "accurate" we mean that the computed answer has relative error less than 1, i.e., has some correct leading digits. We also address efficiency, by which we mean algorithms that run in polynomial time in the size of the input. Our results will depend strongly on the model of arithmetic: Most of our results will use the so-called Traditional Model (TM). We give a set of necessary and sufficient conditions to decide whether a high accuracy algorithm exists in the TM, and describe progress toward a decision procedure that will take any problem and provide either a high accuracy algorithm or a proof that none exists. When no accurate algorithm exists in the TM, it is natural to extend the set of available accurate operations by a library of additional operations, such as $x+y+z$, dot products, or indeed any enumerable set which could then be used to build further accurate algorithms. We show how our accurate algorithms and decision procedure for finding them extend to this case. Finally, we address other models of arithmetic, and the relationship between (im)possibility in the TM and (in)efficient algorithms operating on numbers represented as bit strings.Comment: 49 pages, 6 figures, 1 tabl

Toward accurate polynomial evaluation in rounded arithmetic

Given a multivariate real (or complex) polynomial $p$ and a domain $\cal D$, we would like to decide whether an algorithm exists to evaluate $p(x)$ accurately for all $x \in {\cal D}$ using rounded real (or complex) arithmetic. Here ``accurately'' means with relative error less than 1, i.e., with some correct leading digits. The answer depends on the model of rounded arithmetic: We assume that for any arithmetic operator $op(a,b)$, for example $a+b$ or $a \cdot b$, its computed value is $op(a,b) \cdot (1 + \delta)$, where $| \delta |$ is bounded by some constant $\epsilon$ where $0 < \epsilon \ll 1$, but $\delta$ is otherwise arbitrary. This model is the traditional one used to analyze the accuracy of floating point algorithms.Our ultimate goal is to establish a decision procedure that, for any $p$ and $\cal D$, either exhibits an accurate algorithm or proves that none exists. In contrast to the case where numbers are stored and manipulated as finite bit strings (e.g., as floating point numbers or rational numbers) we show that some polynomials $p$ are impossible to evaluate accurately. The existence of an accurate algorithm will depend not just on $p$ and $\cal D$, but on which arithmetic operators and which constants are are available and whether branching is permitted. Toward this goal, we present necessary conditions on $p$ for it to be accurately evaluable on open real or complex domains ${\cal D}$. We also give sufficient conditions, and describe progress toward a complete decision procedure. We do present a complete decision procedure for homogeneous polynomials $p$ with integer coefficients, {\cal D} = \C^n, and using only the arithmetic operations $+$, $-$ and $\cdot$.Comment: 54 pages, 6 figures; refereed version; to appear in Foundations of Computational Mathematics: Santander 2005, Cambridge University Press, March 200

How Many Subpopulations is Too Many? Exponential Lower Bounds for Inferring Population Histories

Reconstruction of population histories is a central problem in population genetics. Existing coalescent-based methods, like the seminal work of Li and Durbin (Nature, 2011), attempt to solve this problem using sequence data but have no rigorous guarantees. Determining the amount of data needed to correctly reconstruct population histories is a major challenge. Using a variety of tools from information theory, the theory of extremal polynomials, and approximation theory, we prove new sharp information-theoretic lower bounds on the problem of reconstructing population structure -- the history of multiple subpopulations that merge, split and change sizes over time. Our lower bounds are exponential in the number of subpopulations, even when reconstructing recent histories. We demonstrate the sharpness of our lower bounds by providing algorithms for distinguishing and learning population histories with matching dependence on the number of subpopulations. Along the way and of independent interest, we essentially determine the optimal number of samples needed to learn an exponential mixture distribution information-theoretically, proving the upper bound by analyzing natural (and efficient) algorithms for this problem.Comment: 38 pages, Appeared in RECOMB 201

Reduced Order and Surrogate Models for Gravitational Waves

We present an introduction to some of the state of the art in reduced order and surrogate modeling in gravitational wave (GW) science. Approaches that we cover include Principal Component Analysis, Proper Orthogonal Decomposition, the Reduced Basis approach, the Empirical Interpolation Method, Reduced Order Quadratures, and Compressed Likelihood evaluations. We divide the review into three parts: representation/compression of known data, predictive models, and data analysis. The targeted audience is that one of practitioners in GW science, a field in which building predictive models and data analysis tools that are both accurate and fast to evaluate, especially when dealing with large amounts of data and intensive computations, are necessary yet can be challenging. As such, practical presentations and, sometimes, heuristic approaches are here preferred over rigor when the latter is not available. This review aims to be self-contained, within reasonable page limits, with little previous knowledge (at the undergraduate level) requirements in mathematics, scientific computing, and other disciplines. Emphasis is placed on optimality, as well as the curse of dimensionality and approaches that might have the promise of beating it. We also review most of the state of the art of GW surrogates. Some numerical algorithms, conditioning details, scalability, parallelization and other practical points are discussed. The approaches presented are to large extent non-intrusive and data-driven and can therefore be applicable to other disciplines. We close with open challenges in high dimension surrogates, which are not unique to GW science.Comment: Invited article for Living Reviews in Relativity. 93 page

Randomness as a computational strategy : on matrix and tensor decompositions

Matrix and tensor decompositions are fundamental tools for finding structure and data processing. In particular, the efficient computation of low-rank matrix approximations is an ubiquitous problem in the area of machine learning and elsewhere. However, massive data arrays pose a computational challenge for these techniques, placing significant constraints on both memory and processing power. Recently, the fascinating and powerful concept of randomness has been introduced as a strategy to ease the computational load of deterministic matrix and data algorithms. The basic idea of these algorithms is to employ a degree of randomness as part of the logic in order to derive from a high-dimensional input matrix a smaller matrix, which captures the essential information of the original data matrix. Subsequently, the smaller matrix is then used to efficiently compute a near-optimal low-rank approximation. Randomized algorithms have been shown to be robust, highly reliable, and computationally efficient, yet simple to implement. In particular, the development of the randomized singular value decomposition can be seen as a milestone in the era of ‘big data’. Building up on the great success of this probabilistic strategy to compute low-rank matrix decompositions, this thesis introduces a set of new randomized algorithms. Specifically, we present a randomized algorithm to compute the dynamic mode decomposition, which is a modern dimension reduction technique designed to extract dynamic information from dynamical systems. Then, we advocate the randomized dynamic mode decomposition for background modeling of surveillance video feeds. Further, we show that randomized algorithms are embarrassingly parallel by design and that graphics processing units (GPUs) can be utilized to substantially accelerate the computations. Finally, the concept of randomized algorithms is generalized for tensors in order to compute the canonical CANDECOMP/PARAFAC (CP) decomposition

A Tale of Two Hungarians: Tridiagonalizing Random Matrices

The Hungarian physicist Eugene Wigner introduced random matrix models in physics to describe the energy spectra of atomic nuclei. As such, the main goal of Random Matrix Theory (RMT) has been to derive the eigenvalue statistics of matrices drawn from a given distribution. The Wigner approach gives powerful insights into the properties of complex, chaotic systems in thermal equilibrium. Another Hungarian, Cornelius Lanczos, suggested a method of reducing the dynamics of any quantum system to a one-dimensional chain by tridiagonalizing the Hamiltonian relative to a given initial state. In the resulting matrix, the diagonal and off-diagonal Lanczos coefficients control transition amplitudes between elements of a distinguished basis of states. We connect these two approaches to the quantum mechanics of complex systems by deriving analytical formulae relating the potential defining a general RMT, or, equivalently, its density of states, to the Lanczos coefficients and their correlations. In particular, we derive an integral relation between the average Lanczos coefficients and the density of states, and, for polynomial potentials, algebraic equations that determine the Lanczos coefficients from the potential. We obtain these results for generic initial states in the thermodynamic limit. As an application, we compute the time-dependent ``spread complexity'' in Thermo-Field Double states and the spectral form factor for Gaussian and Non-Gaussian RMTs.Comment: 30 pages. References adde

Electronic correlations in inhomogeneous model systems: numerical simulation of spectra and transmission

Many fascinating features in condensed matter systems emerge due to the interaction between electrons. Magnetism is such a paramount consequence, which is explained in terms of the exchange interaction of electrons. Another prime example is the metal-to-Mott-insulator transition, where the energy cost of Coulomb repulsion competes against the kinetic energy, the latter favoring delocalization. While systems of correlated electrons are exciting and show remarkable and technologically promising physical properties, they are difficult to treat theoretically. A single-particle description is insufficient; the quantum many-body problem of interacting electrons has to be solved. In the present thesis, we study physical properties of half-metallic ferromagnets which are used in spintronic devices. Half-metals exhibit a metallic spin channel, while the other spin channel is insulating; they are characterized by a high spin polarization. This thesis contributes to the development of numerical methods and applies them to models of half-metallic ferromagnets. Throughout this work, the single-band Hubbard Hamiltonian is considered, and electronic correlations are treated within dynamical mean-field theory. Instead of directly solving the lattice model, the dynamical mean-field theory amounts to solving a local, effective impurity problem that is determined self-consistently. At finite temperatures, this impurity problem is solved employing continuous-time quantum Monte Carlo algorithms formulated in the action formalism. As these algorithms are formulated in imaginary time, an analytic continuation is required to obtain spectral functions. We formulate a version of the N-point Padé algorithm that calculates the location of the poles in a least-squares sense. To directly obtain spectra for real frequencies, we employ Hamiltonian-based tensor network methods at zero temperature. We also summarize the ideas of the density matrix renormalization group algorithm, and of the time evolution using the time-dependent variational principle, employing a diagrammatic notation. Real materials never display perfect translational symmetry. Thus, realistic models require the inclusion of disorder effects. In this work, we discuss these within a single-site approximation, the coherent potential approximation, and combine it with the dynamical mean-field theory, allowing to treat interacting electrons in multicomponent alloys on a local level. We extend this combined scheme to off-diagonal disorder, that is, disorder in the hopping amplitudes, by employing the Blackman–Esterling–Berk formalism. For this purpose, we illustrate the ideas of this formalism using tensor diagrams and provide an efficient implementation. The structure of the effective medium is discussed, and a concentration scaling is proposed that resolves some of its peculiarities. The limit of vanishing hopping between different components is discussed and solved analytically for the Bethe lattice with a general coordination number. We exemplify the combined algorithm for a Bethe lattice, showing results that exhibit alloy-band-insulator to correlated-metal to Mott-insulator transitions. We study models of half-metallic ferromagnets to elucidate the effects of local electronic correlations on the spectral function. To model half-metallicity, a static spin splitting is used to produce the half-metallic density of states. Applying the Padé analytic continuation to the self-energy instead of the Green’s function produces reliable spectral functions agreeing with the zero-temperature results obtained for real frequencies. To address transport properties, we investigate the interface of a half-metallic layer and a metallic, band insulating, or Mott insulating layer. We observe charge reconstruction which induces metallicity at the interface; quasiparticle states are present in the Mott insulating layer even for a large Hubbard interaction. The transmission through a barrier made of such a single interacting half-metallic layer sandwiched by metallic leads is studied employing the Meir–Wingreen formalism. This allows for a transparent calculation of the transmission in the presence of the Hubbard interaction. For a strong coupling of the central layer to the leads, we identify high intensity bound states which do not contribute to the transmission. For small coupling, on the other hand, we find resonant states which enhance the transmission. In particular, we demonstrate that even for a single half-metallic layer, highly polarized transmissions are achievable

Cleaning large correlation matrices: tools from random matrix theory

This review covers recent results concerning the estimation of large covariance matrices using tools from Random Matrix Theory (RMT). We introduce several RMT methods and analytical techniques, such as the Replica formalism and Free Probability, with an emphasis on the Marchenko-Pastur equation that provides information on the resolvent of multiplicatively corrupted noisy matrices. Special care is devoted to the statistics of the eigenvectors of the empirical correlation matrix, which turn out to be crucial for many applications. We show in particular how these results can be used to build consistent "Rotationally Invariant" estimators (RIE) for large correlation matrices when there is no prior on the structure of the underlying process. The last part of this review is dedicated to some real-world applications within financial markets as a case in point. We establish empirically the efficacy of the RIE framework, which is found to be superior in this case to all previously proposed methods. The case of additively (rather than multiplicatively) corrupted noisy matrices is also dealt with in a special Appendix. Several open problems and interesting technical developments are discussed throughout the paper.Comment: 165 pages, article submitted to Physics Report

Modified Fourier expansions: theory, construction and applications

Modified Fourier expansions present an alternative to more standard algorithms for the approximation of nonperiodic functions in bounded domains. This thesis addresses the theory of such expansions, their effective construction and computation, and their application to the numerical solution of partial differential equations. As the name indicates, modified Fourier expansions are closely related to classical Fourier series. The latter are naturally defined in the d-variate cube, and, in an analogous fashion, we primarily study modified Fourier expansions in this domain. However, whilst Fourier coefficients are commonly computed with the Fast Fourier Transform (FFT), we use modern numerical quadratures instead. In contrast to the FFT, such schemes are adaptive, leading to great potential savings in computational cost. Standard algorithms for the approximation of nonperiodic functions in $d$-variate cubes exhibit complexities that grow exponentially with dimension. The aforementioned quadratures permit the design of approximations based on modified Fourier expansions that do not possess this feature. Consequently, such schemes are increasingly effective in higher dimensions. When applied to the numerical solution of boundary value problems, such savings in computational cost impart benefits over more commonly used polynomial-based methods. Moreover, regardless of the dimensionality of the problem, modified Fourier methods lead to well-conditioned matrices and corresponding linear systems that can be solved cheaply with standard iterative techniques. The theoretical component of this thesis furnishes modified Fourier expansions with a convergence analysis in arbitrary dimensions. In particular, we prove uniform convergence of modified Fourier expansions under rather general conditions. Furthermore, it is known that the notion of modified Fourier expansions can be effectively generalised, resulting in a family of approximation bases sharing many of the features of the modified Fourier case. The purpose of such a generalisation is to obtain both faster rates and higher degrees of convergence. Having detailed the approximation-theoretic properties of modified Fourier expansions, we extend this analysis to the general case and thereby verify this improvement. A central drawback of these expansions is that their convergence rate is both fixed and typically slow. This makes the construction of effective convergence acceleration techniques imperative. In the final part of this thesis, we design and analyse a robust method, applicable in arbitrary numbers of dimensions, for accelerating convergence of modified Fourier expansions. When employed in the approximation of multivariate functions, this culminates in efficient, high-order approximants comprising relatively small numbers of terms