Backward Error Analysis Research Articles

Estimates of the accuracy of numerical solutions using regularization

This article explores the accuracy of numerical solutions, and suggests methods for analyzing accuracy depending on the properties of the problem. Numerical studies of complex expensive objects of technology and physics require that the computational results be obtained with guaranteed accuracy. It also depends on the fact that in the work of technical objects there are large intervals of operation time not observed experimentally. Therefore, there is a need to describe the location of the observed and calculated values, as well as the accuracy with which they are calculated. The effect of strong growth in estimates of error bounds is manifested for a large number of methods used to estimate the error of a numerical solution. This means the lack of correctness of algorithms for evaluating the accuracy of numerical solutions due to the failure of the stability conditions with respect to perturbations of the right-hand side. For many problems, among all the algorithms, the backward analysis of errors turned out to be the most effective method for assessing the accuracy of numerical solutions. The backward analysis of errors consists in the fact that when assessing the accuracy (error) of a numerical solution, the numerical solution is considered as an exact solution to a problem close to the original problem. A backward error analysis was proposed and developed in the algorithms of J Wilkinson in the context of the numerical solution of problems of linear algebra, and in the algorithms of V V Voevodin, who widely distributed it to many areas of numerical analysis. In the framework of the backward error analysis, the regularization of the algorithm for estimating the error of a numerical solution is reasonably applied. This article explores methods for the backward analysis of errors of numerical solutions.

On pole-swapping algorithms for the eigenvalue problem

Pole-swapping algorithms, which are generalizations of the QZ algorithm for the generalized eigenvalue problem, are studied. A new modular (and therefore more flexible) convergence theory that applies to all pole-swapping algorithms is developed. A key component of all such algorithms is a procedure that swaps two adjacent eigenvalues in a triangular pencil. An improved swapping routine is developed, and its superiority over existing methods is demonstrated by a backward error analysis and numerical tests. The modularity of the new convergence theory and the generality of the pole-swapping approach shed new light on bi-directional chasing algorithms, optimally packed shifts, and bulge pencils, and allow the design of novel algorithms.

Uniform-in-time weak error analysis for stochastic gradient descent algorithms via diffusion approximation

Diffusion approximation provides weak approximation for stochastic gradient descent algorithms in a finite time horizon. In this paper, we introduce new tools motivated by the backward error analysis of numerical stochastic differential equations into the theoretical framework of diffusion approximation, extending the validity of the weak approximation from finite to infinite time horizon. The new techniques developed in this paper enable us to characterize the asymptotic behavior of constant-step-size SGD algorithms for strongly convex objective functions, a goal previously unreachable within the diffusion approximation framework. Our analysis builds upon a truncated formal power expansion of the solution of a stochastic modified equation arising from diffusion approximation, where the main technical ingredient is a uniform-in-time weak error bound controlling the long-term behavior of the expansion coefficient functions near the global minimum. We expect these new techniques to greatly expand the range of applicability of diffusion approximation to cover wider and deeper aspects of stochastic optimization algorithms in data science.

Discrete-Time Control Design Based on Symplectic Integration: Linear Systems

We propose a novel approach to design discrete-time state feedback controllers for sampled control systems with guaranteed stability under arbitrary sampling times that fulfill the Nyquist-Shannon condition. The key idea is borrowed from symplectic integration and backward error analysis of Hamiltonian systems: The closed-loop target system is the discretized version of a continuous-time system with appropriately shaped Hamiltonian, where a symplectic scheme is used for discretization. We adopt this argumentation for a systematic discrete-time design procedure for controllable linear systems based on the implicit midpoint rule. We motivate the approach on the example of two basic linear systems under zero order hold sampling, we show the construction of target systems with desired eigenvalues based on the discrete-time controller canonical form, and we illustrate the quality of the main result on a random sixth order system.

Sharper Probabilistic Backward Error Analysis for Basic Linear Algebra Kernels with Random Data

Standard backward error analyses for numerical linear algebra algorithms provide worst-case bounds that can significantly overestimate the backward error. Our recent probabilistic error analysis, which assumes rounding errors to be independent random variables [SIAM J. Sci. Comput., 41 (2019), pp. A2815--A2835], contains smaller constants but its bounds can still be pessimistic. We perform a new probabilistic error analysis that assumes both the data and the rounding errors to be random variables and assumes only mean independence. We prove that for data with zero or small mean we can relax the existing probabilistic bounds of order $\sqrt{n}\mkern1muu$ to much sharper bounds of order $u$, which are independent of $n$. Our fundamental result is for summation and we use it to derive results for inner products, matrix--vector products, and matrix--matrix products. The analysis answers the open question of why random data distributed on $[-1,1]$ leads to smaller error growth for these kernels than random data distributed on [0,1]. We also propose a new algorithm for multiplying two matrices that transforms the rows of the first matrix to have zero mean and we show that it can achieve significantly more accurate results than standard matrix multiplication.

