Accumulation Of Round-off Errors Research Articles

Software for "Exact Matrix Factorization Updates for Nonlinear Programming"

LU and Cholesky matrix factorization algorithms are core subroutines used to solve systems of linear equations (SLEs) encountered when solving an optimization prob-lem. Standard floating-point algorithms are highly efficient but remain susceptible to the accumulation of round-off errors, which can lead solvers to return feasibility and optimal-ity claims that are actually invalid. This paper introduces a novel direct solution approach for solving sequences of closely related SLEs encountered in nonlinear programming effi-ciently and without round-off errors. Specifically, it introduces rank-one update algorithms for the round-off error–free factorization framework, a tool set built on integer-preserving arithmetic that has led to the development and implementation of extremely reliable sub-routines for solving SLEs occurring in linear programming. The formal guarantees of the presented algorithms are established through the derivation of theoretical insights. Their advantages are supported with computational experiments, which demonstrate upward of 75×improvements over exact factorization runtimes on fully dense matrices with more than one million entries. A significant advantage of the featured integer-preserving frame-work is that the length of any matrix coefficient produced by its algorithms is bounded polynomially in the size of the inputs without having to resort to greatest common divisor operations, which are required by and thereby hinder an efficient implementation of exact rational arithmetic approaches.

Exact Matrix Factorization Updates for Nonlinear Programming

