Finite Precision Arithmetic Research Articles

On the accuracy of Prony's method for recovery of exponential sums with closely spaced exponents

In this paper we establish accuracy bounds of Prony's method (PM) for recovery of sparse measures from incomplete and noisy frequency measurements, or the so-called problem of super-resolution, when the minimal separation between the points in the support of the measure may be much smaller than the Rayleigh limit. In particular, we show that PM is optimal with respect to the previously established min-max bound for the problem, in the setting when the measurement bandwidth is constant, with the minimal separation going to zero. Our main technical contribution is an accurate analysis of the inter-relations between the different errors in each step of PM, resulting in previously unnoticed cancellations. We also prove that PM is numerically stable in finite-precision arithmetic. We believe our analysis will pave the way to providing accurate analysis of known algorithms for the super-resolution problem in full generality.

Stability of the Lanczos algorithm on matrices with regular spectral distributions

We study the stability of the Lanczos algorithm run on problems whose eigenvector empirical spectral distribution is near to a reference measure with well-behaved orthogonal polynomials. We give a backwards stability result which can be upgraded to a forward stability result when the reference measure has a density supported on a single interval with square root behavior at the endpoints. Our analysis implies the Lanczos algorithm run on many large random matrix models is in fact forward stable, and hence nearly deterministic, even when computations are carried out in finite precision arithmetic. Since the Lanczos algorithm is not forward stable in general, this provides yet another example of the fact that random matrices are far from “any old matrix”, and care must be taken when using them to test numerical algorithms.

A spectrum adaptive kernel polynomial method.

The kernel polynomial method (KPM) is a powerful numerical method for approximating spectral densities. Typical implementations of the KPM require an a prior estimate for an interval containing the support of the target spectral density, and while such estimates can be obtained by classical techniques, this incurs addition computational costs. We propose a spectrum adaptive KPM based on the Lanczos algorithm without reorthogonalization, which allows the selection of KPM parameters to be deferred to after the expensive computation is finished. Theoretical results from numerical analysis are given to justify the suitability of the Lanczos algorithm for our approach, even in finite precision arithmetic. While conceptually simple, the paradigm of decoupling computation from approximation has a number of practical and pedagogical benefits, which we highlight with numerical examples.

Newton's method revisited: How accurate do we have to be?

SummaryWe analyze the convergence of quasi‐Newton methods in exact and finite precision arithmetic using three different techniques. We derive an upper bound for the stagnation level and we show that any sufficiently exact quasi‐Newton method will converge quadratically until stagnation. In the absence of sufficient accuracy, we are likely to retain rapid linear convergence. We confirm our analysis by computing square roots and solving bond constraint equations in the context of molecular dynamics. In particular, we apply both a symmetric variant and Forsgren's variant of the simplified Newton method. This work has implications for the implementation of quasi‐Newton methods regardless of the scale of the calculation or the machine.

Simple and Robust Boolean Operations for Triangulated Surfaces

Boolean operations on geometric models are important in numerical simulation and serve as essential tools in the fields of computer-aided design and computer graphics. The accuracy of these operations is heavily influenced by finite precision arithmetic, a commonly employed technique in geometric calculations, which introduces numerical approximations. To ensure robustness in Boolean operations, numerical methods relying on rational numbers or geometric predicates have been developed. These methods circumvent the accumulation of rounding errors during computation, thus preserving accuracy. Nonetheless, it is worth noting that these approaches often entail more intricate operation rules and data structures, consequently leading to longer computation times. In this paper, we present a straightforward and robust method for performing Boolean operations on both closed and open triangulated surfaces. Our approach aims to eliminate errors caused by floating-point operations by relying solely on entity indexing operations, without the need for coordinate computation. By doing so, we ensure the robustness required for Boolean operations. Our method consists of two main stages: (1) Firstly, candidate triangle intersection pairs are identified using an octree data structure, and then parallel algorithms are employed to compute the intersection lines for all pairs of triangles. (2) Secondly, closed or open intersection rings, sub-surfaces, and sub-blocks are formed, which is achieved entirely by cleaning and updating the mesh topology without geometric solid coordinate computation. Furthermore, we propose a novel method based on entity indexing to differentiate between the union, subtraction, and intersection of Boolean operation results, rather than relying on inner and outer classification. We validate the effectiveness of our method through various types of Boolean operations on triangulated surfaces.

