Backward Error Research Articles

Backward error analysis and inverse eigenvalue problems for Hankel and Symmetric-Toeplitz structures

This work deals with the study of structured backward error analysis of Hankel and symmetric-Toeplitz matrix pencils. These structured matrix pencils belong to the class of symmetric matrix pencils with some additional properties that a symmetric matrix pencil does not have in general. The perturbation analysis of these two structures is discussed one by one to depict the additional properties explicitly. Present work shows the entrywise structured perturbation of matrix pencils in Frobenius norm such that the specified eigenpairs become exact eigenpairs of an appropriately perturbed matrix pencil. The framework used here maintains the sparsity in the perturbation of the above-structured matrix pencils. Further, the backward error results help for solving a variety of inverse eigenvalue problems.

Backward error analysis of linearizing–balancing strategies for heavily damped quadratic eigenvalue problem

A classical approach for solving quadratic eigenvalue problem (QEP) is via linearization. However, it can suffer from numerical instability when the norms of coefficient matrices vary widely. Two strategies of Betcke’s balancing for heavily damped QEP are considered. One strategy is to balance matrix polynomial before linearizing (called balancing–linearizing for short). The other strategy is first to linearize the matrix polynomial, then balancing (called linearizing–balancing for short). We analyze the backward error of approximate eigenpairs computed by these two strategies, and find that the backward error of approximate eigenpairs by linearizing–balancing is smaller than that by balancing–linearizing under relatively mild conditions. Numerical experiments are presented to demonstrate the advantages of linearizing-balancing.

Solving block low-rank linear systems by LU factorization is numerically stable

Abstract Block low-rank (BLR) matrices possess a blockwise low-rank property that can be exploited to reduce the complexity of numerical linear algebra algorithms. The impact of these low-rank approximations on the numerical stability of the algorithms in floating-point arithmetic has not previously been analysed. We present rounding error analysis for the solution of a linear system by LU factorization of BLR matrices. Assuming that a stable pivoting scheme is used, we prove backward stability: the relative backward error is bounded by a modest constant times $\varepsilon $, where the low-rank threshold $\varepsilon $ is the parameter controlling the accuracy of the blockwise low-rank approximations. In addition to this key result, our analysis offers three new insights into the numerical behaviour of BLR algorithms. First, we compare the use of a global or local low-rank threshold and find that a global one should be preferred. Second, we show that performing intermediate recompressions during the factorization can significantly reduce its cost without compromising numerical stability. Third, we consider different BLR factorization variants and determine the update–compress–factor variant to be the best. Tests on a wide range of matrices from various real-life applications show that the predictions from the analysis are realized in practice.

Chaos Theory and ARTFA: Emerging Tools for Interpreting ECG Signals to Diagnose Cardiac Arrhythmias

Timely detection of cardiac abnormalities from an Electrocardiogram (ECG) signal is very essential. This requires its appropriate and efficient processing. In the literature, most of the researchers focussed on linear techniques that were applied on filtered ECG datasets leaving an ample scope for exploring the use of non-linear techniques in the presence and absence of natural noise-processes. Therefore, there is a need of supplementing the existing research on ECG signal interpretation by using non-linear techniques on noisy ECG data. Non-linear techniques are expected to yield supplementary clues about the non-linearities in the underlying cardiovascular system. One such promising non-linear technique, known as chaos theory (analysis), has been considered here to estimate reliable and robust thresholds for R-peak detection using fractal dimension, Approximate Entropy, Sample Entropy, correlation dimension, and Lyapunov exponent based on time-delay dimension (embedding). Also, time–frequency analysis techniques have shown their effectiveness for analyzing such types of non-linear and non-stationary signals due to simultaneous interpretation of the signal in both time and frequency domain. Among existing TFA techniques such as wavelet transform, short time Fourier transform, Hilbert transform, Auto-regressive Time Frequency Analysis (ARTFA) offers good time–frequency resolution. Therefore, Chaos theory and ARTFA have been considered in this paper. First, raw ECG signal was filtered using Savitzky Golay Digital Filtering (SGDF) because it retains all important signal features after filtering. Second, a novel optimal trajectory detection step was proposed on the basis of phase space reconstruction (attractors) in chaos theory. Third, ARTFA has been used to find the spectral components of the extracted features using chaos theory. Here, ARTFA has been used for finding the autoregressive coefficients in the first step and time–frequency description in the second step. Burg method has been considered for Auto-regressive modeling to fit an ARTFA model for analyzing ECG signal by minimizing (least squares) the forward and backward prediction errors. MIT-BIH arrhythmia database has been considered for validating the present research effort. Some real time signals were also tested to explore direct usage of the proposed technique in practical applications. The obtained results show that the proposed technique enhances sensitivity, positive predictive value and accuracy thereby improving the detection of arrhythmias. Computational cost of the proposed technique is reduced to a great extent by using the chaos theory (analysis) yielding efficient detection performance. All results have been obtained in MATLAB environment R2011a. Improved values viz. 99.96% SE, 99.97% PPV, and 99.93% ACC are obtained.

