Showing papers in "Siam Review in 2016"

Computational Methods for Linear Matrix Equations

Abstract: Given the square matrices $A, B, D, E$ and the matrix $C$ of conforming dimensions, we consider the linear matrix equation $A{\mathbf X} E+D{\mathbf X} B = C$ in the unknown matrix ${\mathbf 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.

Approximating Spectral Densities of Large Matrices

Abstract: In physics, it is sometimes desirable to compute the so-called density of states (DOS), also known as the spectral density, of a real symmetric matrix $A$. The spectral density can be viewed as a probability density distribution that measures the likelihood of finding eigenvalues near some point on the real line. The most straightforward way to obtain this density is to compute all eigenvalues of $A$, but this approach is generally costly and wasteful, especially for matrices of large dimension. There exist alternative methods that allow us to estimate the spectral density function at much lower cost. The major computational cost of these methods is in multiplying $A$ with a number of vectors, which makes them appealing for large-scale problems where products of the matrix $A$ with arbitrary vectors are relatively inexpensive. This article defines the problem of estimating the spectral density carefully and discusses how to measure the accuracy of an approximate spectral density. It then surveys a few kno...

Diffuse Interface Models on Graphs for Classification of High Dimensional Data

Abstract: This paper is a republication of an MMS paper [A. L. Bertozzi and A. Flenner, Multiscale Model. Simul., 10 (2012), pp. 1090--1118] describing a new class of algorithms for classification of high dimensional data. These methods combine ideas from spectral methods on graphs with nonlinear edge/region detection methods traditionally used in the PDE-based imaging community. The algorithms use the Ginzburg--Landau functional, extended to a graphical framework, which has classical PDE connections to total variation minimization. Convex splitting algorithms allow us to quickly find minimizers of the proposed model and take advantage of fast spectral solvers of linear graph-theoretic problems. We review the diverse computational examples presented in the original paper, involving both basic clustering and semisupervised learning for different applications. Case studies include feature identification in images, segmentation in social networks, and segmentation of shapes in high dimensional datasets. Since the pape...

Symplectic Runge--Kutta Schemes for Adjoint Equations, Automatic Differentiation, Optimal Control, and More

Abstract: The study of the sensitivity of the solution of a system of differential equations with respect to changes in the initial conditions leads to the introduction of an adjoint system, whose discretization is related to reverse accumulation in automatic differentiation. Similar adjoint systems arise in optimal control and other areas, including classical mechanics. Adjoint systems are introduced in such a way that they exactly preserve a relevant quadratic invariant (more precisely, an inner product). Symplectic Runge--Kutta and partitioned Runge--Kutta methods are defined through the exact conservation of a differential geometric structure, but may be characterized by the fact that they preserve exactly quadratic invariants of the system being integrated. Therefore, the symplecticness (or lack of symplecticness) of a Runge--Kutta or partitioned Runge--Kutta integrator should be relevant to understanding its performance when applied to the computation of sensitivities, to optimal control problems, and in othe...

Computing Fundamental Matrix Decompositions Accurately via the Matrix Sign Function in Two Iterations: The Power of Zolotarev's Functions

Abstract: The symmetric eigenvalue decomposition and the singular value decomposition (SVD) are fundamental matrix decompositions with many applications. Conventional algorithms for computing these decomposi...

Stochastic Models of Evidence Accumulation in Changing Environments

Abstract: Organisms and ecological groups accumulate evidence to make decisions. Classic experiments and theoretical studies have explored this process when the correct choice is fixed during each trial. However, we live in a constantly changing world. What effect does such impermanence have on classical results about decision making? To address this question we use sequential analysis to derive a tractable model of evidence accumulation when the correct option changes in time. Our analysis shows that ideal observers discount prior evidence at a rate determined by the volatility of the environment, and the dynamics of evidence accumulation is governed by the information gained over an average environmental epoch. A plausible neural implementation of an optimal observer in a changing environment shows that, in contrast to previous models, neural populations representing alternate choices are coupled through excitation. Our work builds a bridge between statistical decision making in volatile environments and stochast...

Canards, folded nodes and mixed-mode oscillations in piecewise-linear slow-fast systems

Abstract: Canard-induced phenomena have been extensively studied in the last three decades, from both the mathematical and the application viewpoints. Canards in slow-fast systems with (at least) two slow variables, especially near folded-node singularities, give an essential generating mechanism for mixed-mode oscillations (MMOs) in the framework of smooth multiple timescale systems. There is a wealth of literature on such slow-fast dynamical systems and many models displaying canard-induced MMOs, particularly in neuroscience. In parallel, since the late 1990s several papers have shown that the canard phenomenon can be faithfully reproduced with piecewise-linear (PWL) systems in two dimensions, although very few results are available in the three-dimensional case. The present paper aims to bridge this gap by analyzing canonical PWL systems that display folded singularities, primary and secondary canards, with a similar control of the maximal winding number as in the smooth case. We also show that the singular phas...

Border-Collision Bifurcations in $\mathbb{R}^N$

Abstract: For piecewise-smooth maps, new dynamics can be created by varying parameters such that a fixed point collides with a surface on which the map is nonsmooth. If the map is continuous and has bounded derivatives near but not on the surface, this is referred to as a border-collision bifurcation. Dynamics near border-collision bifurcations are well approximated by piecewise-linear continuous maps. A lack of differentiability allows for substantial complexity and a wide variety of invariant sets can be created in border-collision bifurcations, such as invariant circles and chaotic sets, and several attractors may be created simultaneously. Yet many calculations that would be intractable for smooth nonlinear maps can be performed exactly for piecewise-linear maps and in recent years several new results have been obtained for border-collision bifurcations. This article reviews border-collision bifurcations with a general emphasis on results that apply to maps of any number of dimensions. The paper covers applicat...

Near-optimal perfectly matched layers for indefinite Helmholtz problems

Abstract: A new construction of an absorbing boundary condition for indefinite Helmholtz problems on unbounded domains is presented. This construction is based on a near-best uniform rational interpolant of the inverse square root function on the union of a negative and positive real interval, designed with the help of a classical result by Zolotarev. Using Krein's interpretation of a Stieltjes continued fraction, this interpolant can be converted into a three-term finite difference discretization of a perfectly matched layer (PML) which converges exponentially fast in the number of grid points. The convergence rate is asymptotically optimal for both propagative and evanescent wave modes. Several numerical experiments and illustrations are included.

Geodesic Transport Barriers in Jupiter's Atmosphere: A Video-Based Analysis

Abstract: Jupiter's zonal jets and Great Red Spot are well known from still images. However, the planet's atmosphere is highly unsteady, which suggests that the actual material transport barriers delineating its main features should be time-dependent. Rare video footage of Jupiter's clouds provides an opportunity to verify this expectation from optically reconstructed velocity fields. Available videos, however, provide short-time and temporally aperiodic velocity fields that defy classical dynamical systems analyses focused on asymptotic features. To this end, we use here the recent theory of geodesic transport barriers to uncover finite-time mixing barriers in the wind field extracted from a video captured by NASA's Cassini space mission. More broadly, the approach described here provides a systematic and frame-invariant way to extract dynamic coherent structures from time-resolved remote observations of unsteady continua.

Enumerating Rigid Sphere Packings

Abstract: Packing problems, which ask how to arrange a collection of objects in space to meet certain criteria, are important in a great many physical and biological systems, where geometrical arrangements at small scales control behavior at larger ones. In many systems there is no single, optimal packing that dominates, but rather one must understand the entire set of possible packings. As a step in this direction we enumerate rigid clusters of identical hard spheres for $n\leq 14$ and clusters with the maximum number of contacts for $n\leq 19$. A rigid cluster is one that cannot be continuously deformed while maintaining all contacts. This is a nonlinear notion that arises naturally because such clusters are the metastable states when the spheres interact with a short-range potential, as is the case in many nano- or microscale systems. We believe that our lists are nearly complete, except for a small number of highly singular clusters (linearly floppy but nonlinearly rigid). The data contains some major geometric...

Computationally Tractable Counterparts of Distributionally Robust Constraints on Risk Measures

Abstract: In optimization problems appearing in fields such as economics, finance, or engineering, it is often important that a risk measure of a decision-dependent random variable stays below a prescribed level. At the same time, the underlying probability distribution determining the risk measure's value is typically known only up to a certain degree and the constraint should hold for a reasonably wide class of probability distributions. In addition, the constraint should be computationally tractable. In this paper we review and generalize results on the derivation of tractable counterparts of such constraints for discrete probability distributions. Using established techniques in robust optimization, we show that the derivation of a tractable robust counterpart can be split into two parts, one corresponding to the risk measure and the other to the uncertainty set. This holds for a wide range of risk measures and uncertainty sets for probability distributions defined using statistical goodness-of-fit tests or pro...

Pairwise Compatibility Graphs: A Survey

Abstract: A graph $G=(V,E)$ is a pairwise compatibility graph (PCG) if there exists an edge-weighted tree $T$ and two nonnegative real numbers $d_{min}$ and $d_{max}$ such that each leaf $u$ of $T$ is a node of $V$ and there is an edge $(u,v) \in E$ if and only if $d_{min} \leq d_T (u, v) \leq d_{max}$, where $d_T (u, v)$ is the sum of weights of the edges on the unique path from $u$ to $v$ in $T$. In this article, we survey the state of the art concerning this class of graphs and some of its subclasses.

An evolutionary model of tumor cell kinetics and the emergence of molecular heterogeneity driving gompertzian growth.

Abstract: We describe a cell-molecular based evolutionary mathematical model of tumor development driven by a stochastic Moran birth-death process. The cells in the tumor carry molecular information in the form of a numerical genome which we represent as a four-digit binary string used to differentiate cells into 16 molecular types. The binary string is able to undergo stochastic point mutations that are passed to a daughter cell after each birth event. The value of the binary string determines the cell fitness, with lower fit cells (e.g. 0000) defined as healthy phenotypes, and higher fit cells (e.g. 1111) defined as malignant phenotypes. At each step of the birth-death process, the two phenotypic sub-populations compete in a prisoner's dilemma evolutionary game with the healthy cells playing the role of cooperators, and the cancer cells playing the role of defectors. Fitness, birth-death rates of the cell populations, and overall tumor fitness are defined via the prisoner's dilemma payoff matrix. Mutation parameters include passenger mutations (mutations conferring no fitness advantage) and driver mutations (mutations which increase cell fitness). The model is used to explore key emergent features associated with tumor development, including tumor growth rates as it relates to intratumor molecular heterogeneity. The tumor growth equation states that the growth rate is proportional to the logarithm of cellular diversity/heterogeneity. The Shannon entropy from information theory is used as a quantitative measure of heterogeneity and tumor complexity based on the distribution of the 4-digit binary sequences produced by the cell population. To track the development of heterogeneity from an initial population of healthy cells (0000), we use dynamic phylogenetic trees which show clonal and sub-clonal expansions of cancer cell sub-populations from an initial malignant cell. We show tumor growth rates are not constant throughout tumor development, and are generally much higher in the subclinical range than in later stages of development, which leads to a Gompertzian growth curve. We explain the early exponential growth of the tumor and the later saturation associated with the Gompertzian curve which results from our evolutionary simulations using simple statistical mechanics principles related to the degree of functional coupling of the cell states. We then compare dosing strategies at early stage development, mid-stage (clinical stage), and late stage development of the tumor. If used early during tumor development in the subclinical stage, well before the cancer cell population is selected for growth, therapy is most effective at disrupting key emergent features of tumor development.

Communicability Angle and the Spatial Efficiency of Networks

Abstract: We introduce the concept of communicability angle between a pair of nodes in a graph. We provide strong analytical and empirical evidence that the average communicability angle for a given network accounts for its spatial efficiency on the basis of the communications among the nodes in a network. We determine characteristics of the spatial efficiency of more than a hundred real-world networks that represent complex systems arising in a diverse set of scenarios. In particular, we find that the communicability angle correlates very well with the experimentally measured value of the relative packing efficiency of proteins that are represented as residue networks. We finally show how we can modulate the spatial efficiency of a network by tuning the weights of the edges of the networks. This allows us to predict the effects of external stresses on the spatial efficiency of a network as well as to design strategies to improve important parameters in real-world systems.

Optimal Mixing Enhancement by Local Perturbation

Abstract: We consider the problem of how to apply local perturbations to optimally enhance the mixing of a (possibly time-dependent) dynamical system. We develop a flexible modeling approach based on the transfer operator of the dynamical system, and pose the problem in the language of convex optimization. The optimal local perturbations can then be efficiently computed, at discrete time instants, by standard convex optimization techniques. The local perturbations satisfy physical constraints, such as preservation of the invariant measure of the dynamics (for example, for incompressible fluid flow, the perturbations preserve volume), and a variety of other physical constraints can also be easily enforced. We show that one can achieve surprisingly large speedups in mixing via optimizing the diffusion, as compared to fixed diffusion protocols. Finally, we indicate how one might alternatively try to use local perturbation to push a mass density toward a particular region of the domain.

Optimization in High Dimensions via Accelerated, Parallel, and Proximal Coordinate Descent

Abstract: We propose a new randomized coordinate descent method for minimizing the sum of convex functions, each of which depends on a small number of coordinates only. Our method (APPROX) is simultaneously Accelerated, Parallel, and PROXimal; this is the first time such a method has been proposed. In the special case when the number of processors is equal to the number of coordinates, the method converges at the rate $2\bar{\omega}\bar{L} R^2/(k+1)^2 $, where $k$ is the iteration counter, $\bar{\omega}$ is a data-weighted average degree of separability of the loss function, $\bar{L}$ is the average of Lipschitz constants associated with the coordinates and individual functions in the sum, and $R$ is the distance of the initial point from the minimizer. We show that the method can be implemented without the need to perform full-dimensional vector operations, which is the major bottleneck of accelerated coordinate descent, rendering it impractical. The fact that the method depends on the average degree of separabili...

Restoring Definiteness via Shrinking, with an Application to Correlation Matrices with a Fixed Block

Abstract: Indefinite approximations of positive semidefinite matrices arise in various data analysis applications involving covariance matrices and correlation matrices. We propose a method for restoring positive semidefiniteness of an indefinite matrix $M_0$ that constructs a convex linear combination $S(\alpha) = \alpha M_1 + (1-\alpha)M_0$ of $M_0$ and a positive semidefinite target matrix $M_1$. In statistics, this construction for improving an estimate $M_0$ by combining it with new information in $M_1$ is known as shrinking. We make no statistical assumptions about $M_0$ and define the optimal shrinking parameter as $\alpha_* = \min \{\alpha \in [0,1] : {$S(\alpha)$ is positive semidefinite}\}$. We describe three \alg s for computing $\alpha_*$. One algorithm is based on the bisection method, with the use of Cholesky factorization to test definiteness; a second employs Newton's method; and a third finds the smallest eigenvalue of a symmetric definite generalized eigenvalue problem. We show that weights that r...

Hyperspherical Sparse Approximation Techniques for High-Dimensional Discontinuity Detection

Abstract: This work proposes a hyperspherical sparse approximation framework for detecting jump discontinuities in functions in high-dimensional spaces. The need for a novel approach results from the theoretical and computational inefficiencies of well-known approaches, such as adaptive sparse grids, for discontinuity detection. Our approach constructs the hyperspherical coordinate representation of the discontinuity surface of a function. Then sparse approximations of the transformed function are built in the hyperspherical coordinate system, with values at each point estimated by solving a one-dimensional discontinuity detection problem. Due to the smoothness of the hypersurface, the new technique can identify jump discontinuities with significantly reduced computational cost, compared to existing methods. Several approaches are used to approximate the transformed discontinuity surface in the hyperspherical system, including adaptive sparse grid and radial basis function interpolation, discrete least squares proj...

Modeling, Simulating, and Parameter Fitting of Biochemical Kinetic Experiments

Abstract: In many chemical and biological applications, systems of differential equations containing unknown parameters are used to explain empirical observations and experimental data. The differential equations are typically nonlinear and difficult to analyze, requiring numerical methods to approximate the solutions. Compounding this difficulty are the unknown parameters in the differential equation system, which must be given specific numerical values in order for simulations to be run. Estrogen receptor protein dimerization is used as an example to demonstrate model construction, reduction, simulation, and parameter estimation. Mathematical, computational, and statistical methods are applied to empirical data to deduce kinetic parameter estimates and guide decisions regarding future experiments and modeling. The process demonstrated serves as a pedagogical example of quantitative methods being used to extract parameter values from biochemical data models.

Variations on a Theme of Euler

Abstract: This paper is an exploration of numerical methods for solving initial-value problems for ordinary differential equations. We look in detail at one “simple” problem, namely, the leaky bucket, and use it to explore the notions of explicit marching methods vs. implicit marching methods, interpolation, backward error, and stiffness. While the leaky bucket example is very well known, and these topics are studied in a great many textbooks, we will here emphasize backward error in a way that might be new and that we hope will be useful for your students. Indeed, the paper is intended to be a resource for students themselves. We will also use two techniques not normally seen in a first course, a new one, namely, “optimal backward error,” and an old one, namely, “the method of modified equations.”

A Probabilistic Analysis of the Game of the Goose

Abstract: We analyze the traditional board game the Game of the Goose. We are particularly interested in the probabilities of the different players winning, and we show that we can determine these probabilities exactly for up to six players and using simulation for any number of players. Our original motivation to investigate this game came from progress in stochastic process theories, which prompted the question of whether such methods are capable of dealing with well-known probabilistic games. As these games have large state spaces, this is not trivial. As a side effect we find that some common wisdom about the game is not true.

A Newton--Galerkin Method for Fluid Flow Exhibiting Uncertain Periodic Dynamics

Abstract: The determination of limit-cycles plays an important role in characterizing complex dynamical systems, such as unsteady fluid flows. In practice, dynamical systems are described by models equations involving parameters which are seldom exactly known, leading to parametric uncertainties. These parameters can be suitably modeled as random variables, so if the system possesses almost surely a stable time periodic solution, limit-cycles become stochastic, too. This paper introduces a novel numerical method for the computation of stable stochastic limit-cycles based on the spectral stochastic finite element method with polynomial chaos (PC) expansions. The method is designed to overcome the limitation of PC expansions related to convergence breakdown for long term integration. First, a stochastic time scaling of the model equations is determined to control the phase-drift of the stochastic trajectories and allowing for accurate low order PC expansions. Second, using the rescaled governing equations, we aim at ...

Analysis of methods used to reconstruct the flight path of Malaysia Airlines flight 370

Abstract: We analyze methods similar to those used by a government-appointed team of engineers to reconstruct possible flight paths of Malaysia Airlines flight MH370, which disappeared in the southern Indian Ocean on March 8th, 2014. These methods make use of data obtained from a series of messages that were exchanged between a satellite communications system and the aircraft over a six hour period after radar contact was lost. We develop a series of three progressively more realistic models to reconstruct the flight path from the satellite data. We validate these models using data from simulated flight paths, and by comparing the results we obtained for MH370 to those found by the government-appointed team of engineers. Our analysis helps to explain how the significant uncertainty in the data has made it difficult to narrow down the priority search area for the aircraft. The exposition is aimed at students with a background in vector calculus, matrix analysis, and numerical analysis. A range of modeling, analytica...