Gerchberg-Saxton Based FIR Filter for Electronic Dispersion Compensation in IM/DD Transmission Part I: Theory and Simulation

Intensity-modulation and direct-detection (IM/DD) transmission over short-reach optical fiber links, require electronic dispersion compensation (EDC) at the transmitter and/or electronic equalization at the receiver. Recently, the iterative Gerchberg-Saxton (GS) algorithm was demonstrated for EDC in IM/DD systems, through treating the amplitude at the transmitter and the phase prior-to the direct detection receiver as a degree of freedom. In Part I of this work, three GS approaches using finite impulse response (FIR) filters for EDC in IM/DD systems are demonstrated. The first two are closely related and rely on a cascaded FIR structure, while the third offers a novel non-iterative EDC solution using a single GS optimized static FIR filter. This is achieved through decoupling pattern dependent aspects of transmission from the GS iterations by targeting a single impulse at the DD receiver. With every successive iteration an impulse response for the GS filter emerges and sets the FIR tap weights. It is also demonstrated that closed-form analytical expressions for the GS filter impulse response can be obtained through small-signal frequency-domain analysis. The FIR filter is simulated using 8-bit finite-precision arithmetic. An adaptive <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$T$</tex-math></inline-formula> -spaced post feed-forward equalizer (FFE) is utilized for mitigating residual chromatic dispersion. It is shown, that a <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$T/2$</tex-math></inline-formula> -spaced pre-EDC FIR filter with 417 taps can support 56 Gb/s non-return-to-zero (NRZ) on-off keying (OOK) transmission over 80 km of single mode fiber (SMF) with a chirp-free Mach-Zehnder modulator (MZM). Part II, presents experimental demonstration of the non-iterative GS FIR filter proposed and simulated in this article.

Franck-Condon spectra of unbound and imaginary-frequency vibrations via correlation functions: A branch-cut free, numerically stable derivation.

Molecular electronic spectra can be represented in the time domain as auto-correlation functions of the initial vibrational wavepacket. We present a derivation of the harmonic vibrational auto-correlation function that is valid for both real and imaginary harmonic frequencies. The derivation rests on Lie algebra techniques that map otherwise complicated exponential operator arithmetic to simpler matrix formulas. The expressions for the zero- and finite-temperature harmonic auto-correlation functions have been carefully structured both to be free of branch-cut discontinuities and to remain numerically stable with finite-precision arithmetic. Simple extensions correct the harmonic Franck-Condon approximation for the lowest-order anharmonic and Herzberg-Teller effects. Quantitative simulations are shown for several examples, including the electronic absorption spectra of F2, HOCl, CH2NH, and NO2.

On the rational approximation of Markov functions, with applications to the computation of Markov functions of Toeplitz matrices

We investigate the problem of approximating the matrix function f(A) by r(A), with f a Markov function, r a rational interpolant of f, and A a symmetric Toeplitz matrix. In a first step, we obtain a new upper bound for the relative interpolation error 1 − r/f on the spectral interval of A. By minimizing this upper bound over all interpolation points, we obtain a new, simple, and sharp a priori bound for the relative interpolation error. We then consider three different approaches of representing and computing the rational interpolant r. Theoretical and numerical evidence is given that any of these methods for a scalar argument allows to achieve high precision, even in the presence of finite precision arithmetic. We finally investigate the problem of efficiently evaluating r(A), where it turns out that the relative error for a matrix argument is only small if we use a partial fraction decomposition for r following Antoulas and Mayo. An important role is played by a new stopping criterion which ensures to automatically find the degree of r leading to a small error, even in presence of finite precision arithmetic.

Explicit Solutions of the Singular Yang–Baxter-like Matrix Equation and Their Numerical Computation

