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

Power-law distributions occur in many situations of scientific interest and have significant consequences for our understanding of natural and man-made phenomena. Unfortunately, the detection and characterization of power laws is complicated by the large fluctuations that occur in the tail of the distribution—the part of the distribution representing large but rare events—and by the difficulty of identifying the range over which power-law behavior holds. Commonly used methods for analyzing power-law data, such as least-squares fitting, can produce substantially inaccurate estimates of parameters for power-law distributions, and even in cases where such methods return accurate answers they are still unsatisfactory because they give no indication of whether the data obey a power law at all. Here we present a principled statistical framework for discerning and quantifying power-law behavior in empirical data. Our approach combines maximum-likelihood fitting methods with goodness-of-fit tests based on the Kolmogorov—Smirnov (KS) statistic and likelihood ratios. We evaluate the effectiveness of the approach with tests on synthetic data and give critical comparisons to previous approaches. We also apply the proposed methods to twenty-four real-world data sets from a range of different disciplines, each of which has been conjectured to follow a power-law distribution. In some cases we find these conjectures to be consistent with the data, while in others the power law is ruled out.

Low-rank matrix approximations, such as the truncated singular value decomposition and the rank-revealing QR decomposition, play a central role in data analysis and scientific computing. This work surveys and extends recent research which demonstrates that randomization offers a powerful tool for performing low-rank matrix approximation. These techniques exploit modern computational architectures more fully than classical methods and open the possibility of dealing with truly massive data sets. This paper presents a modular framework for constructing randomized algorithms that compute partial matrix decompositions. These methods use random sampling to identify a subspace that captures most of the action of a matrix. The input matrix is then compressed—either explicitly or implicitly—to this subspace, and the reduced matrix is manipulated deterministically to obtain the desired low-rank factorization. In many cases, this approach beats its classical competitors in terms of accuracy, robustness, and/or speed. These claims are supported by extensive numerical experiments and a detailed error analysis. The specific benefits of randomized techniques depend on the computational environment. Consider the model problem of finding the k dominant components of the singular value decomposition of an m × n matrix. (i) For a dense input matrix, randomized algorithms require O(mn log(k)) floating-point operations (flops) in contrast to O(mnk) for classical algorithms. (ii) For a sparse input matrix, the flop count matches classical Krylov subspace methods, but the randomized approach is more robust and can easily be reorganized to exploit multi-processor architectures. (iii) For a matrix that is too large to fit in fast memory, the randomized techniques require only a constant number of passes over the data, as opposed to O(k) passes for classical algorithms. In fact, it is sometimes possible to perform matrix approximation with a single pass over the data.

Numerical simulation of large-scale dynamical systems plays a fundamental role in studying a wide range of complex physical phenomena; however, the inherent large-scale nature of the models often leads to unmanageable demands on computational resources. Model reduction aims to reduce this computational burden by generating reduced models that are faster and cheaper to simulate, yet accurately represent the original large-scale system behavior. Model reduction of linear, nonparametric dynamical systems has reached a considerable level of maturity, as reflected by several survey papers and books. However, parametric model reduction has emerged only more recently as an important and vibrant research area, with several recent advances making a survey paper timely. Thus, this paper aims to provide a resource that draws together recent contributions in different communities to survey the state of the art in parametric model reduction methods. Parametric model reduction targets the broad class of problems for which the equations governing the system behavior depend on a set of parameters. Examples include parameterized partial differential equations and large-scale systems of parameterized ordinary differential equations. The goal of parametric model reduction is to generate low-cost but accurate models that characterize system response for different values of the parameters. This paper surveys state-of-the-art methods in projection-based parametric model reduction, describing the different approaches within each class of methods for handling parametric variation and providing a comparative discussion that lends insights to potential advantages and disadvantages in applying each of the methods. We highlight the important role played by parametric model reduction in design, control, optimization, and uncertainty quantification-settings that require repeated model evaluations over different parameter values.

The affine rank minimization problem consists of finding a matrix of minimum rank that satisfies a given system of linear equality constraints. Such problems have appeared in the literature of a diverse set of fields including system identification and control, Euclidean embedding, and collaborative filtering. Although specific instances can often be solved with specialized algorithms, the general affine rank minimization problem is NP-hard because it contains vector cardinality minimization as a special case. In this paper, we show that if a certain restricted isometry property holds for the linear transformation defining the constraints, the minimum-rank solution can be recovered by solving a convex optimization problem, namely, the minimization of the nuclear norm over the given affine space. We present several random ensembles of equations where the restricted isometry property holds with overwhelming probability, provided the codimension of the subspace is sufficiently large. The techniques used in our analysis have strong parallels in the compressed sensing framework. We discuss how affine rank minimization generalizes this preexisting concept and outline a dictionary relating concepts from cardinality minimization to those of rank minimization. We also discuss several algorithmic approaches to minimizing the nuclear norm and illustrate our results with numerical examples.

