"...Common misconceptions regarding the embedding approach are addressed including whether or not it results in an average value control model (no), is necessary to "tweak" the algorithm to get bang-bang solutions (no), requires infinite switching (no), has real-time capability (yes), or reduction to a classical nonlinear optimization problem (a desirable yes)."

"We present a new algorithm to solve polynomial equations, and publish its code, which is 1.6-3 times faster than the ZROOTS subroutine that is commercially available from Numerical Recipes, depending on application. The largest improvement, when compared to naive solvers, comes from a fail-safe procedure that permits us to skip the majority of the calculations in the great majority of cases, without risking catastrophic failure in the few cases that these are actually required. Second, we identify a discriminant that enables a rational choice between Laguerre's Method and Newton's Method (or a new intermediate method) on a case-by-case basis. We briefly review the history of root solving and demonstrate that "Newton's Method" was discovered neither by Newton (1671) nor by Raphson (1690), but only by Simpson (1740). Some of the arguments leading to this conclusion were first given by the British historian of science Nick Kollerstrom in 1992, but these do not appear to have penetrated the astronomical community. Finally, we argue that Numerical Recipes should voluntarily surrender its copyright protection for non-profit applications, despite the fact that, in this particular case, such protection was the major stimulant for developing our improved algorithm."

REBOUND is a new multi-purpose N-body code which is freely available under an open-source license. It was designed for collisional dynamics such as planetary rings but can also solve the classical N-body problem. It is highly modular and can be customized easily to work on a wide variety of different problems in astrophysics and beyond.

"A simplified model of natural convection, similar to the Lorenz (1963) system, is compared to computational fluid dynamics simulations in order to test data assimilation methods and better understand the dynamics of convection. The thermosyphon is represented by a long time flow simulation, which serves as a reference "truth". Forecasts are then made using the Lorenz-like model and synchronized to noisy and limited observations of the truth using data assimilation. The resulting analysis is observed to infer dynamics absent from the model when using short assimilation windows.

Furthermore, chaotic flow reversal occurrence and residency times in each rotational state are forecast using analysis data. Flow reversals have been successfully forecast in the related Lorenz system, as part of a perfect model experiment, but never in the presence of significant model error or unobserved variables. Finally, we provide new details concerning the fluid dynamical processes present in the thermosyphon during these flow reversals."

"…While several practical uses of the method have been proposed in Monte Carlo applications, these require an implementation framework that is flexible, general and efficient. We present such a framework for functions that are either strictly linear, concave, or convex on the unit interval using a series of envelope functions defined through a cascade, and show that this method not only greatly reduces the number of input bits needed in practice compared to other currently proposed solutions for more specific problems, but can easily be coupled to more asymptotically efficient methods to allow for theoretically strong results."

"This paper considers the fundamental problem of recovering a general signal, an image for example, from the magnitude of its Fourier transform. This problem, also known as phase retrieval, arises in many applications and has challenged engineers, physicists, and mathematicians for decades. Its origin comes from the fact that detectors can often times only record the squared modulus of the Fresnel or Fraunhofer diffraction pattern of the radiation that is scattered from an object. In such settings, one cannot measure the phase of the optical wave reaching the detector and, therefore, much information about the scattered object or the optical field is lost since, as is well known, the phase encodes a lot of the structural content of the image we wish to form."

"We present a very fast algorithm for general matrix factorization of a data matrix for use in the statistical analysis of high-dimensional data via latent factors. Such data are prevalent across many application areas and generate an ever-increasing demand for methods of dimension reduction in order to undertake the statistical analysis of interest. Our algorithm uses a gradient-based approach which can be used with an arbitrary loss function provided the latter is differentiable. The speed and effectiveness of our algorithm for dimension reduction is demonstrated in the context of supervised classification of some real high-dimensional data sets from the bioinformatics literature."

"We consider the problem of parameter estimation for a system of ordinary differential equations from noisy observations on a solution of the system. In case the system is nonlinear, as it typically is in practical applications, an analytic solution to it usually does not exist. Consequently, straightforward estimation methods like the ordinary least squares method depend on repetitive use of numerical integration in order to determine the solution of the system for each of the parameter values considered, and to find subsequently the parameter estimate that minimises the objective function. This induces a huge computational load to such estimation methods. We propose an estimator that is defined as a minimiser of an appropriate distance between a nonparametrically estimated derivative of the solution and the right-hand side of the system applied to a nonparametrically estimated solution.…"