On dissipative symplectic integration with applications to gradient-based optimization

Recently, continuous-time dynamical systems have proved useful in providing conceptual and quantitative insights into gradient-based optimization, widely used in modern machine learning and statistics. An important question that arises in this line of work is how to discretize the system in such a way that its stability and rates of convergence are preserved. In this paper we propose a geometric framework in which such discretizations can be realized systematically, enabling the derivation of ‘rate-matching’ algorithms without the need for a discrete convergence analysis. More specifically, we show that a generalization of symplectic integrators to non-conservative and in particular dissipative Hamiltonian systems is able to preserve rates of convergence up to a controlled error. Moreover, such methods preserve a shadow Hamiltonian despite the absence of a conservation law, extending key results of symplectic integrators to non-conservative cases. Our arguments rely on a combination of backward error analysis with fundamental results from symplectic geometry. We stress that although the original motivation for this work was the application to optimization, where dissipative systems play a natural role, they are fully general and not only provide a differential geometric framework for dissipative Hamiltonian systems but also substantially extend the theory of structure-preserving integration.

Structured backward errors in linearizations

A standard approach to calculate the roots of a univariate polynomial is to compute the eigenvalues of an associated confederate matrix instead, such as, for instance, the companion or comrade matrix. The eigenvalues of the confederate matrix can be computed by Francis's QR algorithm. Unfortunately, even though the QR algorithm is provably backward stable, mapping the errors back to the original polynomial coefficients can still lead to huge errors. However, the latter statement assumes the use of a non-structure-exploiting QR algorithm. In [J. L. Aurentz et al., Fast and backward stable computation of roots of polynomials, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 942–973] it was shown that a structure-exploiting QR algorithm for companion matrices leads to a structured backward error in the companion matrix. The proof relied on decomposing the error into two parts: a part related to the recurrence coefficients of the basis (a monomial basis in that case) and a part linked to the coefficients of the original polynomial. In this article we prove that the analysis can be extended to other classes of comrade matrices. We first provide an alternative backward stability proof in the monomial basis using structured QR algorithms; our new point of view shows more explicitly how a structured, decoupled error in the confederate matrix gets mapped to the associated polynomial coefficients. This insight reveals which properties have to be preserved by a structure-exploiting QR algorithm to end up with a backward stable algorithm. We will show that the previously formulated companion analysis fits into this framework, and we analyze in more detail Jacobi polynomials (comrade matrices) and Chebyshev polynomials (colleague matrices).

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.

Stochastic Rounding and Its Probabilistic Backward Error Analysis

Stochastic rounding rounds a real number to the next larger or smaller floating-point number with probabilities 1 minus the relative distances to those numbers. It is gaining attention in deep lear...

Structured Backward Errors for Eigenvalues of Linear Port-Hamiltonian Descriptor Systems

When computing the eigenstructure of matrix pencils associated with the passivity analysis of perturbed port-Hamiltonian descriptor system using a structured generalized eigenvalue method, one should make sure that the computed spectrum satisfies the symmetries that corresponds to this structure and the underlying physical system. We perform a backward error analysis and show that for matrix pencils associated with port-Hamiltonian descriptor systems and a given computed eigenstructure with the correct symmetry structure there always exists a nearby port-Hamiltonian descriptor system with exactly that eigenstructure. We also derive bounds for how near this system is and show that the stability radius of the system plays a role in that bound.

Sharp Error Estimates on a Stochastic Structure-Preserving Scheme in Computing Effective Diffusivity of 3D Chaotic Flows