We derive several explicit formulae for finding infinitely many solutions of the equation \(AXA=XAX\), when A is singular. We start by splitting the equation into a couple of linear matrix equations and then show how the projectors commuting with A can be used to get families containing an infinite number of solutions. Some techniques for determining those projectors are proposed, which use, in particular, the properties of the Drazin inverse, spectral projectors, the matrix sign function, and eigenvalues. We also investigate in detail how well-known similarity transformations like Jordan and Schur decompositions can be used to obtain new representations of the solutions. The computation of solutions by the suggested methods using finite precision arithmetic is also a concern. Difficulties arising in their implementation are identified and ideas to overcome them are discussed. Numerical experiments shed some light on the methods that may be promising for solving numerically the matrix equation.

LMBOPT: a limited memory method for bound-constrained optimization

Recently, Neumaier and Azmi gave a comprehensive convergence theory for a generic algorithm for bound constrained optimization problems with a continuously differentiable objective function. The algorithm combines an active set strategy with a gradient-free line search CLS along a piecewise linear search path defined by directions chosen to reduce zigzagging. This paper describes LMBOPT, an efficient implementation of this scheme. It employs new limited memory techniques for computing the search directions, improves CLS by adding various safeguards relevant when finite precision arithmetic is used, and adds many practical enhancements in other details. The paper compares LMBOPT and several other solvers on the unconstrained and bound constrained problems from the CUTEst collection and makes recommendations on which solver to use and when. Depending on the problem class, the problem dimension, and the precise goal, the best solvers are LMBOPT, ASACG, and LMBFG-EIG-MS.

Sequential method for fast neural population activity reconstruction in the cortex from incomplete noisy measurements

During recent years there has been a growing interest in stochastic dynamic neural fields employed for modeling and predictions in biomedical and technical systems. In this paper, given some incomplete noisy data available from sensors, we propose and explore a state estimation method for fast restorations of membrane potential in the cortex based on such measurements and the Amari equation used for simulations of neural population activity in a stochastic setting. Our novel technique relies upon a Galerkin-type spectral approximation utilized within the conventional state-space approach. Translating a stochastic system into its state-space form creates a straightforward and fruitful way to the data-driven parameter estimation, filtering, prediction and smoothing. The present study is particularly focused on establishing a nonlinear stochastic Galerkin-spectral-approximation-induced system of large size, which is further estimated by the traditional extended Kalman filter (EKF). The efficiency of calculations is the main purpose of our research. That is why the fast filtering solution devised is based on processing the incoming data incrementally, that is, by processing measurements one at a time, rather than handling them as a unified high-dimensional vector. Such sequential filters suit well for dealing with large data sets as well as with real-time on-line computations. Also, their derivation and substantiation is of great interest in the context of neural network training because of large stochastic systems arisen there. In comparison to the batch filtering, our novel algorithm reduces the computational cost of membrane potential reconstructions in terms of the amount of grid nodes N accepted in the underlying spacial discretization, significantly. Apart from its computation efficiency, this sequential method is more robust to round-off errors committed within a computer-based finite precision arithmetics than the classical EKF because of the (N × N)-matrix inversion elimination from such membrane potential calculations. The superior performance of our technique is examined and confirmed in comparison to the batch one on two known scenarios in the dynamic neural field modeling.

Nearly reducible finite Markov chains: Theory and algorithms.

