The purpose of this paper is to provide a comprehensive presentation and interpretation of the Ensemble Kalman Filter (EnKF) and its numerical implementation. The EnKF has a large user group, and numerous publications have discussed applications and theoretical aspects of it. This paper reviews the important results from these studies and also presents new ideas and alternative interpretations which further explain the success of the EnKF. In addition to providing the theoretical framework needed for using the EnKF, there is also a focus on the algorithmic formulation and optimal numerical implementation. A program listing is given for some of the key subroutines. The paper also touches upon specific issues such as the use of nonlinear measurements, in situ profiles of temperature and salinity, and data which are available with high frequency in time. An ensemble based optimal interpolation (EnOI) scheme is presented as a cost-effective approach which may serve as an alternative to the EnKF in some applications. A fairly extensive discussion is devoted to the use of time correlated model errors and the estimation of model bias.

The Ensemble Kalman Filter: Theoretical formulation and practical implementation

A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and variance simply requires the evaluation of the corresponding Gaussian density. In the $n$ -dimensional setting, however, it requires the inversion of an $n \times n$ covariance matrix, $C$ , as well as the evaluation of its determinant, $\det (C)$ . In many cases, such as regression using Gaussian processes, the covariance matrix is of the form $C = \sigma ^2 I + K$ , where $K$ is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix $C$ is typically dense, causing standard direct methods for inversion and determinant evaluation to require $\mathcal {O}(n^3)$ work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix $C$ can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an $\mathcal {O} (n\,\log^2\, n)$ algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant $\det (C)$ , permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining $K$ . Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.

Fast Direct Methods for Gaussian Processes

A three-dimensional imaging system which distributes the optical illumination over the full field-of-view is sought after. Here, the authors demonstrate the capability of reconstructing 128 × 128 pixel resolution three-dimensional scenes to an accuracy of 3 mm as well as real-time video with a frame-rate up to 12 Hz.

Single-pixel three-dimensional imaging with time-based depth resolution

The vast quantity of information brought by big data as well as the evolving computer hardware encourages success stories in the machine learning community. In the meanwhile, it poses challenges for the Gaussian process regression (GPR), a well-known nonparametric, and interpretable Bayesian model, which suffers from cubic complexity to data size. To improve the scalability while retaining desirable prediction quality, a variety of scalable GPs have been presented. However, they have not yet been comprehensively reviewed and analyzed to be well understood by both academia and industry. The review of scalable GPs in the GP community is timely and important due to the explosion of data size. To this end, this article is devoted to reviewing state-of-the-art scalable GPs involving two main categories: global approximations that distillate the entire data and local approximations that divide the data for subspace learning. Particularly, for global approximations, we mainly focus on sparse approximations comprising prior approximations that modify the prior but perform exact inference, posterior approximations that retain exact prior but perform approximate inference, and structured sparse approximations that exploit specific structures in kernel matrix; for local approximations, we highlight the mixture/product of experts that conducts model averaging from multiple local experts to boost predictions. To present a complete review, recent advances for improving the scalability and capability of scalable GPs are reviewed. Finally, the extensions and open issues of scalable GPs in various scenarios are reviewed and discussed to inspire novel ideas for future research avenues.

/pdf/when-gaussian-process-meets-big-data-a-review-of-scalable-489iuo4vst.pdf

When Gaussian Process Meets Big Data: A Review of Scalable GPs

GPflow is a Gaussian process library that uses TensorFlow for its core computations and Python for its front end The distinguishing features of GPflow are that it uses variational inference as the primary approximation method, provides concise code through the use of automatic differentiation, has been engineered with a particular emphasis on software testing and is able to exploit GPU hardware

/pdf/gpflow-a-gaussian-process-library-using-tensorflow-2l6t2dv7rk.pdf

GPflow: a Gaussian process library using tensorflow

. We study the rank of sub-matrices arising out of kernel functions, F ( xxx,yyy ) : R d × R d (cid:55)→ R , where xxx,yyy ∈ R d with F ( xxx,yyy ) is smooth everywhere except along the line yyy = xxx . Such kernel functions are frequently encountered in a wide range of applications such as N body problems, Green’s functions, integral equations, geostatistics, kriging, Gaussian processes, radial basis function interpolation, etc. One of the challenges in dealing with these kernel functions is that the corresponding matrix associated with these kernels is large and dense; thereby, the computational cost of matrix operations is high. In this article, we prove new theorems bounding the numerical rank of sub-matrices arising out of these kernel functions. Under reasonably mild assumptions, we prove that certain sub-matrices are rank-deﬁcient in ﬁnite precision. The rank of these sub-matrices depends on the dimension of the ambient space and also on the type of interaction between the hyper-cubes containing the corresponding set of particles. This rank structure can be leveraged to reduce the computational cost of certain matrix operations such as matrix-vector products, solving linear systems, etc. We also present numerical results on the growth of rank of certain sub-matrices in 1 D, 2 D, 3 D and 4 D, which, not surprisingly, agrees with the theoretical results.