In this paper, we study the problem of computing the effective diffusivity for particles moving in chaotic flows. Instead of solving a convection-diffusion type cell problem in the Eulerian formulation (arising from homogenization theory for parabolic equations), we compute the motion of particles in the Lagrangian formulation, which is modeled by stochastic differential equations (SDEs). A robust numerical integrator based on a splitting method was proposed to solve the SDEs and a rigorous error analysis for the numerical integrator was provided using the backward error analysis (BEA) technique [35]. However, the upper bound in the error estimate is not sharp. To improve our result, we propose a new and uniform in time error analysis for the numerical integrator that allows us to get rid of the exponential growth factor in our previous error estimate. Our new error analysis is based on a probabilistic approach, which interprets the solution process generated by our numerical integrator as a Markov process. By exploring the ergodicity of the solution process, we prove the convergence analysis of our method in computing effective diffusivity over infinite time. We present numerical results to verify the accuracy and efficiency of the proposed method in computing effective diffusivity for several chaotic flows, especially the Arnold-Beltrami-Childress (ABC) flow and Kolmogorov flow in three-dimensional space.

Acceptable Solutions and Backward Errors for Tensor Complementarity Problems

Backward error analysis reveals the numerical stability of algorithms and provides elaborate stopping criteria for iterative methods. Compared with numerical linear algebra problems, the backward error analysis for optimization problems is more rarely conducted in the literature. This paper is devoted to the backward error analysis for several generalizations of tensor complementarity problems. We first present sufficient and necessary conditions for the acceptable solutions for the extended tensor complementarity problem, the vertical tensor complementarity problem, and an extended form of tensor complementarity problem. Next, the backward errors for tensor complementarity problem are also proposed, which can be used to verify the stability of the tensor complementarity problem algorithms. Finally, some numerical examples are reported to illustrate the proposed backward errors for tensor complementarity problems.

Min-max elementwise backward error for roots of polynomials and a corresponding backward stable root finder

A new measure called min-max elementwise backward error is introduced for approximate roots of scalar polynomials p(z). Compared with the elementwise relative backward error, this new measure allows for larger relative perturbations on the coefficients of p(z) that do not participate much in the overall backward error. By how much these coefficients can be perturbed is determined via an associated max-times polynomial and its tropical roots. An algorithm is designed for computing the roots of p(z). It uses a companion linearization C(z)=A−zB of p(z) to which we added an extra zero leading coefficient, and an appropriate two-sided diagonal scaling that balances A and makes B graded in particular when there is variation in the magnitude of the coefficients of p(z). An implementation of the QZ algorithm with a strict deflation criterion for eigenvalues at infinity is then used to obtain approximations to the roots of p(z). Under the assumption that this implementation of the QZ algorithm exhibits a graded backward error when B is graded, we prove that our new algorithm is min-max elementwise backward stable. Several numerical experiments show the superior performance of the new algorithm compared with the MATLAB roots function. Extending the algorithm to polynomial eigenvalue problems leads to a new polynomial eigensolver that exhibits excellent numerical behaviour compared with other existing polynomial eigensolvers, as illustrated by many numerical tests.

New Numerical Algorithm for Deflation of Infinite and Zero Eigenvalues and Full Solution of Quadratic Eigenvalue Problems

This article presents a new method for computing all eigenvalues and eigenvectors of quadratic matrix pencil Q (λ)=λ 2 M + λ C + K . It is an upgrade of the quadeig algorithm by Hammarlinget al., which attempts to reveal and remove by deflation a certain number of zero and infinite eigenvalues before QZ iterations. Proposed modifications of the quadeig framework are designed to enhance backward stability and to make the process of deflating infinite and zero eigenvalues more numerically robust. In particular, careful preprocessing allows scaling invariant/component-wise backward error and thus a better condition number. Further, using an upper triangular version of the Kronecker canonical form enables deflating additional infinite eigenvalues, in addition to those inferred from the rank of M . Theoretical analysis and empirical evidence from thorough testing of the software implementation confirm superior numerical performances of the proposed method.

Verification of a Nondestructive Method for Assessing the Humidity of Saline Brick Walls in Historical Buildings

The paper presents the results of the verification of the neural method for assessing the humidity of saline brick walls. The method was previously developed by the authors and can be useful for the nondestructive assessment of the humidity of walls in historic buildings when destructive intervention during testing is not possible due to conservation restrictions. However, before being implemented in construction practice, this method requires validation by verification on other historic buildings, which to date has not been done. The paper presents the results of such verification, which has never been carried out before, and thus extends the scope of knowledge related to the issue. For experimental verification of the artificial neural network (ANN), the results of moisture tests of two selected historic buildings, other than those used for ANN learning and testing processes, were used. An artificial unidirectional multilayer neural network with backward error propagation and the algorithm for learning conjugate gradient (CG) was found to be useful for this purpose. The obtained satisfactory value of the linear correlation coefficient R of 0.807 and low average absolute error |Δf| of 1.16% confirms this statement. The values of average relative error |RE| of 19.02%, which were obtained in this research, were not very high for an in-situ study. Moreover, the relative error values |RE| were mostly in the range of 15% to 25%.

