Forward Error Analysis Research Articles

Accurate polynomial fitting and evaluation via Arnoldi

Given sample data points $ \{(x_j,f_j)\}_{j = 1}^N $, in [Brubeck, Nakatsukasa, and Trefethen, SIAM Review, 63 (2021), pp. 405-415], an Arnoldi-based procedure is proposed to accurately evaluate the best fitting polynomial, in the least squares sense, at new nodes $ \{s_j\}_{j = 1}^M $, based on the Vandermonde basis. Numerical tests indicated that this procedure can in general achieve high accuracy. The main purpose of this paper is to perform a forward rounding error analysis in finite precision. Our result establishes sensitivity factors regarding the accuracy of the algorithm, and provides a theoretical justification for why the algorithm works. For least-squares approximation on an interval, we propose a variant of this Arnoldi-based evaluation by using the Chebyshev polynomial basis. Numerical tests are reported to demonstrate our forward rounding error analysis.

Sensitivity Analysis of the Data Assimilation-Driven Decomposition in Space and Time to Solve PDE-Constrained Optimization Problems

This paper is presented in the context of sensitivity analysis (SA) of large-scale data assimilation (DA) models. We studied consistency, convergence, stability and roundoff error propagation of the reduced-space optimization technique arising in parallel 4D Variational DA problems. The results are helpful to understand the reliability of DA, to assess what confidence one can have that the simulation results are correct and to determine its configuration in any application. The main contributions of the present work are as follows. By using forward error analysis, we derived the number of conditions of the parallel approach. We found that the parallel approach reduces the number of conditions, revealing that it is more appropriate than the standard approach usually implemented in most operative software. As the background values are used as initial conditions of local PDE models, we analyzed stability with respect to time direction. Finally, we proved consistency of the proposed approach by analyzing local truncation errors of each computational kernel.

Two Taylor Algorithms for Computing the Action of the Matrix Exponential on a Vector

The action of the matrix exponential on a vector eAtv, A∈Cn×n, v∈Cn, appears in problems that arise in mathematics, physics, and engineering, such as the solution of systems of linear ordinary differential equations with constant coefficients. Nowadays, several state-of-the-art approximations are available for estimating this type of action. In this work, two Taylor algorithms are proposed for computing eAv, which make use of the scaling and recovering technique based on a backward or forward error analysis. A battery of highly heterogeneous test matrices has been used in the different experiments performed to compare the numerical and computational properties of these algorithms, implemented in the MATLAB language. In general, both of them improve on those already existing in the literature, in terms of accuracy and response time. Moreover, a high-performance computing version that is able to take advantage of the computational power of a GPU platform has been developed, making it possible to tackle high dimension problems at an execution time significantly reduced.

An Improved Taylor Algorithm for Computing the Matrix Logarithm

The most popular method for computing the matrix logarithm is a combination of the inverse scaling and squaring method in conjunction with a Padé approximation, sometimes accompanied by the Schur decomposition. In this work, we present a Taylor series algorithm, based on the free-transformation approach of the inverse scaling and squaring technique, that uses recent matrix polynomial formulas for evaluating the Taylor approximation of the matrix logarithm more efficiently than the Paterson–Stockmeyer method. Two MATLAB implementations of this algorithm, related to relative forward or backward error analysis, were developed and compared with different state-of-the art MATLAB functions. Numerical tests showed that the new implementations are generally more accurate than the previously available codes, with an intermediate execution time among all the codes in comparison.

Towards Lower Precision Adaptive Filters: Facts From Backward Error Analysis of RLS