Finite Markov chains, memoryless random walks on complex networks, appear commonly as models for stochastic dynamics in condensed matter physics, biophysics, ecology, epidemiology, economics, and elsewhere. Here, we review exact numerical methods for the analysis of arbitrary discrete- and continuous-time Markovian networks. We focus on numerically stable methods that are required to treat nearly reducible Markov chains, which exhibit a separation of characteristic timescales and are therefore ill-conditioned. In this metastable regime, dense linear algebra methods are afflicted by propagation of error in the finite precision arithmetic, and the kinetic Monte Carlo algorithm to simulate paths is unfeasibly inefficient. Furthermore, iterative eigendecomposition methods fail to converge without the use of nontrivial and system-specific preconditioning techniques. An alternative approach is provided by state reduction procedures, which do not require additional a priori knowledge of the Markov chain. Macroscopic dynamical quantities, such as moments of the first passage time distribution for a transition to an absorbing state, and microscopic properties, such as the stationary, committor, and visitation probabilities for nodes, can be computed robustly using state reduction algorithms. The related kinetic path sampling algorithm allows for efficient sampling of trajectories on a nearly reducible Markov chain. Thus, all of the information required to determine the kinetically relevant transition mechanisms, and to identify the states that have a dominant effect on the global dynamics, can be computed reliably even for computationally challenging models. Rare events are a ubiquitous feature of realistic dynamical systems, and so the methods described herein are valuable in many practical applications.

Data-dependent orthogonal polynomials on generalized circles: A unified approach applied to [formula omitted]-domain identification

The performance of algorithms in system identification and control, depends on their implementation in finite-precision arithmetic. The aim of this paper is to develop a unified approach for numerically reliable system identification that combines the numerical advantages of data-dependent orthogonal polynomials and the discrete-time δ-domain parametrization. In this paper, earlier results for discrete orthogonal polynomials on the real-line and the unit-circle are generalized to obtain an approach for the construction of orthogonal polynomials on generalized circles in the complex plane. This enables the formulation of a unified framework for the numerically reliable identification of systems expressed in the δ-domain, as well as in the traditional Laplace and Z-domains. An example is presented which shows the significant numerical advantages of the δ-domain approach for the identification of fast-sampled systems.

Computing the partition function of the Sherrington–Kirkpatrick model is hard on average

We establish the average-case hardness of the algorithmic problem of exact computation of the partition function associated with the Sherrington–Kirkpatrick model of spin glasses with Gaussian couplings and random external field. In particular, we establish that unless P=#P, there does not exist a polynomial-time algorithm to exactly compute the partition function on average. This is done by showing that if there exists a polynomial time algorithm, which exactly computes the partition function for inverse polynomial fraction (1/nO(1)) of all inputs, then there is a polynomial time algorithm, which exactly computes the partition function for all inputs, with high probability, yielding P=#P. The computational model that we adopt is finite-precision arithmetic, where the algorithmic inputs are truncated first to a certain level N of digital precision. The ingredients of our proof include the random and downward self-reducibility of the partition function with random external field; an argument of Cai et al. (In STACS 99 (Trier) (1999) 90–99 Springer) for establishing the average-case hardness of computing the permanent of a matrix; a list-decoding algorithm of Sudan (In 37th Annual Symposium on Foundations of Computer Science (Burlington, VT, 1996) (1996) 164–172 IEEE Comput. Soc. Press), for reconstructing polynomials intersecting a given list of numbers at sufficiently many points; and near-uniformity of the log-normal distribution, modulo a large prime p. To the best of our knowledge, our result is the first one establishing a provable hardness of a model arising in the field of spin glasses. Furthermore, we extend our result to the same problem under a different real-valued computational model, for example, using a Blum–Shub–Smale machine (In [Proceedings 1988] 29th Annual Symposium on Foundations of Computer Science (1988) 387–397 IEEE) operating over real-valued inputs. We establish that, if there exists a polynomial time algorithm which exactly computes the partition function for 34+1nO(1) fraction of all inputs, then there exists a polynomial time algorithm, which exactly computes the partition function for all inputs, with high probability, yielding P=#P. Our proof uses the random self-reducibility of the partition function, together with a control over the total variation distance for log-normal random variables in presence of a convex perturbation, and the Berlekamp–Welch algorithm.

Exact linesearch limited-memory quasi-Newton methods for minimizing a quadratic function