"In this paper, a parametric level set method for reconstruction of obstacles in general inverse problems is considered. General evolution equations for the reconstruction of unknown obstacles are derived in terms of the underlying level set parameters. We show that using the appropriate form of parameterizing the level set function results a significantly lower dimensional problem, which bypasses many difficulties with traditional level set methods, such as regularization, re-initialization and use of signed distance function.…"

"We propose a degree-based coarse graining approach that not just accelerates the evaluation of dynamics on complex networks, but also satisfies the consistency conditions for both equilibrium statistical distributions and nonequilibrium dynamical flows. For the Ising model and susceptible-infected-susceptible epidemic model, we introduce these required conditions explicitly and further prove that they are satisfied by our coarse-grained network construction within the annealed network approximation. Finally, we numerically show that the phase transitions and fluctuations on the coarse-grained network are all in good agreements with those on the original one."

"We introduced an effective algorithm for k-nearest neighbor problem which works on multiple GPUs. By an experiment, we have shown that it runs more than 330 times faster than an implementation on a single core of an up-to-date CPU. We have also shown that the algorithm is effective from the viewpoint of parallelism of GPUs. That is because 1) there is no synchronization between GPUs until the very end of the process and 2) the workload is well balanced."

"Many computerized methods for RNA-RNA interaction structure prediction have been developed. Recently, $O(N^6)$ time and $O(N^4)$ space dynamic programming algorithms have become available that compute the partition function of RNA-RNA interaction complexes. However, few of these methods incorporate the knowledge concerning related sequences, thus relevant evolutionary information is often neglected from the structure determination. Therefore, it is of considerable practical interest to introduce a method taking into consideration both thermodynamic stability and sequence covariation. We present the \emph{a priori} folding algorithm \texttt{ripalign}, whose input consists of two (given) multiple sequence alignments (MSA). \texttt{ripalign} outputs (1) the partition function, (2) base-pairing probabilities, (3) hybrid probabilities and (4) a set of Boltzmann-sampled suboptimal structures consisting of canonical joint structures that are compatible to the alignments.…"

"A new floating-point arithmetic called precision arithmetic is developed to track precision for arithmetic calculations. It uses a novel rounding scheme to avoid excessive rounding error propagation of conventional floating-point arithmetic. Unlike interval arithmetic, its uncertainty tracking is based on statistics and its bounding range is much tighter. Generic standards and systematic methods for validating uncertainty-bearing arithmetics are discussed. The precision arithmetic is found to be better than interval arithmetic in uncertainty-tracking for linear algorithms. "

"The well known algorithm of VOLKER STRASSEN for matrix multiplication can only be used for $(m2^k \times m2^k)$ matrices. For arbitrary $(n \times n)$ matrices one has to add zero rows and columns to the given matrices to use STRASSEN's algorithm. … The aim of this work is to give a detailed analysis of the number of additional zero rows and columns and the additional work caused by STRASSEN's bad parameters. STRASSEN used the parameters $m$ and $k$ to show that his matrix multiplication algorithm needs less than $4.7n^{\log_2 7}$ flops. We can show in this paper, that these parameters cause an additional work of approx. 20 % in the worst case in comparison to the optimal strategy for the worst case. This is the main reason for the search for better parameters."

"We present a novel method for the calculation of free energy landscapes. Our approach involves a history dependent bias potential which is evaluated on a grid. The corresponding free energy landscape is constructed via a histogram reweighting procedure a posteriori. Due to the presence of the bias potential, our method can also be used to accelerate rare events. In addition, the calculated free energy landscape is not restricted to the actual choice of collective variables and can in principle be extended to all variables of interest without further numerical effort. We present numerical results for the alanine dipeptide and the Met-Enkephalin in explicit solution to illustrate our approach."