A two-step symmetric method for charged-particle dynamics in a normal or strong magnetic field

The study of the long time conservation for numerical methods poses interesting and challenging questions from the point of view of geometric integration. In this paper, we analyze the long time energy and magnetic moment conservations of two-step symmetric methods for charged-particle dynamics. A two-step symmetric method is proposed and its long time behaviour is shown not only in a normal magnetic field but also in a strong magnetic field. The approaches to dealing with these two cases are based on the backward error analysis and modulated Fourier expansion, respectively. It is obtained from the analysis that the method has better long time conservations than the variational method which was researched recently in the literature.

Estimating the backward error for the least-squares problem with multiple right-hand sides

Let A and B be m×n and m×d matrices, and let X˜ be an approximate solution to the problem minX⁡‖AX−B‖F. In 1996, Sun found an explicit expression for the optimal backward error—the size of the smallest perturbation to A (and possibly B) such that X˜ is an exact solution to the perturbed problem. The expression requires finding the difference of two potentially close numbers, and so its numerical evaluation can be unstable. We offer an estimate of the backward error that can be evaluated stably and when d=1 is identical to the Karlson-Waldén estimate of 1997. We prove that this estimate always approximates the optimal backward error to within a factor of 2.

Backward error measures for roots of polynomials

We analyze different measures for the backward error of a set of numerical approximations for the roots of a polynomial. We focus mainly on the element-wise mixed backward error introduced by Mastronardi and Van Dooren, and the tropical backward error introduced by Tisseur and Van Barel. We show that these measures are equivalent under suitable assumptions. We also show relations between these measures and the classical element-wise and norm-wise backward error measures.

A Sparse Spike Deconvolution Algorithm Based on a Recurrent Neural Network and the Iterative Shrinkage-Thresholding Algorithm

Conventional sparse spike deconvolution algorithms that are based on the iterative shrinkage-thresholding algorithm (ISTA) are widely used. The aim of this type of algorithm is to obtain accurate seismic wavelets. When this is not fulfilled, the processing stops being optimum. Using a recurrent neural network (RNN) as deep learning method and applying backpropagation to ISTA, we have developed an RNN-like ISTA as an alternative sparse spike deconvolution algorithm. The algorithm is tested with both synthetic and real seismic data. The algorithm first builds a training dataset from existing well-logs seismic data and then extracts wavelets from those seismic data for further processing. Based on the extracted wavelets, the new method uses ISTA to calculate the reflection coefficients. Next, inspired by the backpropagation through time (BPTT) algorithm, backward error correction is performed on the wavelets while using the errors between the calculated reflection coefficients and the reflection coefficients corresponding to the training dataset. Finally, after performing backward correction over multiple iterations, a set of acceptable seismic wavelets is obtained, which is then used to deduce the sequence of reflection coefficients of the real data. The new algorithm improves the accuracy of the deconvolution results by reducing the effect of wrong seismic wavelets that are given by conventional ISTA. In this study, we account for the mechanism and the derivation of the proposed algorithm, and verify its effectiveness through experimentation using theoretical and real data.

AIR: Iterative refinement acceleration using arbitrary dynamic precision

The increased degree of concurrent operations by lower precision arithmetic enables high performance for iterative refinement. Most of related work present statically defined mixed precision arithmetic approaches, while adapting a level of arithmetic precision dynamically in a loop with one-bit granularity can further improve the performance. This paper presents Arbitrary Dynamic Precision Iterative Refinement algorithm (AIR) that minimizes the total significand bit-width to solve iterative refinement. AIR detects the number of cancellation bits dynamically per iteration and uses the information to provide the least sufficient significand bit-width for the next iteration. We prove that AIR is a backward stable algorithm and can bring up to 2−3× speedups over a mixed precision iterative refinement depending on the characteristics of hardware. Our software demonstration shows that AIR requires only 83% of the significand bits required by mixed precision iterative refinement that solve linear systems for double precision accuracy for backward error with 32 × 32 standard normally distributed matrices.

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.