Lower precision arithmetic can improve the throughput of adaptive filters, while requiring less hardware resources and less power. Such benefits are crucial for adaptive filters, especially for IoT and wearable applications. In order to apply lower precision arithmetic to adaptive filters, a clear rounding error analysis framework is required, since lower precision arithmetic can degrade the filter performance. Previously, rounding error analyses of adaptive filters were based on forward error analysis. This limited the descriptiveness of rounding error impact on adaptive filter performance in relation to other external variables such as measurement noise, regularisation, and numerical stability of an algorithm. To overcome such limitations, we first present a new backward error analysis framework for adaptive Recursive Least Squares (RLS) filters. Our framework transforms finite precision arithmetic adaptive filters into exact arithmetic adaptive filters with the input data corrupted by rounding error noise that is additive to measurement noise. Findings throughout our backward error analysis framework can provide a guide on how to apply lower precision arithmetic to adaptive filters: (i) the magnitudes of the rounding error noise depend on the numerical stability of the implementation algorithm, arithmetic precision, and regularisation, (ii) the rounding error noise is independently additive to measurement noise, (iii) a higher regularisation is recommended for lower precision arithmetic adaptive filters, and (iv) adaptive filters using lower precision arithmetic have equivalent filter performance to those using higher precision if the magnitudes of rounding error noise are lower than measurement noise.

High-Performance Computation of the Exponential of a Large Sparse Matrix

Computation of the large sparse matrix exponential has been an important topic in many fields, such as network and finite-element analysis. The existing scaling and squaring algorithm (SSA) is not suitable for the computation of the large sparse matrix exponential as it requires greater memories and computational cost than is actually needed. By introducing two novel concepts, i.e., real bandwidth and bandwidth, to measure the sparsity of the matrix, the sparsity of the matrix exponential is analyzed. It is found that for every matrix computed in the squaring phase of the SSA, a corresponding sparse approximate matrix exists. To obtain the sparse approximate matrix, a new filtering technique in terms of forward error analysis is proposed. Combining the filtering technique with the idea of keeping track of the incremental part, a competitive algorithm is developed for the large sparse matrix exponential. The proposed method can primarily alleviate the over-scaling problem due to the filtering technique. Three sets of numerical experiments, including one large matrix with a dimension larger than 2e6 , are conducted. The numerical experiments show that, compared with the expm function in MATLAB, the proposed algorithm can provide higher accuracy at lower computational cost and with less memory.

Compensated Evaluation of Tensor Product Surfaces in CAGD

In computer-aided geometric design, a polynomial surface is usually represented in Bézier form. The usual form of evaluating such a surface is by using an extension of the de Casteljau algorithm. Using error-free transformations, a compensated version of this algorithm is presented, which improves the usual algorithm in terms of accuracy. A forward error analysis illustrating this fact is developed.

Simulation of harmonic oscillators on the lattice

This work deals with the simulation of a two‐dimensional ideal lattice having simple tetragonal geometry. The harmonic character of the oscillators give rise to a system of second‐order linear differential equations, which can be recast into matrix form. The explicit solutions which govern the dynamics of this system can be expressed in terms of matrix trigonometric functions. For the derivation we employ the Lagrangian formalism to determine the correct solutions, which extremize the underlying action of the system. In the numerical evaluation we develop diverse state‐of‐the‐art algorithms which efficiently tackle equations with matrix sine and cosine functions. For this purpose, we introduce two special series related to trigonometric functions. They provide approximate solutions of the system through a suitable combination. For the final computation an algorithm based on Taylor expansion with forward and backward error analysis for computing those series had to be devised. We also implement several MATLAB programs which simulate and visualize the two‐dimensional lattice and check its energy conservation.

Accurate solutions of weighted least squares problems associated with rank-structured matrices

In this paper, we consider how to accurately solve the weighted least squares (WLS) problem associated with a class of rank-structured matrices admitting bidiagonal representations (BRs) from weighted polynomial regression. We develop a new algorithm to solve the structured WLS problem, provided that BRs are available. A mechanism is exploited to guarantee that all the solution components are computed with a desirable high accuracy. Forward error analysis and numerical experiments are performed to confirm the high accuracy. In particular, our random examples demonstrate the high relative accuracy of each computed solution component, independently of the ill-conditioning of weighted matrices.

Compensated de Casteljau algorithm in K times the working precision

In computer aided geometric design a polynomial is usually represented in Bernstein form. This paper presents a family of compensated algorithms to accurately evaluate a polynomial in Bernstein form with floating point coefficients. The principle is to apply error-free transformations to improve the traditional de Casteljau algorithm. At each stage of computation, round-off error is passed on to first order errors, then to second order errors, and so on. After the computation has been “filtered” (K−1) times via this process, the resulting output is as accurate as the de Casteljau algorithm performed in K times the working precision. Forward error analysis and numerical experiments illustrate the accuracy of this family of algorithms.