"We present a framework for simulating signal propagation in geometric networks (i.e. networks that can be mapped to geometric graphs in some space) and for developing algorithms that estimate (i.e. map) the state and functional topology of complex dynamic geometric net- works. Within the framework we define the key features typically present in such networks and of particular relevance to biological cellular neural networks: Dynamics, signaling, observation, and control. The framework is particularly well-suited for estimating functional connectivity in cellular neural networks from experimentally observable data, and has been implemented using graphics processing unit (GPU) high performance computing. Computationally, the framework can simulate cellular network signaling close to or faster than real time. We further propose a standard test set of networks to measure performance and compare different mapping algorithms."

Also: how about an inverse Doomsday algorithm. "We propose a simplification of a key component in the Doomsday Algorithm for calculating the day-of-the-week of any given date.…"

"The main goal of the paper is to introduce methods which compute B\'ezier curves faster than Casteljau's method does. These methods are based on the spectral factorization of a $n\times n$ Bernstein matrix, $B^e_n(s)= P_nG_n(s)P_n^{-1}$, where $P_n$ is the $n\times n$ lower triangular Pascal matrix. So we first calculate the exact optimum positive value $t$ in order to transform $P_n$ in a scaled Toeplitz matrix, which is a problem that was partially solved by X. Wang and J. Zhou (2006). Then fast Pascal matrix-vector multiplications and strategies of polynomial evaluation are put together to compute B\'ezier curves. Nevertheless, when $n$ increases, more precise Pascal matrix-vector multiplications allied to affine transformations of the vectors of coordinates of the control points of the curve are then necessary to stabilize all the computation."

"We propose a novel distributed algorithm to cluster graphs. The algorithm recovers the solution obtained from spectral clustering without the need for expensive eigenvalue/vector computations. We prove that, by propagating waves through the graph, a local fast Fourier transform yields the local component of every eigenvector of the Laplacian matrix, which are used to cluster graphs. For large graphs, the proposed algorithm is orders of magnitude faster than random walk based approaches. We prove the equivalence of the proposed algorithm to spectral clustering and derive convergence rates. We also demonstrate the benefit of using this decentralized clustering algorithm to accelerate distributed estimation for sensor networks and for efficient computation of distributed multi-agent search strategies."

"Recovery of the sparsity pattern (or support) of a sparse vector from a small number of noisy linear samples is a common problem that arises in signal processing and statistics. In the high dimensional setting, it is known that recovery with a vanishing fraction of errors is impossible if the sampling rate and per-sample signal-to-noise ratio (SNR) are finite constants independent of the length of the vector. In this paper, it is shown that recovery with an arbitrarily small but constant fraction of errors is, however, possible, and that in some cases a computationally simple thresholding estimator is near-optimal.…"

"We consider the problem of doing fast and reliable estimation of the number of non-zero entries in a sparse boolean matrix product. This problem has applications in databases and computer algebra. Let n denote the total number of non-zero entries in the input matrices. We show how to compute a 1 +- epsilon approximation (with small probability of error) in expected time O(n) for any epsilon > 4/\sqrt[4]{n}. The previously best estimation algorithm, due to Cohen (JCSS 1997), uses time O(n/epsilon^2). We also present a variant using O(sort(n)) I/Os in expectation in the cache-oblivious model. In contrast to these results, the currently best algorithms for computing a sparse boolean matrix product use time omega(n^{4/3}) (resp. omega(n^{4/3}/B) I/Os), even if the result matrix has only z=O(n) nonzero entries.…"

"We present a novel algorithm for segmentation of natural images that harnesses the principle of minimum description length (MDL). Our method is based on observations that a homogeneously textured region of a natural image can be well modeled by a Gaussian distribution and the region boundary can be effectively coded by an adaptive chain code. The optimal segmentation of an image is the one that gives the shortest coding length for encoding all textures and boundaries in the image, and is obtained via an agglomerative clustering process applied to a hierarchy of decreasing window sizes as multi-scale texture features. The optimal segmentation also provides an accurate estimate of the overall coding length and hence the true entropy of the image. We test our algorithm on the publicly available Berkeley Segmentation Dataset. It achieves state-of-the-art segmentation results compared to other existing methods."