Quasi‐shuffle algebras and renormalisation of rough differential equations

The objective of this work is to compare several approaches to the process of renormalisation in the context of rough differential equations using the substitution bialgebra on rooted trees known from backward error analysis of $B$-series. For this purpose, we present a so-called arborification of the Hoffman--Ihara theory of quasi-shuffle algebra automorphisms. The latter are induced by formal power series, which can be seen to be special cases of the cointeraction of two Hopf algebra structures on rooted forests. In particular, the arborification of Hoffman's exponential map, which defines a Hopf algebra isomorphism between the shuffle and quasi-shuffle Hopf algebra, leads to a canonical renormalisation that coincides with Marcus' canonical extension for semimartingale driving signals. This is contrasted with the canonical geometric rough path of Hairer and Kelly by means of a recursive formula defined in terms of the coaction of the substitution bialgebra.

Contact variational integrators

We present geometric numerical integrators for contact flows that stem from a discretization of Herglotz’ variational principle. First we show that the resulting discrete map is a contact transformation and that any contact map can be derived from a variational principle. Then we discuss the backward error analysis of our variational integrators, including the construction of a modified Lagrangian. Throughout the paper we use the damped harmonic oscillator as a benchmark example to compare our integrators to their symplectic analogues.

Perturbation analysis for matrix joint block diagonalization

The matrix joint block-diagonalization problem (jbdp) of a given matrix set A={Ai}i=1m is about finding a nonsingular matrix W such that all WTAiW are block-diagonal. It includes the matrix joint diagonalization problem (jdp) as a special case for which all WTAiW are required diagonal. Generically, such a matrix W may not exist, but there are practical applications such as multidimensional independent component analysis (MICA) for which it does exist under the ideal situation, i.e., no noise is presented. However, in practice noises do get in and, as a consequence, the matrix set is only approximately block-diagonalizable, i.e., one can only make all W˜TAiW˜ nearly block-diagonal at best, where W˜ is an approximation to W, obtained usually by computation. The main goal of this paper is to develop a perturbation theory for jbdp to address, among others, the question: how accurate this W˜ is. Previously such a theory for jdp has been discussed, but no effort had been attempted for jbdp because, in large part, there is no quantitative way to describe solution uniqueness of jbdp until 2017 when Cai and Liu (2017) [9] successfully obtained a necessary and sufficient uniqueness condition. Based on the condition, in this article, we will establish a perturbation theory for jbdp. Our main contributions include an error bound for the approximate block-diagonalizer W˜ and a backward error analysis for jbdp. Numerical tests are presented to validate the theoretical results.

Error estimates for semi-Lagrangian finite difference methods applied to Burgers' equation in one dimension

We make an analytic study of the diffusive, dispersive and overall errors, which arise when using semi-implicit semi-Lagrangian (SISL) finite difference methods to approximate those travelling wave solutions of the one-dimensional Burgers' equation with small diffusion, which develop sharp fronts. For the case of a fixed uniform spatial mesh, with piecewise linear interpolation, a backward error analysis approach is used to construct a precise formal analytic description of the front profile of the numerical approximation to this solution. From this description it is possible to obtain precise estimates of the front width and the front speed in terms of the spatial and temporal step size and to express the overall solution error in terms of these. These formal estimates agree closely with numerical calculations, both qualitatively and quantitatively, and display a roughly periodic behaviour as the number Nx of mesh points increases, and the CFL number passes through integer values. In particular, they show that despite the otherwise poor resolution of the method, the front width is closely approximated when the CFL number is close to an integer, and the front speed is closely approximated when it is close to a half integer. The overall L2 error also shows super-convergence for certain values of Nx. This possibly motivates doing two calculations with different Nx when using the SISL method on such problems to separately minimise the diffusive and dispersive errors. Similar errors in the front width and speed are observed for a number of different interpolation schemes with and without flux limiters.

Backward error analysis for linearizations in heavily damped quadratic eigenvalue problem