Fast Taylor polynomial evaluation for the computation of the matrix cosine

In this work we introduce a new method to compute the matrix cosine. It is based on recent new matrix polynomial evaluation methods for the Taylor approximation and a mixed forward and backward error analysis. The matrix polynomial evaluation methods allow to evaluate the Taylor polynomial approximation of the matrix cosine function more efficiently than using Paterson–Stockmeyer method. A sequential Matlab implementation of the new algorithm is provided, giving better efficiency and accuracy than state-of-the-art algorithms. Moreover, we provide an implementation in Matlab that can use NVIDIA GPUs easily and efficiently.

A Truncated Taylor Series Algorithm for Computing the Action of Trigonometric and Hyperbolic Matrix Functions

A new algorithm is derived for computing the action $f(tA)B$, where $A$ is an $n\times n$ matrix, $B$ is $n\times n_0$ with $n_0 \ll n$, and $f$ is cosine, sinc, sine, hyperbolic cosine, hyperbolic sinc, or hyperbolic sine function. In the case of $f$ being even, the computation of $f(tA^{1/2})B$ is possible without explicitly computing $A^{1/2}$, where $A^{1/2}$ denotes any matrix square root of $A$. The algorithm offers six independent output options given $t$, $A$, $B$, and a tolerance. For each option, actions of a pair of trigonometric or hyperbolic matrix functions are simultaneously computed. The algorithm scales the matrix $A$ down by a positive integer $s$, approximates $f(s^{-1}tA^\sigma)B$, where $\sigma$ is either 1 or 1/2, by using a truncated Taylor series, and finally uses the recurrences of the Chebyshev polynomials of the first and second kind to recover $f(tA^\sigma)B$. The selection of the scaling parameter and the degree of Taylor polynomial is based on a forward error analysis and a sequence of the form $\|A^k\|^{1/k}$ in such a way that the overall computational cost of the algorithm is minimized. Shifting is used where applicable as a preprocessing step to reduce the scaling parameter. The algorithm works for any matrix $A$, and its computational cost is dominated by the formation of products of $A$ with $n\times n_0$ matrices that could take advantage of the implementation of level-3 BLAS. Our numerical experiments show that the new algorithm behaves in a forward stable fashion and in most problems outperforms the existing algorithms in terms of CPU time, computational cost, and accuracy.

Two algorithms for computing the matrix cosine function

The computation of matrix trigonometric functions has received remarkable attention in the last decades due to its usefulness in the solution of systems of second order linear differential equations. Several state-of-the-art algorithms have been provided recently for computing these matrix functions. In this work, we present two efficient algorithms based on Taylor series with forward and backward error analysis for computing the matrix cosine. A MATLAB implementation of the algorithms is compared to state-of-the-art algorithms, with excellent performance in both accuracy and cost.

Retrieval of ice cloud properties using an optimal estimation algorithm and MODIS infrared observations. Part I: Forward model, error analysis, and information content.

An optimal estimation (OE) retrieval method is developed to infer three ice cloud properties simultaneously: optical thickness (τ), effective radius (reff ), and cloud-top height (h). This method is based on a fast radiative transfer (RT) model and infrared (IR) observations from the MODerate resolution Imaging Spectroradiometer (MODIS). This study conducts thorough error and information content analyses to understand the error propagation and performance of retrievals from various MODIS band combinations under different cloud/atmosphere states. Specifically, the algorithm takes into account four error sources: measurement uncertainty, fast RT model uncertainty, uncertainties in ancillary datasets (e.g., atmospheric state), and assumed ice crystal habit uncertainties. It is found that the ancillary and ice crystal habit error sources dominate the MODIS IR retrieval uncertainty and cannot be ignored. The information content analysis shows that, for a given ice cloud, the use of four MODIS IR observations is sufficient to retrieve the three cloud properties. However, the selection of MODIS IR bands that provide the most information and their order of importance varies with both the ice cloud properties and the ambient atmospheric and the surface states. As a result, this study suggests the inclusion of all MODIS IR bands in practice since little a priori information is available.