LU and Cholesky matrix factorization algorithms are core subroutines used to solve systems of linear equations (SLEs) encountered when solving an optimization problem. Standard floating-point algorithms are highly efficient but remain susceptible to the accumulation of round-off errors, which can lead solvers to return feasibility and optimality claims that are actually invalid. This paper introduces a novel direct solution approach for solving sequences of closely related SLEs encountered in nonlinear programming efficiently and without round-off errors. Specifically, it introduces rank-one update algorithms for the round-off error–free factorization framework, a tool set built on integer-preserving arithmetic that has led to the development and implementation of extremely reliable subroutines for solving SLEs occurring in linear programming. The formal guarantees of the presented algorithms are established through the derivation of theoretical insights. Their advantages are supported with computational experiments, which demonstrate upward of 75× improvements over exact factorization runtimes on fully dense matrices with more than one million entries. A significant advantage of the featured integer-preserving framework is that the length of any matrix coefficient produced by its algorithms is bounded polynomially in the size of the inputs without having to resort to greatest common divisor operations, which are required by and thereby hinder an efficient implementation of exact rational arithmetic approaches. History: Accepted by Antonio Frangioni, Area Editor for Design & Analysis of Algorithms–Continuous. Supplemental Material: The software that supports the findings of this study is available within the paper and its Supplemental Information ( https://pubsonline.informs.org/doi/suppl/10.1287/ijoc.2021.0331 ) as well as from the IJOC GitHub software repository ( https://github.com/INFORMSJoC/2021.0331 ). The complete IJOC Software and Data Repository is available at https://informsjoc.github.io/ .

A reduced order with data assimilation model: Theory and practice

Numerical simulations are extensively used as a predictive tool to better understand complex air flows and pollution transport on the scale of individual buildings, city blocks and entire cities. Fast-running Non-Intrusive Reduced Order Model (NIROM) for predicting the turbulent air flows has been proved to be an efficient method to provide numerical forecasting results. However, due to the reduced space on which the model operates, the solution includes uncertainties. Additionally, any computational methodology contributes to uncertainty due to finite precision and the consequent accumulation and amplification of round-off errors. Taking into account these uncertainties is essential for the acceptance of any numerical simulation. In this paper we combine the NIROM method with Data Assimilation (DA), the main question is how to incorporate data (e.g. from physical measurements) in models in a suitable way, in order to improve model predictions and quantify prediction uncertainty. Here, the focus is on the prediction of nonlinear dynamical systems (the classical application example being weather forecasting). DA is an uncertainty quantification technique used to incorporate observed data into a prediction model in order to improve numerical forecasted results. The Reduced Order Data Assimilation (RODA) model we propose in this paper achieves both efficiency and accuracy including Variational DA into NIROM. The model we present is applied to the pollutant dispersion within an urban environment.

Spherical radial basis function approximation of some physical quantities measured

The article treats the spherical interpolation and approximation, i.e., a way of measured data processing in case when the scalar data sampling is performed at nodes on the unit 2D sphere surface in the 3D Euclidean space (in general, on the (d−1)-dimensional sphere surface in the d-dimensional space). We use the spherical radial basis function interpolation with an inverse multiquadric basis function and a second degree polynomial in Cartesian coordinates as a trend.The formulae of this type can be useful in the interpretation of various physical measurements and have wide applications in geosciences, e.g., in the treatment of anisotropy of magnetic susceptibility data measured. We present the advantages of the formula chosen on numerical examples. However, in practical computation with high number of sampling nodes, the matrix of the system for determining interpolation coefficients may be ill-conditioned. Then the standard double precision LU factorization does not give a reliable solution due to the effect of round-off error accumulation and some more sophisticated means of solving the system are to be used.

Optimized integrating factor technique for Schrödinger-like equations

The integrating factor technique is widely used to solve numerically (in particular) the Schrödinger equation in the context of spectral methods. Here, we present an improvement of this method exploiting the freedom provided by the gauge condition of the potential. Optimal gauge conditions are derived considering the equation and the temporal numerical resolution with an adaptive embedded scheme of arbitrary order. We illustrate this approach with the nonlinear Schrödinger (NLS) and with the Schrödinger–Newton (SN) equations. We show that this optimization increases significantly the overall computational speed, sometimes by a factor five or more. This gain is crucial for long time simulations, as, with larger time steps, less computations are performed and the overall accumulation of round-off errors is reduced.

Existence of the first magic angle for the chiral model of bilayer graphene

We consider the chiral model of twisted bilayer graphene introduced by Tarnopolsky, Kruchkov, and Vishwanath (TKV). TKV proved that for inverse twist angles α such that the effective Fermi velocity at the moiré K point vanishes, the chiral model has a perfectly flat band at zero energy over the whole Brillouin zone. By a formal expansion, TKV found that the Fermi velocity vanishes at α ≈ 0.586. In this work, we give a proof that the Fermi velocity vanishes for at least one α between 0.57 and 0.61 by rigorously justifying TKV’s formal expansion of the Fermi velocity over a sufficiently large interval of α values. The idea of the proof is to project the TKV Hamiltonian onto a finite-dimensional subspace and then expand the Fermi velocity in terms of explicitly computable linear combinations of modes in the subspace while controlling the error. The proof relies on two propositions whose proofs are computer-assisted, i.e., numerical computation together with worst-case estimates on the accumulation of round-off error, which show that round-off error cannot possibly change the conclusion of computation. The propositions give a bound below on the spectral gap of the projected Hamiltonian, an Hermitian 80 × 80 matrix whose spectrum is symmetric about 0, and verify that two real 18-th order polynomials, which approximate the numerator of the Fermi velocity, take values with a definite sign when evaluated at specific values of α. Together with TKV’s work, our result proves the existence of at least one perfectly flat band of the chiral model.

An Adaptive Method for Numerical Differentiation of Difficult-to-Compute Functions

An adaptive approach to the numerical differentiation of difficult-to-compute functions is considered. Complex dependencies, which are the result of multiple superpositions of functions or the product of various algorithmic processes, are knowingly difficult to study directly. To establish the nature of the behavior of such dependencies, one has to resort to numerical analysis. One of the important characteristics of functions is a derivative, which indicates the direction and rate of change of a dependence. However, with difficult-to-compute functions, the available a priori information is not always sufficient to achieve the appropriate accuracy of the solution by known means. The loss of accuracy occurs due to the accumulation of round-off errors that grow in proportion to the number of calculated values of a function. In this case, it is necessary to pass on to the posterior approach in order to determine the behavior of the function and move away from the scheme of equidistant nodes, relying on an adaptive way of studying the local situation in the domain of the function. This paper implements an adaptive method for finding derivatives of a function with a minimum of restrictive requirements for the class of functions and the form of their assignment. Due to this, the costs of calculating the function have been significantly reduced with the result that their number has been brought to almost the optimal level. At the same time, the amount of RAM used has sharply decreased. There is no need for a preliminary analysis of the problem of establishing the class of the function under study, in the involvement of special functions or transformation of initial conditions for using standard tables of weight coefficients, etc. For research, it is enough to assign a continuous and bounded function on a fixed segment and a minimum step, which is indirectly responsible for ensuring the required accuracy of differentiation. The effectiveness of the proposed method is demonstrated on a number of test examples. The developed method can be used in more complex problems, for example, in solving some types of differential and integral equations, as well as for a wide range of optimization problems in a wide variety of areas of applied analysis and synthesis.

Practical considerations in the implementation of time-domain acoustic full waveform inversion

Full waveform inversion (FWI) plays a major role in the oil and gas industry as a state-of-the-art technique that produces quantitative subsurface structures with high-fidelity images. Various FWI studies have been conducted, and these suggest that FWI is a promising inversion method. Recently, many attempts have been made toward three-dimensional (3D) and four-dimensional FWI applications (which were difficult to perform in the past) because of the progress made in computer science and the growth of computer resources. To manage the very large data requirement of 3D problems, a time-domain FWI that is relatively efficient in terms of memory demands must be implemented. However, it could encounter practical issues, leading to failure in its convergence. In this paper, we introduce these practical issues and several alternative methods for mitigating them. The first issue is the bandpass filtering of the observed seismograms. We suggest that the frequency-domain filter based on a reference wavelet would be optimal in terms of both bandpass filtering and source wavelet estimation. The second issue is related to acoustic approximation. We show that a simple density model comprising only water and solid layers is a reasonable option to address seafloor reflectivity properly. The last issue is the accumulation of round-off errors due to the massive computation of the objective function. We demonstrate that a simple modification of the error calculation can resolve this round-off error problem.

Extended Kalman Filter Using Orthogonal Polynomials

This paper reports a new extended Kalman filter where the underlying nonlinear functions are linearized using a Gaussian orthogonal basis of a weighted $\mathcal {L}_{2}$ space. As we are interested in computing the states’ mean and covariance with respect to Gaussian measure, it would be better to use a linearization, that is optimal with respect to the same measure. The resulting first-order polynomial coefficients are approximately calculated by evaluating the integrals using (i) third-order Taylor series expansion (ii) cubature rule of integration. Compared to direct integration-based filters, the proposed filter is far less susceptible to the accumulation of round-off errors leading to loss of positive definiteness. The proposed algorithms are applied to four nonlinear state estimation problems. We show that our proposed filter consistently outperforms the traditional extended Kalman filter and achieves a competitive accuracy to an integration-based square root filter, at a significantly reduced computing cost.

A transient interval reliability analysis for software rejuvenation models with phase expansion

The phenomenon of software aging refers to the continuing degradation of software system performance with the operation time and is usually caused by the aging-related bugs such as the memory leak and the accumulation of round-off errors. Software rejuvenation acts as one of the proactive fault management techniques against the software aging. In this paper, we evaluate the interval reliability for two basic stochastic models with periodic software rejuvenation by Garg et al. (1995) and Suzuki et al. (2002a, b). The interval reliability is one of the most generalized dependability measures that involve commonly used reliability function and steady-state availability as the special cases, and is helpful for the system design during a fixed mission period. From the mathematical point of view, the interval reliability for the software rejuvenation models leads to the transient analysis. We focus on the phase expansion approach for solving the transient solutions for the basic software rejuvenation models. The phase expansion is an approximate technique that replaces arbitrary probability distributions by the phase-type (PH) distributions. Benefiting from the phase expansion, we can numerically derive the transient interval reliability for two software rejuvenation models. In numerical examples, we discuss the accuracy of the phase expansion and also reveal quantitative properties of the interval reliability measures.

Vector Penalty-Projection Methods for Open Boundary Conditions with Optimal Second-Order Accuracy

Recently, a new family of splitting methods, the so-called vector penalty-projection methods (VPP) were introduced by Angot et al. [2, 3] to compute the solution of unsteady incompressible fluid flows and to overcome most of the drawbacks of the usual incremental projection methods. Two different parameters are related to the VPP methods: the augmentation parameter r ≥ 0 and the penalty parameter 0 < e ≤ 1. In this paper, we deal with the time-dependent incompressible Navier-Stokes equations with outflow boundary conditions using the VPP methods. The spatial dis-cretization is based on the finite volume scheme on a Marker and Cells (MAC) staggered grid. Furthermore, two different second-order time discretization schemes are investigated: the second-order Backward Difference Formula (BDF2) known also as Gear's scheme and the Crank-Nicolson scheme. We show that the VPP methods provide a second-order convergence rate for both velocity and pressure in space and time even in the presence of open boundary conditions with small values of the augmentation parameter r typically 0 ≤ r ≤ 1 and a penalty parameter e small enough typically e = 10 −10. The resulting constraint on the discrete divergence of velocity is not exactly equal to zero but is satisfied approximately as O(e δt) where e is the penalty parameter (taken as small as desired) and δt is the time step. The choice r = 0 requires a special attention to avoid the accumulation of round-off errors for very small values of e. Indeed, it is important in this case to directly correct the pressure gradient by taking account of the velocity correction issued from the vector penalty-projection step. Finally, the efficiency and the second-order accuracy of the method are illustrated by several numerical test cases including homogeneous or non-homogeneous given traction on the boundary.

Simple but accurate periodic solutions for the nonlinear pendulum equation

Abstract Despite its elementary structure, the simple pendulum oscillations are described by a nonlinear differential equation whose exact solution for the angular displacement from vertical as a function of time cannot be expressed in terms of an elementary function, so either a numerical treatment or some analytical approximation is ultimately demanded. Such solutions have been thoroughly investigated due to the abundance of distinct pendular systems in nature and, more recently, due to the availability of automatic data acquisition systems in undergraduate laboratories. However, it is well-known that numerical solutions to differential equations usually loose accuracy (due to accumulation of roundoff errors) and polynomial approximations diverge after long time intervals. In this work, I take a few terms of the Fourier series expansion of the elliptic function sn (u;k) as a source of accurate periodic solutions for the pendulum equation. Interestingly, these approximations remain accurate for arbitrarily long time intervals, even for large amplitudes, which shows its adequacy for the analysis of experimental data gathered in classical mechanics classes.

Fast Conjugate Heat Transfer Simulation of Long Transient Flexible Operations Using Adaptive Time Stepping

The emerging renewable energy market calls for more advanced prediction tools for turbine transient operations in fast startup/shutdown cycles. Reliable numerical analysis of such transient cycles is complicated by the disparity in time scales of the thermal responses in fluid and solid domains. Obtaining fully coupled time-accurate unsteady conjugate heat transfer (CHT) results under these conditions would require to march in both domains using the time-step dictated by the fluid domain: typically, several orders of magnitude smaller than the one required by the solid. This requirement has strong impact on the computational cost of the simulation as well as being potentially detrimental to the accuracy of the solution due to accumulation of round-off errors in the solid. A novel loosely coupled CHT methodology has been recently proposed, and successfully applied to both natural and forced convection cases that remove these requirements through a source-term based modeling (STM) approach of the physical time derivative terms in the relevant equations. The method has been shown to be numerically stable for very large time steps with adequate accuracy. The present effort is aimed at further exploiting the potential of the methodology through a new adaptive time stepping approach. The proposed method allows for automatic time-step adjustment based on estimating the magnitude of the truncation error of the time discretization. The developed automatic time stepping strategy is applied to natural convection cases under long (2000 s) transients: relevant to the prediction of turbine thermal loads during fast startups/shutdowns. The results of the method are compared with fully coupled unsteady simulations showing comparable accuracy with a significant reduction of the computational costs.

ANALYTICAL INVERSE FOR THE SYMMETRIC CIRCULANT TRIDIAGONAL MATRIX

Finding the inverse of a matrix is an open problem especially when it comes to engineering problems due to their complexity and running time (cost) of matrix inversion algorithms. An optimum strategy to invert a matrix is, first, to reduce the matrix to a simple form, only then beginning a mathematical procedure. For symmetric matrices, the preferred simple form is tridiagonal. This makes tridiagonal matrices of high interest in applied mathematics and engineering problems. This study presents a time efficient, exact analytical approach for finding the inverse, decomposition, and solving linear systems of equations where symmetric circulant matrix appears. This matrix appears in many researches and it is different from ordinary tridiagonal matrices as there are two corner elements. For finding the inverse matrix, a set of matrices are introduced that any symmetric circulant matrix could be decomposed into them. After that, the exact analytical inverse of this set is found which gives the inverse of circulant matrix. Moreover, solving linear equations can be carried out using implemented decomposition where this matrix appears as a coefficient matrix. The methods principal strength is that it is as stable as any other direct methods (i.e., execute in a predictable number of operations). It is straightforward, understandable, solid as a rock and an exceptionally good psychological backup for those times when something is going wrong and you think it might be your linear equation solver. The downside of present method, and every other direct method, is the accumulation of round off errors.

Asymptotics and Well-Posedness of the Derived Distribution Density in a Study of Biovariability

In our recent work (Wang, Burgei, and Zhou, 2018) we studied the hearing loss injury among subjects in a crowd with a wide spectrum of heterogeneous individual injury susceptibility due to biovariability. The injury risk of a crowd is defined as the average fraction of injured. We examined mathematically the injury risk of a crowd vs the number of acoustic impulses the crowd is exposed to, under the assumption that all impulses act independently in causing injury regardless of whether one is preceded by another. We concluded that the observed dose-response relation can be explained solely on the basis of biovariability in the form of heterogeneous susceptibility. We derived an analytical solution for the distribution density of injury susceptibility, as a power series expansion in terms of scaled log individual non-injury probability. While theoretically the power series converges for all argument values, in practical computations with IEEE double precision, at large argument values, the numerical accuracy of the power series summation is completely wiped out by the accumulation of round-off errors. In this study, we derive a general asymptotic approximation at large argument values, for the distribution density. The combination of the power series and the asymptotics provides a practical numerical tool for computing the distribution density. We then use this tool to verify numerically that the distribution obtained in our previous theoretical study is indeed a proper density. In addition, we will also develop a very efficient and accurate Pade approximation for the distribution density.

Analytical solution for the dynamic model of tumbling mills

Optimisation of grinding circuits is invariably dependent on sound process models together with process simulators that can solve the process models accurately. Most of the process models are solved numerically because analytical solutions are not available, which can lead to errors in the results due to the numerical approximation of mathematical equations. Whiten [1], and Valery Jnr & Morrell [2,3] have developed a dynamic model with numerical simulation for autogenous and semi-autogenous mills, and validated the model with dynamic response of mills in terms of power draw, grinding charge level, slurry level and product size distribution to changes in feed rate, feed size, feed hardness and water addition [2,3]. In this work, an analytical solution for their dynamic model of tumbling mills has been developed based on the knowledge of solutions to the first-order nonhomogeneous linear differential equations. Two algorithms, Direct Single Time method (DST) and Direct Multiple Time method (DMT), were applied to obtain the analytical solutions respectively. It was found that analytical solutions are more accurate than the traditional finite difference numerical methods. However, the DST analytical method has a drawback of numerical instability due to the accumulation of round-off errors which are amplified by exponential functions, whilst the DMT method can provide stable solutions. To test the DMT analytical method, two cases of SAG mill dynamic operation were studied with both the traditional numerical method and the newly developed analytical method, further proving the robustness and feasibility of the analytical solutions.

Validated Computation of Heteroclinic Sets

In this work we develop a method for computing mathematically rigorous enclosures of some one dimensional manifolds of heteroclinic orbits for nonlinear maps. Our method exploits a rigorous curve following an argument built on high order Taylor approximation of the local stable/unstable manifolds. The curve following argument is a uniform interval Newton method applied on short line segments. The definition of the heteroclinic sets involve compositions of the map, and we use a Lohner-type representation to overcome the accumulation of roundoff errors. Our argument requires precise control over the local unstable and stable manifolds so that we must first obtain validated a posteriori error bounds on the truncation errors associated with the manifold approximations. We illustrate the utility of our method by proving some computer assisted theorems about heteroclinic invariant sets for a volume preserving map of $\mathbb{R}^3$.

Implicit predictor–corrector central finite difference scheme for the equations of magnetohydrodynamic simulations

Nowadays most of supercomputers are based on the frame of PC cluster; therefore, the efficiency of parallel computing is of importance especially with the increasing computing scale. This paper proposes a high-order implicit predictor–corrector central finite difference (iPCCFD) scheme and demonstrates its high efficiency in parallel computing. Of special interests are the large scale numerical studies such as the magnetohydrodynamic (MHD) simulations in the planetary magnetosphere. An iPCCFD scheme is developed based on fifth-order central finite difference method and fourth-order implicit predictor–corrector method in combination with elimination-of-the-round-off-errors (ERE) technique. We examine several numerical studies such as one-dimensional Brio–Wu shock tube problem, two-dimensional Orszag–Tang vortex system, vortex type K–H instability, kink type K–H instability, field loop advection, and blast wave. All the simulation results are consistent with many literatures. iPCCFD can minimize the numerical instabilities and noises along with the additional diffusion terms. All of our studies present relatively small numerical errors without employing any divergence-free reconstruction. In particular, we obtain fairly stable results in the two-dimensional Brio–Wu shock tube problem which well conserves ∇⋅B=0 throughout the simulation. The ERE technique removes the accumulation of roundoff errors in the uniform or non-disturbed system. We have also shown that iPCCFD is characterized by the high order of accuracy and the low numerical dissipation in the circularly polarized Alfvén wave tests. The proposed iPCCFD scheme is a parallel-efficient and high precision numerical scheme for solving the MHD equations in hyperbolic conservation systems.

Combined high-order algorithms in robust least-squares estimation with harmonic regressor and strictly diagonally dominant information matrix

This article describes new high-order algorithms in the least-squares problem with harmonic regressor and strictly diagonally dominant information matrix. Estimation accuracy and the number of steps to achieve this accuracy are controllable in these algorithms. Simplified forms of the high-order matrix inversion algorithms and the high-order algorithms of direct calculation of the parameter vector are found. The algorithms are presented as recursive procedures driven by estimation errors multiplied by the gain matrices, which can be seen as preconditioners. A simple and recursive (with respect to order) algorithm for update of the gain matrix, which is associated with Neumann series, is found. It is shown that the limiting form of the algorithm (algorithm of infinite order) provides perfect estimation. A new form of the gain matrix is also a basis for unification method of high-order algorithms. New combined and fast convergent high-order algorithms of recursive matrix inversion and algorithms of direct calculation of the parameter vector are presented. The stability of algorithms is proved and explicit transient bound on estimation error is calculated. New algorithms are simple, fast and robust with respect to round-off error accumulation.

An efficient implementation of the ensemble Kalman filter based on an iterative Sherman–Morrison formula

We present a practical implementation of the ensemble Kalman filter (EnKF) based on an iterative Sherman–Morrison formula. The new direct method exploits the special structure of the ensemble-estimated error covariance matrices in order to efficiently solve the linear systems involved in the analysis step of the EnKF. The computational complexity of the proposed implementation is equivalent to that of the best EnKF implementations available in the literature when the number of observations is much larger than the number of ensemble members, as typically is case in practice. Moreover, the proposed method provides the best theoretical complexity when it is compared to generic formulations of matrix inversion based on the Sherman–Morrison formula. The stability analysis of the proposed method is carried out and a pivoting strategy is discussed in order to reduce the accumulation of round-off errors without increasing the computational effort. A parallel implementation is discussed as well. Computational experiments carried out using an oceanic quasi-geostrophic model reveal that the proposed algorithm yields the same accuracy as other EnKF implementations, but scales better with regard to the number of observations.