SummaryHeavily damped quadratic eigenvalue problem (QEP) is a special type of QEPs. It has a large gap between small and large eigenvalues in absolute value. One common way for solving QEP is to linearize the original problem via linearizations. Previous work on the accuracy of eigenpairs of not heavily damped QEP focuses on analyzing the backward error of eigenpairs relative to linearizations. The objective of this paper is to explain why different linearizations lead to different errors when computing small and large eigenpairs. To obtain this goal, we bound the backward error of eigenpairs relative to the linearization methods. Using these bounds, we build upper bounds of growth factors for the backward error. We present results of numerical experiments that support the predictions of the proposed methods.

The e-MoM approach for approximating matrix functionals

The approximation of matrix functionals appears in many applications arising from the fields of mathematics, statistics, mechanics, networks, machine learning and physics. In this paper, we estimate matrix functionals of the form XTf(A)Y, where A∈Rp×p is a given diagonalizable matrix, X, Y∈Rp×k are skinny “block vectors” with k≪p columns and f is a smooth function defined on the spectrum of the matrix A. We apply a direct approach based on the extrapolation of the moments of the given matrix, for estimating this kind of matrix functionals. This approach avoids the application of the polarization identity, is fairly inexpensive and leads to a stable procedure. Moreover, we develop a detailed backward error analysis for the derived estimates. Several numerical results illustrating the effectiveness of the direct method are presented and concrete classes of matrices, suitable for the extrapolation estimates, are proposed.

Refinement of thermostated molecular dynamics using backward error analysis.

Kinetic energy equipartition is a premise for many deterministic and stochastic molecular dynamics methods that aim at sampling a canonical ensemble. While this is expected for real systems, discretization errors introduced by the numerical integration may lead to deviations from equipartition. Fortunately, backward error analysis allows us to obtain a higher-order estimate of the quantity that is actually subject to equipartition. This is related to a shadow Hamiltonian, which converges to the specified Hamiltonian only when the time-step size approaches zero. This paper deals with discretization effects in a straightforward way. With a small computational overhead, we obtain refined versions of the kinetic and potential energies, whose sum is a suitable estimator of the shadow Hamiltonian. Then, we tune the thermostatting procedure by employing the refined kinetic energy instead of the conventional one. This procedure is shown to reproduce a canonical ensemble compatible with the refined system, as opposed to the original one, but canonical averages regarding the latter can easily be recovered by reweighting. Water, modeled as a rigid body, is an excellent test case for our proposal because its numerical stability extends up to time steps large enough to yield pronounced discretization errors in Verlet-type integrators. By applying our new approach, we were able to mitigate discretization effects in equilibrium properties of liquid water for time-step sizes up to 5 fs.

Weak backward error analysis for stochastic Hamiltonian Systems

This paper presents a study of the approximation error corresponding to a symplectic scheme of weak order one for a stochastic autonomous Hamiltonian system. A backward error analysis is done at the level of the Kolmogorov equation associated with the initial stochastic Hamiltonian system. An expansion of the weak error and expansions of the ergodic averages and of the invariant measures associated with the numerical scheme are obtained in terms of powers of the discretization step size and the solutions of the modified Kolmogorov equation.

Vector Spaces of Generalized Linearizations for Rectangular Matrix Polynomials