I wonder to what extent this sort of thing might be useful in heuristics for Pareto-ranking

"A key problem in multiobjective linear programming is to find the set of all efficient extreme points in objective space. In this paper we introduce oriented projective geometry as an efficient and effective framework for solving this problem. The key advantage of oriented projective geometry is that we can work with an "optimally simple" but unbounded efficiency-equivalent polyhedron, yet apply to it the familiar theory and algorithms that are traditionally restricted to bounded polytopes. We apply these techniques to Benson's outer approximation algorithm, using oriented projective geometry to remove an exponentially large complexity from the algorithm and thereby remove a significant burden from the running time."

"Sparse modeling is a powerful framework for data analysis and processing. Traditionally, encoding in this framework is performed by solving an L1-regularized linear regression problem, commonly referred to as Lasso or Basis Pursuit. In this work we combine the sparsity-inducing property of the Lasso model at the individual feature level, with the block-sparsity property of the Group Lasso model, where sparse groups of features are jointly encoded, obtaining a sparsity pattern hierarchically structured. This results in the Hierarchical Lasso (HiLasso), which shows important practical modeling advantages.…"

"… In this paper, we provide a formal introduction to riffled independence and present algorithms for using riffled independence within Fourier-theoretic frameworks which have been explored by a number of recent papers. Additionally, we propose an automated method for discovering sets of items which are riffle independent from a training set of rankings. We show that our clustering-like algorithms can be used to discover meaningful latent coalitions from real preference ranking datasets and to learn the structure of hierarchically decomposable models based on riffled independence."

"Monte Carlo methods play an important role in scientific computation, especially when problems have a vast phase space. In this chapter an introduction to the Monte Carlo method is given. Concepts such as Markov chains, detailed balance, critical slowing down, and ergodicity, as well as the Metropolis algorithm are explained. The Monte Carlo method is illustrated by numerically studying the critical behavior of the two-dimensional Ising ferromagnet using finite-size scaling methods. Furthermore, advanced Monte Carlo methods are described (parallel tempering Monte Carlo) and illustrated with nontrivial models from the physics of glassy systems."

"The beam-splitting method for fine-heterogeneity cor- rection will inevitably multiply beams to transport and thus will slow down dose calculation. With the GDS al- gorithm, the dose convolution is made only once after all the beams have been transported, which minimizes the impact of the beam multiplication on computing time. In fact, for the beams individually split into several tens, the calculation time increased only by several times with the GDS. This algorithmic framework will thus enable fast and accurate treatment planning of heavy charged particle ra- diotherapy in the presence of density heterogeneity finer than the size of intrinsic beam blurring."

"This paper presents the design and analysis of parallel approximation algorithms for facility-location problems, including $\NC$ and $\RNC$ algorithms for (metric) facility location, $k$-center, $k$-median, and $k$-means. These problems have received considerable attention during the past decades from the approximation algorithms community, concentrating primarily on improving the approximation guarantees. In this paper, we ask, is it possible to parallelize some of the beautiful results from the sequential setting?"

I want to make this sort of thing less newsworthy.

awesome

"We overview our recently introduced theory of n-fold integer programming which enables the polynomial time solution of fundamental linear and nonlinear integer programming problems in variable dimension. We demonstrate its power by obtaining the first polynomial time algorithms in several application areas including multicommodity flows and privacy in statistical databases."

"We consider the problem of developing an efficient multi-threaded implementation of the matrix-vector multiplication algorithm for sparse matrices with structural symmetry. Matrices are stored using the compressed sparse row-column format (CSRC), designed for profiting from the symmetric non-zero pattern observed in global finite element matrices. Unlike classical compressed storage formats, performing the sparse matrix-vector product using the CSRC requires thread-safe access to the destination vector. To avoid race conditions, we have implemented two partitioning strategies. In the first one, each thread allocates an array for storing its contributions, which are later combined in an accumulation step. We analyze how to perform this accumulation in four different ways.…"