The main focus in this paper is exact linesearch methods for minimizing a quadratic function whose Hessian is positive definite. We give a class of limited-memory quasi-Newton Hessian approximations which generate search directions parallel to those of the BFGS method, or equivalently, to those of the method of preconditioned conjugate gradients. In the setting of reduced Hessians, the class provides a dynamical framework for the construction of limited-memory quasi-Newton methods. These methods attain finite termination on quadratic optimization problems in exact arithmetic. We show performance of the methods within this framework in finite precision arithmetic by numerical simulations on sequences of related systems of linear equations, which originate from the CUTEst test collection. In addition, we give a compact representation of the Hessian approximations in the full Broyden class for the general unconstrained optimization problem. This representation consists of explicit matrices and gradients only as vector components.

Secure multiparty computations in floating-point arithmetic

Abstract Secure multiparty computations enable the distribution of so-called shares of sensitive data to multiple parties such that the multiple parties can effectively process the data while being unable to glean much information about the data (at least not without collusion among all parties to put back together all the shares). Thus, the parties may conspire to send all their processed results to a trusted third party (perhaps the data providers) at the conclusion of the computations, with only the trusted third party being able to view the final results. Secure multiparty computations for privacy-preserving machine-learning turn out to be possible using solely standard floating-point arithmetic, at least with a carefully controlled leakage of information less than the loss of accuracy due to roundoff, all backed by rigorous mathematical proofs of worst-case bounds on information loss and numerical stability in finite-precision arithmetic. Numerical examples illustrate the high performance attained on commodity off-the-shelf hardware for generalized linear models, including ordinary linear least-squares regression, binary and multinomial logistic regression, probit regression and Poisson regression.

On the Convergence Rate of Variants of the Conjugate Gradient Algorithm in Finite Precision Arithmetic

On the Convergence Rate of Variants of the Conjugate Gradient Algorithm in Finite Precision Arithmetic

On choices of formulations of computing the generalized singular value decomposition of a large matrix pair

For the computation of the generalized singular value decomposition (GSVD) of a large matrix pair $(A,B)$ of full column rank, the GSVD is commonly formulated as two mathematically equivalent generalized eigenvalue problems, so that a generalized eigensolver can be applied to one of them and the desired GSVD components are then recovered from the computed generalized eigenpairs. Our concern in this paper is, in finite precision arithmetic, which generalized eigenvalue formulation is numerically preferable to compute the desired GSVD components more accurately. We make a detailed perturbation analysis on the two formulations and show how to make a suitable choice between them. Numerical experiments illustrate the results obtained.

An Error-Free Transformation for Matrix Multiplication with Reproducible Algorithms and Divide and Conquer Methods

This paper discusses accurate numerical algorithms for matrix multiplication. Matrix multiplication is a basic and important problem in numerical linear algebra. Numerical computations using floating-point arithmetic can be quickly performed on existing computers. However, the accumulation of rounding errors due to finite precision arithmetic is a critical problem. An error-free transformation for matrix multiplication is reviewed in this paper. Such a transformation is extremely useful for developing accurate numerical algorithms for matrix multiplication. One advantage of the transformation is that it exploits Basic Linear Algebra Subprograms (BLAS). We provide a rounding error analysis of reproducible algorithms for matrix multiplication based on the error-free transformation. In addition, we propose an error-free transformation for matrix multiplication that can be utilized with the divide and conquer methods.

The J-Orthogonal Square-Root MATLAB-Based Continuous-Discrete Unscented Kalman Filtering Method

The paper suggests a general solution to the square-rooting problem existed for the Unscented Kalman Filter (UKF) since its appearance in the late 1990s. As properly noted in engineering literature, the previously suggested Cholesky-based UKF implementations are, in fact, the ‘pseudo’ square-root versions. Their key feature is the utilization of one-rank Cholesky update procedure required at each filtering step because of the possibly negative sigma points’ weights. In a finite precision arithmetic, the resulting downdated matrix might be not a positive definite matrix. This yields a failure of the UKF estimator in practice. We resolve this problem by suggesting a novel square-root approach based on the J-orthogonal matrix utilization for updating the required Cholesky factors. Additionally, we explain how the MATLAB language with built-in numerical integration schemes developed for solving ordinary differential equations can be easily and effectively used for accurate calculations when implementing the continuous-discrete UKF time update stage.