Global error analysis and inertial manifold reduction

Four types of global error for initial value problems are considered in a common framework. They include classical forward error analysis and shadowing error analysis together with extensions of both to include rescaling of time. To determine the amplification of the local error that bounds the global error we present a linear analysis similar in spirit to condition number estimation for linear systems of equations. We combine these ideas with techniques for dimension reduction of differential equations via a boundary value formulation of numerical inertial manifold reduction. These global error concepts are exercised to illustrate their utility on the Lorenz equations and inertial manifold reductions of the Kuramoto–Sivashinsky equation.

A Note on the Perturbation of arithmetic expressions

In this paper we present the theoretical foundation of forward error analysis of numerical algorithms under;• Approximations in "built-in" functions.• Rounding errors in arithmetic floating-point operations.• Perturbations of data.The error analysis is based on linearization method. The fundamental tools of the forward error analysis are system of linear absolute and relative a prior and a posteriori error equations and associated condition numbers constituting optimal of possible cumulative round – off errors. The condition numbers enable simple general, quantitative bounds definitions of numerical stability. The theoretical results have been applied a Gaussian elimination, and have proved to be very effective means of both a priori and a posteriori error analysis.

A Note on the Perturbation of arithmetic expressions

In this paper we present the theoretical foundation of forward error analysis of numerical algorithms under;• Approximations in "built-in" functions.• Rounding errors in arithmetic floating-point operations.• Perturbations of data.The error analysis is based on linearization method. The fundamental tools of the forward error analysis are system of linear absolute and relative a prior and a posteriori error equations and associated condition numbers constituting optimal of possible cumulative round – off errors. The condition numbers enable simple general, quantitative bounds definitions of numerical stability. The theoretical results have been applied a Gaussian elimination, and have proved to be very effective means of both a priori and a posteriori error analysis.

Explicit MPC: Hard constraint satisfaction under low precision arithmetic

MPC is becoming increasingly implemented on embedded systems, where low precision computation is preferred either to reduce costs, speedup execution or reduce power consumption. However, in a low precision implementation, constraint satisfaction cannot be guaranteed. To enforce constraint satisfaction under numerical errors, we adopt tools from forward error analysis to compute an error bound on the output of the embedded controller. We treat this error as a state disturbance and use it to inform the design of a constraint-tightening robust controller. The technique is validated via a practical implementation on an FPGA evaluation board.

Automated backward error analysis for numerical code

Numerical code uses floating-point arithmetic and necessarily suffers from roundoff and truncation errors. Error analysis is the process to quantify such uncertainty in the solution to a problem. Forward error analysis and backward error analysis are two popular paradigms of error analysis. Forward error analysis is more intuitive and has been explored and automated by the programming languages (PL) community. In contrast, although backward error analysis is more preferred by numerical analysts and the foundation for numerical stability, it is less known and unexplored by the PL community. To fill the gap, this paper presents an automated backward error analysis for numerical code to empower both numerical analysts and application developers. In addition, we use the computed backward error results to also compute the condition number, an important quantity recognized by numerical analysts for measuring how sensitive a function is to changes or errors in the input. Experimental results on Intel X87 FPU functions and widely-used GNU C Library functions demonstrate that our analysis is effective at analyzing the accuracy of floating-point programs.

Recent Patents on Rotating Double-Prism Scanning Mechanism

In recent years, many double-prism mechanism patents have become promising solutions for precision engineering applications, which can perform superior optical pointing and tracking functions. The paper focuses on surveying recent patents related to rotating double-prism scanner used in the fields of free-space laser communications, laser radar, optical switches, laser scanning, and other precision engineering. A general overview of various rotating double-prism mechanisms is provided, and various aspects of rotating double-prism researches are introduced, including forward problem, inverse problem, drive mechanism, mechanical analysis, error analysis and control system. Patents related to rotating double-prism mechanisms, issued since 2010, are selected for the survey. With a comparison of these patents, some valuable conclusions are drawn to predict the future research and development in rotating double-prism scanner. Keywords: Beam scanning, control system, mechanical structure, rotating double prisms, scanning mechanism.