"We consider methods for finding high-precision approximations to simple zeros of smooth functions. As an application, we give fast methods for evaluating the elementary functions log(x), exp(x), sin(x) etc. to high precision. For example, if x is a positive floating-point number with an n-bit fraction, then (under rather weak assumptions) an n-bit approximation to log(x) or exp(x) may be computed in time asymptotically equal to 13M(n)lg(n), where M(n) is the time required to multiply floating-point numbers with n-bit fractions. Similar results are given for the other elementary functions.…"

"In this paper, the performances of the quasi-Newton BFGS algorithm, the NEWUOA derivative free optimizer, the Covariance Matrix Adaptation Evolution Strategy (CMA-ES), the Differential Evolution (DE) algorithm and Particle Swarm Optimizers (PSO) are compared experimentally on benchmark functions reflecting important challenges encountered in real-world optimization problems. Dependence of the performances in the conditioning of the problem and rotational invariance of the algorithms are in particular investigated."

"…Since the produced mesh is applied to the meshless-type methods, the connectivity of the points is not used and only the grid points are important, though the grid lines are utilized in the adapting process. The performance of the adaptive points is exam- ined by considering a collocation meshless method which is based on interpolation in terms of a set of radial basis functions. … Some experimental results will be presented to illustrate the effectiveness of the proposed method."

"Our numerical studies suggest that the proposed CP-WOPT approach is accurate and scalable. CP-WOPT can recover the underlying factors successfully with large amounts of missing data, e.g., 90% missing entries for tensors of size 50 × 40 × 30. We have also studied how CP-WOPT can scale to problems of larger sizes, e.g., 1000 × 1000 × 1000, and recover CP factors from large, sparse tensors with 99.5% missing data.…"

"Although care needs to be taken in the implementation of normal random number generators like fastnorm, and the end-user should be aware of the small but unavoidable defects discussed in §§5.6-5.7, these generators have such a performance advantage over more conventional generators that they can not be ignored in applications where the speed of generation of pseudo- random numbers is critical."

"An exact, one-to-one transform is presented that not only allows digital circular convolutions, but is free from multiplications and quantisation errors for transform lengths of arbitrary powers of two. The transform is analogous to the Discrete Fourier Transform, with the canonical harmonics replaced by a set of cyclic integers computed using only bit-shifts and additions modulo a prime number. The prime number may be selected to occupy contemporary word sizes or to be very large for cryptographic or data hiding applications. The transform is an extension of the Rader Transforms via Carmichael's Theorem. These properties allow for exact convolutions that are impervious to numerical overflow and to utilise Fast Fourier Transform algorithms."

"…This approach to analyze functions is hence more suitable for computers (in order to save computation time), and the simplicity of the definition allows further research in other areas of Classical Analysis."

"We describe von Neumann's elegant idea for sampling from the exponential distribution, Forsythe's generalization for sampling from a probability distribution whose density has the form exp(-G(x)), where G(x) is easy to compute (e.g. a polynomial), and my refinement of these ideas to give an efficient algorithm for generating pseudo-random numbers with a normal distribution. Later developments are also mentioned."

"We present an efficient solution to the following problem, of relevance in a numerical optimization scheme: calculation of integrals of the type \[\iint_{T \cap \{f\ge0\}} \phi_1\phi_2 \, dx\,dy\] for quadratic polynomials $f,\phi_1,\phi_2$ on a plane triangle $T$. The naive approach would involve consideration of the many possible shapes of $T\cap\{f\geq0\}$ (possibly after a convenient transformation) and parameterizing its border, in order to integrate the variables separately. Our solution involves partitioning the triangle into smaller triangles on which integration is much simpler."

"NArray is an Numerical N-dimensional Array class. Supported element types are 1/2/4-byte Integer, single/double-precision Real/Complex, and Ruby Object. This extension library incorporates fast calculation and easy manipulation of large numerical arrays into the Ruby language. NArray has features similar to NumPy, but NArray has vector and matrix subclasses."

"I recently got back reviews of a paper in which I used automatic differentiation. Therein, a reviewer clearly thought I was using finite difference, or “numerical” differentiation. This has led me to wondering: Why don’t machine learning people use automatic differentiation more? Why don’t they use it…constantly? Before recklessly speculating on the answer, let me briefly review what automatic differentiation (henceforth “autodiff”) is. Specifically, I will be talking about reverse-mode autodiff."