In this paper we survey the primary research, both theoretical and applied, in the area of robust optimization (RO). Our focus is on the computational attractiveness of RO approaches, as well as the modeling power and broad applicability of the methodology. In addition to surveying prominent theoretical results of RO, we also present some recent results linking RO to adaptable models for multistage decision-making problems. Finally, we highlight applications of RO across a wide spectrum of domains, including finance, statistics, learning, and various areas of engineering.

This survey provides an overview of higher-order tensor decompositions, their applications, and available software. A tensor is a multidimensional or N-way array. Decompositions of higher-order tensors (i.e., N-way arrays with N ≥ 3) have applications in psychometrics, chemometrics, signal processing, numerical linear algebra, computer vision, numerical analysis, data mining, neuroscience, graph analysis, and elsewhere. Two particular tensor decompositions can be considered to be higher-order extensions of the matrix singular value decomposition: CANDECOMP/PARAFAC (CP) decomposes a tensor as a sum of rank-one tensors, and the Tucker decomposition is a higher-order form of principal component analysis. There are many other tensor decompositions, including INDSCAL, PARAFAC2, CANDELINC, DEDICOM, and PARATUCK2 as well as nonnegative variants of all of the above. The N-way Toolbox, Tensor Toolbox, and Multilinear Engine are examples of software packages for working with tensors.

A full-rank matrix ${\bf A}\in {\Bbb R}^{n\times m}$ with n < m generates an underdetermined system of linear equations Ax = b having infinitely many solutions. Suppose we seek the sparsest solution, i.e., the one with the fewest nonzero entries. Can it ever be unique? If so, when? As optimization of sparsity is combinatorial in nature, are there efficient methods for finding the sparsest solution? These questions have been answered positively and constructively in recent years, exposing a wide variety of surprising phenomena, in particular the existence of easily verifiable conditions under which optimally sparse solutions can be found by concrete, effective computational methods. Such theoretical results inspire a bold perspective on some important practical problems in signal and image processing. Several well-known signal and image processing problems can be cast as demanding solutions of undetermined systems of equations. Such problems have previously seemed, to many, intractable, but there is considerable evidence that these problems often have sparse solutions. Hence, advances in finding sparse solutions to undetermined systems have energized research on such signal and image processing problemsto striking effect. In this paper we review the theoretical results on sparse solutions of linear systems, empirical results on sparse modeling of signals and images, and recent applications in inverse problems and compression in image processing. This work lies at the intersection of signal processing and applied mathematics, and arose initially from the wavelets and harmonic analysis research communities. The aim of this paper is to introduce a few key notions and applications connected to sparsity, targeting newcomers interested in either the mathematical aspects of this area or its applications.

This paper develops a novel framework for phase retrieval, a problem which arises in X-ray crystallography, diffraction imaging, astronomical imaging, and many other applications. Our approach, called PhaseLift, combines multiple structured illuminations together with ideas from convex programming to recover the phase from intensity measurements, typically from the modulus of the diffracted wave. We demonstrate empirically that a complex-valued object can be recovered from the knowledge of the magnitude of just a few diffracted patterns by solving a simple convex optimization problem inspired by the recent literature on matrix completion. More importantly, we also demonstrate that our noise-aware algorithms are stable in the sense that the reconstruction degrades gracefully as the signal-to-noise ratio decreases. Finally, we introduce some theory showing that one can design very simple structured illumination patterns such that three diffracted figures uniquely determine the phase of the object we wish to recover.