Numerical rank of kernel functions

Ignition delay time (IDT) is an important global combustion property that affects the thermal efficiency of the engine and emissions (particularly NO ). IDT is generally measured by performing experiments using Shock-tube and Rapid Compression Machine (RCM). The numerical calculation of IDT is a computationally expensive and time-consuming process. Arrhenius type empirical correlations offer an inexpensive alternative to obtain IDT. However, such correlations have limitations as these typically involve ad-hoc parameters and are valid only for a specific fuel and particular range of temperature/pressure conditions. This study aims to formulate a data-driven scientific way to obtain IDT for new fuels without performing detailed numerical calculations or relying on ad-hoc empirical correlations. We propose a rigorous, well-validated data-driven study to obtain IDT for new fuels using a regression-based clustering algorithm. In our model, we bring in the fuel structure through the overall activation energy ( ) by expressing it in terms of the different bonds present in the molecule. Gaussian Mixture Model (GMM) is used for the formation of clusters, and regression is applied over each cluster to generate models. The proposed algorithm is used on the ignition delay dataset of straight-chain alkanes (C H for n = 4 to 16). The high level of accuracy achieved demonstrates the applicability of the proposed algorithm over IDT data. The algorithm and framework discussed in this article are implemented in python and made available at https://doi.org/10.5281/zenodo.5774617.

A data-driven framework to predict ignition delays of straight-chain alkanes

. This article looks at the eigenvalues of magicsquares generated by the MATLAB’s magic( n ) function. The magic() function constructs doubly even ( n = 4 k ) magic squares, singly even ( n = 4 k + 2) magic squares and odd ( n = 2 k + 1) magic squares using diﬀerent algorithms. The doubly even magic squares are constructed by a criss-cross method that involves reﬂecting the entries of a simple square about the center. The odd magic squares are constructed using the Siamese method. The singly even magic squares are constructed using a lower-order odd magic square (Strachey method). We obtain approximations of eigenvalues of odd and singly even magic squares with error bounds. Further, the eigenvalues of doubly even magic squares are exactly obtained. The approximation of the spectra involves some interesting connections with the spectrum of g-circulant matrices and the use of Bauer-Fike theorem.

Spectrum of MATLABs magic squares

This article discusses a more general and numerically stable Rybicki Press algorithm, which enables inverting and computing determinants of covariance matrices, whose elements are sums of exponentials. The algorithm is true in exact arithmetic and relies on introducing new variables and corresponding equations, thereby converting the matrix into a banded matrix of larger size. Linear complexity banded algorithms for solving linear systems and computing determinants on the larger matrix enable linear complexity algorithms for the initial semi-separable matrix as well. Benchmarks provided illustrate the linear scaling of the algorithm.

Generalized Rybicki Press algorithm

Using accurate multi-component diffusion treatment in numerical combustion studies remains formidable due to the computational cost associated with solving for diffusion velocities. To obtain the diffusion velocities, for low density gases, one needs to solve the Stefan-Maxwell equations along with the zero diffusion flux criteria, which scales as $\mathcal{O}(N^3)$, when solved exactly. In this article, we propose an accurate, fast, direct and robust algorithm to compute multi-component diffusion velocities. To our knowledge, this is the first provably accurate algorithm (the solution can be obtained up to an arbitrary degree of precision) scaling at a computational complexity of $\mathcal{O}(N)$ in finite precision. The key idea involves leveraging the fact that the matrix of the reciprocal of the binary diffusivities, $V$, is low rank, with its rank being independent of the number of species involved. The low rank representation of matrix $V$ is computed in a fast manner at a computational complexity of $\mathcal{O}(N)$ and the Sherman-Morrison-Woodbury formula is used to solve for the diffusion velocities at a computational complexity of $\mathcal{O}(N)$. Rigorous proofs and numerical benchmarks illustrate the low rank property of the matrix $V$ and scaling of the algorithm.

Sivaram Ambikasaran

Papers

Numerical rank of kernel functions

A data-driven framework to predict ignition delays of straight-chain alkanes

Spectrum of MATLABs magic squares

Generalized Rybicki Press algorithm

An accurate, fast, mathematically robust, universal, non-iterative algorithm for computing multi-component diffusion velocities