The complete eigenvalue problem associated with a rectangular matrix polynomial is typically solved via the technique of linearization. This work introduces the concept of generalized linearizations of rectangular matrix polynomials. For a given rectangular matrix polynomial, it also proposes vector spaces of rectangular matrix pencils with the property that almost every pencil is a generalized linearization of the matrix polynomial which can then be used to solve the complete eigenvalue problem associated with the polynomial. The properties of these vector spaces are similar to those introduced in the literature for square matrix polynomials and in fact coincide with them when the matrix polynomial is square. Further, almost every pencil in these spaces can be `trimmed' to form many smaller pencils that are strong linearizations of the matrix polynomial which readily yield solutions of the complete eigenvalue problem for the polynomial. These linearizations are easier to construct and are often smaller than the Fiedler linearizations introduced in the literature for rectangular matrix polynomials. Additionally, a global backward error analysis applied to these linearizations shows that they provide a wide choice of linearizations with respect to which the complete polynomial eigenvalue problem can be solved in a globally backward stable manner.

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.

Perturbation analysis for palindromic and anti-palindromic nonlinear eigenvalue problems

A structured backward error analysis for an approximate eigenpair of structured nonlinear matrix equations with T-palindromic, H-palindromic, T-anti-palindromic, and H-anti-palindromic structures is conducted. We construct a minimal structured perturbation in the Frobenius norm such that an approximate eigenpair becomes an exact eigenpair of an appropriately perturbed nonlinear matrix equation. The present work shows that our general framework extends existing results in the literature on the perturbation theory of matrix polynomials.

Backward error analysis of polynomial approximations for computing the action of the matrix exponential

We describe how to perform the backward error analysis for the approximation of $$\exp (A)v$$ by $$p(s^{-1}A)^sv$$ , for any given polynomial p(x). The result of this analysis is an optimal choice of the scaling parameter s which assures a bound on the backward error, i.e. the equivalence of the approximation with the exponential of a slightly perturbed matrix. Thanks to the SageMath package expbea we have developed, one can optimize the performance of the given polynomial approximation. On the other hand, we employ the package for the analysis of polynomials interpolating the exponential function at so called Leja–Hermite points. The resulting method for the action of the matrix exponential can be considered an extension of both Taylor series approximation and Leja point interpolation. We illustrate the behavior of the new approximation with several numerical examples.

On-the-fly backward error estimate for matrix exponential approximation by Taylor algorithm

In this paper we show that it is possible to estimate the backward error for the approximation of the matrix exponential on-the-fly , without the need to precompute in high precision quantities related to specific accuracies. In this way, the scaling parameter s and the degree m of the truncated Taylor series (the underlying method) are adapted to the input matrix and not to a class of matrices sharing some spectral properties, such as with the classical backward error analysis. The result is a very flexible method which can approximate the matrix exponential at any desired accuracy. Moreover, the risk of overscaling, that is the possibility to select a value s larger than necessary, is mitigated by a new choice as the sum of two powers of two. Finally, several numerical experiments in MATLAB with data in double and variable precision and with different tolerances confirm that the method is accurate and often faster than available good alternatives.

Structured backward error analysis of linearized structured polynomial eigenvalue problems

We introduce a new class of structured matrix polynomials, namely, the class of M_A-structured matrix polynomials, to provide a common framework for many classes of structured matrix polynomials that are important in applications: the classes of (skew-)symmetric, (anti-)palindromic, and alternating matrix polynomials. Then, we introduce the families of M_A-structured strong block minimal bases pencils and of M_A-structured block Kronecker pencils,, and show that any M_A-structured odd-degree matrix polynomial can be strongly linearized via an M_A-structured block Kronecker pencil. Finally, for the classes of (skew-)symmetric, (anti-)palindromic, and alternating odd-degree matrix polynomials, the M_A-structured framework allows us to perform a global and structured backward stability analysis of complete structured polynomial eigenproblems, regular or singular, solved by applying to an M_A-structured block Kronecker pencil a structurally backward stable algorithm that computes its complete eigenstructure, like the palindromic-QR algorithm or the structured versions of the staircase algorithm. This analysis allows us to identify those M_A-structured block Kronecker pencils that yield a computed complete eigenstructure which is the exact one of a slightly perturbed structured matrix polynomial.These pencils include (modulo permutations) the well-known block-tridiagonal and block-antitridiagonal structure-preserving linearizations. Our analysis incorporates structure to the recent (unstructured) backward error analysis performed for block Kronecker linearizations by Dopico, Lawrence, P\'erez and Van Dooren, and share with it its key features, namely, it is a rigorous analysis valid for finite perturbations, i.e., it is not a first order analysis, it provides precise bounds, and it is valid simultaneously for a large class of structure-preserving strong linearizations.

Block Kronecker linearizations of matrix polynomials and their backward errors

We introduce a new family of linearizations of matrix polynomials---which we call Kronecker pencils---and perform a backward stability analysis of complete polynomial eigenproblems. These problems are solved by applying any backward stable algorithm to a block Kronecker pencil, such as the staircase algorithm for singular pencils or the QZ algorithm for regular pencils. This analysis allows us to identify those block Kronecker pencils that yield a computed complete eigenstructure which is exactly that of a slightly perturbed matrix polynomial. These favorable pencils include the famous Fiedler linearizations, which are just a very particular case of block Kronecker pencils. Thus, our analysis offers the first proof available in the literature of global backward stability for Fiedler pencils. In addition, the theory developed for block Kronecker pencils is much simpler than the theory available for Fiedler especially in the case of rectangular matrix polynomials. The global backward error analysis in this work presents for the first time the following key properties: it is a rigorous analysis valid for finite perturbations (i.e., it is not a first order analysis), it provides precise bounds, it is valid simultaneously for a large class of linearizations, and it establishes a framework that may be generalized to other classes of linearizations. These features are related to the fact that block Kronecker pencils are a particular case of the new family of strong block minimal bases pencils, which include certain perturbations of block Kronecker pencils; this allows us to extend the results in this paper to other contexts.