A recently developed nonlocal vector calculus is exploited to provide a variational analysis for a general class of nonlocal diffusion problems described by a linear integral equation on bounded domains in ℝn. The nonlocal vector calculus also enables striking analogies to be drawn between the nonlocal model and classical models for diffusion, including a notion of nonlocal flux. The ubiquity of the nonlocal operator in applications is illustrated by a number of examples ranging from continuum mechanics to graph theory. In particular, it is shown that fractional Laplacian and fractional derivative models for anomalous diffusion are special cases of the nonlocal model for diffusion that we consider. The numerous applications elucidate different interpretations of the operator and the associated governing equations. For example, a probabilistic perspective explains that the nonlocal spatial operator appearing in our model corresponds to the infinitesimal generator for a symmetric jump process. Sufficient conditions on the kernel of the nonlocal operator and the notion of volume constraints are shown to lead to a well-posed problem. Volume constraints are a proxy for boundary conditions that may not be defined for a given kernel. In particular, we demonstrate for a general class of kernels that the nonlocal operator is a mapping between a volume constrained subspace of a fractional Sobolev subspace and its dual. We also demonstrate for other particular kernels that the inverse of the operator does not smooth but does correspond to diffusion. The impact of our results is that both a continuum analysis and a numerical method for the modeling of anomalous diffusion on bounded domains in ℝn are provided. The analytical framework allows us to consider finite-dimensional approximations using discontinuous and continuous Galerkin methods, both of which are conforming for the nonlocal diffusion equation we consider; error and condition number estimates are derived.

When uncontrollable resources fluctuate, optimal power flow (OPF), routinely used by the electric power industry to redispatch hourly controllable generation (coal, gas, and hydro plants) over control areas of transmission networks, can result in grid instability and, potentially, cascading outages. This risk arises because OPF dispatch is computed without awareness of major uncertainty, in particular fluctuations in renewable output. As a result, grid operation under OPF with renewable variability can lead to frequent conditions where power line flow ratings are significantly exceeded. Such a condition, which is borne by our simulations of real grids, is considered undesirable in power engineering practice. Possibly, it can lead to a risky outcome that compromises grid stability-line tripping. Smart grid goals include a commitment to large penetration of highly fluctuating renewables, thus calling to reconsider current practices, in particular the use of standard OPF. Our chance-constrained (CC) OPF corrects the problem and mitigates dangerous renewable fluctuations with minimal changes in the current operational procedure. Assuming availability of a reliable wind forecast parameterizing the distribution function of the uncertain generation, our CC-OPF satisfies all the constraints with high probability while simultaneously minimizing the cost of economic redispatch. CC-OPF allows efficient implementation, e. g., solving a typical instance over the 2746-bus Polish network in 20 seconds on a standard laptop.

This paper develops a novel framework for phase retrieval, a problem which arises in X-ray crystallography, diffraction imaging, astronomical imaging, and many other applications. Our approach, called PhaseLift, combines multiple structured illuminations together with ideas from convex programming to recover the phase from intensity measurements, typically from the modulus of the diffracted wave. We demonstrate empirically that a complex-valued object can be recovered from the knowledge of the magnitude of just a few diffracted patterns by solving a simple convex optimization problem inspired by the recent literature on matrix completion. More importantly, we also demonstrate that our noise-aware algorithms are stable in the sense that the reconstruction degrades gracefully as the signal-to-noise ratio decreases. Finally, we introduce some theory showing that one can design very simple structured illumination patterns such that three diffracted figures uniquely determine the phase of the object we wish to recover.

We review a general class of models for self-organized dynamics based on alignment. The dynamics of such systems is governed solely by interactions among individuals or "agents," with the tendency to adjust to their "environmental averages." This, in turn, leads to the formation of clusters, e. g., colonies of ants, flocks of birds, parties of people, rendezvous in mobile networks, etc. A natural question which arises in this context is to ask when and how clusters emerge through the self-alignment of agents, and what types of "rules of engagement" influence the formation of such clusters. Of particular interest to us are cases in which the self-organized behavior tends to concentrate into one cluster, reflecting a consensus of opinions, flocking of birds, fish, or cells, rendezvous of mobile agents, and, in general, concentration of other traits intrinsic to the dynamics. Many standard models for self-organized dynamics in social, biological, and physical sciences assume that the intensity of alignment increases as agents get closer, reflecting a common tendency to align with those who think or act alike. Moreover, "similarity breeds connection" reflects our intuition that increasing the intensity of alignment as the difference of positions decreases is more likely to lead to a consensus. We argue here that the converse is true: when the dynamics is driven by local interactions, it is more likely to approach a consensus when the interactions among agents increase as a function of their difference in position. Heterophily, the tendency to bond more with those who are different rather than with those who are similar, plays a decisive role in the process of clustering. We point out that the number of clusters in heterophilious dynamics decreases as the heterophily dependence among agents increases. In particular, sufficiently strong heterophilious interactions enhance consensus.

The study of the stability properties of switched and hybrid systems gives rise to a number of interesting and challenging mathematical problems. The objective of this paper is to outline some of these problems, to review progress made in solving them in a number of diverse communities, and to review some problems that remain open. An important contribution of our work is to bring together material from several areas of research and to present results in a unified manner. We begin our review by relating the stability problem for switched linear systems and a class of linear differential inclusions. Closely related to the concept of stability are the notions of exponential growth rates and converse Lyapunov theorems, both of which are discussed in detail. In particular, results on common quadratic Lyapunov functions and piecewise linear Lyapunov functions are presented, as they represent constructive methods for proving stability and also represent problems in which significant progress has been made. We also comment on the inherent difficulty in determining stability of switched systems in general, which is exemplified by NP-hardness and undecidability results. We then proceed by considering the stability of switched systems in which there are constraints on the switching rules, through both dwell-time requirements and state-dependent switching laws. Also in this case the theory of Lyapunov functions and the existence of converse theorems are reviewed. We briefly comment on the classical Lur'e problem and on the theory of stability radii, both of which contain many of the features of switched systems and are rich sources of practical results on the topic. Finally we present a list of questions and open problems which provide motivation for continued research in this area.

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

Google's PageRank method was developed to evaluate the importance of web-pages via their link structure. The mathematics of PageRank, however, are entirely general and apply to any graph or network in any domain. Thus, PageRank is now regularly used in bibliometrics, social and information network analysis, and for link prediction and recommendation. It's even used for systems analysis of road networks, as well as biology, chemistry, neuroscience, and physics. We'll see the mathematics and ideas that unite these diverse applications.

Mixed-mode oscillations (MMOs) are trajectories of a dynamical system in which there is an alternation between oscillations of distinct large and small amplitudes. MMOs have been observed and studied for over thirty years in chemical, physical, and biological systems. Few attempts have been made thus far to classify different patterns of MMOs, in contrast to the classification of the related phenomena of bursting oscillations. This paper gives a survey of different types of MMOs, concentrating its analysis on MMOs whose small-amplitude oscillations are produced by a local, multiple-time-scale "mechanism." Recent work gives substantially improved insight into the mathematical properties of these mechanisms. In this survey, we unify diverse observations about MMOs and establish a systematic framework for studying their properties. Numerical methods for computing different types of invariant manifolds and their intersections are an important aspect of the analysis described in this paper.

It is well known that the trapezoidal rule converges geometrically when applied to analytic functions on periodic intervals or the real line. The mathematics and history of this phenomenon are reviewed, and it is shown that far from being a curiosity, it is linked with computational methods all across scientific computing, including algorithms related to inverse Laplace transforms, special functions, complex analysis, rational approximation, integral equations, and the computation of functions and eigenvalues of matrices and operators.

This paper presents a review and critical analysis of the mathematical literature concerning the modeling of vehicular traffic and crowd phenomena. The survey of models deals with the representation scales and the mathematical frameworks that are used for the modeling approach. The paper also considers the challenging objective of modeling complex systems consisting of large systems of individuals interacting in a nonlinear manner, where one of the modeling difficulties is the fact that these systems are difficult to model at a global level when based only on the description of the dynamics of individual elements. The review is concluded with a critical analysis focused on research perspectives that consider the development of a unified modeling strategy.

High order accurate weighted essentially nonoscillatory (WENO) schemes are relatively new but have gained rapid popularity in numerical solutions of hyperbolic partial differential equations (PDEs) and other convection dominated problems. The main advantage of such schemes is their capability to achieve arbitrarily high order formal accuracy in smooth regions while maintaining stable, nonoscillatory, and sharp discontinuity transitions. The schemes are thus especially suitable for problems containing both strong discontinuities and complex smooth solution features. WENO schemes are robust and do not require the user to tune parameters. At the heart of the WENO schemes is actually an approximation procedure not directly related to PDEs, hence the WENO procedure can also be used in many non-PDE applications. In this paper we review the history and basic formulation of WENO schemes, outline the main ideas in using WENO schemes to solve various hyperbolic PDEs and other convection dominated problems, and present a collection of applications in areas including computational fluid dynamics, computational astronomy and astrophysics, semiconductor device simulation, traffic flow models, computational biology, and some non-PDE applications. Finally, we mention a few topics concerning WENO schemes that are currently under investigation.