work that allows the probabilistic error analysis of algorithms under SR. In this context, we decompose the error into a martingale plus a drift. We show that the drift is zero for algorithms with multi-linear errors, while the probabilistic analysis of the martingale term leads to probabilistic error bounds in O( √ n). We show that the drift is negligible at the first order compared to the martingale term for the variance computation, and we prove probabilistic error bounds in O( √ n).

I express my sincere gratitude to Devan Sohier, Pablo de Oliveira Castro, and Eric Petit for entrusting me with the opportunity to pursue my PhD under their guidance at the Laboratoire d'informatique Parallélisme Réseaux Algorithmes Distribués (LI-PaRAD). Working with them for three years was an invaluable experience. I greatly appreciated the inspiring discussions with them, which helped me understand things better and greatly impacted my academic journey. I learned a lot from their research experience and was honored to contribute to an exciting research project. I am grateful for their guidance, which was vital to my academic and personal development. Furthermore, I am infinitely grateful for the patience you have afforded me during the more challenging moments, particularly health problems. Thank you for that.

I would like to express my gratitude to Stef Graillat and Paul Zimmermann for agreeing to be my thesis referees. I am deeply appreciative of their time, expertise, and dedication, which played a key role in ensuring the academic rigor and excellence of my research. I am thankful for their comments and questions that have enhanced the content. I also extend warm thanks to Sylvie Boldo, Alexandre d'Aspremont, and Ilse C.F Ipsen for agreeing to be part of my jury. Special thanks to Ilse C.F Ipsen, who has provided me with constructive and valuable feedback, which has helped me a lot to improve the manuscript.

I extend my sincere gratitude to the members of the Interflop project for welcoming me into their team. It is truly an honor to participate in this innovative and dynamic project. Special thanks to David Defour and Bruno Lathuiliere for their warm welcome and support. I express my sincere gratitude to all members of LI-PaRAD who have played a significant role in fostering a pleasant work atmosphere. I am particularly grateful to Thomas Dufaud for his invaluable support over three years.

My daily life at the LI-Parad would be impossible without Max Hoffer and Aurelien Delval, with whom I shared my office. I thank them wholeheartedly for their help and support. Special thanks to Mohamed El Ibrahimi and Yassine Ennabo, my apartment mates for more than a year; their support has been a source of motivation throughout this period. I would also like to express my gratitude to Yassin Bouhafs, Anas Bouali, Omar Belkhouad, Mohammed Mansouri Mjahed, and Salah-Eddine EL Bouzidi for their encouragement despite our physical distance. Now, let us address the most sensitive part. My deepest and most heartfelt thoughts go to my parents. Their unwavering support and unconditional love have been the cornerstone of my journey. I am grateful to them for their sacrifices, encouragement, and constant kindness. Special thanks to my brothers and my betrothed Soukaina Doumagui, for their unconditional support. What I have achieved would not have been possible without them. Thank you from the bottom of my heart for being the source of my strength and inspiration.

-Introduction

. General Context

The problem of the representation of numbers has ancient roots. Around 2000 BC, the Babylonians used a number system based on the number 60 [START_REF] Neugebauer | The Exact Sciences in Antiquity[END_REF].

That system allowed the Babylonians to perform arithmetic operations relatively efficiently. Archimedes is attributed to inventing a system for naming extremely large numbers in his famous treatise Arenarius (The Sand Reckoner, c. 287-212 BC) [START_REF] Hirshfeld | Eureka man: The life and legacy of Archimedes[END_REF]. While it does not directly correspond to our modern concept of exponents, it prefigures a form of exponential representation. Over time, various civilizations introduced distinct numerical notations, including the "Hindu-Arabic" numerals developed first in India within the Hindu culture, followed by its adoption and refinement in the Arabic numeral system.

Not long after this, Al-Khwarizmi wrote his book "Al-Kitab al-Mukhtasar fi Hisab al-Jabr wal-Muqabala," which translates to "The Compendious Book on Calculation by Completion and Balancing." This work laid the foundation for the field of algebra and introduced the "algorithm" notion. It also played a crucial role in transferring the Greek and Hindu-Arabic numeral system to Europe at the beginning of the 12th century, which was the groundwork for modern mathematics and science development.

Over time, these developments set the stage for the scientific revolution in the 16th and 17th centuries by luminaries like Johannes Kepler and Isaac Newton. In particular, Newton significantly contributed to the development of calculus, a fundamental stepping stone for solving mathematical physics problems such as ordinary differential equations (ODEs). By introducing the concept of derivatives and integrals, Newton provided a powerful mathematical framework for understanding change and accumulation. Interestingly, it took about 2000 years of migration of astronomical knowledge from Mesopotamia via Greeks, Hindus, and Arabs to arrive at a truly numerical system. For more detail, we refer to [START_REF] Neugebauer | The Exact Sciences in Antiquity[END_REF][START_REF] Donald | The art of computer programming[END_REF].

In the 20th century, scientific development underwent a qualitative leap thanks to Hilbert's questions. In the domain of scientific computing, where there was a growing need to automate solutions for scientific problems, researchers confronted the challenge of developing a computer to expedite calculations. Turing made important contributions to this development, marking a significant milestone in the history of science and technology. Turing addressed Hilbert's decidability question using the abstract device [START_REF] Turing | On computable numbers, with an application to the entscheidungsproblem[END_REF], now known as a Turing machine. Furthermore, he participated in constructing the first programmable general-purpose electronic digital computer, ENIAC, during World War II in 1945, which was a groundbreaking achievement. ENIAC was used to calculate artillery firing tables and perform other complex mathematical operations.

ENIAC was the first computer to incorporate hardware support for floatingpoint arithmetic, a method of representing and performing arithmetic operations on real numbers in a computer. Floating-point arithmetic is commonly used in various scientific and engineering applications, including simulations and scientific computing. It provides a compromise between precision and range by allowing the representation of various types and sizes of numbers while maintaining reasonable accuracy. This flexibility makes it essential for solving complex scientific problems in fields such as physics, chemistry, biology, and engineering. However, floating-point arithmetic is not exact and can introduce rounding errors.

Rounding error is a phenomenon that can occur during the representation of numbers or elementary operations. It arises because not all real numbers can be accurately represented in floating-point format due to the limited number of bits available. As a result, these numbers must be rounded to representable floating-point numbers, leading to errors.

The limited precision may cause various types of numerical errors, such as catastrophic cancellation, which occurs when subtracting two nearly equal numbers (i.e., two numbers whose difference has a smaller exponent than either number). Absorption occurs when adding a small value and a much larger value (i.e., the exponent of the larger value is significantly greater than that of the smaller value). Stagnation, an extreme case of absorption, occurs when the result of a floating-point calculation is repeatedly rounded to the same value, and then the information in the updates is lost.

Rounding errors can propagate through repeated operations using inaccurate numbers. Therefore, the way we perform calculations can affect how errors add up. This becomes especially important when we have long sequences of operations, as these errors can accumulate and make our final outcome less accurate [60, p. 8]. This is relevant in various domains, such as scientific and numerical computing, where precision and accuracy are crucial for obtaining reliable results.

In the history of scientific computing, stochastic rounding [START_REF] Parker | Monte Carlo Arithmetic: exploiting randomness in floatingpoint arithmetic[END_REF], a computing paradigm developed as a model for performing computations using precise or imprecise data, has attracted a lot of attention. It can be used to estimate empirically the numerical error of computer programs. Since the 1940s, John von Neumann and Goldstine [START_REF] Neumann | Numerical inverting of matrices of high order[END_REF] first proposed modeling rounding errors as random errors to obtain probabilistic error bounds. They suggested a probabilistic analysis of rounding errors in elementary operations such as addition or division. Although they were unaware of the nature of these errors, they treated them as random variables with a known maximum size, leading to probabilistic investigations.

After World War II, there has been a rapid increase in computing resources and simulation complexity. During the 1950s, 1960s, and 1970s, the field of computing experienced significant growth and the development of various computer systems. However, before the 1980s, the representation of floatingpoint numbers was a complex and challenging issue for programmers working with different computers and systems, as described in William Kahan's paper [START_REF] Kahan | Why do we need a floating-point arithmetic standard[END_REF]. There was a need for a unique format to represent these numbers.

Therefore, the publication of the IEEE-754 standard [START_REF] Institute | group of the Microprocessor Standards Subcommittee and A[END_REF] in 1985 established a unified approach for floating-point arithmetic, making computations more reliable and consistent across different computers. A significant revision was published in 2008 [START_REF] Zuras | Ieee standard for floating-point arithmetic[END_REF], and a minor revision was released in 2019 [START_REF]IEEE standard for floating-point arithmetic, IEEE Std[END_REF].

The IEEE-754 norm [START_REF]IEEE standard for floating-point arithmetic, IEEE Std[END_REF] defines five rounding modes for floating-point arithmetic: round to nearest ties to even (the default rounding mode and RN in the following), round to nearest ties away, round to zero, round to +∞, and round to -∞. These modes are characterized by their deterministic nature, in which the rounded value of a number x is determined by an exact value (depending on x), and an arbitrary sequence of rounded elementary operations will always produce the same result.

Stochastic rounding [START_REF] Croci | Stochastic rounding: implementation, error analysis and applications[END_REF] can also be used as a replacement for the default deterministic rounding mode in numerical computations. SR is a probabilistic rounding mode: an inexact computation is rounded to the next smaller or larger floating-point number with probability depending on the distances to those numbers. To the best of our knowledge, Forsythe first presented the use of SR to reduce the accumulation of rounding errors. In 1949, at the fifty-second meeting of the American Mathematical Society, he proposed [START_REF] Forsythe | Round-off errors in numerical integration on automatic machinery. preliminary report[END_REF], a rounding mode called "random round-off " in the context of solving simple ordinary differential equations. This rounding mode involves rounding up or down with a certain probability.

. Problematic

The quality of a numerical computation is usually measured by its accuracy, and in general, finding good error bounds for numerical algorithms is difficult. Probabilistic error bounds can be an alternative proposition to measure the numerical error of algorithms. The probabilistic analysis takes into account the distribution of inputs and analyzes algorithm behavior on average. It provides error bounds valid from certain probability. Note that considering rounding errors as random variables to calculate probabilistic error bounds is not new. It dates back to the work of von Neumann and Goldstine [START_REF] Neumann | Numerical inverting of matrices of high order[END_REF], Henrici [START_REF] Henrici | Discrete variable methods in ordinary differential equations[END_REF][START_REF]Elements of numerical analysis[END_REF][START_REF]Test of probabilistic models for the propagation of roundoff errors[END_REF], among others.

In the worst case, the error bound of a computation with n operations is proportional to nu, where u is the maximal error of each operation. This sit-uation is attainable; for instance, consider the summation of n real numbers, each affected by an error of the same sign with a maximal error. However, Wilkinson [71, sec 1.33] had the intuition that the roundoff error accumulated in n operations is typically proportional to √ nu rather than nu. Von Neumann and Goldstine [START_REF] Neumann | Numerical inverting of matrices of high order[END_REF] observed that assuming rounding errors as independent random variables uniformly distributed in [-u; u], the dispersion of the sum is given by

√ nu/ √ 3,
where n is the number of random variables. Furthermore, the Central Limit Theorem (CLT) [START_REF] Dacunha-Castelle | Probability and Statistics[END_REF] aligns with this intuition in the context of a sum of independent random variables. This creates a strong temptation to use probabilistic analysis instead of deterministic analysis. Several numerical analyses [START_REF] Croci | Stochastic rounding: implementation, error analysis and applications[END_REF] have demonstrated higher accuracy when using SR as a replacement for RN. However, the mathematical studies in this area remain limited and do not align with the observed numerical advantages of SR. In this dissertation, our primary focus is on conducting a theoretical analysis of stochastic rounding as a rounding mode. Specifically, we address the following questions:

• What precise meaning can we give to Wilkinson's intuition in SR, and to what extent does it hold? For some algorithms, under the independence assumption, the rule of thumb that one can replace a nu error bound with √ nu has been proven using probabilistic bounds. However, using SR and without additional assumptions, demonstrating the validity of this property can be challenging.

• What can we infer about the variance of an algorithm under SR? The variance analysis of a SR computation has yet to attract any attention in the literature despite allowing the use of several probabilistic properties.

• How can we enhance existing probabilistic error bounds, and what is the behavior of SR in low-precision? Current probabilistic bounds are based on the Azuma-Hoeffding inequality, but other concentration inequalities can be employed to ensure tight error bounds.

• What impact does SR have on complex algorithms? For algorithms with multi-linear error, SR is unbiased and provides tight probabilistic error bounds. However, it is essential to investigate non-linear algorithms under SR to confirm the extension of SR benefits to these problems.

. Motivation

Let us illustrate how SR works through a simple example: let z = 1.2•10 0 + 4.8 • 10 -1 = 1.2 + 0.48 = 1.68. With RN, keeping two significant digits ẑ = 1.7. While for SR, ẑ = 1.7 with probability 0.8 and ẑ = 1.6 with probability 0.2. The expected result is the exact value: 0.8 × 1.7 + 1.6 × 0.2 = 1.68. SR results are concentrated around the exact result because its randomness breaks the bias of RN [START_REF] Parker | Monte Carlo Arithmetic: exploiting randomness in floatingpoint arithmetic[END_REF], which in this case always rounds upwards.

Stochastic rounding has two main applications [START_REF] Croci | Stochastic rounding: implementation, error analysis and applications[END_REF]. First, it can be used to estimate empirically the numerical error of computer programs; SR introduces a random noise in each floating-point operation, and then a statistical analysis of the set of sampled outputs can be applied to estimate the effect of rounding errors. To make this simulation available, various tools such as verificarlo [START_REF] Denis | Verificarlo: Checking floating point accuracy through Monte Carlo arithmetic[END_REF], Verrou [START_REF] Févotte | Verrou: a cestac evaluation without recompilation[END_REF], and Cadna [START_REF] Jézéquel | Cadna: a library for estimating roundoff error propagation[END_REF] have been developed. Second, stochastic rounding can also be used as a replacement for the default deterministic rounding mode in numerical computations. It has been demonstrated that in multiple domains such as neural networks, ODEs, and PDEs [START_REF] Croci | Stochastic rounding: implementation, error analysis and applications[END_REF], and in low-precision, SR provides better results compared to RN.

Studying algorithms under SR in low-precision, especially bfloat-16, is becoming increasingly attractive due to its higher speed and lower energy consumption. In this regard, Artificial Intelligence (AI) has motivated research on stochastic rounding due to its higher accuracy in various applications such as deep learning and optimization. In AI, particularly deep learning, neural network training involves the accumulation of gradients during backpropagation. The deterministic rounding modes can introduce biases into these accumulated gradients. At the same time, Gupta et al. [START_REF] Gupta | Deep learning with limited numerical precision[END_REF] have shown that stochastic rounding can be used as an alternative rounding mode, leveraging its ability to reduce gradient biases in gradient updates.

Connolly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] have shown that SR successfully prevents the phenomenon of stagnation that takes place in various applications such as neural networks, ODEs, and PDEs. In particular, Gupta et al. show in [START_REF] Gupta | Deep learning with limited numerical precision[END_REF] that deep neural networks are prone to stagnation during the training phase. For PDEs, solved via Runge-Kutta finite difference methods in low precision, SR avoids stagnation in the computations of the heat equation solution as proved in [START_REF] Croci | Effects of round-to-nearest and stochastic rounding in the numerical solution of the heat equation in low precision[END_REF].

Despite its potential advantages, hardware units implementing SR have yet to be widely available in most computer systems. However, this rounding mode has been successfully integrated into diverse specialized processors such as Graphcore IPUs [START_REF] Limited | c ai-float™-mixed precision arithmetic for AI: A hardware perspective[END_REF], which supports SR for 32 bits floating point, bi-nary32, and 16 bits floating point, binary16, or Intel neuromorphic chip Loihi [START_REF] Davies | Loihi: A neuromorphic manycore processor with on-chip learning[END_REF] to improve the accuracy of biological neuron and synapse models. Also, AMD [START_REF] Loh | Stochastic rounding logic[END_REF], NVIDIA [START_REF] Alben | Stochastic rounding of numerical values[END_REF], IBM [START_REF] Bradbury | Reproducible stochastic rounding for out of order processors[END_REF][START_REF]Stochastic rounding floating-point multiply instruction using entropy from a register[END_REF], and other computing companies [START_REF] Henry | Processor with memory array operable as either cache memory or neural network unit memory[END_REF][START_REF] Kanter | Apparatus and methods for hardware-efficient unbiased rounding[END_REF][START_REF] Lifsches | In-memory stochastic rounder[END_REF] own several related patents. These developments support the idea of hardware implementations using SR becoming more available in the future.

. Theoretical Framework of Stochastic Rounding

Numerical accuracy is usually measured through the error between the actual computation and a reference value, such as the exact mathematical solution or a measure obtained by experimentation (often computed in high precision compared to the precision used). In the literature and to investigate the SR effect on algorithms, several probabilistic tools have been used, such as the independence of random variables, the Central Limit Theorem [START_REF] Dacunha-Castelle | Probability and Statistics[END_REF], and concentration inequalities [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF].

Relying on the independence assumption of random errors, Higham and Mary [START_REF] Higham | A new approach to probabilistic rounding error analysis[END_REF] have shown that in various linear algebra computations, such as the inner product, a probabilistic bound of the error proportional to n ln (n)u can be achieved rather than the deterministic bound in O(nu). Their approach uses the Azuma-Hoeffding inequality for independent random variables, an inequality that provides a bound on how a sum of n independent random variables deviates from its expected value [START_REF] Boucheron | Concentration Inequalities: A Nonasymptotic Theory of Independence[END_REF]. Moreover, by employing the unbiased nature of SR [START_REF] Parker | Monte Carlo Arithmetic: exploiting randomness in floatingpoint arithmetic[END_REF], the expected value coincides with the exact value, which allows obtaining a probabilistic bound on the absolute error.

Although the independence assumption is not always true, errors can accumulate and propagate throughout the computation process, leading to correlated errors. Factors such as iterative algorithms and numerical approximations contribute to the dependence on random errors. This result was the first that shows a probabilistic error bound valid for any n, unlike results obtained by applying the central limit theorem, which applies only as n → +∞.

In collaboration with Connolly [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF], they have demonstrated the same bounds by proving the following result: since the mean of each rounding error is zero regardless of the previous computation and assuming a perfect random generator, the random errors under SR-nearness satisfy the mean independence, a weaker property than the independence of random errors. Consequently, using this tool, they constructed a martingale [START_REF] Dacunha-Castelle | Probability and Statistics[END_REF]. This development is particularly significant as it enables the retrieval of various probabilistic properties, including concentration inequalities and CLT.

A martingale [START_REF] Dacunha-Castelle | Probability and Statistics[END_REF] is a sequence of random variables where the expected value of the next variable, given all the previous ones, is equal to the current value. In the context of the probabilistic error analysis with SR, various methods exist to construct the martingale, and the technique used impacts the quality of the final result. For instance, in the case of the inner product, Ipsen, and Zhou [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] form a martingale in a different way than Connelly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF]. Their method demonstrates that the probabilistic bound of the forward error is proportional to √ nu rather than nu when nu ≪ 1.

The martingale central limit theorem also implies that under certain conditions, the error converges in distribution to a normal distribution that is characterized by its mean and variance [START_REF] Dacunha-Castelle | Probability and Statistics[END_REF]. This behavior is often observed in practice. In this case, the number of significant digits can be estimated by -log( σ |µ| ) where σ is the standard deviation (the square root of the variance) and µ is the expected value [START_REF] Parker | Monte Carlo Arithmetic: exploiting randomness in floatingpoint arithmetic[END_REF]. Sohier et al. [START_REF] Sohier | Confidence intervals for stochastic arithmetic[END_REF] have shown that this propo-sition is valid but under normality assumption with a probability 0.68.

. Contributions

The recent theoretical developments presented in last section have shown that when using SR, the probabilistic error bound of algorithms such as the inner product is proportional to

√

nu instead of nu. The overall goal of this thesis was to pursue this intuition and generalize this benefit to other complex algorithms.

The first research question was a study of the bias and a comparison of two stochastic rounding modes (SR-nearness and SR-up-or-down) and RN-nearest32 on rectangular integration, which is at the basis of Euler's Method for ODE. In Chapter 3, which is an extension of the paper "The Positive Effects of Stochastic Rounding in Numerical Algorithms" In 29th IEEE Symposium on Computer Arithmetic ARITH 2022, E-M. El Arar, D. Sohier, P. de Oliveira Castro, and E. Petit (see [START_REF] Arar | The positive effects of stochastic rounding in numerical algorithms[END_REF]), we demonstrate through two examples (the constant function and the cosine function) that bias can result in a loss of accuracy.

For the constant function, an exact expression and an estimation of the bias are given for SR-up-or-down. We show how the accumulation of errors with both SR-up-or-down and RN-nearest32 can significantly impact the accuracy of computations, even on simple algorithms, and that SR-nearness can remain unbiased and provide the full expected precision on them.

We assume an error-free cosine function and only focus on the errors accumulated through elementary operations. In addition to the previous results for the constant function, we also present an expression for the method error of the cosine function. We conclude with numerical experiments that validate the theoretical research discussed in this chapter.

Since SR-nearness is unbiased and satisfies the mean independent property, we thus decided to focus on this rounding mode and explore more findings and insights that can contribute to enhancing SR as a rounding mode. We analyzed several algorithms and demonstrated probabilistic bounds on the error in O( √ nu) instead of O(nu) for the deterministic bounds.

In chapter 4, which is an extension of the paper "Stochastic Rounding Variance and Probabilistic Bounds: A New Approach" SIAM Journal on Scientific Computing 2022, E-M. El Arar, D. Sohier, P. de Oliveira Castro, and E. Petit (see [START_REF]Stochastic rounding variance and probabilistic bounds: A new approach[END_REF]), we apply the approach based on Azuma-Hoeffding inequality [53, p. 303] and martingales [53, p. 295] (AH method in the following) to the sequential summation (Sub-section 4.1.1) in which we explain how this method can be used to establish probabilistic error bounds in O( √ nu) for algorithms with exact inputs. We also extend this method to derive a new probabilistic bound on the forward error of the Horner algorithm in O( √ nu) (Sub-section 4.1.2).

The analysis of this algorithm differs from that of the summation or the inner product due to the presence of an error affecting one of the multiplication operands. This extends the approach to a whole new class of algorithms, and is a step towards its application to all numerical schemes.

Hallman and Ipsen [START_REF] Hallman | Precision-aware deterministic and probabilistic error bounds for floating point summation[END_REF] have studied pairwise summation in the context of SR, showing that the forward error for a sum of n values has a probabilistic bound in O( log(n)u) instead O(log(n)u) for RN. In Sub-section 4.1.3, we propose a more straightforward technique to show the martingale presence of this algorithm, which improves Hallman and Ipsen pairwise summation error bound [START_REF] Hallman | Precision-aware deterministic and probabilistic error bounds for floating point summation[END_REF].

To our knowledge, the variance analysis of a SR computation has not attracted any attention in the literature. In Section 4.2, we have also introduced a novel approach referred to as the BC method in the following. It exclusively relies on error variance information and exhibits enhanced accuracy for a large problem size n. This new method uses the Bienaymé-Chebyshev inequality [8, p. 19] to establish a probabilistic error bound. Based on the mean independence property, we have presented Lemma 11, a general framework applicable to a wide class of algorithms that allows to compute deterministic bounds of the variance.

We have illustrated the applicability of this method to the previous algorithms studied with the AH method. We demonstrated that the use of Bienaymé-Chebyshev inequality combined with the previous variance bound leads to probabilistic bound also in O( √ nu). Note that both methods yield to obtain probabilistic bounds depending on three parameters: the precision u, the problem size n, and the probability λ that a SR-nearness computation has an error greater than the bound.

In Section 4.3, we analyze these probabilistic bounds, and we show that the one obtained by the BC method is tighter in many cases. We also demonstrate superior accuracy with BC method in low-precision formats, which are becoming critical in high-performance computing to reduce computation and storage costs. In conclusion, in this chapter, we have demonstrated that any multi-linear transformation of errors under SR-nearness forms a martingale. We have shown that using the AH method ensures a tight probabilistic bound in probability, while the BC method guarantees a tight probabilistic bound for a large problem size n.

In previous theoretical investigations concerning SR error bounds, as discussed in Chapter 4, the focus has been on algorithms where the resulting error is a multi-linear function of each operation rounding error. The martingale directly stems from this error model without additional terms. In Chapter 5, which is an extension of the paper "Bounds on Non-Linear Errors for Variance Computation with Stochastic Rounding" preprint (submitted 2023), E-M. El Arar, D. Sohier, P. de Oliveira Castro, and E. Petit (see [START_REF] Arar | Bounds on non-linear errors for variance computation with stochastic rounding[END_REF]), we focus on the probabilistic error analysis of algorithms under SR whatever the nature of errors. Using the Doob-Meyer decomposition [START_REF] Doob | Stochastic Processes, Probability and Statistics Series[END_REF], we introduce a general framework that enables computing probabilistic error bounds for algorithms involving both linear and non-linear sources of error under SR. To the best of our knowledge, this represents the first theoretical exploration of non-linear problems with SR.

Section 5.1 shows that under SR, the error of an algorithm can be decomposed into a martingale plus a drift. We demonstrate that the drift is zero for algorithms with multi-linear error while using the AH or BC method to the martingale term leads to probabilistic error bounds in O( √ nu). Furthermore, we have shown that in the general case, the drift is dominated by the martingale at the first order, which allows us to derive tight probabilistic error bounds.

We apply this general framework to two algorithms that compute the variance: textbook-variance in Section 5.2, and two-pass-variance in Section 5.3. In the case of the textbook-variance algorithm, we exploit the fact that one part of the error constructs directly a martingale, and we use the generalization to the remaining part of the error. We use the generalization to the entire error for the two-pass-variance algorithm, and for both algorithms, we show probabilistic error bounds in O( √ nu). The applicability of this generalization to these algorithms illustrates its flexibility and adaptability across various situations.

We also proposed an alternative approach based on similar techniques to handle the variance computation problem. BC and AH methods in Sec- No study to date has theoretically examined the division case under SR. In this thesis, we have investigated the division problem through two examples: the computation of Møller-Plesset Perturbation Theory (MP2), a method used to estimate the correlation energy of molecules. This problem was proposed by Anthony Scemama from the "Laboratoire de Chimie et Physique Quantiques (LCPQ)." The second example is the computation of the inverse of the sum of n real numbers, for which we proposed a computation model to estimate the error. For both problems, we demonstrated probabilistic error bounds in O( √ nu).

In conclusion, in this thesis we have demonstrated that using stochastic rounding allows us to obtain probabilistic error bounds in O( √ nu). Using the BC method, we have established tight probabilistic bounds for large problem size n, which gives an interest to use SR in low-precision. Moreover, the generalization presented in Section 5.1 provides a functional framework for the probabilistic analysis of algorithms under SR.

15

-Floating-point Arithmetic and Stochastic Rounding

This chapter comprehensively overviews the floating-point background used in this thesis. Section 2.4 presents a state of the art regarding error analysis with stochastic rounding and describes two stochastic rounding modes: SR-nearness and SR-up-or-down. Sub-section 2.4.2 shows that SR-nearness is unbiased and satisfies the mean independence property, an assumption weaker than independence yet powerful enough to yield significant results by martingale theory. We review the works proposed by Connelly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] and Ipsen, and Zhou [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] that show probabilistic bounds on the forward error for the inner product proportional respectively to n ln(n)u and √ nu rather than to the deterministic bound that is proportional to nu. We point out the differences in these two works and illustrate the advantages of each method. • e ∈ Z is the exponent.

. Floating-Point Representation

The precision p in the representation (2.1) is the number of significant digits or bits that can be used to represent the fractional and integral parts of x. The exponent e belongs to a subset of Z defined by a specific range. This range is determined by the standard or specification of the floating-point format being used.

Modern computers use the IEEE-754 standard [START_REF]IEEE standard for floating-point arithmetic, IEEE Std[END_REF] for implementing floatingpoint operations that define different formats, such as binary-32 (single precision) and binary-64 (double precision). These formats ensure consistent representation and interoperability across different computer systems. However, some recent hardware has proposed other formats to improve numerical precision, reduce memory storage, and consume less energy. One such format is bfloat-16, originally proposed by Google and formalized by Intel [START_REF]Bfloat16-hardware numerics definition[END_REF]. Because the representation of numbers is finite, e min and e max are also finite, and their values depend on the normalization used for the mantissa. Note that decimal formats also exist, providing alternative representations for real numbers in computer systems. Some values require special encoding and cannot be expressed using the representation (2.1) such as +0, -0, +∞, and -∞. We also have • NaN (not a number): Any invalid operation will return a NaN. For instance √ -2 or 0/0.

• subnormal numbers (also called denormalized numbers) are represented by setting the exponent bits to all 0s, i.e., e = e min , and using a non-zero mantissa.

In this dissertation, we don't take into account special FP values such as underflow, overflow, denormals, and NaNs. Detailed information on the floatingpoint format most generally in use in current computer systems is defined in the IEEE-754 standard. For a comprehensive understanding of floating-point arithmetic, we highly recommend referring to [START_REF] Muller | Handbook of floatingpoint arithmetic[END_REF] and [START_REF] Boldo | Floating-point arithmetic[END_REF] for further details.

. Rounding Error

Certain inputs are not exact representations in FP numbers. Additionally, the results of elementary floating-point operations are typically not precisely represented as FP numbers. This leads to a fundamental aspect of floatingpoint arithmetic known as rounding error. In order to specify a floating-point arithmetic system, it is necessary to establish a protocol for rounding the result of an operation to a representable FP number.

Denote F the set of normal floating-point numbers. For a real number x, upward rounding ⌈x⌉ and downward rounding ⌊x⌋ are defined by: ⌈x⌉ = min{y ∈ F : y ≥ x}, ⌊x⌋ = max{y ∈ F : y ≤ x}, and by definition, ⌊x⌋ ≤ x ≤ ⌈x⌉, with equalities if and only if x ∈ F. The IEEE-754 standard defines five rounding modes.

• Round to nearest, ties to even: rounds x to the nearest FP number. If x falls exactly halfway between two consecutive FP numbers, it is rounded to the FP number whose least significant digit is even. This is the default rounding mode in the IEEE-754 standard.

• Round to nearest, ties away from zero: rounds x to the nearest FP number. If x falls exactly halfway between two consecutive FP numbers, it is rounded to the FP number with the larger magnitude.

• Round toward +∞: round x to ⌈x⌉.

• Round toward -∞: round x to ⌊x⌋.

• Round toward 0: rounding to the representable number closest to zero.

For a given rounding mode, denote fl(x) the floating-point approximation of a real number x ̸ = 0 that is one of ⌊x⌋ or ⌈x⌉. The relative error δ of this approximation is given by:

δ = fl(x) -x x .
(2.2)

Thus fl(x) = x(1 + δ) ∈ F.
Definition 1. The unit roundoff u of a radix β, precision p, is defined as

u = β 1-p .
Theorem 1. The relative error δ in the Equation (2.2) satisfies |δ| ≤ u.

While the IEEE-754 mode RN (round to nearest, ties to even) has the stronger property that |δ| ≤ 1 2 u.

Proof. Since x = m × β e-p where β p-1 ≤ |m| < β p , we have β e-1 ≤ |x| < β e . If |x| = β e-1 then x is a FP number and thus δ = 0. Otherwise, 1 |x| < β 1-e and the distance between two consecutive FP numbers is β e-p . So, for round to nearest

|fl(x) -x| ≤ 1 2 β e-p , Otherwise, |fl(x) -x| ≤ β e-p .
Hence,

|δ| = fl(x) -x x ≤ u.
Let x, y ∈ F and op ∈ {+, -, * , /}. For IEEE-754 rounding modes [START_REF]IEEE standard for floating-point arithmetic, IEEE Std[END_REF] and stochastic rounding [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF], the standard model defines the approximation fl(x op y) as:

fl(x op y) = (x op y)(1 + δ).

(2.3)

Assume that x is a real number that is not representable: x ∈ R \ F. The machine-epsilon or the distance between the two FP numbers enclosing x is ϵ(x) = ⌈x⌉ -⌊x⌋ = β e-p . Since β p-1 ≤ |m| < β p , then β e-1 ≤ |x| < β e and |ϵ(x)| = β e-1 u ≤ |x|u.

(2.4)

The fraction of ϵ(x) rounded away, as shown in figure 2

.1, is p(x) = x-⌊x⌋ ⌈x⌉-⌊x⌋ ⌊x⌋ ⌈x⌉ x 1 2 ϵ(x) p(x)ϵ(x) Figure 2.1: p(x)
is the fraction of ϵ(x) to be rounded away.

We note x , the greatest integer less than or equal to x. The following lemma gives an important property of downward rounding. Lemma 1. Let x ∈ R \ F. β p-e ⌊x⌋ = β p-e x , where e is the exponent.

Proof. We know that β p-e ⌊x⌋, β p-e ⌈x⌉ ∈ Z, and ⌊x⌋ < x < ⌈x⌉, then β p-e ⌊x⌋ < β p-e x < β p-e ⌈x⌉. We thus have β p-e ⌊x⌋ ≤ β p-e x < β p-e ⌈x⌉.

Since ⌈x⌉ -⌊x⌋ = β e-p , then β p-e ⌈x⌉ -β p-e ⌊x⌋ = 1 and β p-e ⌊x⌋ ≤ β p-e x < β p-e ⌊x⌋ + 1.

. Deterministic Error Analysis

As stated previously in Section 2.2, rounding errors arise when representing a non-representable number as the nearest FP number. The accumulation of these errors can significantly reduce the accuracy of the computed result. In this section, we present some of the deterministic methods proposed to investigate rounding errors.

. Forward and Backward Error

To the best of our knowledge, the earliest proposal of backward error analysis was by Turing [START_REF] Turing | Rounding-off errors in matrix processes[END_REF], and Von Neumann and Goldstine [START_REF] Neumann | Numerical inverting of matrices of high order[END_REF]. As summarized by Higham [START_REF]What is backward error?[END_REF]:

Backward error is a measure of error associated with an approximate solution to a problem. Whereas the forward error is the distance between the approximate and true solutions, the backward error is how much the data must be perturbed to produce the approximate solution.

In other words, forward error is the error between a computed result ŷ and the exact computation y, while backward error is the error in inputs that allows to compute ŷ. The condition number seems first to have been used by Turing [START_REF] Turing | Rounding-off errors in matrix processes[END_REF] and is determined by the maximum value of the quotient: relative forward error relative backward error , it measures the sensitivity of the solution to small perturbations in the input data [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF].

One approach to evaluate computation accuracy is establishing an upper bound on the forward error. Let us consider the computation fl(x op y) = (x op y)(1 + δ), where op ∈ {+, -, * , /}. So, the forward error δ satisfies |δ| ≤ u. In the following, we study the summation algorithm as an example, and we bound its relative error. For real numbers a 1 , . . . , a n , denote s = n i=1 a i . The computed ŝ can be expressed as

ŝ = (((a 1 + a 2 )(1 + δ 2 ) + a 3 )(1 + δ 3 ) . . . + a n ) (1 + δ n ) = n i=1 a i n k=max (2,i) (1 + δ k ).
The forward error satisfies

|ŝ -s| = n i=1 a i   n k=max (2,i) (1 + δ k ) -1   ≤ n i=1 |a i | n k=max (2,i) (1 + δ k ) -1 by triangle inequality ≤ n i=1 |a i | (1 + u) n-1 -1 .
Finally, the relative error satisfies

|ŝ -s| |s| ≤ n i=1 |a i | | n i=1 a i | γ n-1 (u), (2.5)
where

γ n (u) = (1 + u) n -1 = nu + O(u 2 ) for nu ≪ 1.
The forward error of a summation of n floating point numbers is proportional to nu. The quantity

n i=1 |a i | | n i=1 a i|
is the condition number of n i=1 a i using the 1-norm.

In order to evaluate the accuracy of a calculation, we can also use interval arithmetic that provides a systematic approach to handle uncertainty and errors in numerical computations. In the next sub-section, we briefly explain how interval arithmetic works.

. Interval Arithmetic

Interval arithmetic, as introduced by Moore in 1963 [START_REF] Moore | Interval arithmetic and automatic error analysis in digital computing[END_REF] and developed in [START_REF] Neumaier | Interval methods for systems of equations[END_REF][START_REF] Rump | Verification methods: Rigorous results using floating-point arithmetic[END_REF], offers a systematic way to handle uncertainty and imprecision in computations. It operates by representing a quantity as an interval, with lower and upper bounds containing all the possible computation values. The elementary operations are redefined to handle intervals operands and guarantee that the resulting interval provides rigorous bounds on the computation.

For instance, let x, y be such that x belongs to the interval [⌊x⌋; ⌈x⌉] and y belongs to the interval [⌊y⌋; ⌈y⌉]. Then, the addition [⌊x⌋; ⌈x⌉] + [⌊y⌋; ⌈y⌉] can be implemented as follows

[⌊x⌋; ⌈x⌉] + [⌊y⌋; ⌈y⌉] = [⌊z⌋; ⌈z⌉],
where ⌊z⌋ = ⌊x⌋ + ⌊y⌋ and ⌈z⌉ = ⌈x⌉ + ⌈y⌉ are obtained with rounding toward -∞ and +∞, respectively.

Interval arithmetic is a powerful tool that provides an interval (often with the minimal range) that encompasses all possible exact results [START_REF] Boldo | Floating-point arithmetic[END_REF]. It is possible to refine the analysis by considering a sophisticated object such as zonotopes [START_REF] Goubault | Static analysis of numerical algorithms[END_REF]. However, its practical application is costly, as it requires abstract approximation methods. Moreover, due to the conservative nature of intervals, these methods tend to return intervals too large when the algorithm is complex. The error bound can be achievable in these cases but often is not significant and misrepresents experimental results. In these cases, statistic techniques based on Monte Carlo arithmetic can help understand and optimize complex HPC programs. Several extensions have been proposed to the original model, including affine arithmetic [START_REF] De Oliveira Castro | High Performance Computing code optimizations: Tuning performance and accuracy[END_REF] and Taylor models [START_REF] Nedialkov | Interval arithmetic, affine arithmetic, taylor series methods: why, what next?[END_REF].

. State of the Art of Stochastic Rounding

Forsythe has presented the use of stochastic rounding to reduce the accumulation of rounding errors [START_REF] Forsythe | Round-off errors in numerical integration on automatic machinery. preliminary report[END_REF]. John von Neumann and Goldstine first proposed modeling rounding errors as random errors to obtain probabilistic error bounds in the 1940s [START_REF] Neumann | Numerical inverting of matrices of high order[END_REF]. They treated rounding errors as random variables with known averages and maximum values of β -s /2, where β is the basis of the digital representation, and 2s is the number of places used. They observed that if we assume m independent random variables equi-distributed in [-β -s /2; β -s /2], the maximum of the sum of these random variables is bounded by mβ -s /2, while the dispersion is given by √ mβ -s / √ 12. This encouraged the consideration of probabilistic analysis instead of deterministic analysis. This section provides a comprehensive overview of the research done on estimating numerical errors using stochastic arithmetic.

. Stochastic Arithmetic

The Stochastic Arithmetic field proposes automatic methods for estimating the number of significant digits for complex programs. Two main methods have been proposed: MCA "Monte Carlo Arithmetic" by Parker [START_REF] Parker | Monte Carlo Arithmetic: exploiting randomness in floatingpoint arithmetic[END_REF] (equivalent to SR-nearness in the following), CESTAC "Contrôle et Estimation Stochastique des Arrondis de Calculs" by Vignes [START_REF] Vignes | Discrete stochastic arithmetic for validating results of numerical software[END_REF] (equivalent to SR-up-or-down in the following). The idea is to substitute the error term δ within each operation with a random variable that simulates the rounding error and to run the computation several times while storing the various outcomes. Then, a statistical analysis is applied to the set of samples in order to assess the quality of the result.

Throughout this dissertation, x = fl(x) is the approximation of the real number x under stochastic rounding.

SR-up-or-down involves rounding a floating-point number x up or down with probability 1 2 .

x = ⌈x⌉ with probability This mode can be expressed [60, p. 34] in terms of p(x): since the two outcomes of SR-up-or-down mode are equiprobable, we have E(x) = ⌈x⌉+⌊x⌋ 2 , which allow us to write the bias as

E(x -x) = ⌈x⌉ + ⌊x⌋ 2 -x, because p(x) = x-⌊x⌋ ⌈x⌉-⌊x⌋ , E(x -x) = (⌈x⌉ -⌊x⌋)( 1 2 -p(x)) = ϵ(x)( 1 2 -p(x)).
Thus, we conclude that SR-up-or-down is biased and the expected value depends on p(x) and ϵ(x). SR-nearness involves rounding a FP number x up or down with probability depending on p(x) = x-⌊x⌋ ⌈x⌉-⌊x⌋ .

x = ⌈x⌉ with probability p(x), ⌊x⌋ with probability 1 -p(x), ⌊x⌋ ⌈x⌉ x 1 -p(x) p(x) This definition does not include exceptional cases such as overflow, underflow, and rounding of infinities and NaNs. A definition addressing these limit cases can be found in [13, sec 5.a]. However, our analysis does not take these special cases into account.

The rounding SR-nearness mode is unbiased

E(x) = p(x)⌈x⌉ + (1 -p(x))⌊x⌋ = p(x)(⌈x⌉ -⌊x⌋) + ⌊x⌋ = x.
We highlight that there exists an alternative expression to present stochastic rounding. Consider a real number x such that x / ∈ F. Let x be the random variable of the distribution of results after random rounding of x. Then x = random_round(x) = round p (x + β e-p ξ), (2.6) where ξ is a random variable that can be discrete or continuous and round p is the default IEEE-754 rounding mode to the nearest with p precision (for more details, we refer to [17, sec 5.4.3]).

One way to evaluate the numerical accuracy of computation is the number of significant digits which measures the relative error. The key idea is to count the number of accurate digits in the floating-point mantissa against a reference. Stott Parker [START_REF] Parker | Monte Carlo Arithmetic: exploiting randomness in floatingpoint arithmetic[END_REF] proposed that the number of significant digits can be defined as -log( σ µ ), where µ represents the mean and σ represents the standard deviation. Sohier et al. [START_REF] Sohier | Confidence intervals for stochastic arithmetic[END_REF] have demonstrated that this proposition holds true but under the normality assumption with a probability 0.68. They also introduce a new quantity of interest: the number of digits contributing to the accuracy of the final result.

As mentioned in the introduction of this chapter, SR can be used as a replacement for RN in numerical computation. The growing interest in SR arises from its unbiased property and its positive effect in various domains, such as neural networks and ODEs, especially in low-precision. In the next sub-section, we will review the theoretical investigations proposed to analyze errors with SR. We show that rounding errors under SR satisfy the mean independence property, an assumption weaker than independence yet powerful enough to yield important results by martingale theory.

. Error Analysis via Martingales

Wilkinson [71, sec 1.33] had the intuition that the roundoff error accumulated in n operations is typically proportional to √ nu rather than nu. In order to validate this intuition, several works have addressed the issue of estimating the error of computation under stochastic rounding.

Motivated by the unbiased nature of SR-nearness and assuming independent random errors, Higham and Mary [START_REF] Higham | A new approach to probabilistic rounding error analysis[END_REF] have shown that for the inner product, a probabilistic bound of the error proportional to n ln (n)u can be achieved rather than the deterministic bound in nu. However, in general, and under SR-nearness, the error terms in algorithms appear as a sequence of random variables such that the independence property does not hold. However, a weaker yet fruitful assumption, called mean independence, does.

Definition 2. A random variable Y is said to be mean independent from ran- dom variable X if its conditional mean E[Y /X] = E(Y ). The random sequence (X 1 , X 2 , , . . . , ) is mean independent if E[X k /X 1 , ..., X k-1 ] = E(X k ) for all k.
Proposition 1. Let X and Y be two real random variables:

1. If X and Y are independent then X is mean independent from Y .

2. If X is mean independent from Y then X and Y are uncorrelated.

The reciprocals of these two implications are false.

Let a, b ∈ F and c ← a op b the exact result of an elementary operation op ∈ {+, -, * , /}. Under SR-nearness, the relative error δ, such that ĉ = (a op b)(1 + δ), is a random variable satisfying E(δ) = 0 and |δ| ≤ u.

In the following, we recall [12, Lem 5.2] (and its proof by adapting the same notations) which shows that SR-nearness errors satisfy the mean independence property. Moreover, (⌈c⌉-c)/c and (⌊c⌋-c)/c are themselves random variables that are entirely determined by δ 1 , • • • , δ k-1 , and so the conditional expectation of each given

Lemma 2. Consider a sequence of elementary operations

c k ← a k op k b k , with δ k the error of their kth operation, that is to say, ĉk = (â k op k bk )(1 + δ k ). Using SR- nearness, the δ k are random variables such that E[δ k /δ 1 , • • • , δ k-1 ] = E(δ k ) = 0.
δ 1 , • • • , δ k-1 is itself. Therefore,
we obtain

E(δ k |δ 1 , • • • , δ k-1 ) =pE ⌈c⌉ -c c δ 1 , • • • , δ k-1 + (1 -p)E ⌊c⌋ -c c δ 1 , • • • , δ k-1 = p ⌈c⌉ -c c + (1 -p) ⌊c⌋ -c c = 0.
Lemma 2 (that has been proven in [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF]Lem 5.2]) will be the fundamental tool for the theoretical analysis of algorithms under SR. It substitutes the independence assumption used for the error analysis under SR by a weaker property satisfied by SR-nearness. The mean independence property is sufficient to improve the error analysis of algorithms with SR-nearness. It leads to obtain a martingale (Definition 3), which is a sequence of random variables such that the expected value of the next value in the sequence, given all the past values, is equal to the current value. Using Azuma-Hoeffding inequality [53, p. 303], allows to obtain probabilistic bounds on the error in O( √ nu).

The full process is detailed after Lemma 3 for the inner product.

Definition 3. A sequence of random variables M 1 , • • • , M n is a martingale with respect to the sequence X 1 , • • • , X n if, for all k, • M k is a function of X 1 , • • • , X k , • E(|M k |) < ∞, and
• E[M k /X 1 , • • • , X k-1 ] = M k-1 . If E[M k /X 1 , . . . , X k-1 ] ≥ M k-1 , M 1 , . . . , M n is called sub-martingale. Lemma 3. (Azuma-Hoeffding inequality). Let M 0 , • • • , M n be a martingale with respect to a sequence X 1 , • • • , X n . We assume that there exist a k < b k such that a k ≤ M k -M k-1 ≤ b k for k = 1 : n. Then, for any A > 0 P(|M n -M 0 | ≥ A) ≤ 2 exp - 2A 2 n k=1 (b k -a k ) 2 .
In the particular case

a k = -b k and λ = 2 exp -A 2 2 n k=1 b 2 k we have P   |M n -M 0 | ≤ n k=1 b 2 k 2 ln(2/λ)   ≥ 1 -λ,
where 0 < λ < 1.

Remark 1. In Lemma 3, if b k is equal to a constant b, n k=1 b 2 k = n k=1 b 2 = |b| √ n.
Interestingly, this inequality shows a √ n in the final bound, similar to the results obtained by applying the Central Limit Theorem (CLT). While the CLT is only applicable as n → ∞, this inequality establishes bounds valid for all finite values of n, which is of great interest since, in numerical computation, the problem size is always a finite real number.

Under SR-nearness, the inner product

y = a ⊤ b, where a, b ∈ R n is de- fined as ŷ = n i=1 a i b i (1 + δ 2(i-1) ) n k=i (1 + δ 2k-1
). Since each |δ k | ≤ u, the worst case of the forward error of the computed ŷ is in O(nu). Based on the mean independence of errors established in Lemma 2, Connelly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] and Ipsen, and Zhou [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] have investigated this problem for SR-nearness. Both works build on the mean independence property of SR-nearness. This allows them to form a martingale, and then to apply the Azuma-Hoeffding concentration inequality. The difference between these two works is in the way they form the martingale. In [12, sec 3], the martingale is built using the errors accumulated in the whole process

ψ i = (1 + δ 2(i-1) ) n k=i (1 + δ 2k-1 )
for all 1 ≤ i ≤ n. This approach uses the inclusion-exclusion principle to generalize the bound to the summation which results in a pessimistic n in the probability.

They prove

|ŷ -y| |y| ≤ Kγ n (λ), with probability at least 1 -2n exp -λ 2 2 , where K = n i=1 |a i b i | | n i=1 a i b i| and γn (λ) = exp λ √ nu+nu 2 1-u -1.
The factor n in the probability disrupts the √ nu property.

δ = 2n exp -λ 2 2 implies that λ = 2 ln (2n/δ) and |ŷ -y| |y| ≤ Kγ n 2 ln (2n/δ) , (AH1-IP) with probability at least 1 -δ. When nu ≪ 1, we have γn ( 2 ln (2n/δ)) = exp 2n ln (2n/δ)u + nu 2 1 -u -1 = u 2n ln (2n/δ) + O(u 2 ) = u √ 2n ln 2n -2n ln δ + O(u 2 ) = O(u n ln n).
The probabilistic bound in O(u n ln n) is better than the deterministic bound in O(nu). But, it is possible to obtain a probabilistic bound in O( √ nu).

In [46, sec 4], the martingale is formed by following step-by-step how the error accumulates in the recursive summation of the inner product. In particular, the authors distinguish between the multiplications and additions computed at each step and carefully monitor their mean independences. This approach leads to the following probabilistic bound

|ŷ -y| |y| ≤ K uγ 2n (u) ln(2/δ), (AH2-IP)
with probability at least 1 -δ. This technique avoids the inclusion-exclusion principle and when nu ≪ 1, it leads to

uγ 2n (u) ln(2/δ) = u √ 2n ln 2 -2n ln δ + O(u 2 ).
Note that when nu ≪ 1, (AH1-IP) and (AH2-IP) differ only in the factor √ ln n that appears in (AH1-IP) due to the use of the martingale property on each partial sum necessitating to use the inclusion-exclusion principle. The effect of this factor on the bound behavior is comprehensively illustrated in Section 4.4

through numerical experiments. All in all, (AH2-IP) is proportional to u √ n, while (AH1-IP) is proportional to u √ n ln n.

. Stochastic Rounding and Applications

Recent theoretical/numerical developments have revealed that SR provides better results than the IEEE-754 default rounding mode in multiple domains [START_REF] Croci | Stochastic rounding: implementation, error analysis and applications[END_REF]. Connolly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] have shown that numerous numerical linear algebra algorithms, including the inner product and the triangular system solution, are unbiased when using SR. Moreover, the forward error of these algorithms has a probabilistic bound in O( √ nu) instead O(nu) for the deterministic bound.

The positive effect of SR also extends to calculating the solution of ordinary differential equations in low-precision. Several studies [START_REF] Hopkins | Stochastic rounding and reduced-precision fixed-point arithmetic for solving neural ordinary differential equations[END_REF][START_REF] Fasi | Algorithms for stochastically rounded elementary arithmetic operations in ieee 754 floating-point arithmetic[END_REF] have shown that incorporating SR into numerical integration schemes for ODE solvers yields notable benefits. By introducing randomness into the rounding process, SR reduces the accumulation of round-off errors and mitigates the numerical instability that often arises in long-term simulations. Hopkins et al. [START_REF] Hopkins | Stochastic rounding and reduced-precision fixed-point arithmetic for solving neural ordinary differential equations[END_REF] have demonstrated that fixed-point ODE solvers exhibit greater robustness when using SR than other rounding algorithms. Furthermore, through experimental analysis, they have established that fixed-point ODE solvers with stochastic rounding achieve higher accuracy compared to single-precision floatingpoint ODE solvers. Their study also indicated that utilizing just 6 bits in the residual is sufficient for a high-performance stochastic rounding algorithm within an ODE solver.

For partial differential equations (PDEs), Croci and Giles [START_REF] Croci | Effects of round-to-nearest and stochastic rounding in the numerical solution of the heat equation in low precision[END_REF] have conducted a study on the accumulation of rounding errors in the solution of the heat equation in low-precision using Runge-Kutta finite difference methods with RN and SR. They demonstrate the implementation of a scheme that effectively reduces rounding errors and derives a priori estimates for both local and global rounding errors. While the worst-case scenario for local errors is O(u) with respect to the discretization parameters (mesh size and timestep), the RN solution always stagnates for small enough ∆t. Until stagnation occurs, the global error grows at a rate of O(u∆t -1 ). In contrast, stagnation and the accumulation of rounding errors can be avoided with SR. They prove that the global rounding errors are only O(u∆t -1/4 ) in one dimension and essentially bounded (up to logarithmic factors) in higher dimensions.

In the domain of neural networks, the utilization of SR is not a new concept. It was initially introduced as probabilistic rounding by Höhfeld and Fahlman [START_REF] Höhfeld | Probabilistic rounding in neural network learning with limited precision[END_REF][START_REF] Hoehfeld | Learning with limited numerical precision using the cascade-correlation algorithm[END_REF] in 1992, specifically within the context of the Cascade-Correlation algorithm. Gupta et al. [START_REF] Gupta | Deep learning with limited numerical precision[END_REF] subsequently demonstrated that deep networks could be successfully trained in half-precision fixed-point arithmetic by using SR, with little or no degradation in the classification accuracy. Su et al. [START_REF] Su | Towards high performance low bitwidth training for deep neural networks[END_REF] have shown that SR can be used for training deep neural networks in 8-bit fixedpoint arithmetic using SR and analyzed the success of SR in this context through simulation and experiments. Wang et al. [START_REF] Wang | Training deep neural networks with 8-bit floating point numbers[END_REF] also employed SR for training neural networks using 8-bit floating-point arithmetic, resulting in a reduction of bit-precision for additions down to 16 bits instead 36 bits as well as speed improvements of approximately 2 -4 times compared to 32-bit training.

Paxton et al. [START_REF] Paxton | Climate modeling in low precision: Effects of both deterministic and stochastic rounding[END_REF] study the effectiveness of low-precision arithmetic for climate simulations, especially the effects of RN and SR in chaotic ODE and PDE systems related to climate modeling: the Lorenz system, heat diffusion, a nonlinear shallow water approximation for flow over a ridge, and a coarse resolution global atmospheric model with simplified parametrizations. They find that SR can effectively mitigate rounding errors across various applications, and the results also provide evidence that SR could be relevant to nextgeneration climate models.

Kimpson et al. [START_REF] Kimpson | Climate-change modelling at reduced floating-point precision with stochastic rounding[END_REF] have demonstrated the possibility of modeling a changing climate system using SR in conjunction with reduced-precision FP numbers. Their study reveals that employing SR significantly improves the performance of half-precision computations, making them comparable to the solutions achieved with single-precision computations. Various applications, including quantum mechanics and quantum computing, utilize SR to enhance their outcomes. For a comprehensive overview of several applications using SR, we recommend [13, sec 7].

. Simulation of SR in Software

One approach to implementing a mathematical operator with stochastic rounding is to follow three key steps. First, the operator must be evaluated using high-precision floating-point arithmetic to obtain a more accurate result. Second, drawing a pseudo-random number from some uniform distribution to determine the rounding direction. Third, the high-precision result is rounded to the desired precision level, ensuring that the final value aligns with the intended level of precision required for the application. The availability of arithmetic operations beyond the working precision simplifies the implementation of this approach. This can be achieved through hardware methods, such as emulating binary32 using binary64 format or through software employing arbitrary precision libraries like the GNU Multiprecision Library (GNU MPFR) [START_REF] Fousse | Mpfr: A multiple-precision binary floating-point library with correct rounding[END_REF].

Once the high-precision result is available, the rounding step can be performed in several ways. In the case of the MATLAB function chop [START_REF] Higham | Simulating low precision floating-point arithmetic[END_REF], the random numbers are compared to a threshold value drawn from the uniform distribution over ]0; 1[, and the rounding operation is performed based on this comparison. For instance, let x = x(1 + δ) and p(x) = x-⌊x⌋ ⌈x⌉-⌊x⌋ . If p(x) is greater than the threshold, x is rounded up; otherwise, it is rounded down.

In the case of verificarlo [START_REF] Denis | Verificarlo: Checking floating point accuracy through Monte Carlo arithmetic[END_REF] that is built upon the LLVM compiler [START_REF] Lattner | Llvm: A compilation framework for lifelong program analysis & transformation[END_REF], the random numbers are rounded by adding a random noise (SR-nearness: the random variable ξ in the Equation (2.6) is uniformly distributed on ]-1 2 ; 1 2 [) and using the default rounding mode in the IEEE-754 standard. Three expressions are possible: Random Rounding (RR) which introduces perturbation only on the output, Precision Bounding (PB) which introduces perturbation only on the input, and Full MCA (MCA) which introduces perturbation on operand(s) and the result. We recommend referring to [17, sec 6] for a comprehensive exploration and in-depth explication of the work process of SR-nearness in Verificarlo.

Verrou [START_REF] Févotte | Verrou: a cestac evaluation without recompilation[END_REF] is built upon Valgrind [START_REF] Nethercote | Valgrind: a framework for heavyweight dynamic binary instrumentation[END_REF]. Verrou intercepts floating point operations at runtime and replaces them with their random rounding counterparts. It uses both SR-up-or-down (the random variable ξ in the Equation (2.6) is uniformly distributed on ] -p(x); 1 -p(x)[) and SR-nearness. The Interflop project establishes a shared interface between Verrou and Verificarlo, facilitating tool interoperability. In both tools, the end-user performs a statistical analysis to conclude the numerical quality of the results.

There are other tools available for simulating SR in software. Cadna [START_REF] Jézéquel | Cadna: a library for estimating roundoff error propagation[END_REF] computes three times the result with SR-up-or-down and estimates statistically the numerical error. Fasi and Mikaitis in [START_REF] Fasi | Algorithms for stochastically rounded elementary arithmetic operations in ieee 754 floating-point arithmetic[END_REF] have proposed two algorithms for emulating stochastic rounding for both square root and the elementary operations {+, -, * , /}. They showcase the value of these algorithms through diverse applications where stochastic rounding is favorable. Klöwer's Julia software package StochasticRounding.jl [START_REF] Paxton | Climate modeling in low precision: Effects of both deterministic and stochastic rounding[END_REF] uses integer operations to perform the stochastic rounding. It exports three floating-point formats with SR, Float32sr, Float16sr, and BFloat16sr.

Each of these tools uses a random number generator for simulation purposes. The choice of the generator directly impacts the quality of the simulation. Therefore, opting for a generator with a sufficiently large period that is evenly distributed relative to the characteristics of floating-point numbers is essential.

In this dissertation, we use Verificarlo [START_REF] Denis | Verificarlo: Checking floating point accuracy through Monte Carlo arithmetic[END_REF] to simulate the stochastic rounding errors. In the next chapter, we will focus on SR as a rounding mode. We analyze the biases of the two stochastic rounding modes presented in Subsection 2.4.1: SR-nearness and SR-up-or-down. We demonstrate that IEEE-754 default rounding modes and SR-up-or-down accumulate rounding errors across iterations, while SR-nearness, being unbiased, does not.

-Rectangular Rule with Stochastic Rounding

This chapter presents our first contribution: we compare SR-nearness and SR-up-or-down in the computation of integrals using rectangular integration, which is the basis of Euler's explicit method for ordinary differential equations (ODEs). Our investigation focuses on the constant function and the cosine function. We show that SR-nearness remains unbiased in these two examples. However, an exact expression and an estimation of the bias are given for SR-up-or-down. We show how the accumulation of errors with both SR-upor-down and IEEE-754 modes leads to results significantly less accurate than with SR-nearness. Additionally, we provide an expression for the method error of the cosine function, which supports the numerical observations.

. Integrating a Constant Function

Rectangular integration rule is a classic approximation for performing numerical integration: the area under a curve is approximated by a sum of

N rectangle areas. b a f (t)dt ≈ N -1 k=0 hf (a + kh)
where h = b-a N . In particular, using Euler's forward method, the rectangular rule is one of the resolution techniques for ODE.

Verrou's tutorial [START_REF] Févotte | Etudier la qualité numérique d'un code avec Verrou[END_REF] integrates the cosine function with the rectangular rule; with deterministic round to nearest or SR-up-or-down modes, the solution is biased. When the number of integration steps grows, this bias can become high and degrade the quality of the solution. In this section, we show why deterministic and SR-up-or-down modes can be biased with the rectangular rule and how the accumulation of errors with both previous modes leads to results significantly less accurate than the unbiased mode SR-nearness.

We perform the analysis on a constant function f (t) = 1 for all t ∈ [0; 1]. With f constant, the evaluation error is zero, making it clear how the numerical error accumulates on the summation.

Denote x = 1 = N -1
k=0 h, where h = 1/N . The distribution x is produced by summing N times the integration step h. We note ŝk the random variable for the partial sum at step 0 ≤ k ≤ N -1 and s k the exact expected result, with ŝN-1 = x.

SR-up-or-down: As shown before, for each ŝk we introduce a bias corre- sponding to

E(ŝ k -s k ) = ϵ(s k )( 1 2 -p(s k )), from the definition of p(s k ), we have 0 < p(s k ) < 1, then -1 2 < 1 2 -p(s k ) < 1 2 and |E(ŝ k -s k )| < 1 2 ϵ(s k ).
Table 3.1 shows these different values for N = 20. Interestingly, in this table, we note that p(s k ) is constant between two successive powers of the base except for the first value. For example for 9 < k < 20, s k stays within [2 -1 ; 2 0 ) and both p(s k ) and E(ŝ k -s k ) are constant. In the following, we show why that is always the case. Suppose s k , s k+1 ∈ [β e ; β e+1 ). Then ϵ(s k ) = β e-p . At each step, the next partial sum is computed as, s k+1 = ŝk + h. In that case, using the lemma 1, we have

k s k p(s k ) E(ŝ k -s k ) ϵ(s k ) 2
p(s k+1 ) = β p-e (s k+1 -⌊s k+1 ⌋) = β p-e s k+1 -β p-e s k+1 = β p-e ŝk + β p-e h -β p-e ŝk + β p-e h .
Since ŝk ∈ F, we have β p-e ŝk ∈ Z and β p-e ŝk + β p-e h = β p-e ŝk + β p-e h . Finally p(s k+1 ) = β p-e h -β p-e h . Thus p(s k+1 ) depends only on h and e. Recursively for all l > 0 satisfying

s k+l ∈ [β e ; β e+1 ), p(s k+l ) = β p-e h -β p-e h is constant. The bias E(ŝ k+l -s k+l ) = ϵ(s k+l )( 1 2 -p(s k+l )) = β e-p ( 1 2 -β p-e h -β p-e h ),
is also constant in this interval.

Between two successive powers of the base, p(s) remains constant, as well as the bias. Because the bias is constant (and, consequently, its sign too), it accumulates across iterations. The total bias can be written as

E(x -x) = N -1 k=0 E(ŝ k -s k ) = N -1 k=0 ϵ(s k )( 1 2 -p(s k )),
and

|E(x -x)| < 1 2 N -1 k=0 ϵ(s k ).
For a large N , we can neglect the effect of the first partial sum in each power-of-the-base interval. As we can observe in Table 3.1 the last power-ofthe-base interval, [ 1 2 ; 1[ contains more summation terms and has a larger ϵ.

Then, its bias usually dominates in the final result.

Numerical experiment. The computations were done thanks to Python

BigFloat arbitrary precision library. SR and RN values were computed in bi-nary32 and the reference was computed in binary64. We have verified numerically that the above expression for the bias closely predicts the bias measured with SR-up-or-down.

We consider a fixed number of iterations N. We ran one time the C program in Listing 1 with each of the two previously defined stochastic rounding modes as well as round to nearest.

The following figure plots the three distributions over N . On the other hand, SR-up-or-down mode and RN-binary-32 samples have a bias, which confirms the previous results for SR-up-or-down mode. The maximal amplitude of the bias for both SR-up-or-down and RN-binary-32 increases with N because of errors accumulation. The bias is reproducible and constant across different runs.

. Integrating the Cosine Function

Another example that illustrates the effect of rounding errors in numerical computations is the integral of the cosine function. In the following we consider the evaluation of float dx = (pi/2)/N; float s = 0.0; float x = dx/2; for (int i=0 ; i < N ; i++) { s += dx*cos(x);

x += dx; } return s;

Listing 2: Rectangle integration of a cosine function.

It is important to emphasize that here we neglect the error of the cosine function and only study the error accumulated during summation and multiplication. Figure 3.2 illustrates the two distributions across the variable N . The C program in Listing 2 produces two errors. The first error, denoted by e 1 , results from the instruction x += dx, while the second error, denoted by e 2 , results from the statement s += dx*cos(x). Thus, for all k ∈ 0 ; N -1 , it can be stated that:

ŝk = ŝk-1 + dx cos(x k + e k 1 ) + e k 2 .
The mean value theorem implies that there exist

z k ∈ [x k ; x k + e k 1 ] such that: -sin(z k ) = cos(x k + e k 1 ) -cos(x k ) e k 1 .
It follows that

ŝk = ŝk-1 + dx (cos(x k ) -e k 1 sin(z k )) + e k 2 = ŝk-1 + dx cos(x k ) -dx e k 1 sin(z k ) + e k 2 .
The Sub-section 3.1 demonstrates that the errors e k 1 all have the same sign in each interval [β e ; β e+1 ) except the first, implying that the term e k 1 dx sin(z k ) also shares the same sign. By replacing x += dx in Listing 2 with x = dx*k, we observe numerically that the error is negligible, implying that the source of error is e k 1 . This suggests that the bias of e 2 is insignificant compared to e 1 .

These arguments explain the bias accumulation observed in Figure 3.2.

One key aspect in evaluating the accuracy of numerical computations is the numerical error. It encompasses both the rounding error, arising from the rounding function and computational limitations, and the method error, arising from the inherent approximations made when applying a specific numerical technique. The following lemma calculates the method error of the previous program 2.

Lemma 4.

The method error at each step k of this algorithm is given by:

ŝk -s k = dx 2 sin (dx(k + 1)) sin( dx 2 ) -s k .
In particular, for

k = N -1, ŝN-1 -s N -1 = π 4N 1 sin( π 4N ) -1.
Proof. At each step k, we have

ŝk = k p=0 dx cos((p + 1 2 )dx) = dx 2 k p=0 (e idx ) p+ 1 2 + (e -idx ) p+ 1 2 = dx 2 e i dx 2 1 -(e idx ) k+1 1 -e idx + e -i dx 2 1 -(e -idx ) k+1 1 -e -idx = dx 2 e idx(k+1) -1 e i dx 2 -e -i dx 2 + 1 -e -idx(k+1) e i dx 2 -e -i dx 2 = dx 2 e idx(k+1) -e -idx(k+1) e i dx 2 -e -i dx 2 = dx 2 sin (dx(k + 1)) sin( dx 2 )
.

In the particular case k = N -1, s N -1 = 1 and since dx = π 2N we obtain :

ŝN-1 -s N -1 = π 4N sin(N dx) sin( π 4N ) -1 = π 4N 1 sin( π 4N ) -1.
Remark 2. Since for large N , π 4N 1 sin( π 4N ) ≈ 1, the method error converge to 0. As N increases, the method error decreases because the method becomes more accurate with a larger rectangle number, which confirms Remark 2. On the other hand, the numerical error increases due to the accumulation of rounding errors in numerical calculations. For large values of N , the rounding error dominates the method error, representing the overall numerical error in the calculation.

In conclusion, the examples presented in this section demonstrate that rounding errors can accumulate even through simple computations, such as summation, and significantly impact the final result. The stochastic rounding mode SR-nearness is unbiased not only for one elementary operation but, even in other numerical methods such as rectangular integration, it is much closer to the expected value than SR-up-or-down or RN-binary32, in particular for large N . The remainder of this dissertation focuses on SR-nearness and its potential applications in numerical computation. We analyze other algorithms and demonstrate probabilistic bounds on the error in O( √ nu) instead of O(nu) for the deterministic bounds.

-Error Analysis for Algorithms with Multi-Linear Error

Stochastic arithmetic has been used to empirically estimate the numerical error in complex programs. Stochastic rounding (SR) introduces a random perturbation in each floating-point operation, followed by a statistical analysis of the sampled output set to estimate the impact of rounding errors. Methods such as MCA [START_REF] Parker | Monte Carlo Arithmetic: exploiting randomness in floatingpoint arithmetic[END_REF] or CESTAC [START_REF] Vignes | Discrete stochastic arithmetic for validating results of numerical software[END_REF] have been introduced to simulate this effect. However, recent theoretical findings like [START_REF] Arar | The positive effects of stochastic rounding in numerical algorithms[END_REF][START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] and numerical simulations such as [START_REF] Paxton | Climate modeling in low precision: Effects of both deterministic and stochastic rounding[END_REF] demonstrate superior outcomes when utilizing SR as a replacement for RN in numerical computations. Connolly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] have shown that SR avoids the phenomenon of stagnation in sums, where small values are ignored by RN when they are too small relative to the sum. This phenomenon arises in various domains, such as neural networks and ODEs. In the remainder of this thesis, we investigate the probabilistic properties of SR as a rounding mode. This chapter is structured around the theoretical analysis of algorithms with multi-linear errors under SR, i.e, algorithms linear separately in each random error. By investigating this class of algorithms, we aim to understand the effects and implications of SR on their overall accuracy. The theoretical analysis demonstrates that they remain unbiased under SR. Furthermore, we prove probabilistic error bounds in O( √ nu): given a fixed probability q ∈]0; 1[, one can obtain bounds (depending on q) in O( √ nu).

We specifically focus on SR-nearness as a substitute for RN in numerical computations. As shown in Chapter 3, this stochastic rounding mode is unbiased for summation. Drawing upon the mean independence property (Lemma 2), we prove that the unbiased characteristic extends to algorithms with multi-linear errors. We also introduce and discuss two primary approaches that aim to provide bounds for the forward error of computations under SR-nearness. To demonstrate the practical applicability of these methods, we illustrate their benefits in algorithms such as sequential and pairwise summations, inner product, and the Horner algorithm. In particular, their forward errors have probabilistic bounds in O( √ nu) versus deterministic bounds in O(nu), which allows to apply them more easily in low-precision formats: we only need to suppose that nu 2 ≪ 1 to apply these results.

The first method AH is based on martingales and Azuma-Hoeffding inequality. To apply this technique, multiple approaches can be employed to construct the martingale [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF][START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF], thereby influencing the accuracy of the forward error bound [START_REF]Stochastic rounding variance and probabilistic bounds: A new approach[END_REF]. We use the same method described in [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] for the sequential summation and Horner algorithm. The analysis of the latest al-gorithm differs from that of the inner product discussed in Sub-section 2.4.2 due to the presence of an error affecting one of the multiplication operands. This extends the approach to a whole new class of algorithms, and is a step towards its application to various numerical schemes. For pairwise summation, unlike the technique proposed in [START_REF] Hallman | Precision-aware deterministic and probabilistic error bounds for floating point summation[END_REF], we introduce a novel technique to construct the martingale from the summation tree such that each level correspond to a term of the martingale sequence. We show that both probabilistic bounds of the error are in O( √ nu) with an advantage to our approach.

In Section 4.1, we demonstrate that any computation tree can accumulate a stochastic process given by ψ K = k∈K (1 + δ k ) in each of its inputs, accompanied by a natural filtration F where K ⊂ N, and certain δ k can be zero.

Additionally, various methods can be used to construct a martingale from this stochastic process. We demonstrate in Sub-section 4.1.4 that the error of any multi-linear algorithm achieves a probabilistic bound in O( √ nu), where n is the number of nodes. For a fixed probability, the AH method exhibits rapid growth for large values of n, particularly the approach proposed by Ipsen and Zhou in [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] for the inner product, which provides a probabilistic bound equivalent to O(exp(nu)) when nu ≫ 1. Additionally, the approach proposed by Connolly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] for the inner product yields a probabilistic bound equivalent to O(exp(nu)) when n log(n)u 2 ≫ 1. Furthermore, applying the martingale theory may pose challenges for certain algorithms, as demonstrated in the case of Horner's polynomial evaluation.

The second method BC is based on a bound of the error variance and Bienaymé-Chebyshev inequality. Based on the mean independence property, we have presented Lemma 11, a general framework applicable to a wide class of algorithms that allows to compute deterministic bound of the variance. This new approach is simple since it requires only information on the error variance. It demonstrates superior accuracy in low-precision formats, where SR has shown favorable outcomes across various domains such as climate modeling [START_REF] Kimpson | Climate-change modelling at reduced floating-point precision with stochastic rounding[END_REF], deep neural networks [START_REF] Gupta | Deep learning with limited numerical precision[END_REF], and PDEs [START_REF] Croci | Effects of round-to-nearest and stochastic rounding in the numerical solution of the heat equation in low precision[END_REF].

We use the Bienaymé-Chebyshev inequality to establish a probabilistic error bound. Moreover, this method improves the accuracy of the probabilistic bounds for large values of n, particularly when nu 2 ≪ 1. Theoretical analysis (Section 4.3) and numerical experiments (Section 4.4) illustrate the previous results and compare the probabilistic bounds to deterministic ones through two algorithms: the inner product and the Horner algorithm.

. AH Method

John von Neumann and Goldstine [START_REF] Neumann | Numerical inverting of matrices of high order[END_REF] have demonstrated that for m independent random variables equi-distributed in [-β -s /2; β -s /2], the maximum of the sum of these random variables is bounded by mβ -s /2, while the dispersion is given by √ mβ -s / √ 12. Also the Central Limit Theorem (CLT) [START_REF] Dacunha-Castelle | Probability and Statistics[END_REF], which states that the summation or average of independent and identically distributed random variables converges to a Gaussian (normal) distribution, ensuring favorable asymptotic results. This aligns with Wilkinson's intuition [START_REF] Wilkinson | Rounding errors in algebraic processes[END_REF],

which states that the accumulated roundoff error in n operations is usually proportional to √ nu rather than nu.

Based on the assumption of independent rounding errors, Higham and

Mary [START_REF] Higham | A new approach to probabilistic rounding error analysis[END_REF] have proposed probabilistic error bounds in O( n log(n)u) for various linear algebra computations. However, achieving the characteristic of independence is frequently unattainable. To overcome this hypothesis, with Connelly, they show in [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] that this bound always holds for SR due to mean independence (Lemma 2), a property that lies between independence and uncorrelatedness.

The mean independence of random errors is not an additional assumption but a property satisfied by SR. The idea is to construct a martingale (Definition 3) from this property and then use Azuma-Hoeffding inequality for martingales [53, p. 303] to establish probabilistic error bounds. This approach will be referred to as the AH method throughout the remainder of this dissertation.

. Sequential Summation

As discussed in Sub-section 2.3.1, the forward error of a summation of n floating point numbers is proportional to nu. We investigate this error under SR-nearness in the sequel using the method proposed by Ipsen, and Zhou in [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] for the inner product. We show O( √ nu) probabilistic bound on the forward error.

Consider s = n i=1 a i , we have

Stochastic rounding

Exact computation

ŝ1 = a 1 s 1 = a 1 ŝ2 = (ŝ 1 + a 2 )(1 + δ 2 ) s 2 = s 1 + a 2 ŝk = (ŝ k-1 + a k )(1 + δ k ) s k = s k-1 + a k ŝn = ŝ s n = s It follows that ŝ = n i=1 a i n k=max (i,2) (1 + δ k ) = n k=1 a i ψ i , (4.1)
where

ψ i = n k=max (i,2) (1 + δ k ) for all 1 ≤ i ≤ n.
Theorem 2. For all 0 < λ < 1, the computed ŝ satisfies under SR-nearness

|ŝ -s| |s| ≤ K uγ 2(n-1) (u) ln(2/λ), (4.2
)

with probability at least 1 -λ, where K = n i=1 |a i | | n i=1 a i|
is the condition number of n k=1 a i using the 1-norm and γ 2(n-1

) (u) = (1 + u) 2(n-1) -1. Proof. Denote F k = {δ 1 , , . . . , δ k } and Z k = ŝk -s k for all 1 ≤ k ≤ n, with Z 1 = 0 and Z n = ŝ -s. Let us show that Z 1 , . . . , Z n form a martingale with respect to F n-1 . It is straightforward that Z k is a function of δ 1 , . . . , δ k . E(|Z k |) is finite because a k are finite and |δ k | ≤ u for all 1 ≤ k ≤ n. Regarding the third assumption in definition 3 (page 26) we have 1 ≤ k ≤ n Z k = ŝk -s k = (ŝ k-1 + a k )(1 + δ k ) -s k-1 -a k = ŝk-1 -s k-1 + (ŝ k-1 + a k )δ k = Z k-1 + (ŝ k-1 + a k )δ k . It follows that E[Z k /F k-1 ] = E[(Z k-1 + (ŝ k-1 + a k )δ k )/F k-1 ] = Z k-1 + (ŝ k-1 + a k )E[δ k /F k-1 ] = Z k-1 by Lemma 2.
Therefore, the sequence Z 1 , . . . , Z n form a martingale with respect to F n-1 .

Moreover,

|Z k -Z k-1 | = |(ŝ k-1 + a k )δ k | ≤ u (|a k | + |(ŝ k-2 + a k-1 )(1 + δ k-1 )|) ≤ u (|a k | + (1 + u)(|ŝ k-2 | + |a k-1 |)) .

By induction we have

|Z k -Z k-1 | ≤ u |a 1 | (1 + u) k-2 + k i=2 |a i | (1 + u) k-i = uC k , where C k = |a 1 | (1 + u) k-2 + k i=2 |a i | (1 + u) k-i . Since Z n = ŝ -s and Z 1 = 0,
Azuma-Hoeffding inequality (Lemma 3) yields

P   |ŝ -s| ≤ n k=2 u 2 C 2 k 2 ln(2/λ)   ≥ 1 -λ,
where 0 < λ < 1. We have

C 2 k ≤ (1 + u) 2(k-2) k i=1 |a i | 2 ≤ (1 + u) 2(k-2) n i=1 |a i | 2 ,
and

n k=2 C 2 k ≤ n i=1 |a i | 2 n k=2 (1 + u) 2(k-2) = n i=1 |a i | 2 (1 + u) 2(n-1) -1 (1 + u) 2 -1 = n i=1 |a i | 2 γ 2(n-1) (u) u 2 + 2u . Since, u u+2 ≤ u 2 n k=2 u 2 C 2 k ≤ n i=1 |a i | 2 uγ 2(n-1) (u) 2 .
Il follows that

|ŝ -s| ≤ n i=1 |a i | uγ 2(n-1) (u) 2 2 ln(2/λ).
with probability at least 1 -λ. Finally

|ŝ -s| |s| ≤ K uγ 2(n-1) (u) ln(2/λ),
with probability at least 1 -λ.

This example demonstrates the direct applicability of the AH method to algorithms with exact inputs. Furthermore, this method is also valid for the inner product [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF][START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF], where multiplication operands are exact (in the sense of multiplying two inputs), and the errors are accumulated through summations. In the following, we analyze the Horner algorithm and show that the AH method can still be used, even if the multiplication operand is affected by an error.

. Horner Algorithm

Horner algorithm is an efficient way of evaluating polynomials. When performed in floating-point arithmetic, this algorithm may suffer from catastrophic cancellations and yield a computed value less accurate than expected.

Let P (x) = n i=0 a i x i , Horner rule consists in writing this polynomial as

P (x) = (((a n x + a n-1 )x + a n-2 )x • • • + a 1 )x + a 0 .
We define by induction the following sequence Operation Stochastic rounding Exact computation

r0 = a n r 0 = a n * r2k-1 = r2k-2 x(1 + δ 2k-1 ) r 2k-1 = r 2k-2 x + r2k = (r 2k-1 + a n-k )(1 + δ 2k ) r 2k = r 2k-1 + a n-k Output r2n = P (x) r 2n = P (x)
for all 1 ≤ k ≤ n, with δ 2k-1 and δ 2k the rounding errors from the products and the additions respectively. Let δ 0 = 0, we have

r2n = n i=0 a i x i 2n k=2(n-i) (1 + δ k ). (4.3)
Let us denote Z i := ri -r i for all 0 ≤ i ≤ 2n. The total forward error is

|Z 2n | = |r 2n -r 2n | = | P (x) -P (x)| and | P (x) -P (x)| = n i=0 a i x i   2n k=2(n-i) (1 + δ k ) -1   ≤ n i=0 |a i x i |γ 2n (u).
Finally,

| P (x) -P (x)| |P (x)| ≤ Kγ 2n (u), (4.4)
where

K = n i=0 |aix i | |P (x)|
is the condition number of the polynomial evaluation.

The deterministic bound is proportional to nu. In the following, we prove a probabilistic bound in O( √ nu). The partial sum forward errors satisfy

Z 2k-1 = r2k-1 -r 2k-1 = r2k-2 x(1 + δ 2k-1 ) -r 2k-2 x = xZ 2k-2 + r2k-2 xδ 2k-1 , Z 2k = r2k -r 2k = (r 2k-1 + a n-k )(1 + δ 2k ) -r 2k-1 -a n-k = Z 2k-1 + (r 2k-1 + a n-k )δ 2k ,
for all 1 ≤ k ≤ n. The sequence Z 0 , . . . , Z 2n does not form a martingale with respect to δ 1 , . . . , δ 2n due to the multiplication in odd steps,

E[Z 2k-1 /δ 1 , . . . , δ 2k-2 ] = xZ 2k-2 .
In order to form a martingale and use the Azuma-Hoeffding inequality, we define the following variable change

Y i = Z i x (i+1)/2 ,
where (i + 1)/2 is the greatest integer less than or equal to (i + 1)/2, we thus have

Y 2k-1 = Y 2k-2 + 1 x k-1 r2k-2 δ 2k-1 , Y 2k = Y 2k-1 + 1 x k (r 2k-1 + a n-k )δ 2k , (4.5)
for all 1 ≤ k ≤ n with Y 0 = 0.

Theorem 3. The sequence of random variables Y 0 , • • • , Y 2n is a martingale with respect to δ 1 , • • • , δ 2n .
Proof. We check that the three conditions of Definition 3 are satisfied. Throughout the proof, we note the set

F k = {δ 1 , • • • , δ k }. • The recursion of Z k shows that Y i is a function of δ 1 , • • • , δ i for all 1 ≤ i ≤ 2n. • E(|Y i |) is finite because x and a k are finite for all n -i ≤ k ≤ n and |δ j | ≤ u for all 1 ≤ j ≤ i. • We prove that E[Y i /F i-1 ] = Y i-1
by distinguishing the even and odd cases. Firstly, using the mean independence of δ 1 , • • • , δ 2k-1 and Equation (4.5) we obtain

E[Y 2k-1 /F 2k-2 ] = E[Y 2k-2 /F 2k-2 ] + E[ 1 x k-1 r2k-2 δ 2k-1 /F 2k-2 ] = Y 2k-2 + 1 x k-1 r2k-2 E[δ 2k-1 /F 2k-2 ] = Y 2k-2 .
Secondly, using the mean independence of δ 1 , • • • , δ 2k and Equation (4.5)

we obtain

E[Y 2k /F 2k-1 ] = E[Y 2k-1 /F 2k-1 ] + E[ 1 x k (r 2k-1 + a n-k )δ 2k /F 2k-1 ] = Y 2k-1 + 1 x k (r 2k-1 + a n-k )E[δ 2k /F 2k-1 ] = Y 2k-1 .
The martingale does not manifest explicitly, we need a change of variable to exhibit it. To apply Azuma-Hoeffding inequality, we need to bound the martingale steps. The following lemma presents bounds on the martingale

increments |Y i -Y i-1 | for all 1 ≤ i ≤ 2n. Lemma 5. The above martingale Y 0 , • • • , Y 2n satisfies |Y i -Y i-1 | ≤ C i u, for all 1 ≤ i ≤ 2n, where C 2k-1 = |a n |(1 + u) 2k-2 + k-1 j=1 |a n-j ||x| -j (1 + u) 2(k-j)-1 , C 2k = |a n |(1 + u) 2k-1 + k j=1 |a n-j ||x| -j (1 + u) 2(k-j) , for all 1 ≤ k ≤ n. Proof. Note that Y 0 = 0, then |Y 1 -Y 0 | = |Y 1 | ≤ |a n |u
and the equality holds for C 1 . Using Equation (4.5)

|Y 2k-1 -Y 2k-2 | ≤ 1 |x| k-1 |r 2k-2 |u.
Moreover,

|r 2k-2 | ≤ |r 2k-3 |(1 + u) + |a n-k+1 |(1 + u) ≤ |r 2k-4 ||x|(1 + u) 2 + |a n-k+1 |(1 + u).

By induction we obtain

|r 2k-2 | ≤ |a n ||x| k-1 (1 + u) 2k-2 + k-1 j=1 |a n-j ||x| k-j-1 (1 + u) 2(k-j)-1 . (4.6)
This completes the proof for C 2k-1 for all 1 ≤ k ≤ n. For C 2k for all 1 ≤ k ≤ n, using Equation (4.5)

|Y 2k -Y 2k-1 | ≤ 1 |x| k |r 2k-1 + a n-k |u.
Moreover,

|r 2k-1 + a n-k | ≤ |r 2k-1 | + |a n-k | ≤ |r 2k-2 | |x| (1 + u) + |a n-k | .
The inequality (4.6) implies

|x| (1 + u) |r 2k-2 | ≤   |a n ||x| k-1 (1 + u) 2k-2 + k-1 j=1 |a n-j ||x| k-j-1 (1 + u) 2(k-j)-1   |x| (1 + u) = |a n ||x| k (1 + u) 2k-1 + k-1 j=1 |a n-j ||x| k-j (1 + u) 2(k-j) .
Finally,

|Y 2k -Y 2k-1 | ≤ u |x| k |r 2k-1 + a n-k | ≤ u |x| k   |a n ||x| k (1 + u) 2k-1 + k-1 j=1 |a n-j ||x| k-j (1 + u) 2(k-j) + |a n-k |   = u   |a n |(1 + u) 2k-1 + k j=1 |a n-j ||x| -j (1 + u) 2(k-j)   .
This completes the proof for C 2k for all 1 ≤ k ≤ n.

We now have all the tools to state and demonstrate the main result of this sub-section: Theorem 4. For all 0 < λ < 1, the computed P (x) satisfies under SR-nearness

| P (x) -P (x)| |P (x)| ≤ K uγ 4n (u) ln(2/λ), (4.7
)

with probability at least 1 -λ, where K = n i=0 |a i x i | |P (x)|
is the condition number of the polynomial evaluation.

Proof. Recall that |r 2n -r 2n | = |Z 2n | = |x n ||Y 2n |.
Therefore, Y 0 , . . . , Y 2n is a martingale with respect to δ 1 , ..., δ 2n and Lemma 5 implies

|Y i -Y i-1 | ≤ C i u for all 1 ≤ i ≤ 2n.
Using the Azuma-Hoeffding inequality yields

P   |Y 2n | ≤ u 2n i=1 C 2 i 2 ln(2/λ)   ≥ 1 -λ, it follows that |Z 2n | ≤ u 2n i=1 (|x| n C i ) 2 2 ln(2/λ),
with probability at least 1 -λ, where

|x| n C 2k = |a n ||x| n (1 + u) 2k-1 + k j=1 |a n-j x n-j |(1 + u) 2(k-j) ≤ (1 + u) 2k-1 k j=0 |a n-j x n-j | ≤ (1 + u) 2k-1 n j=0 |a j x j |, for all 1 ≤ k ≤ n. Hence, (|x| n C 2k ) 2 ≤ (1 + u) 2(2k-1)   n j=0 |a j x j |   2 .
In a similar way,

(|x| n C 2k-1 ) 2 ≤ (1 + u) 2(2k-2)   n j=0 |a j x j |   2 . Thus, 2n i=1 (|x| n C i ) 2 ≤   n j=0 |a j x j |   2 2n-1 i=0 ((1 + u) 2 ) i =   n j=0 |a j x j |   2 ((1 + u) 2 ) 2n -1 (1 + u) 2 -1 =   n j=0 |a j x j |   2 γ 4n (u) u 2 + 2u .
As a result,

| P (x) -P (x)| = |Z 2n | ≤ n j=0 |a j x j | uγ 4n (u) 2 + u 2 ln(2/λ),
with probability at least 1 -λ. Finally,

| P (x) -P (x)| |P (x)| ≤ K uγ 4n (u) ln(2/λ),
with probability at least 1 -λ.

We have extended the method proposed in [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] to derive a new probabilistic bound on the forward error of the Horner algorithm. This illustrates how the AH method can be applied (with some work) to any algorithm based on a fixed sequence of sum and products. We now turn to the pairwise summation and show the benefits of AH method through this algorithm.

. Pairwise Summation

It is known [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF] that the accumulator implementation of a sum of n numbers s = n i=1 x i using a binary tree leads to a deterministic error bound in O(log 2 (n)u). Hallman and Ipsen in [START_REF] Hallman | Precision-aware deterministic and probabilistic error bounds for floating point summation[END_REF] have studied this problem under stochastic rounding. Using the AH method (with an alternative form of Azuma-Hoeffding inequality), they have shown probabilistic bound in O( log 2 (n)u).

The key idea of their approach involved constructing a martingale through a systematic recurrence on the computational tree, incorporating two probability parameters.

In the following, we investigate the forward error made by the pairwise summation under SR-nearness using the AH method. This approach's feature lies in its flexibility for building the martingale. We construct a martingale straight from the tree levels for the pairwise summation and then use Azuma-Hoeffding inequality. This technique has the advantage of building a martingale from the entire tree. We compare our probabilistic bound to the bound proposed by Hallman and Ipsen in [START_REF] Hallman | Precision-aware deterministic and probabilistic error bounds for floating point summation[END_REF].

Let us assume the tree structure as follows:

S h 1 S h-1 1 . . . S 1 1 x 1 x 2 S 1 2 . . . S h-1 2 . . . . . . S 1 2 h-2 S 1 2 h-1 x 2 h-1 x 2 h
Considering h the height of the summation tree, if 2 h-1 < n < 2 h , we set the absent 2 h -n inputs to zero. Without loss of generality, let us then assume that n = 2 h . Denote S 0 i = x i and

S k i = S k-1 2i-1 + S k-1 2i for all 1 ≤ i ≤ 2 h-k and 1 ≤ k ≤ h. We have S k l = l2 k i=(l-1)2 k +1 x i and S h 1 = 2 h i=1 x i = s. Let Ŝ0 i = S 0 i and Ŝk i = ( Ŝk-1 2i-1 + Ŝk-1 2i )(1+δ k i ) for all 1 ≤ i ≤ 2 h-k and 1 ≤ k ≤ h.
We have

Ŝk l = l2 k i=(l-1)2 k +1 x i k j=1 (1 + δ j ⌈ i 2 j ⌉ ).
In particular

Ŝh 1 = 2 h i=1 x i h j=1 (1 + δ j ⌈ i 2 j ⌉ ) = 2 h i=1 x i ψ i , (4.8)
where

ψ i = h j=1 (1 + δ j ⌈ i 2 j ⌉ ) for all 1 ≤ i ≤ 2 h . Theorem 5. For all k ≥ 1, let F k = {δ j i , 1 ≤ i ≤ 2 h-j , 1 ≤ j ≤ k}. Denote M 0 = 0 and for k > 0, M k = 2 h-k i=1 Ŝk i -S k i .
Therefore, M h = Ŝh 1 -S h 1 and M 0 , . . . , M h form a martingale with respect to

F h-1 = {δ j i , 1 ≤ i ≤ 2 h-j , 1 ≤ j ≤ h -1}.
Proof. The recursion of Ŝk i shows that M k is a function of F k-1 for all 1 ≤ k ≤ h. Moreover, E(|M k |) is finite because x i are finite for all 1 ≤ i ≤ 2 h and the relative errors |δ i j | ≤ u. Let us prove the third point:

M k = 2 h-k i=1 Ŝk i -S k i = 2 h-k i=1 ( Ŝk-1 2i-1 + Ŝk-1 2i )(1 + δ k i ) -(S k-1 2i-1 + S k-1 2i ) = M k-1 + 2 h-k i=1 ( Ŝk-1 2i-1 + Ŝk-1 2i )δ k i .
By definition of F k , we have for all

1 ≤ k ≤ h E[M k /F k-1 ] = E   M k-1 + 2 h-k i=1 ( Ŝk-1 2i-1 + Ŝk-1 2i )δ k i /F k-1   = E[M k-1 /F k-1 ] + E   2 h-k i=1 ( Ŝk-1 2i-1 + Ŝk-1 2i )δ k i /F k-1   by linearity = M k-1 + 2 h-k i=1 ( Ŝk-1 2i-1 + Ŝk-1 2i )E[δ k i /F k-1 ] = M k-1
by mean independence.

Therefore, M 0 , . . . , M h form a martingale with respect to F h-1 .

We now have all the tools to state and demonstrate the main result of this sub-section: Theorem 6. For all 0 < λ < 1, the computed Ŝh

1 satisfies under SR-nearness Ŝh 1 -S h 1 S h 1 ≤ κ uγ 2⌈log 2 (n)⌉ (u) ln(2/λ), (4.9
)

with probability at least 1 -λ, where κ = n i=1 |x i | | n i=1 x i|
is the condition number of the summation of the x i .

Proof. We need firstly to bound the martingale steps. Equation (4.8) yields

|M k -M k-1 | ≤ 2 h-k i=1 ( Ŝk-1 2i-1 + Ŝk-1 2i )δ k i ≤ u 2 h-k i=1 Ŝk-1 2i-1 + Ŝk-1 2i ≤ u(1 + u) k-1 2 h-k i=1 2 k-1 (2i-1) m=2 k-1 (2i-2)+1 x m + 2 k-1 (2i) m=2 k-1 (2i-1)+1 x m ≤ u(1 + u) k-1 2 h-k i=1 2 k i m=2 k (i-1)+1 |x m | = u(1 + u) k-1 2 h i=1 |x m | = u(1 + u) k-1 ∥x∥ 1 . Denote C k = u(1 + u) k-1 ∥x∥ 1 , Azuma-Hoeffding inequality implies |M h | ≤ h k=1 C 2 k 2 ln(2/λ), with probability at least 1 -λ. Now h k=1 C 2 k = u 2 ∥x∥ 2 1 h k=1 (1 + u) 2(k-1) = u 2 ∥x∥ 2 1 (1 + u) 2h -1 (1 + u) 2 -1 = u ∥x∥ 2 1 γ 2h (u) u + 2 .
Since, u u+2 ≤ u 2 and h = ⌈log 2 (n)⌉, we have

|M h | ≤ ∥x∥ 1 u γ 2⌈log 2 (n)⌉ (u) 2 2 ln(2/λ),
with probability at least 1 -λ. Finally

Ŝh 1 -S h 1 S h 1 ≤ κ uγ 2⌈log 2 (n)⌉ (u) ln(2/λ),
with probability at least 1 -λ.

Comparison with Hallman and Ipsen pairwise bound [START_REF] Hallman | Precision-aware deterministic and probabilistic error bounds for floating point summation[END_REF]. The probabilistic bound proposed in [33, cor, 2.10] to the pairwise summation forward error is

Ŝh 1 -S h 1 S h 1 ≤ κu √ h 2 ln(2/δ)(1 + ϕ n,h,η ), (4.10)
with probability at least 1 -(η + δ), where h is the height of the computational tree, κ is the same condition number, and ϕ n,h,η ≡ λ n,η √ 2hu exp λ 2 n,η hu 2 with λ n,η ≡ 2 ln(2n/η). The figure 4.1 compares the probabilistic bounds (4.9) and (4.10) with probability 0.9, using identical values for the two probabilistic parameters of the bound (4.10). Although both are in O( log 2 (n)u), the figure clearly shows that in this case, the bound obtained from the AH method consistently yields lower values for all n and u. Note that we observe numerically convergence of the bound (4.10) to the bound (4.9) when η approaches 0 and δ approaches λ. Nevertheless, the AH method remains simple and intuitive.

. Generalization

This section illustrates that for a computation tree with additions/subtractions, the error has a probabilistic bound in O( √ h) under SR-nearness, where h = ⌈log 2 (n)⌉ is the tree height and n is the number of nodes. In the case of multiplication of n numbers, this result remains valid but with a probabilistic error bound in O( √ n), where n is the number of nodes. This result can be generalized to any complete computation tree with a sequence of elementary operations {+, -, * } and multi-linear errors. Let A, A x and A y algorithms based on elementary operations {+, -, * }. Suppose that last operation in A is: z ← x op y x y • y corresponds to the outcome of the algorithm A y (the right sub-tree computation).

• op ∈ {+, -, * }

Note that the errors of A, A x and A y are multi-linear, i.e, in the multiplication case, A x and A y do not share any instructions. In other words, the degree of each error is one throughout the entire computation.

Our goal is to build by induction a martingale from the computation of A by taking into account the last operation in use. Suppose that the last operation in A is an addition, i.e, z ← x + y. Consider the relative errors Φ, X and Ψ associated respectively to x, y, and z. Note x, y, and z their respective theoretical values, and x, ŷ and ẑ their computed values. Therefore:

     x = x(1 + Φ), ŷ = y(1 + X), ẑ = z(1 + Ψ).
We have z = x + y, then, there exists δ such that ẑ = (x + ŷ)(1 + δ). Hence,

Ψ = ẑ -z z = x + ŷ x + y (1 + δ) -1 = 1 + x x + y Φ + y x + y X (1 + δ) -1.
Since we assume the inputs are exact, the relative error associated with an input is zero. The(0) 0 i=0 forms a martingale. Suppose by induction that there exist constants K x ≥ 1 and K y ≥ 1, and martingales

(Φ i ) k-1 i=0 and (X i ) l-1 i=0 with their ith step |Φ i -Φ i-1 | and |X i -X i-1 | bounded respectively by K x u(1 + u) i-1 and K y u(1 + u) i-1 , such that Φ 0 = 0, X 0 = 0, Φ = Φ k-1 and X = X l-1 .
Note that Φ 0 = 0, X 0 = 0 are martingales. Lemma 6. Let m = max{k, l} + 1. The stochastic process (Ψ i ) m-1 i=0 such that Ψ m-1 = Ψ, and for all 0 ≤ i < m -1,

Ψ i = x x + y Φ min{i,k} + y x + y X min{i,l} ,
forms a martingale.

Proof. Without loss of generality, let us assume that k ≤ l. Then, m = l + 1

and

Ψ i = x x+y Φ i + y x+y X i for all 0 ≤ i ≤ k -2 Ψ i = x x+y Φ + y x+y X i for all k -1 ≤ i ≤ l -1.
Note that

(Φ i ) m-2 i=0 with Φ i = Φ for all k -1 ≤ i ≤ m -2 and (X i ) m-2
i=0 are martingales by induction hypothesis. Since the martingale set is a vector space, as a linear combination of them, (Ψ i ) m-2 i=0 is a martingale. Moreover, by mean independence of δ from Φ and X we have

E[Ψ m-1 /Ψ m-2 ] = E 1 + x x + y Φ + y x + y X (1 + δ) -1/Ψ m-2 = 1 + x x + y Φ + y x + y X E[(1 + δ)/Ψ m-2 ] -1 = Ψ m-2 .
Thus, (Ψ) m-1 i=0 is a martingale and Ψ m-1 = Ψ.

In this lemma, we have shown a martingale by induction when the last operation of the algorithm A is an addition. In order to use Azuma-Hoeffding inequality (Lemma 3), we have to bound the martingale increments.

Lemma 7. Let K z = |x| |x+y| K x + |y| |x+y| K y . The martingale (Ψ i ) m-1 i=0 satisfies |Ψ i -Ψ i-1 | ≤ uC i ,
where

C i = K z (1 + u) i-1 for all 1 ≤ i ≤ m -1.
Proof. For all 0 ≤ i < m -1, by induction hypothesis we have

|Ψ i -Ψ i-1 | = x x + y (Φ min{i,k-1} -Φ min{i-1,k-1} ) + y x + y (X min{i,l-1} -X min{i-1,l-1} ) ≤ x x + y Φ min{i,k-1} -Φ min{i-1,k-1} + y x + y X min{i,l-1} -X min{i-1,l-1} ≤ |x| |x + y| K x u(1 + u) min{i-1,k-1} + |y| |x + y| K y u(1 + u) min{i-1,l-1} ≤ K z u(1 + u) i-1 .

56

Moreover,

|Ψ m-1 -Ψ m-2 | = 1 + x x + y Φ + y x + y X δ ≤ u x x + y (Φ + 1) + y x + y (X + 1) .
Since Φ 0 = 0,

|Φ + 1| = |Φ k-1 + 1| = 1 + k-1 j=1 (Φ j -Φ j-1 ) ≤ 1 + k-1 j=1 |Φ j -Φ j-1 | ≤ 1 + k-1 j=1 uK x (1 + u) j-1 = 1 + uK x (1 + u) k-1 -1 u = 1 + K x (1 + u) k-1 -K x ≤ K x (1 + u) k-1 .

The same method shows that |X

+ 1| ≤ K y (1 + u) l-1 . It follows that |Ψ m-1 -Ψ m-2 | ≤ u |x| |x + y| K x (1 + u) k-1 + |y| |x + y| K y (1 + u) l-1 ≤ K z u(1 + u) m-2 .
Corollary 1. For all 0 < λ < 1, the computed ẑ satisfies under SR-nearness

|ẑ -z| |z| ≤ K z uγ 2(m-1) (u) ln(2/λ), (4.11)
with probability at least 1 -λ.

Proof. Using Azuma-Hoeffding inequality we have

|ẑ -z| |z| = |Ψ m-1 | ≤ m-1 i=1 u 2 C 2 i 2 ln(2/λ), with probability at least 1 -λ. Moreover, m-1 i=1 u 2 C 2 i = u 2 K 2 z m-1 i=1 (1 + u) 2(i-1) = u 2 K 2 z (1 + u) 2(m-1) -1 u 2 + 2u ≤ uK 2 z γ 2(m-1) (u) 2 .
Finally,

|ẑ -z| |z| ≤ K z uγ 2(m-1) (u) ln(2/λ),
with probability at least 1 -λ.

Suppose now that the last operation in A is a multiplication, i.e, z ← x × y. Consider the relative errors Φ, X, and Ψ associated respectively to x, y, and z. Note x, y, and z their respective theoretical values, and x, ŷ and ẑ their computed values, with x = x(1 + Φ), ŷ = y(1 + X), and ẑ = z(1 + Ψ).

We have z = x × y, then there exists δ such that ẑ = x × ŷ(1 + δ). Hence,

Ψ = ẑ -z z = x × ŷ x × y (1 + δ) -1 = (1 + Φ)(1 + X)(1 + δ) -1.
We know by induction that Φ and X are the last terms of two martingales.

At the same time, the multiplication of two martingales is not necessary a martingale. Consequently, in contrast to the addition case, we have to establish an order of operations in the construction of the martingale Ψ. Note that they are all equivalent and leads asymptotically to obtain the same result. In figure 4.2, we assume that the left sub-tree is computed before the right sub tree, which means that in the computation of x, we assume that we don't have any operation on y. This choice, leads to built two martingales

(Φ i ) k-1 i=0
and (X i ) m-2 i=0 such that Φ 0 = 0, Φ = Φ k-1 , X j = 0 for all 0 ≤ j ≤ k -1, X = X m-2 , and random errors in Φ are different from those of X (thanks to the multi-linearity of errors in the computation of z). The following lemma shows that Ψ is the last term of a martingale built from (Φ i ) k-1 i=0 and (X i ) m-2 i=0 .

Lemma 8. The stochastic process (Ψ i ) m-1 i=0 such that

Ψ i =      Φ i = (1 + Φ i )(1 + 0) -1 for all 0 ≤ i ≤ k -1 (1 + Φ)(1 + X i ) -1 for all k ≤ i ≤ m -2 (1 + Φ)(1 + X)(1 + δ) -1 for i = m -1,
forms a martingale.

Proof. For all 0 ≤ i ≤ k -1, by construction of Ψ, we have

Ψ i = Φ i . Since (Φ i ) k-1
i=0 is a martingale, we have

E[Ψ i /Ψ i-1 ] = E[Φ i /Ψ i-1 ] = Φ i-1 = Ψ i-1 .

58

Moreover,

E[Ψ k /Ψ k-1 ] = E[(1 + Φ)(1 + X k ) -1/Ψ k-1 ] = (1 + Φ)E[(1 + X k )/Ψ k-1 ] -1 = (1 + Φ) -1 by Lemma 2 = Φ. Since (X i ) m-2 i=0 is a martingale, for all k < i ≤ m -2, E[Ψ i /Ψ i-1 ] = E[(1 + Φ)(1 + X i ) -1/Ψ i-1 ] = (1 + Φ)E[(1 + X i )/Ψ i-1 ] -1 = (1 + Φ)(1 + X i-1 ) -1 = Ψ i-1 .
By the mean independence of δ and Ψ m-2 ,

E[Ψ m-1 /Ψ m-2 ] = E[(1 + Φ)(1 + X m-2 )(1 + δ) -1/Ψ m-2 ] = (1 + Φ)(1 + X m-2 )E[(1 + δ)/Ψ m-2 ] -1 = Ψ m-2 .
In order to use Azuma-Hoeffding inequality (Lemma 3), we have to bound the martingale increments. We can show by induction that there exist constants K x ≥ 1 and K y ≥ 1, such that the ith step |Φ i -Φ i-1 | and |X i -X i-1 | are bounded respectively by K x u(1+u) i-1 and K y u(1+u) i-k (because X j = 0 for all 0 ≤ j ≤ k -1).

Lemma 9. Let

K z = K x K y . The martingale (Ψ i ) m-1 i=0 satisfies |Ψ i -Ψ i-1 | ≤ uC i , where C i = K z (1 + u) i-1 for all 1 ≤ i ≤ m -1. Proof. For all 1 ≤ i ≤ k -1, |Ψ i -Ψ i-1 | = |Φ i -Φ i-1 | ≤ K x u(1 + u) i-1 . Moreover, for all k ≤ i ≤ m -2, |Ψ i -Ψ i-1 | = |(1 + Φ)(1 + X i ) -(1 + Φ)(1 + X i-1 )| = |(1 + Φ)(X i -X i-1 )| ≤ |1 + Φ i | K y u(1 + u) i-k .
As for the summation case,

|1 + Φ i | ≤ K x (1+u) k-1 . Then, for all k ≤ i ≤ m-2, |Ψ i -Ψ i-1 | ≤ K x (1 + u) k-1 K y u(1 + u) i-k = K z u(1 + u) i-1 . Finally, |Ψ m-1 -Ψ m-2 | = |(1 + Φ)(1 + X)δ| ≤ uK x (1 + u) k-1 K y (1 + u) m-k-1 ≤ uK z (1 + u) m-2 .
. Corollary 2. For all 0 < λ < 1, the computed ẑ satisfies under SR-nearness

|ẑ -z| |z| ≤ K z uγ 2(m-1) (u) ln(2/λ), (4.12)
with probability at least 1 -λ.

Proof. Using Azuma-Hoeffding inequality we have

|ẑ -z| |z| = |Ψ m-1 | ≤ m-1 i=1 u 2 C 2 i 2 ln(2/λ), with probability at least 1 -λ. Moreover, m-1 i=1 u 2 C 2 i = u 2 K 2 z m-1 i=1 (1 + u) 2(i-1) = u 2 K 2 z (1 + u) 2(m-1) -1 u 2 + 2u ≤ uK 2 z γ 2(m-1) (u) 2 .
Finally,

|ẑ -z| |z| ≤ K z uγ 2(m-1) (u) ln(2/λ),
with probability at least 1 -λ.

Interestingly, Corollaries 1 and 2 show that the error of any algorithm A with elementary operations {+, -, * } and multi-linear errors has a probabilistic bound in O( √ nu), where n is the number of operations.

Discussion: The Azuma-Hoeffding inequality is commonly used in mathematical analysis to establish an upper bound on the deviation of the sum of a sequence of independent and bounded random variables, martingales in the previous studies. This inequality is particularly useful in the AH method, which provides more precise bounds for higher probabilities. However, Azuma-Hoeffding inequality assumes that the martingale steps are bounded by their worst-case scenario, which may lead to less accuracy in certain cases.

In numerical analysis, where problem sizes are often large in domains like numerical simulation and data analysis, ensuring high-quality results can be challenging under stochastic rounding, especially for tight probabilistic bounds of the error with a fixed probability and for a large n. To overcome this problem, in the next section, we present a new method based on a bound of the error variance and Bienaymé-Chebyshev inequality. The latter provides a general upper bound on the probability that a random variable deviates from its mean by a certain amount.

. BC Method

Built upon the mean independence property, we propose a new approach to characterize SR errors based on the variance computation. We pinpoint common error patterns in a set of numerical algorithms and introduce Lemma 11, which presents a general framework applicable to a wide range of algorithms. This framework enables the computation of a deterministic upper bound on the variance of the accumulated error in a computation under SR-nearness. This new method employs the Bienaymé-Chebyshev inequality to establish a probabilistic error bound.

Interestingly, this approach also gives probabilistic error bounds in O( √ nu)

and leads to better asymptotic behavior than the AH method in several situations. Our method has the advantage of providing a tighter probabilistic bound for all algorithms fitting our model. We show that this method is always better with probability at most 0.7678. Otherwise, for a fixed probability and with respect to the unit roundoff u, the bound proposed by this method remains tight from a certain problem size n. This approach will be referred to as the BC method throughout the remainder of this dissertation. Now, let us recall the Bienaymé-Chebyshev inequality [8, p. 19].

Lemma 10. (Bienaymé-Chebyshev inequality) Let X be a random variable with finite expected value E(X) and finite non-zero variance V (X). For any real number α > 0,

P |X -E(X)| ≤ α V (X) ≥ 1 - 1 α 2 .
We now turn to bound the variance of the error in computation. If x = x(1 + δ) is the result of an elementary operation rounded with SR-nearness, then E(x) = x and

V (x) = E(x 2 ) -x 2 = ⌈x⌉ 2 p(x) + ⌊x⌋ 2 (1 -p(x)) -x 2 = p(x)(⌈x⌉ 2 -⌊x⌋ 2 ) -(x 2 -⌊x⌋ 2 ) = p(x)ϵ(x)(⌈x⌉ + ⌊x⌋) -p(x)ϵ(x)(x + ⌊x⌋) = p(x)ϵ(x)(⌈x⌉ -x) = ϵ(x) 2 p(x)(1 -p(x)). 61 Using (2.4) leads to V (x) ≤ x 2 u 2
4 , in particular V (x) ≤ x 2 u 2 . Lemma 11 below allows to estimate the variance of the accumulated errors in a sequence of additions and multiplications. This accumulation manifests as a product of errors, a phenomenon that naturally arises when defining the relative error in relation to the standard model (2.3).

Let K a subset of N of cardinality n. Assume that δ 1 , δ 2 , ... in that order are random errors on elementary operations obtained from SR-nearness. Let us denote

ψ K = k∈K (1 + δ k ). Since |δ k | ≤ u for all k ∈ K we have |ψ K | ≤ (1 + u) n . Let K△K ′ = (K ∪ K ′ ) \ (K ∩ K ′ ).
The following lemma gives some properties of ψ that allow to bound the variance of errors in an algorithm consisting of a fixed sequence of sums and products.

Lemma 11. Under SR-nearness ψ

K satisfies 1. E(ψ K ) = 1. 2. Let K ′ ⊂ N such that |K∩K ′ | = m, under the assumption that ∀ j ∈ K△ K ′ , k ∈ K ∩ K ′ , with j < k we have 0 ≤ Cov(ψ K , ψ K ′ ) ≤ γ m (u 2 ). 3. V (ψ K ) ≤ γ n (u 2 ),
where

γ n (u 2 ) = (1 + u 2 ) n -1 ≈ exp (nu 2 ) -1 = nu 2 + O(u 4 ).
Remark 3. To compute a bound on Cov(ψ K , ψ K ′ ), we assume that common errors in inputs tend to appear toward the end of the computation. This happens when errors accumulate while evaluating algorithms based on elementary operations. However, for the bound on the variance, we have K = K ′ ; therefore, this assumption remains irrelevant.

Proof. The first point is an immediate consequence of [12, lem 6.1]. The third point is a particular case of the second with K = K ′ . Let us prove point 2.

Cov(ψ K , ψ K ′ ) = E(ψ K ψ K ′ ) -E(ψ K )E(ψ K ′ ) = E(ψ K ψ K ′ ) -1. Assume that K ∩ K ′ = {k 1 , . . . , k m }. Let us denote Q m := ψ K ψ K ′ = j∈K△K ′ (1 + δ j ) km l=k 1 (1 + δ l ) 2 ,
such that j < k i for all j ∈ K△K ′ and i ∈ {1, . . . , m}.

We prove by induction over

m that 1 ≤ E(Q m ) ≤ (1 + u 2 ) m . For m = 0, we have K ∩ K ′ = ∅ and Q 0 = j∈K△K ′ (1 + δ j ), from the first point E(Q 0 ) = 1.
Assume that the inequality holds for Q m-1 .

Q m = (1 + δ km ) 2 k m-1 l=k 1 (1 + δ l ) 2 j∈K△K ′ (1 + δ j ) = (1 + δ km ) 2 Q m-1 .
Let us denote S K△K ′ = {δ j , j ∈ K△K ′ }, using the law of total expectation E(X) = E(E[X/Y ]) and lemma 2 we have

E(Q m ) = E (1 + δ km ) 2 Q m-1 = E E[(1 + δ km ) 2 Q m-1 /S K△K ′ , δ k 1 , . . . , δ k m-1 ] = E Q m-1 E[(1 + δ km ) 2 /S K△K ′ , δ k 1 , . . . , δ k m-1 ] = E Q m-1 E[1 + δ 2 km /S K△K ′ , δ k 1 , . . . , δ k m-1 ] .
Since |δ km | ≤ u, we have by induction

1 ≤ E(Q m-1 ) ≤ E Q m-1 E[1 + δ 2 km /S K△K ′ , δ k 1 , . . . , δ k m-1 ] ≤ E Q m-1 (1 + u 2 ) . Thus, 1 ≤ E (Q m ) ≤ (1 + u 2 ) m . Finally, by induction, the claim is proven 0 ≤ E (Q m ) -1 = Cov(ψ K , ψ K ′ ) ≤ γ m (u 2 ).
Under SR-nearness, Lemma 11 can now be used to derive a variance bound for many algorithms, such as summation, inner products, matrix-vector and matrix-matrix products, solutions of triangular systems, and the Horner algorithm.

. Sequential Summation

Recall that from Sub-section 4.1.1, the computed ŝ satisfies:

ŝ = n i=1 a i n k=max (i,2) (1 + δ k ) = n k=1 a i ψ K i ,
where K i = {max (i, 2), . . . , n} and ψ Proof. By expectation of linearity, E(ŝ) = n i=1 a i E(ψ K i ). Lemma 11 shows that for all 1 ≤ i ≤ n, E(ψ K i ) = 1 and V (ψ K i ) ≤ γ n i (u 2 ), where n i is the cardinality of K i . It follows that E(ŝ) = s and

K i = n k=max (i,2) (1+δ k ) for all 1 ≤ i ≤ n.
V (ŝ) = V n i=1 a i ψ K i = σ n i=1 a i ψ K i 2 ≤ n i=1 σ(a i ψ K i ) 2 since σ(X + Y ) ≤ σ(X) + σ(Y ) = n i=1 |a i | V (ψ K i ) 2 ≤ n i=1 |a i | γ n i (u 2 ) 2 by Lemma 11. Since γ n i (u 2 ) ≤ γ n-1 (u 2 ), V (ŝ) ≤ γ n-1 (u 2 ) n i=1 |a i | 2 . (4.14) Bienaymé-Chebyshev inequality implies P |ŝ -E(ŝ)| ≤ V (ŝ)/λ ≥ 1 -λ. Thus |ŝ -s| |s| ≤ 1 |s| V (ŝ)/λ ≤ n i=1 |a i | |s| γ n-1 (u 2 )/λ = K γ n-1 (u 2 )/λ,
with probability at least 1 -λ.

Remark 4. Because E(ŝ) = s, under a normality assumption of ŝ, the number of significant bits can be lower-bounded by

-log 2 σ(ŝ) |E(ŝ)| ≥ -log 2 K γ n-1 (u 2 ) ≈ -log 2 (K) -log 2 (u) - 1 2 log 2 (n -1).
With K is the condition number and σ(ŝ) is the standard deviation of ŝ.

. Inner Product

In this sub-section, we bound the forward error of the inner product using the BC method. Consider the inner product

s n = y = a 1 b 1 + • • • + a n b n , evaluated from left to right, i.e, s i = s i-1 + a i b i , starting with s 1 = a 1 b 1 . Let δ 0 = 0, the computed ŝi satisfies ŝ1 = a 1 b 1 (1 + δ 1 ) and ŝi = (ŝ i-1 + a i b i (1 + δ 2i-2 ))(1 + δ 2i-1 ), |δ 2i-2 |, |δ 2i-1 | ≤ u,
for all 2 ≤ i ≤ n, where δ 2i-2 and δ 2i-1 represent the rounding errors from the products and additions, respectively. We thus have

ŷ = ŝn = n i=1 a i b i (1 + δ 2i-2 ) n k=i (1 + δ 2k-1 ).
Theorem 8. For all 0 < λ < 1, the computed ŷ satisfies under SR-nearness

E(ŷ) = y and |ŷ -y| |y| ≤ K γ n (u 2 )/λ, (4.15
)

with probability at least 1-λ, where K = n i=1 |a i b i | | n i=1 a i b i | is the condition number using the 1-norm for the computed y = n i=1 a i b i . Proof. For all 1 ≤ i ≤ n, we have ŷ = n i=1 a i b i (1 + δ 2i-2 ) n k=i (1 + δ 2k-1 ) = n i=1 a i b i ψ K i ,
with K 1 = {1, 3, . . . , 2n -1} and K i = {2i -2, 2i -1, 2i + 1, . . . , 2n -1} for all 2 ≤ i ≤ n. Lemma 11 shows that for all 1 ≤ i ≤ n, E(ψ

K i ) = 1 and V (ψ K i ) ≤ γ n i (u 2 )
, where n i is the cardinality of K i . Hence

E(ŷ) = E n i=1 a i b i ψ K i = n i=1 a i b i E(ψ K i ) = y. And V (ŷ) = V n i=1 a i b i ψ K i ≤ n i=1 |a i b i | V (ψ K i ) 2 since σ(X + Y ) ≤ σ(X) + σ(Y ) ≤ n i=1 |a i b i | γ n i (u 2 ) 2 by Lemma 11 ≤ γ n (u 2 )( n i=1 |a i b i |) 2 since γ n i (u 2 ) ≤ γ n (u 2 ) = y 2 K 2 γ n (u 2 ).
Thus, Bienaymé-Chebyshev inequality implies

|ŷ -y| |y| ≤ V (ŷ)/λ |y| ≤ K γ n (u 2 )/λ,
with probability at least 1 -λ.

Remark 5. Because E(ŷ) = y, under a normality assumption of ŷ, the number of significant bits can be lower-bounded by

-log 2 σ(ŷ) |E(ŷ)| ≥ -log 2 K γ n (u 2 ) ≈ -log 2 (K) -log 2 (u) - 1 2 log 2 (n).
With K is the condition number and σ(ŷ) is the standard deviation of ŷ.

. Horner Algorithm

In this sub-section, we bound the forward error of the Horner algorithm using the BC method. Recall that for the Horner algorithm , the computed r2n

satisfies r2n = n i=0 a i x i 2n k=2(n-i) (1 + δ k ),
Theorem 9. For all 0 < λ < 1, the computed r2n in the Equation (4.3) satisfies under SR-nearness E(r 2n ) = r 2n and

|r 2n -r 2n | |r 2n | ≤ K γ 2n (u 2 )/λ, (4.16
)

with probability at least 1-λ, where K = n i=0 |a i x i | | n i=0 a i x i |
is the condition number using the 1-norm for the computed P (x) = n i=0 a i x i .

Proof. We have

r2n = n i=0 a i x i 2n k=2(n-i) (1 + δ k ) = n i=0 a i x i ψ K i , with δ 0 = 0, K i = {2(n -i), 2(n -i) + 1, . . . , 2n}
for all 0 ≤ i ≤ n -1, and K n = {1, . . . , 2n}. Lemma 11 implies that for all 0 ≤ i ≤ n, E(ψ

K i ) = 1 and V (ψ K i ) ≤ γ n i (u 2 )
, where n i is the cardinality of K i . Then

E(r 2n ) = E n i=0 a i x i ψ K i = n i=0 a i x i E(ψ K i ) = r 2n .
And

V (r 2n ) = V n i=0 a i x i ψ K i ≤ n i=0 |a i x i | V (ψ K i ) 2 since σ(X + Y ) ≤ σ(X) + σ(Y ) ≤ n i=0 |a i x i | γ n i (u 2 )
2 by Lemma 11

≤ γ 2n (u 2 )( n i=0 |a i x i |) 2 since γ n i (u 2 ) ≤ γ 2n (u 2 ) = r 2 2n K 2 γ 2n (u 2 ).
Thus, Bienaymé-Chebyshev inequality implies

|r 2n -r 2n | |r 2n | ≤ V (r 2n )/λ |r 2n | ≤ K γ 2n (u 2 )/λ,
with probability at least 1 -λ.

Remark 6. Because E(r 2n ) = r 2n , under a normality assumption of r2n , the number of significant bits can be lower-bounded by

-log 2 σ(r 2n ) |E(r 2n )| ≥ -log 2 K γ 2n (u 2 ) ≈ -log 2 (K) -log 2 (u) - 1 2 log 2 (2n).

. Pairwise Summation

In this sub-section, we bound the forward error of the pairwise summation using the BC method. Recall that for the pairwise summation, the computed

Ŝh 1 satisfies Ŝh 1 = 2 h i=1 x i h j=1 (1 + δ j ⌈ i 2 j ⌉ ) = 2 h i=1 x i ψ K i ,
where

ψ K i = h j=1 (1 + δ j ⌈ i 2 j ⌉
) and the cardinality of K i is h for all 1 ≤ i ≤ 2 h . Theorem 10. For all 0 < λ < 1, the computed Ŝh

1 in the Equation (4.8) satisfies under SR-nearness E( Ŝh 1 ) = S h 1 and Ŝh 1 -S h 1 S h 1 ≤ κ γ log 2 (n) (u 2 )/λ, (4.17
)

with probability at least 1 -λ, where κ = n i=1 |x i | | n i=1 x i|
is the condition number of the summation of the x i .

Proof. By expectation linearity, E( Ŝh 1 ) =

2 h i=1 x i E(ψ K i ). Lemma 11 shows that for all 1 ≤ i ≤ 2 h , E(ψ i ) = 1 and V (ψ K i ) ≤ γ h (u 2 ). It follows that E( Ŝh 1 ) = S h 1 and V ( Ŝh 1 ) = V   2 h i=1 x i ψ K i   ≤   2 h i=1 |x i | V (ψ K i )   2 since σ(X + Y ) ≤ σ(X) + σ(Y ) ≤   2 h i=1 |x i | γ h (u 2 )   2 by Lemma 11 ≤ γ h (u 2 )( 2 h i=1 |x i |) 2 = ∥x∥ 2 1 γ h (u 2 ).
Bienaymé-Chebyshev inequality implies

P Ŝh 1 -E( Ŝh 1 ) ≤ V ( Ŝh 1 )/λ ≥ 1 -λ. Thus Ŝh 1 -S h 1 S h 1 ≤ 1 S h 1 V ( Ŝh 1 )/λ ≤ ∥x∥ 1 S h 1 γ h (u 2 )/λ = κ γ h (u 2 )/λ, with probability at least 1 -λ. Since h = log 2 (n) Ŝh 1 -S h 1 S h 1 ≤ κ γ log 2 (n) (u 2 )/λ,
with probability at least 1 -λ.

Remark 7. Because E( Ŝh 1 ) = S h 1 , under a normality assumption of Ŝh 1 , the number of significant bits can be lower-bounded by

-log 2 σ( Ŝh 1 ) |E( Ŝh 1 )| ≥ -log 2 K γ h (u 2 ) ≈ -log 2 (K) -log 2 (u) - 1 2 log 2 (h).

. Generalization

This method can also be generalized to any complete computation tree with a multi-linear sequence of elementary operations {+, -, * }. Using the same notations of Sub-section 4.1.4, in the case of the last operation of A is z ← x + y. We have z = x + y, then there exists δ such that ẑ = (x + ŷ)(1 + δ), with x = x(1 + Φ) and ŷ = y(1 + X).

Lemma 11 shows that we can assume by induction that

V (x(1 + δ)) ≤ K 2 x x 2 γ k (u 2 ), and
V (ŷ(1 + δ)) ≤ K 2 y y 2 γ l (u 2 ).
The following two theorems allow us to compute a probabilistic error bound in O( √ nu) of the algorithm A using the BC method.

Theorem 11. For all 0 < λ < 1, the computed ẑ satisfies under SR-nearness

E(ẑ) = z, and |ẑ -z| |z| ≤ K z γ m-1 (u 2 )/λ, (4.18
)

with probability at least 1 -λ, where

K z = K x |x| |x+y| + K y |y| |x+y| .
Proof. By construction of A, Lemma 11 shows that E(ẑ) = z. Note that

V (ẑ) = V ((x + ŷ)(1 + δ)) = V (x(1 + δ) + ŷ(1 + δ)) ≤ V (x(1 + δ)) + V (ŷ(1 + δ)) 2 .
Lemma 11 shows that

V (x(1 + δ)) ≤ K 2 x x 2 γ k (u 2 ), and
V (ŷ(1 + δ)) ≤ K 2 y y 2 γ l (u 2 ). Then, for m = max{k, l} V (ẑ) ≤ K x |x| γ k (u 2 ) + K y |y| γ l (u 2 ) 2 ≤ γ m-1 (u 2 )(K x |x| + K y |y|) 2 .

Bienaymé-Chebyshev inequality implies

P |ẑ -E(ẑ)| ≤ V (ẑ) /λ ≥ 1 -λ. Thus, |ẑ -z| |z| ≤ 1 |z| V (ẑ)/λ ≤ K x |x| + K y |y| |x + y| γ m-1 (u 2 )/λ = K z γ m-1 (u 2 )/λ,
with probability at least 1 -λ.

In the case of the last operation of A is z ← x × y. We have z = xy, then there exists δ such that ẑ = xŷ(1 + δ), with x = x(1 + Φ) and ŷ = y(1 + X).

Theorem 12. For all 0 < λ < 1, the computed ẑ satisfies under SR-nearness E(ẑ) = z, and

|ẑ -z| |z| ≤ K z γ m-1 (u 2 )/λ, (4.19
)

with probability at least 1 -λ, where

K z = K x K y .
Proof. By construction of A, Lemma 11 shows that E(ẑ) = z. Note that

V (ẑ) = V (xŷ(1 + δ)) = V (xy(1 + Φ)(1 + X)(1 + δ)) ≤ K 2 x x 2 K 2 y y 2 γ m-1 (u 2
) by Lemma 11. Because the worst error accumulation is of degree m-1. Bienaymé-Chebyshev inequality implies

|ẑ -z| |z| ≤ 1 |z| V (ẑ)/λ ≤ K x |x| K y |y| |xy| γ m-1 (u 2 )/λ = K z γ m-1 (u 2 )/λ,
with probability at least 1 -λ.

Discussion:

The BC method is a less restrictive result, as it requires few assumptions and is simple to apply. The Bienaymé-Chebyshev inequality yields an upper bound (variance) on the probability of deviation between the mean of a distribution and a specific value. In the case of SR-nearness, due to the unbiased nature of this stochastic rounding mode, a bound of the variance is sufficient to bound the error. Although this method is less accurate asymptotically in probability, in the next section, we showcase that the skillful handling of the variance bound outweighs this limitation and effectively governs the accuracy of the ultimate bound.

. Error Bounds Analysis

In this section, on the one hand, we compare the deterministic bound versus the probabilistic bounds. On the other hand, we undertake a comparative analysis of the two preceding methods, evaluating the tightness of their probabilistic bounds based on the precision used, the target probability, and the number of operations. We illustrate this analysis through the inner product as an example, but it also applies to the other algorithms. The advantage of the BC method in low-precision is clearly presented in Table 4.1.

. Inner Product

In the beginning, let us recall all bounds obtained for the inner product All bounds have the same condition number K, but differ in the others factor: γ n (u) for (Det-IP), γn ( 2 ln (2n/λ)) for (AH1-IP), uγ 2n (u) ln(2/λ) for (AH2-IP), and γ n (u 2 ) 1/λ for (BC-IP). In the following, for a constant λ, we investigate three cases: nu ≪ 1, nu ≫ 1 and n ln(n)u 2 ≪ 1, and n ln(n)u 2 ≫ 1 and nu 2 ≪ 1.

For n and u such that nu ≪ 1, we have

exp 2n ln (2n/λ)u + nu 2 1 -u -1 = 2n ln (2n/λ)u + O(u 2 ).
Moreover, [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]Lemma 3.1] implies

γ n (u) ≤ nu 1 -nu , it follows that for 2nu < 1, uγ 2n (u) ≤ 2nu 2 1 -2nu = u √ n 2 1 -2nu ,
and for nu 2 < 1

γ n (u 2 ) ≤ nu 2 1 -nu 2 = u √ n 1 √ 1 -nu 2 .
Interestingly, for the inner product, at any fixed probability, when nu ≪ 1, (AH2-IP) and (BC-IP) bounds are proportional to √ nu unlike √ n ln nu for the (AH1-IP) bound. Note that the deterministic bound is in O(nu).

For n and u such that nu ≫ 1 and n ln(n)u 2 ≪ 1, we have

exp 2n ln (2n/λ)u + nu 2 1 -u -1 = 2n ln (2n/λ)u + O(u 2 ).
Furthermore,

uγ 2n (u) ≈ u (exp (2nu) -1) ≈ √ u exp (nu), (4.20)
and

γ n (u 2 ) ≈ exp (nu 2 ) -1 = √ nu + O(u 2 ).
For n and u such that n ln(n)u 2 ≫ 1 and nu 2 ≪ 1, we have show that asymptotically (BC-IP) ≤ (AH1-IP) ≤ (AH2-IP) when nu ≫ 1 and nu 2 ≪ 1. In conclusion, these probabilistic bounds show that the roundoff error accumulated in n operations is proportional to √ nu rather than nu. In the next sub-section, we analyze these two probabilistic methods.

exp 2n ln (2n/λ)u + nu 2 1 -u -1 ≈ exp 2n ln (2n/λ)u + nu 2 1 -u ≈ exp 2n ln (2n/λ)u ,

. BC Method vs AH Method

In the following, we compare the three probabilistic bounds (AH1-IP), (AH2-IP) and (BC-IP) on the inner product forward error (similar reasoning can be applied to the other algorithms with the same result). When nu ≪ 1, at any fixed probability, (AH2-IP) and (BC-IP) are proportional to O( √ nu). First, we focus on this case. The two probabilistic bounds have the same condition number K. Thus, it is enough to compare u 2 γ 2n (u) 2 ln(2/λ) and γ n (u 2 ) 1/λ.

These two bounds depend on n and λ. Firstly, using the binomial theorem, we have

u 2 γ 2n (u) -γ n (u 2 ) = u 2 (1 + u 2 + 2u) n -1 -(1 + u 2 ) n -1 = u 2 n k=1 n k (u 2 + 2u) k - n k=1 n k (u 2 ) k = n k=1 n k u 2 (u 2 + 2u) k -(u 2 ) k ≥ 2 k=1 n k u 2 (u 2 + 2u) k -(u 2 ) k ≥ n(n - 1 2 )u 3 .
We can conclude that and BC bounds. Regarding BC bound, the variance calculation and the mean independence allow to bound the error terms (1+δ) 2 by (1+u 2 ) and avoid all δ terms of degree one because E(δ) = 0. In contrast, the AH1 and AH2 methods require bounded increments leading to terms (1 + u) 2 . As n increases, the advantage of Azuma-Hoeffding inequality for a probability near 1 becomes negligible.

γ n (u 2 ) ≤ u 2 γ 2n (u) for all n ≥ 1.
For all asymptotic comparisons between the bounds in this chapter, we have chosen to work with u → 0, n → ∞ and fixed probability λ, which we think adapted to many if not most current practical use cases. A situation with λ → 0 and fixed n gives the advantage to the Azuma-Hoeffding bounds over the Bienaymé-Chebyshev one.

Table 4.1 illustrates how BC bound is tighter than AH2 bound when n grows. The n threshold above which BC bound is preferable to AH2 bound depends on the format precision. The lower the precision, the lower the threshold becomes. Using SR in low-precision is of high interest in the areas of machine learning [START_REF] Gupta | Deep learning with limited numerical precision[END_REF], PDEs [START_REF] Croci | Effects of round-to-nearest and stochastic rounding in the numerical solution of the heat equation in low precision[END_REF], and ODEs [START_REF] Hopkins | Stochastic rounding and reduced-precision fixed-point arithmetic for solving neural ordinary differential equations[END_REF], motivating the use of our improved BC method.

. Numerical Experiments

This section presents numerical experiments that support and complete the theory presented previously. The various bounds are compared on two Table 4.1: The smallest n such that BC method gives a tighter probabilistic bound than AH2 method for the inner product.

numerical applications: the inner product and the evaluation of the Chebyshev polynomial.

We show that the probabilistic bounds are tighter than the deterministic bound and faithfully capture the behavior of SR-nearness forward error. For an inner product of large vectors, we show that BC bound is smaller than AH1 and AH2 bounds. All SR computations are repeated 30 times with verificarlo [START_REF] Denis | Verificarlo: Checking floating point accuracy through Monte Carlo arithmetic[END_REF]; we plot all samples and the forward error of the average of the 30 SR instances.

. Horner Algorithm

Let us firstly recall the previous error bounds obtained for this algorithm under SR-nearness:

| P (x) -P (x)| |P (x)| ≤ Kγ 2n (u), (Det-H) | P (x) -P (x)| |P (x)| ≤ K uγ 4n (u) ln 2 λ with probability ≥ 1 -λ, (AH-H) | P (x) -P (x)| |P (x)| ≤ K γ 2n (u 2 ) 1 λ with probability ≥ 1 -λ.
(BC-H)

Note that (Det-H) is the bound from Inequality (4.4), (AH-H) is the bound in Theorem 4, and (BC-H) is the bound in Theorem 9. We use Horner's method to evaluate the polynomial P (x) = T N (x) = N 2 i=0 a i (x 2 ) i where T N is the Chebyshev polynomial of even degree N = 2n. Both probabilistic bounds are tighter than the deterministic bound. For N = 20, there is no significant difference between SR-nearness and RN-binary32.

However, as expected, the average of the SR-nearness computations is more precise than the nearest round evaluation for almost all values of x.

= n i=0 |a i x i | | n i=0 a i x i | .
The evaluation in x = 24/26 ≈ 0.923 is plotted for various polynomial degrees N. As expected, when N increases, the deterministic bound grows faster than the probabilistic bounds. The right plot shows that Azuma-Hoeffding bound is tighter for a high probability and a small n. Overall, Chebyshev polynomial numerical experiment illustrates the advantage of the probabilistic error bounds over the deterministic error bound. However, for most of the evaluations in this experiment, RN-binary32 is more accurate than one instance of SR-nearness. This result is unsurprising because the degree n is small. To illustrate the behavior of these errors with a large n, we now turn to the inner product.

. Inner Product

To showcase the advantage of using BC method for large n, we present a numerical application of the inner product for vectors with positive floatingpoint numbers chosen uniformly at random between 0 and 1. For small n, AH1, AH2, and BC bounds are comparable with a slight advantage for (AH2-IP). However, as shown in Sub-section 4.3.1, when nu ≫ 1, the AH2 bound grows exponentially faster than AH1 and BC bounds. Asymptotically, the AH1 and BC bounds are therefore much tighter.

Interestingly, when n increases, a single instance of SR-nearness in bi-nary32 precision is more accurate than RN-binary32. This is because the summation terms are chosen uniformly between 0 and 1. The terms closest to zero are absorbed. With RN-binary32 the absorption errors are biased and will add up, while SR avoids stagnation and mitigates absorption errors. If we choose the terms in [-1; 1], SR and RN-binary32 have the same behavior. In this case, the absorption errors for RN-binary32 compensate because positive and negative errors are uniformly distributed. If we choose the terms in [1/2; 1], no absorption occurs for n < 2 23 , and on this domain, SR and RN-binary32 behave similarly. For many applications, SR results in a smaller accumulated error, for example, by avoiding stagnation effects. It satisfies the mean independence property, allowing tight probabilistic error bounds derived from our variance bound or the martingale theory. In the AH method, the martingale is readily apparent for computations such as summation or inner product. However, revealing the presence of the martingale can be difficult in some situations, as we have shown in our analysis of the Horner algorithm. This method uses the Azuma-Hoeffding inequality, which establishes a probabilistic bound and provides exponential tail bounds, making it a powerful tool for achieving excellent results asymptotically in probability.

For the BC method, the utilization of Lemma 11 dramatically simplifies the application of this approach. The Bienaymé-Chebyshev inequality often offers looser bounds on the tails of random variables, making it less accurate asymptotically in probability. Nevertheless, the quality of our bound on the variance dominates, leading to a final bound in O( √ nu). In several domains, the problem size n of computation is generally very large, as well as probabilistic investigations with a probability of 0.95 or 0.99 are frequently accepted. Table 4.1 shows the advantage of the BC method in low-precision with these probabilities.

Thanks to stochastic rounding and the previous two methods, we have been able to establish probabilistic bounds in O( √ nu) for algorithms with multi-linear errors based on {+, -, * }. It has been observed that any multilinear transformation of errors under SR-nearness forms a martingale: when errors result from elementary operations on exact inputs or affine functions of other errors or previous computations, they form a stochastic process given by ( k∈K (1 + δ k ) -1) with a natural filtration F, where K ⊂ N and certain δ k can be zero.

In order to effectively manage this stochastic process and achieve the desired F-martingale (martingale with respect to F), several techniques have been proposed depending on the algorithm employed. The Sub-section 2.4.2 summarizes the techniques proposed by Connelly et al. [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] and Ipsen, and Zhou [START_REF] Ipsen | Probabilistic error analysis for inner products[END_REF] for the inner product. The Sub-section 4.1.3 compares our technique to the Hallman and Ipsen [START_REF] Hallman | Precision-aware deterministic and probabilistic error bounds for floating point summation[END_REF] method for the pairwise summation. In these two problems, we have shown that the method used to construct the martingale is crucial in enhancing the quality of the probabilistic error analysis.

The shared property of the algorithms studied in this section is that the martingale emerges from the previous stochastic process without supplementary terms (bias or drift). This characteristic allows the direct application of probabilistic inequalities, such as the Azuma-Hoeffding inequality with bounded steps, enabling tight probabilistic bounds for higher probabilities. We have also demonstrated an alternative method based on a bound of the variance and Bienaymé-Chebyshev that ensures tight probabilistic bounds as the problem size increases and for a fixed probability.

Motivated by the applicability of this general framework to any numerical scheme featuring multi-linear errors, our goal in the next chapter is to extend its scope to encompass more complex situations: algorithms with errors that are not necessarily linear. The presence of non-linear errors will inevitably alter the form of the previous stochastic process. Hence, it is crucial to carefully handle this alteration and present a strategy to show the existence of the martingale. Moreover, it is necessary to control the effect of the nonlinear nature on the entire computation process and the quality of the final probabilistic bound. The main objective is to establish a general framework covering the probabilistic error analysis of a wide range of algorithms.

The scripts for reproducing the numerical experiments of this chapter are published in the repository: https://github.com/verificarlo/sr-variance-bounds/.

-Error Analysis for Algorithms with Non-linear Errors

Previous theoretical studies of SR error bounds have only considered algorithms based on sums, and sometimes products of uncorrelated random errors (i.e, ) in which the resulting error is a multi-linear function of each operation rounding error. As shown in Chapter 4, SR error analysis of these algorithms demonstrates the existence of a stochastic process ( k∈K (1+δ k )-1), associated with a natural filtration F, where K ⊂ N and certain δ k can be zero.

Two main methods have been proposed to investigate this stochastic process and bound the forward error of algorithms: AH and BC. For the AH method, the technique used to construct the martingale impacts the quality of the final bound. For the BC method, effectively managing the bound of the variance allows us to obtain tight probabilistic error bounds.

To the best of our knowledge, there is no previous theoretical research on studying non-linear problems with SR. This chapter introduces a general framework that allows the probabilistic error analysis of algorithms with linear and non-linear errors under SR. In Section 5.1, and using the Doob-Meyer decomposition [START_REF] Doob | Stochastic Processes, Probability and Statistics Series[END_REF], we demonstrate how the error of an algorithm can be decomposed in a martingale and a drift. Whatever the studied algorithm, this decomposition illustrates how to separate the errors that can be compensated with SR, corresponding to the martingale term, from the remaining errors, corresponding to the drift.

Our analysis shows that the drift of any computation tree (algorithms with multi-linear errors) is zero, indicating that its average error contribution remains consistent over time. Conversely, the probabilistic error analysis of algorithms with non-linear errors demonstrates a non-zero drift, suggesting that their error average contributions are subject to variations over time.

The study of this decomposition consists in probabilistically bounding the martingale term using the BC or AH method. However, a deterministic drift analysis reveals that it is negligible at the first order. Consequently, the martingale term dominates the entire investigation, and effective management of this component allows the derivation of tight probabilistic bounds. In the following, this approach will be called DM bound.

We show how this general framework can be applied to give SR error bounds for algorithms that compute the variance, called: "textbook-variance" and "two-pass-variance". In 1983, Chan, Golub, and LeVeque proved deterministic error bounds [START_REF] Chan | Algorithms for computing the sample variance: Analysis and recommendations[END_REF] for different algorithms computing the variance of a sample of n data points. These algorithms have non-linear errors due to the presence of squaring in the computation.

Our probabilistic error analysis demonstrates a non-zero drift for both algorithms and proves their SR forward error bounds in the recursive and pairwise cases. For the textbook-variance algorithm (Section 5.2), a part of the error forms a martingale directly, so we apply this decomposition to the remaining part of the error. Since the set of martingales forms a vector space, we can present the error as a martingale and a drift. For the two-pass-variance algorithm (Section 5.3), we apply the decomposition to the entire error. We illustrate how this decomposition separates the errors that can be compensated with SR, which corresponds to the martingale term, from the remainder. For both algorithms, we show probabilistic bounds in O( √ nu) instead of O(nu) for the deterministic bounds.

Although the DM method provides probabilistic bounds in O( √ nu), we have proposed an approach for the variance computation without the utilization of Doob-Meyer decomposition. We have introduced analogous concepts aligned with the general framework. Using simple yet effective techniques, we have successfully managed the non-linearity present in these algorithms, resulting in the use of BC and AH methods to obtain probabilistic bounds in O( √ nu).

. General Framework

The goal of this section is to present a general framework for probabilistic error analysis with stochastic rounding. Given a dataset X, we represent the application of an algorithm to this data by H(X, 0), and the SR computation is represented by H(X, δ). For instance, in the case of a summation of n values x 1 , . . . , x n , the exact computation is given by H(X, 0) = n i=1 x i , and the SR computation is given by

H(X, δ) = n i=1 x i n k=max i,2 (1 + δ k ).
To establish a comprehensive framework, we use the Doob-Meyer decomposition, a fundamental theorem in stochastic analysis [20, p. 296]. The Doob-Meyer decomposition separates a stochastic process into two distinct components: a martingale part and a predictable process. This decomposition provides valuable insights into the behavior and characteristics of the process under consideration. Let us first define a predictable stochastic process [p. 65] [START_REF] Dacunha-Castelle | Probability and Statistics[END_REF].

Definition 4. Given a filtration (F n ) n≥0 , a stochastic process X n is predictable if X 0 is F 0 -measurable, and X n is F n-1 -measurable for all n ≥ 1.
This means that the value of X n is known at the previous time step. Now, let us state the Doob-Meyer decomposition.

Theorem 13 (Doob-Meyer decomposition). Let (F k ) 0≤k≤n and X 0 , . . . , X n an adapted stochastic process locally integrable, meaning that E(|X k |) < ∞ for all 0 ≤ k ≤ n. There exists a martingale M 0 , . . . , M n and a predictable integrable sequence A 0 , . . . , A n starting with A 0 = 0 for which we have:

     X n = M n + A n , A n = n k=1 E[X k -X k-1 /F k-1 ], E(M n ) = 0.
This decomposition is almost surely unique.

Proof. Let us show that A n is predictable:

E[A n /F n-1 ] = E n k=1 E[X k -X k-1 /F k-1 ]/F n-1 = n k=1 E [E[X k -X k-1 /F k-1 ]/F n-1 ] by linearity = n k=1 E[X k -X k-1 /F k-1 ] = A n . The last equality is because E[X k -X k-1 /F k-1 ] is F k-1 -measurable, then
F n-1 -measurable. Let us show that M n forms a martingale. The first two points of the definition 3 are evident. Let us check the point three:

E[M n /F n-1 ] = E[X n -A n /F n-1 ] = E[X n /F n-1 ] -E[A n /F n-1 ] = E[X n /F n-1 ] -A n since A n is predictable = E[X n /F n-1 ] - n k=1 E[X k -X k-1 /F k-1 ] = E[X n /F n-1 ] -E[X n -X n-1 /F n-1 ] -A n-1 = E[X n-1 /F n-1 ] -A n-1 = X n-1 -A n-1 = M n-1 .
Thus, M n forms a martingale with respect to (F k ) 0≤k≤n , and the decomposition is valid. Moreover, using the expectation linearity and the law of total expectation, we have

E(M n ) = E (X n -A n ) = E(X n ) -E n k=1 E[X k -X k-1 /F k-1 ] = E(X n ) - n k=1 E (E[X k -X k-1 /F k-1 ]) = E(X n ) - n k=1 E(X k -X k-1 ) = E(X n ) -E(X n ) since E(X 0 ) = 0 = 0.
Remark 8. X 0 , . . . , X n is a sub-martingale if and only if the predictable process A 0 , . . . , A n is an increasing process.

The martingale component does not exhibit any systematic drift or trend.

It captures the random behavior of the stochastic process X n and reflects the information available up to time n. While the sequence A n represents the cumulative effect of the predictable part of the stochastic process X n . It can be interpreted as the drift of X n .

Our investigation aims to examine the error |H(X, δ) -H(X, 0)| under SRnearness using Doob-Meyer decomposition (DM method). The key idea is to represent any error under SR-nearness as a decomposition of a martingale and a drift:

H(X, δ) -H(X, 0) = M + A. (5.1)
Since each random error δ satisfies |δ| ≤ u, the assumption of locally integrable in Theorem 13 is always valid in the context of the error analysis under SR-nearness. Therefore, this decomposition is valid whatever the nature of the algorithm H and the dataset X.

The martingale component in Equation (5.1) represents the inherent stochastic behavior captured by the decomposition; in other words, the errors that can be compensated with SR. Hence, we conduct a probabilistic analysis of this component using AH or BC methods. Our analysis shows probabilistic bounds in O( √ nu) for this decomposition part instead of the deterministic bounds in O(nu). If the drift is non-zero, we deterministically bound its influence and establish its negligible effect at the first order, which shows that its impact is insignificant compared to the martingale factor.

Interestingly, the use of Equation (5.1) is found to be straightforward when applied to the problems discussed in Chapter 4. In the case of multi-linear errors, the drift component becomes zero owing to the unbiased nature of SR-nearness, resulting in the direct acquisition of the martingale. Theorem 14. Let H an algorithm with multi-linear error and X a dataset. The Equation (5.1) is given by:

H(X, δ) -H(X, 0) = M,
where M is the martingale.

Proof. We associate the proof to the sequential summation, but is remains valid for all algorithms with multi-linear error. We have H(X, δ) = ŝ and H(X, 0) = s. In the proof of Theorem 2, we have shown that ŝ -s = Z n , where Z n is a martingale. Then

H(X, δ) -H(X, 0) = ŝ -s = Z n = M,
and A = 0. This observation highlights using Equation (5.1) in the context of algorithms with multi-linear error, where the decomposition simplifies to a single term.

In the following, we examine SR for non-linear computations via two algorithms that compute the variance: textbook-variance and two-pass-variance. Through the previous general framework, especially, Equation (5.1), we introduce a new approach to establish probabilistic bounds on the forward error (DM method). We also use the BC and AH methods discussed in Chapter 4, and we show probabilistic error bounds in O( √ nu).

. Error Analysis for Textbook-Variance Algorithm

For x ∈ R n , let s = n i=1 x i . For the textbook-variance algorithm, we have the following:

H(X, 0) = n i=1 x 2 i - 1 n s 2 .
The condition number using the 2-norm for the variance computation is defined in [START_REF] Chan | Algorithms for computing the sample variance: Analysis and recommendations[END_REF] as

K 2 = ∥x∥ 2 H(X, 0) .
We define the condition number using the 1-norm by

K 1 = ∥x∥ 1 nH(X, 0) .
Using Cauchy-Schwarz inequality,

K 1 ≤ K 2 .
K 1 can be lower than 1 (for instance, consider n = 4 and x 1 = 1/2, x 2 = 1/4, x 3 = -x 1 and x 4 = -x 2 ). On the one hand, we consider that the computation of s and ŝ is as follows:

Stochastic rounding Exact computation

ŝ1 = x 1 s 1 = x 1 ŝ2 = (ŝ 1 + x 2 )(1 + δ 1 ) s 2 = s 1 + x 2 ŝk = (ŝ k-1 + x k )(1 + δ k-1 ) s k = s k-1 + x k ŝn = ŝ s n = s It follows that ŝn = n i=1 x i n k=max(2,i) (1 + δ k-1 ) = n i=1 x i ϕ i , with ϕ i = n k=max(2,i) (1 + δ k-1 ) for all 1 ≤ i ≤ n. On the other hand, H(X, δ) = n i=1 x 2 i ψ i - 1 n s 2 ψ n+1 ,
where

ψ i = (1 + ϵ i ) n+1 k=max(2,i) (1 + η k ) and ψ n+1 = (1 + ϵ n+1 )(1 + η n+1 )(1 + θ).
For all 1 ≤ i ≤ n, ϵ i and η i represent the rounding errors from the products and additions, respectively. ϵ n+1 represents the error of the square of s, η n+1 represents the error of the subtraction, and θ represents the error of the division of s 2 by n.

Denote

Z k = s k -s k = Z k-1 + ( s k-1 + x k )δ k-1 . Then, Z n = s n -s n .
As we have demonstrated in Sub-section 4.1.1 Z 1 , . . . , Z n form a martingale with respect to δ 1 , . . . , δ n-1 . Since the set of martingales forms a vector space, Z 1 + s, . . . , Z n + s is also a martingale. Denote:

• F k = {δ 1 , . . . , δ k }. • Y k-1 = Z k -Z k-1 = ( s k-1 + x k )δ k-1 for all 2 ≤ k ≤ n. Then Z n = n k=2 Y k-1 . • σ 2 k-1 = E[Y 2 k-1 /F k-2 ]. • A n = n k=2 σ 2 k-1 with A 1 = 0.
Lemma 12. The stochastic process A n is predictable.

Proof.

E[A n /F n-1 ] = E n k=2 σ 2 k-1 /F n-1 = E n k=2 E Y 2 k-1 /F k-2 /F n-1 = n k=2 E[E[Y 2 k-1 /F k-2 ]/F n-1 ]. Since E[Y 2 k-1 /F k-2 ] is F k-2 -measurable, so it is F n-1 -measurable, and for all 2 ≤ k ≤ n, we have E[E[Y 2 k-1 /F k-2 ]/F n-1 ] = E[Y 2 k-1 /F k-2 ] . Then E[A n /F n-1 ] = n k=2 E[Y 2 k-1 /F k-2 ] = A n .
Lemma 13. The stochastic process

X n = (Z n + s) 2 -A n -s 2 forms a martingale with respect to F n .
Proof.

E[X n /F n-1 ] = E[(Z n + s) 2 -A n -s 2 /F n-1 ] = E[(Z n-1 + s + Y n-1 ) 2 /F n-1 ] -A n -s 2 = (Z n-1 + s) 2 + 2(Z n-1 + s)E[Y n-1 /F n-1 ] + E[Y 2 n-1 /F n-1 ] -A n -s 2 = X n-1 because E[Y n-1 /F n-1 ] = 0.
Then, X n forms a martingale with respect to F n .

Note that (Z n + s) 2 is a sub-martingale, and the expression of

(Z n + s) 2 = X n + s 2 + A n is a Doob-Meyer decomposition. Moreover, since Z n = s n -s n H(X, δ) -H(X, 0) = n i=1 x 2 i ψ i - 1 n s 2 ψ n+1 - n i=1 x 2 i + 1 n s 2 = n i=1 x 2 i (ψ i -1) - 1 n (Z n + s) 2 ψ n+1 + 1 n s 2 = n i=1 x 2 i (ψ i -1) - 1 n ψ n+1 (X n + s 2 + A n ) + 1 n s 2 = n i=1 x 2 i (ψ i -1) - 1 n ψ n+1 X n - 1 n s 2 (ψ n+1 -1) - 1 n ψ n+1 A n = M + A.
Finally,

H(X, δ) -H(X, 0) = M + A, (5.2)
where

M = n i=1 x 2 i (ψ i -1) - 1 n ψ n+1 X n - 1 n s 2 (ψ n+1 -1), and A = - 1 n ψ n+1 A n .
Interestingly, M constructs a martingale as the summation of three martingales. In fact, we have shown a martingale from the inner product, which is the case for n i=1 x 2 i (ψ i -1). Clearly, 1 n s 2 (ψ n+1 -1) forms a martingale. Since the set of martingales is a vector space, M is a martingale.

. Bias Analysis

The unbiased nature of SR-nearness extends to various algorithms such as the inner product [START_REF] Connolly | Stochastic rounding and its probabilistic backward error analysis[END_REF] and Horner's rule [START_REF] Arar | The positive effects of stochastic rounding in numerical algorithms[END_REF]. However, it fails to hold in the general case. In the sequel, we compute the bias of the textbook-variance algorithm.

Theorem 15. For the textbook-variance algorithm, the bias is given by

E (H(X, δ) -H(X, 0)) = - 1 n V (ŝ).
Proof. We know that

E (H(X, δ) -H(X, 0)) = E(M + A) = E(M ) + E(A).
Since X n is a martingale, we have E(X n ) = E(X 0 ) and

X 0 = (Z 0 + s) 2 -s 2 -A 0 = s 2 -s 2 = 0.
Then E(X n ) = 0. Lemma 11 shows that E(ψ i ) = 1 for all 1 ≤ i ≤ n + 1. Then,

E n i=1
x 2 i (ψ i -1) = 0, and E

1 n s 2 (ψ n+1 -1) = 0.
It follows that E(M ) = 0. Moreover, the law of total expectation yields

E(A) = - 1 n E(ψ n+1 A n ) = - 1 n E (E[ψ n+1 A n /F n-1 ]) = - 1 n E (A n E[ψ n+1 /F n-1 ]) since A n is predictable = - 1 n E(A n ) by Lemma 2 = - 1 n E n k=2 E[Y 2 k-1 /F k-2 ] = - 1 n n k=2 E E[Y 2 k-1 /F k-2 ] = - 1 n n k=2 E Y 2 k-1 . Recall that Y k-1 = Z k -Z k-1 . Then E Y 2 k-1 = E (Z k -Z k-1 ) 2 = E(Z 2 k ) + E(Z 2 k-1 ) -2E(Z k Z k-1 ) = E(Z 2 k ) + E(Z 2 k-1 ) -2E (E[Z k Z k-1 /F k-1 ]) = E(Z 2 k ) + E(Z 2 k-1 ) -2E (Z k-1 E[Z k /F k-1 ]) = E(Z 2 k ) + E(Z 2 k-1 ) -2E(Z 2 k-1 ) since Z k is a martingale = E(Z 2 k ) -E(Z 2 k-1 ).
Because Z 1 = 0 and E(Z n ) = E( s -s) = 0, we have

E(A) = - 1 n n k=2 E Y 2 k-1 = - 1 n n k=2 E(Z 2 k ) -E(Z 2 k-1 ) = - 1 n E(Z 2 n ) = - 1 n V (Z n ) = - 1 n V (ŝ -s) = - 1 n V (ŝ).
Finally,

E (H(X, δ) -H(X, 0)) = E(M ) + E(A) = - 1 n V (ŝ). Lemma 15. The martingale X 1 , . . . , X n satisfies |X k -X k-1 | ≤ uC k , for all 2 ≤ k ≤ n, where C k = ∥x∥ 2 1 (1 + u) 2(k-2) (2 + u).
Proof. Note that by definition of F k-2

σ 2 k-1 = E[Y 2 k-1 /F k-2 ] = E[( s k-1 + x k ) 2 δ 2 k-1 /F k-2 ] = ( s k-1 + x k ) 2 E[δ 2 k-1 /F k-2 ]. Because Z k = Z k-1 + ( s k-1 + x k )δ k-1 and A k = A k-1 + σ 2 k-1 , we have X k -X k-1 = (Z k + s) 2 -A k -(Z k-1 + s) 2 + A k-1 = (Z k-1 + s + ( s k-1 + x k )δ k-1 ) 2 -A k -(Z k-1 + s) 2 + A k-1 = 2(Z k-1 + s)( s k-1 + x k )δ k-1 + ( s k-1 + x k ) 2 δ 2 k-1 -σ 2 k-1 = 2(Z k-1 + s)( s k-1 + x k )δ k-1 + ( s k-1 + x k ) 2 δ 2 k-1 -E[δ 2 k-1 /F k-2 ] . Since |δ k-1 | ≤ u and 0 ≤ δ 2 k-1 ≤ u 2 , we have δ 2 k-1 -E[δ 2 k-1 /F k-2 ] ≤ u 2 and | s k-1 + x k | ≤ (1 + u) k-2 k i=1 |x i | ≤ (1 + u) k-2 ∥x∥ 1 .
It follows that

|Z k-1 + s| ≤ |Z k-1 | + |s| ≤ ((1 + u) k-2 -1) k-1 i=1 |x i | + |s| ≤ ∥x∥ 1 (1 + u) k-2 . Thus |X k -X k-1 | = 2δ k-1 (Z k-1 + s)( s k-1 + x k ) + ( s k-1 + x k ) 2 δ 2 k-1 -E[δ 2 k-1 /F k-2 ] ≤ 2u |Z k-1 + s| | s k-1 + x k | + u 2 | s k-1 + x k | 2 ≤ 2u(1 + u) 2(k-2) ∥x∥ 2 1 + u 2 (1 + u) 2(k-2) ∥x∥ 2 1 = u ∥x∥ 2 1 (1 + u) 2(k-2) (2 + u).
Theorem 16. For 0 < λ < 1, the martingale X 1 , . . . , X n satisfies under SRnearness

|X n | ≤ ∥x∥ 2 1 2uγ 4(n-1) (u) ln(2/λ), (5.3)
with probability at least 1 -λ.

Proof. Since X 1 = 0, Lemma 3 and Lemma 15 yield

|X n | ≤ n k=2 u 2 C 2 k 2 ln(2/λ),
with probability at least 1 -λ. Furthermore

n k=2 u 2 C 2 k = u 2 n k=2 ∥x∥ 4 1 (1 + u) 4(k-2) (2 + u) 2 = u 2 ∥x∥ 4 1 (2 + u) 2 γ 4(n-1) (u) (1 + u) 4 -1 = u ∥x∥ 4 1 4 + 4u + u 2 4 + 6u + 4u 2 + u 3 γ 4(n-1) (u) ≤ u ∥x∥ 4 1 γ 4(n-1) (u).
Finally,

|X n | ≤ ∥x∥ 2 1 2uγ 4(n-1) (u) ln(2/λ).
We are now in a position to state the main result of this sub-section.

Theorem 17. For all 0 < λ < 1, the computed H(X, δ) satisfies under SRnearness

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤K 2 1 (1 + u) 3 2uγ 4(n-1) (u) ln(4/λ) + u γ 2(n-1) (u) 2 + 1 + K 2 2 uγ 2(n+1) (u) ln(4/λ) -K 2 1 ,
with probability at least 1 -λ, where

K 1 = ∥x∥ 1 √ nH(X,0)
and

K 2 = ∥x∥ 2 √ H(X,0) . Proof. Recall that |H(X, δ) -H(X, 0)| = |M + A| ≤ |M | + |A| , where M = n i=1 x 2 i (ψ i -1) - 1 n ψ n+1 X n - 1 n s 2 (ψ n+1 -1) and A = - 1 n ψ n+1 A n .
Because |ψ n+1 | ≤ (1 + u) 3 and |s| ≤ ∥x∥ 1 , we firstly deduce that

|M | ≤ n i=1 x 2 i (ψ i -1) + 1 n (1 + u) 3 |X n | + 1 n ∥x∥ 2 1 γ 3 (u).
Theorem 16 states that

|X n | ≤ ∥x∥
(u) = ∥x∥ 2 2 uγ 2(n+1) (u) ln(4/λ) + 1 n (1 + u) 3 ∥x∥ 2 1 2uγ 4(n-1) (u) ln(4/λ) + 1 - 1 n ∥x∥ 2 1 ,
with probability at least 1 -λ.

Secondly, A n = n k=2 E[Y 2 k-1 /F k-2 ] = n k=2 ( s k-1 + x k ) 2 E[δ 2 k-1 /F k-2 ],
then

|A n | ≤ u 2 n k=2 | s k-1 + x k | 2 ≤ u 2 n k=2 (1 + u) k-2 k i=1 |x i | 2 ≤ u 2 ∥x∥ 2 1 n k=2 (1 + u) 2(k-2) ≤ u 2 ∥x∥ 2 1 γ 2(n-1) (u) 2u + u 2 ≤ u ∥x∥ 2 1 γ 2(n-1) (u) 2 .
Finally,

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤ |M | |H(X, 0)| + |A| |H(X, 0)| ≤ ∥x∥ 2 1 n |H(X, 0)| (1 + u) 3 2uγ 4(n-1) (u) ln(4/λ) + u γ 2(n-1) (u) 2 + 1 + ∥x∥ 2 2 |H(X, 0)| uγ 2(n+1) (u) ln(4/λ) - ∥x∥ 2 1 n |H(X, 0)| =K 2 1 (1 + u) 3 2uγ 4(n-1) (u) ln(4/λ) + u γ 2(n-1) (u) 2 + 1 + K 2 2 uγ 2(n+1) (u) ln(4/λ) -K 2 1 ,
with probability at least 1 -λ.

Discussion

The DM method contributes to the probabilistic error analysis of algorithms under SR. It allows the decomposition of the error of any algorithm into two parts: a martingale and a drift. The study of this decomposition shows probabilistic error bounds in O( √ nu) instead of O(nu) for the deterministic bounds. However, we will propose similar ideas to the general framework with simple techniques for the variance computation.

. BC Method

This sub-section uses the BC method and Lemma 14 to provide a probabilistic bound on the forward error of the textbook-variance algorithm under SR-nearness. Let us compute H(X, δ) -H(X, 0) from another angle:

|H(X, δ) -H(X, 0)| = n i=1 x 2 i (ψ i -1) - 1 n ( s 2 ψ n+1 -s 2 ) ≤ n i=1 x 2 i (ψ i -1) + 1 n s 2 ψ n+1 -s 2 = n i=1 x 2 i (ψ i -1) + 1 n (( s -s) + s) 2 ψ n+1 -s 2 ≤ n i=1 x 2 i (ψ i -1) + 1 n B,
where

B = ( s -s) 2 ψ n+1 + 2 |s( s -s)ψ n+1 | + s 2 (ψ n+1 -1)
. The following inequality will be used in the proofs of the textbook-variance forward errors by BC and AH methods:

|H(X, δ) -H(X, 0)| ≤ n i=1 x 2 i (ψ i -1) + 1 n B.
(5.4)

Remark 10. To handle the non-linearity of errors, the key idea of this approach is to isolate terms of order 1 in errors and then use the previous results on the inner product or summation. Other decompositions could be used. For instance,

1 n ( s 2 ψ n+1 -s 2 ) = 1 n ( s 2 ψ n+1 -ss + ss -s 2 ) = 1 n ( s( sψ n+1 -s) + s( s -s)) .
Then, we can apply the same properties on ( sψ n+1 -s) and ( s-s). The bounds are different but asymptotically equivalent when nu ≪ 1.

The rounding errors accumulated in the whole process of this algorithm ϕ i and ψ i satisfy for all 1 ≤ i ≤ n, (2,i) and |ψ n+1 | ≤ (1 + u) 3 .

|ϕ i | ≤ (1 + u) n+1-max(2,i) , |ψ i | ≤ (1 + u) n+3-max
Let us compute the deterministic bound of this algorithm. We have n i=1

x 2 i (ψ i -1) ≤ ∥x∥ 2 2 γ n+1 (u). Since |s| ≤ ∥x∥ 1 and | s -s| = | n i=1 x i (ϕ i -1)| ≤ ∥x∥ 1 γ n-1 (u), B ≤ (1 + u) 3 ∥x∥ 2 1 γ 2 n-1 (u) + 2γ n-1 (u) + ∥x∥ 2 1 ((1 + u) 3 -1) = (1 + u) 3 ∥x∥ 2 1 γ 2 n-1 (u) + 2γ n-1 (u) + 1 -∥x∥ 2 1 = (1 + u) 3 ∥x∥ 2 1 (γ n-1 (u) + 1) 2 -∥x∥ 2 1 = ∥x∥ 2 1 (1 + u) 2n+1 -∥x∥ 2 1 = ∥x∥ 2 1 γ 2n+1 (u). Finally |H(X, δ) -H(X, 0)| |H(X, 0)| ≤ K 2 2 γ n+1 (u) + K 2 1 γ 2n+1 (u), (5.5)
where

K 1 = ∥x∥ 1 √ nH(X,0)
and

K 2 = ∥x∥ 2 √ H(X,0)
.

The following theorem presents a probabilistic bound of the forward error of this algorithm through the BC method.

Theorem 18. For all 0 < λ < 1, the computed H(X, δ) satisfies under SRnearness

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤ K 2 2 2γ n+1 (u 2 )/λ+K 2 1 (1 + u) 3 2γ n-1 (u 2 )/λ + 1 2 -1 ,
with probability at least 1 -λ, where

K 1 = ∥x∥ 1 √ nH(X,0)
and

K 2 =
x 2 i (ψ i -1) ≤ ∥x∥ 2 2 2γ n+1 (u 2 )/λ with probability at least 1 - λ 2 , | s -s| ≤ ∥x∥ 1 2γ n-1 (u 2 )/λ with probability at least 1 - λ 2 .
Since,

|ψ n+1 | ≤ (1 + u) 3 and |s| ≤ ∥x∥ 1 , B ≤ (1 + u) 3 ∥x∥ 2 1 2γ n-1 (u 2 )/λ + 2 2γ n-1 (u 2 )/λ + ∥x∥ 2 1 (1 + u) 3 -1 = (1 + u) 3 ∥x∥ 2 1 2γ n-1 (u 2 )/λ + 2 2γ n-1 (u 2 )/λ + 1 -∥x∥ 2 1 = (1 + u) 3 ∥x∥ 2 1 2γ n-1 (u 2 )/λ + 1 2 -∥x∥ 2 1 ,

95

with probability at least 1 -λ 2 . Finally, Lemma 14 shows that

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤ 1 |H(X, 0)| n i=1 x 2 i (ψ i -1) + 1 n |H(X, 0)| B ≤ K 2 2 2γ n+1 (u 2 )/λ + K 2 1 (1 + u) 3 2γ n-1 (u 2 )/λ + 1 2 -1 ,
with probability at least 1 -λ.

. AH Method

This sub-section uses the AH method and Lemma 14 to provide a probabilistic bound of the forward error of the textbook-variance algorithm under SR-nearness.

Theorem 19. For all 0 < λ < 1, the computed H(X, δ) satisfies under SRnearness

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤K 2 2 uγ 2(n+1) (u) ln(4/λ) + K 2 1 (1 + u) 3 uγ 2(n-1) (u) ln(4/λ) + 1 2 -1 , with probability at least 1 -λ. Where K 1 = ∥x∥ 1 √ nH(X,0)
and

K 2 = ∥x∥ 2 √ H(X,0)
. Proof. Equation (5.4) states that

|H(X, δ) -H(X, 0)| ≤ n i=1 x 2 i (ψ i -1) + 1 n B.
Moreover, [46, cor 4.7] shows that n i=1

x 2 i (ψ i -1) ≤ ∥x∥ 2 2 uγ 2(n+1) (u) ln(4/λ) with probability at least 1 - λ 2 , | s -s| ≤ ∥x∥ 1 uγ 2(n-1) (u) ln(4/λ) with probability at least 1 - λ 2 .
Since,

|ψ n+1 | ≤ (1 + u) 3 and |s| ≤ ∥x∥ 1 , B ≤(1 + u) 3 ∥x∥ 2 1 uγ 2(n-1) (u) ln(4/λ) + 2 uγ 2(n-1) (u) ln(4/λ) + ∥x∥ 2 1 (1 + u) 3 -1 =(1 + u) 3 ∥x∥ 2 1 uγ 2(n-1) (u) ln(4/λ) + 2 uγ 2(n-1) (u) ln(4/λ) + 1 -∥x∥ 2 1 =(1 + u) 3 ∥x∥ 2 1 uγ 2(n-1) (u) ln(4/λ) + 1 2 -∥x∥ 2 1 , |H(X, δ) -H(X, 0)| |H(X, 0)| ≤K 2 2 uγ 2(n+1) (u) ln(4/λ) + K 2 1 (1 + u) 3 uγ 2(n-1) (u) ln(4/λ) + 1 2 -1 ,
with probability at least 1 -λ.

. Error Analysis for the Two-pass-Variance Algorithm

This section aims to analyze the error of the two-pass-variance algorithm under SR-nearness. We employ the DM, BC, and AH methods to probabilistically bound the forward error of this algorithm. In Section 5.1, we demonstrated that for any computational tree with multi-linear errors, Equation (5.1) results in a single term: martingale. Additionally, in Section 5.2, we showed how this Equation can be applied to a part of the error for the textbookvariance algorithm while the remaining part forms a martingale. However, we apply this Equation to the entire error and show probabilistic error bounds of O( √ nu) for the two-pass-variance algorithm. For x ∈ R n , let m = 1 n n i=1 x i . For the two-pass-variance algorithm we have:

H(X, 0) = n i=1 (x i -m) 2 .
As for the computation of the summation s in Equation (4.1), the computed m satisfies m = 1

n n i=1 x i n+1 k=max(2,i) (1 + δ k-1 )
, where δ n represents the rounding error of the division by n. Hence, H(X, δ) satisfy

H(X, δ) = n i=1 (x i -m) 2 ψ i ,
where

ψ i = (1 + ϵ i ) 2 (1 + η i ) n k=max(2,i) (1 + θ k ).
For all 1 ≤ i ≤ n, ϵ i , η i and θ i represent the rounding errors of subtraction, square, and addition, respectively. In the following table, we don't take into account the errors accumulated in the computation of m, and we consider the following table:

Stochastic rounding

Exact computation

ŝ1 = (x 1 -m) 2 (1 + ϵ 1 ) 2 (1 + η 1 ) s 1 = (x 1 -m) 2 ŝ2 = ŝ1 + (x 2 -m) 2 (1 + ϵ 2 ) 2 (1 + η 2 ) s 2 = s 1 + (x 2 -m) 2 ŝ3 = ŝ2 (1 + θ 2 ) s 3 = s 2 ŝ2k-2 = ŝ2k-3 + (x k -m) 2 (1 + ϵ k ) 2 (1 + η k ) s 2k-2 = s 2k-3 + (x k -m) 2 ŝ2k-1 = ŝ2k-2 (1 + θ k ) s 2k-1 = s 2k-2 ŝ2n-1 = H(X, δ) s 2n-1 = H(X, 0) Let us denote • X k = ŝk -s k for all 0 ≤ k ≤ 2n -1 with X 0 = 0. Then, H(X, δ) - H(X, 0) = ŝ2n-1 -s 2n-1 = X 2n-1 . • F 0 = {δ i , ϵ 1 , i ∈ [1; n]}. • F 2k-3 = {δ i , ϵ j , η l , θ l , i ∈ [1; n], j ∈ [1; k], and l ∈ [1; k -1]}. • F 2k-2 = {δ i , ϵ j , η j , θ l , i ∈ [1; n], j ∈ [1; k], and l ∈ [1; k -1]}. • A 2k-1 = 2k-1 i=1 E[X i -X i-1 /F k-1 ] for all 1 ≤ k ≤ 2n -1 with A 0 = 0. Note that X 1 -X 0 = (x 1 -m) 2 (1 + ϵ 1 ) 2 (1 + η 1 ) -(x 1 -m) 2 and for all 2 ≤ k ≤ n X 2k-1 = X 2k-2 + ŝ2k-2 θ k , X 2k-2 = X 2k-3 + (x k -m) 2 (1 + ϵ k ) 2 (1 + η k ) -(x k -m) 2 .
(5.6) By construction and from Theorem 13, A k is a predictable process and

M k = X k -A k forms a martingale for all 0 ≤ k ≤ 2n -1. Moreover, H(X, δ) -H(X, 0) = X 2n-1 = M 2n-1 + A 2n-1 = M 2n-1 + 2n-1 k=1 E[X k -X k-1 /F k-1 ] = M 2n-1 + E[X 2n-1 /F 2n-2 ].
Finally, H(X, δ) -H(X, 0) = M + A, (5.7) where

M = X 2n-1 -A 2n-1 , A = E[X 2n-1 /F 2n-2 ].

. Bias Analysis

In the following we compute the bias of the two-pass-variance algorithm using the Equation (5.7).

Theorem 20. For the two-pass-variance algorithm, we have

E (H(X, δ) -H(X, 0)) = 1 n V ( s) + O(u 2 )
,

where 1 n s = m.
Proof. From Theorem 13, we know that E(M ) = 0. Then

E (H(X, δ) -H(X, 0)) = E(M + A) = E(M ) + E(A) = E(A).
Since X 0 = 0 we have

E(A) = E 2n-1 k=1 E[X k -X k-1 /F k-1 ] = E (E[X 1 /F 0 ]) + n k=2 E (E[X 2k-1 -X 2k-2 /F 2k-2 ]) + E (E[X 2k-2 -X 2k-3 /F 2k-3 ]) = E (E[X 1 /F 0 ]) + n k=2 E (E[θ k ŝ2k-2 /F 2k-2 ]) + E (E[X 2k-2 -X 2k-3 /F 2k-3 ]) .
Therefore, by mean independence property and law of total expectation, we have

E (E[X 1 /F 0 ]) = E E[(x 1 -m) 2 (1 + ϵ 1 ) 2 (1 + η 1 )/F 0 ] -(x 1 -m) 2 = E (x 1 -m) 2 (1 + ϵ 1 ) 2 E[(1 + η 1 )/F 0 ] -(x 1 -m) 2 = E (x 1 -m) 2 (1 + ϵ 1 ) 2 -(x 1 -m) 2 = E (x 1 -m) 2 (1 + ϵ 2 1 ) -(x 1 -m) 2 .
While the mean independence property implies

n k=2 E (E[θ k ŝ2k-2 /F 2k-2 ]) = n k=2 E (ŝ 2k-2 E[θ k /F 2k-2 ]) = 0.
For all 2 ≤ k ≤ n, we have

E (E[X 2k-2 -X 2k-3 /F 2k-3 ]) =E E[(x k -m) 2 (1 + ϵ k ) 2 (1 + η k )/F 2k-3 ] -(x k -m) 2 =E (x k -m) 2 (1 + ϵ 2 k ) -(x k -m) 2 .
Because E( m) = m, we have

E(A) = n k=1 E (x k -m) 2 (1 + ϵ 2 k ) -(x k -m) 2 = n k=1 E (x k -m) 2 -(x k -m) 2 + E (x k -m) 2 ϵ 2 k = n k=1 V ( m) + E (x k -m) 2 ϵ 2 k = nV ( m) + n k=1 E (x k -m) 2 ϵ 2 k . Since m = 1 n (1 + δ n ) s and ϵ 2 k ≤ u 2 for all 1 ≤ k ≤ n, n k=1 E (x k -m) 2 ϵ 2 k = O(u 2 ).
Moreover, δ 2 n ≤ u 2 , then

V ( m) = 1 n 2 V ( s(1 + δ n )) = 1 n 2 E s 2 (1 + δ n ) 2 - 1 n 2 E ( s(1 + δ n )) 2 = 1 n 2 E s 2 (1 + δ 2 n ) - 1 n 2 E ( s) 2 by Lemma 2 = 1 n 2 V ( s) + 1 n 2 E s 2 δ 2 n = 1 n 2 V ( s) + O(u 2 ) Therefore E(A) = 1 n V ( s) + O(u 2 ). Finally, E (H(X, δ) -H(X, 0)) = E(A) = 1 n V ( s) + O(u 2 ).
Interestingly, the textbook-variance and two-pass-variance algorithms under SR have an opposed bias at the first order over u.

Remark 11. Lemma 11 shows that

V ( m) = V   1 n n i=1 x i n+1 k=max(2,i) (1 + δ k-1 )   ≤ 1 n 2 ∥x∥ 2 1 γ n (u 2 ).
Therefore,

E(H(X, δ)) = H(X, 0) + nV ( m) + n k=1 E (x k -m) 2 ϵ 2 k ≤ H(X, 0) + nV ( m) + u 2 n k=1 E (x k -m) 2 = H(X, 0) + nV ( m) + u 2 H(X, 0) + n k=1 V (x k -m) = H(X, 0) + nV ( m) + u 2 (H(X, 0) + nV ( m)) ≤ (1 + u 2 ) H(X, 0) + 1 n ∥x∥ 2 1 γ n (u 2 ) = (1 + u 2 )H(X, 0) 1 + K 2 1 γ n (u 2 ) .
where

K 1 = ∥x∥ 1 √ nH(X,0)
.

We now turn to bound the forward error of this algorithm. We use the same three methods used for the textbook-variance algorithm to establish probabilistic bounds on the error: DM, BC and AH methods.

. DM Method

This sub-section uses DM method to provide a probabilistic bound on the forward error of the two-pass-variance algorithm under SR-nearness. We have demonstrated in Equation (5.7) that

H(X, δ) -H(X, 0) = M + A, where M = M 2n-1 = X 2n-1 -A 2n-1 , A = A 2n-1 = E[X 2n-1 /F 2n-2 ].
i -M i-1 | ≤ C i u, for all 1 ≤ i ≤ 2n -1, where C 2k-1 = k i=1 (x i -m) 2 (1 + u) k+1 , for all 1 ≤ k ≤ n C 2k-2 = (x k -m) 2 (1 + ϵ k ) 2 , for all 2 ≤ k ≤ n. Proof. Firstly, |M 2k-1 -M 2k-2 | = |X 2k-1 -X 2k-2 -(A 2k-1 -A 2k-2 )| = |X 2k-1 -X 2k-2 -E[X 2k-1 -X 2k-2 /F 2k-2 ]| = |ŝ 2k-2 θ k -E[ŝ 2k-2 θ k /F 2k-2 ]| = |ŝ 2k-2 θ k -ŝ2k-2 E[θ k /F 2k-2 ]| = |ŝ 2k-2 θ k | by Lemma 2 ≤ u k i=1 (x i -m) 2 (1 + u) k+1 = uC 2k-1 . Secondly, |M 2k-2 -M 2k-3 | = |X 2k-2 -X 2k-3 -(A 2k-2 -A 2k-3 )| = |X 2k-2 -X 2k-3 -E[X 2k-2 -X 2k-3 /F 2k-3 ]| = (x k -m) 2 (1 + ϵ k ) 2 (1 + η k ) -E[(x k -m) 2 (1 + ϵ k ) 2 (1 + η k )/F 2k-3 ] = (x k -m) 2 (1 + ϵ k ) 2 (1 + η k ) -(x k -m) 2 (1 + ϵ k ) 2 E[(1 + η k )/F 2k-3 ] = (x k -m) 2 (1 + ϵ k ) 2 (1 + η k ) -(x k -m) 2 (1 + ϵ k ) 2 by Lemma 2 = (x k -m) 2 (1 + ϵ k ) 2 η k ≤ u(x k -m) 2 (1 + ϵ k ) 2 = uC 2k-2 . Theorem 21. The martingale M 0 , • • • , M 2n-1 from System (5.8) satisfies |M | ≤ (1 + u) 2 uγ 2n (u) ln(2/λ) H(X, 0) + ∥x∥ 2 1 n γ n (u) 2 ,
(5.9)

with probability at least 1 -λ.

Proof. Lemma 16 implies that |M k -M k-1 | ≤ C k u.
Then, using Azuma-Hoeffding inequality yields

P   |M | ≤ u 2 2n-1 k=1 C 2 k 2 ln(2/λ)   ≥ 1 -λ. Partition 2n-1 k=1 C 2 k = n k=1 C 2 2k-1 + n k=2 C 2 2k-2 = C 2 1 + n k=2 C 2 2k-1 + n k=2 C 2 2k-2 .
We have

C 2 1 = (x 1 -m) 4 (1 + u) 4 and n k=2 C 2 2k-1 = n k=2 k i=1 (x i -m) 2 2 (1 + u) 2(k+1) ≤ n i=1 (x i -m) 2 2 n k=2 (1 + u) 2(k+1) .
We also have

n k=2 C 2 2k-2 = n k=2 (x k -m) 4 (1 + ϵ k ) 4 ≤ (1 + u) 4 n k=2 (x k -m) 4 . Furthermore, since n i=1 (x i -m) 4 ≤ n i=1 (x i -m) 2 2 we have 2n-1 k=1 C 2 k = n k=2 C 2 2k-1 + C 2 1 + n k=2 C 2 2k-2 ≤ n i=1 (x i -m) 2 2 n k=2 (1 + u) 2(k+1) + (x 1 -m) 4 (1 + u) 4 + (1 + u) 4 n i=2 (x i -m) 4 = n i=1 (x i -m) 2 2 n k=2 (1 + u) 2(k+1) + (1 + u) 4 n i=1 (x i -m) 4 ≤ n i=1 (x i -m) 2 2 n k=1 (1 + u) 2(k+1) = n i=1 (x i -m) 2 2 (1 + u) 4 (1 + u) 2n -1 (1 + u) 2 -1 = n i=1 (x i -m) 2 2 (1 + u) 4 γ 2n (u) u 2 + 2u . It follows that u 2 2n-1 k=1 C 2 k ≤ u 2 n i=1 (x i -m) 2 2 (1 + u) 4 γ 2n (u) u 2 + 2u = n i=1 (x i -m) 2 (1 + u) 2 u 2 γ 2n (u) u 2 + 2u ≤ n i=1 (x i -m) 2 (1 + u) 2 uγ 2n (u) 2 because u 2 + u ≤ u 2 . Since (x k -m) = (x k -m) + (m -m) and n k=1 (x k -m) = 0, we have n k=1 (x k -m) 2 = n k=1 (x k -m) 2 + n(m -m) 2 = H(X, 0) + n(m -m) 2 .
Note that

| m -m| = 1 n n i=1 x i   n+1 k=max(2,i) (1 + δ k-1 ) -1   ≤ ∥x∥ 1 n γ n (u). Then, n k=1 (x k -m) 2 ≤ H(X, 0) + ∥x∥ 2 1 n γ n (u) 2 .
(5.10) Finally,

|M | ≤ (1 + u) 2 uγ 2n (u) ln(2/λ) H(X, 0) + ∥x∥ 2 1 n γ n (u) 2 ,
with probability at least 1 -λ.

Let us focus now on A.

Theorem 22. The drift A satisfies

|A| ≤ (1 + u) 2 ∥x∥ 2 1 n γ n (u) 2 + H(X, 0)γ 2 (u).
(5.11)

Proof. Since E[θ k /F 2k-2 ] = 0 for all 1 ≤ k ≤ n, and E[η k /F 2k-3 ] = 0 for all 2 ≤ k ≤ n, we have A = 2n-1 k=1 E[X k -X k-1 /F k-1 ] = n k=1 E[X 2k-1 -X 2k-2 /F 2k-2 ] + n k=2 E[X 2k-2 -X 2k-3 /F 2k-3 ] = n k=1 E[ŝ 2k-2 θ k /F 2k-2 ] + n k=2 E[(x k -m) 2 (1 + ϵ k ) 2 (1 + η k )/F 2k-3 ] -(x k -m) 2 = n k=1 ŝ2k-2 E[θ k /F 2k-2 ] + n k=2 (x k -m) 2 (1 + ϵ k ) 2 E[(1 + η k )/F 2k-3 ] -(x k -m) 2 = n k=1 (x k -m) 2 (1 + ϵ k ) 2 -(x k -m) 2 .
It follows that

|A| = n k=1 (x k -m) 2 (1 + ϵ k ) 2 -(x k -m) 2 ≤ n k=1 (x k -m) 2 -(x k -m) 2 + n k=1 (x k -m) 2 (2ϵ k + ϵ 2 k ) ≤ n k=1 (x k -m) 2 -(x k -m) 2 + (2u + u 2 ) n k=1 (x k -m) 2 .
The inequality (5.10) implies that n

k=1 (x k -m) 2 ≤ H(X, 0) + ∥x∥ 2 1 n γ n (u) 2 , then n k=1 (x k -m) 2 -(x k -m) 2 ≤ ∥x∥ 2 1 n γ n (u) 2 .
Theorem 24. For all 0 < λ < 1, the computed H(X, δ) satisfies under SRnearness

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤(1 + u) 4γ n+1 (u 2 ) λ + 4γ n+1 (u 2 ) λ 2K 1 + K 2 1 4γ n+1 (u 2 ) λ + 1 + u, with probability at least 1 -λ. Where K 1 = ∥x∥ 1 √ nH(X,0)
. Proof. Equation (5.12) states that

|H(X, δ) -H(X, 0)| ≤ (1 + u)C + u |H(X, 0)| . Note that | n i=1 φ i | ≤ | n i=1 (φ i -1)| + n. The following quantities        n i=1 (x i -m) 2 (φ i -1) , | m -m| , | n i=1 (x i -m)(φ i -1)| , | n i=1 (φ i -1)| represent the absolute errors of the inner product n i=1 (x i -m) 2 , the average m = 1 n n i=1 x i , the summations s = n i=1 (x i -m) and n i=1 1 respectively. Then [23, sec 5.1] proves that n i=1 (x i -m) 2 (φ i -1) ≤ |H(X, 0)| 4γ n+1 (u 2 ) λ with probability at least 1 - λ 4 , | m -m| ≤ 1 n ∥x∥ 1 4γ n (u 2 ) λ with probability at least 1 - λ 4 , n i=1 (x i -m)(φ i -1) ≤ n i=1 |x i -m| 4γ n+1 (u 2 ) λ with probability at least 1 - λ 4 , n i=1 (φ i -1) ≤ n 4γ n+1 (u 2 ) λ with probability at least 1 - λ 4 .
Using Cauchy-Schwarz inequality, we obtain

n i=1 |x i -m| ≤ n n i=1 (x i -m) 2 = nH(X, 0). Since γ n (u 2 ) ≤ γ n+1 (u 2 ), Lemma 14 implies C ≤ |H(X, 0)| 4γ n+1 (u 2 ) λ + 2 ∥x∥ 1 n 4γ n+1 (u 2 ) λ nH(X, 0) + ∥x∥ 2 1 n 4γ n+1 (u 2 ) λ 4γ n+1 (u 2 ) λ + 1 = |H(X, 0)| 4γ n+1 (u 2 ) λ + 4γ n+1 (u 2 ) λ 2 |H(X, 0)| ∥x∥ 1 nH(X, 0) + ∥x∥ 2 1 n 4γ n+1 (u 2 ) λ + 1 , with probability at least 1 -λ. Finally |H(X, δ) -H(X, 0)| |H(X, 0)| ≤(1 + u) 4γ n+1 (u 2 ) λ + 4γ n+1 (u 2 ) λ 2K 1 + K 2 1 4γ n+1 (u 2 ) λ + 1 + u,
with probability at least 1 -λ.

. AH Method

Theorem 25. For all 0 < λ < 1, the computed H(X, δ) satisfies under SRnearness

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤(1 + u) uγ 2(n+1) (u) ln(8/λ) + uγ 2(n+1) (u) ln(8/λ) 2K 1 + K 2 1 uγ 2(n+1) (u) ln(8/λ) + 1 + u, with probability at least 1 -λ. Where K 1 = ∥x∥ 1 √ nH(X,0)
. Proof. Equation (5.12) states that and [46, cor 4.7] shows that each of the following inequalities holds with probability at least 1 -λ 4 :

|H(X, δ) -H(X, 0)| ≤ (1 + u)C + u |H(X, 0)| . We have | n i=1 φ i | ≤ | n i=1 (φ i -1)| + n,
n i=1 (x i -m) 2 (φ i -1) ≤ |H(X, 0)| uγ 2(n+1) (u) ln(8/λ), | m -m| ≤ 1 n ∥x∥ 1 uγ 2n (u) ln(8/λ), n i=1 (x i -m)(φ i -1) ≤ n i=1 |x i -m| uγ 2(n+1) (u) ln(8/λ), n i=1 (φ i -1) ≤ n uγ 2(n+1) (u) ln(8/λ).
By Cauchy-Schwarz inequality,

n i=1 |x i -m| ≤ n n i=1 (x i -m) 2 = nH(X, 0). Since γ 2n (u) ≤ γ 2(n+1) (u), Lemma 14 implies C ≤ |H(X, 0)| uγ 2(n+1) (u) ln(8/λ) + 2 ∥x∥ 1 n uγ 2(n+1) (u) ln(8/λ) nH(X, 0) + ∥x∥ 2 1 n 2 uγ 2(n+1) (u) ln(8/λ) n uγ 2(n+1) (u) ln(8/λ) + n = |H(X, 0)| uγ 2(n+1) (u) ln(8/λ) + uγ 2(n+1) (u) ln(8/λ) 2 |H(X, 0)| ∥x∥ 1 nH(X, 0) + ∥x∥ 2 1 n uγ 2(n+1) (u) ln(8/λ) + 1 , with probability at least 1 -λ. Finally |H(X, δ) -H(X, 0)| |H(X, 0)| ≤(1 + u) uγ 2(n+1) (u) ln(8/λ) + uγ 2(n+1) (u) ln(8/λ) 2K 1 + K 2 1 uγ 2(n+1) (u) ln(8/λ) + 1 + u,
with probability at least 1 -λ.

. Pairwise Textbook-Variance and Pairwise Two-pass-Variance

In this section, we illustrate the continued applicability of SR results on the forward error of the pairwise summation to the forward error of both two-pass-variance and textbook-variance algorithms. We refer to "pairwise two-pass-variance," the algorithm that computes the variance using a pairwise method to compute the summation and the average m, and "pairwise textbook-variance," the algorithm that computes the variance using a pairwise method to compute the summations. The following theorem derives a probabilistic bound for the pairwise textbook-variance using the BC method.

Theorem 26. For the pairwise textbook-variance, for all 0 < λ < 1, the computed H(X, δ) satisfies under SR-nearness

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤ K 2 2 2γ log(n)+1 (u 2 )/λ+K 2 1 (1 + u) 3 2γ log(n) (u 2 )/λ + 1 2 -1 , with probability at least 1 -λ, where K 1 = ∥x∥ 1 √ nH(X,0) and K 2 = ∥x∥ 2 √ H(X,0) . Proof. Equation (5.4) states that |H(X, δ) -H(X, 0)| ≤ n i=1 x 2 i (ψ i -1) + 1 n B,
where B = ( s -s) 2 ψ n+1 +2 |s( s -s)ψ n+1 |+ s 2 (ψ n+1 -1) . Moreover, using a pairwise method to compute s and n i=1 x 2 i , Sub-section 4.2.4 implies

   ŝ = n i=1 x i log 2 (n) j=1 (1 + δ j ⌈ i 2 j ⌉ ) ψ i = (1 + ϵ i ) j∈K i (1 + η j ) for all 1 ≤ i ≤ n,
where the cardinality of |K i | = log 2 (n) for all 1 ≤ i ≤ n. Note that the square s 2 , the division s 2 n , and the subtraction 2 h i=1 x 2 i -1 n s 2 are computed in the standard case (without a pairwise method). Then, ψ n+1 ≤ (1 + u) 3 as in Equation (5.4). Lemma 11 implies n i=1

x 2 i (ψ i -1) ≤ ∥x∥ 2 2 2γ log 2 (n)+1 (u 2 )/λ with probability at least 1 - λ 2 , | s -s| ≤ ∥x∥ 1 2γ log 2 (n) (u 2 )/λ with probability at least 1 - λ 2 .
Since, |ψ n+1 | ≤ (1 + u) 3 and |s| ≤ ∥x∥ 1 , we have

B ≤ (1 + u) 3 ∥x∥ 2 1 2γ log(n) (u 2 )/λ + 2 2γ log(n) (u 2 )/λ + ∥x∥ 2 1 (1 + u) 3 -1 = (1 + u) 3 ∥x∥ 2 1 2γ log(n) (u 2 )/λ + 2 2γ log(n) (u 2 )/λ + 1 -∥x∥ 2 1 = (1 + u) 3 ∥x∥ 2 1 2γ log(n) (u 2 )/λ + 1 2 -∥x∥ 2 1 .
with probability at least 1 -λ 2 . Finally, Lemma 14 implies

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤ 1 |H(X, 0)| n i=1 x 2 i (ψ i -1) + 1 n |H(X, 0)| B ≤ K 2 2 2γ log(n)+1 (u 2 )/λ + K 2 1 (1 + u) 3 2γ log(n) (u 2 )/λ + 1 2 -1 ,
with probability at least 1 -λ.

The following theorem shows the probabilistic bound for the pairwise textbookvariance using the AH method. Theorem 27. For the pairwise textbook-variance, for all 0 < λ < 1, the computed H(X, δ) satisfies under SR-nearness

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤K 2 2 uγ 2(log(n)+1) (u) ln(4/λ) + K 2 1 (1 + u) 3 uγ 2 log(n) (u) ln(4/λ) + 1 2 -1 ,
with probability at least 1 -λ, where

K 1 = ∥x∥ 1 √ nH(X,0) and K 2 = ∥x∥ 2 √ H(X,0) . 110 Proof. Equation (5.4) states that |H(X, δ) -H(X, 0)| ≤ n i=1 x 2 i (ψ i -1) + 1 n B,
where B = ( s -s) 2 ψ n+1 +2 |s( s -s)ψ n+1 |+ s 2 (ψ n+1 -1) . Moreover, using a pairwise method to compute s and n i=1 x 2 i , Sub-section 4.2.4 implies

   ŝ = n i=1 x i log 2 (n) j=1 (1 + δ j ⌈ i 2 j ⌉ ) ψ i = (1 + ϵ i ) j∈K i (1 + η j ) for all 1 ≤ i ≤ n,
where the cardinality of |K i | = log 2 (n) for all 1 ≤ i ≤ n. Note that the square s 2 , the division s 2 n , and the subtraction 2 h i=1 x 2 i -1 n s 2 are computed in the standard case (without a pairwise method). Then, ψ n+1 ≤ (1 + u) 3 as in Equation (5.4). Sub-section 4.1.3 shows

n i=1 x 2 i (ψ i -1) ≤ ∥x∥ 2 2 uγ 2(log(n)+1) (u) ln(4/λ) with probability at least 1 - λ 2 , | s -s| ≤ ∥x∥ 1 uγ 2 log(n) (u) ln(4/λ) with probability at least 1 - λ 2 .
As the previous proof, we can show that

B ≤ (1 + u) 3 ∥x∥ 2 1 uγ 2 log(n) (u) ln(4/λ) + 1 2 -∥x∥ 2 1 ,
with probability at least 1 -λ 2 , Finally, Lemma 14 implies

|H(X, δ) -H(X, 0)| |H(X, 0)| ≤K 2 2 uγ 2(log(n)+1) (u) ln(4/λ) + K 2 1 (1 + u) 3 uγ 2 log(n) (u) ln(4/λ) + 1 2 -1 ,
with probability at least 1 -λ.

Similar bounds are reached for the pairwise two-pass-variance algorithm using the same methods.

. Error Bound Analysis

Table 5.1 shows the asymptotic forward error bounds for the textbookvariance. Higher order terms in u have been dropped when nu ≪ 1 and uniquely for the BC when nu ≫ 1 and nu 2 ≪ 1, and only dominant terms are shown. For H(X, 0

) = n i=1 x 2 i -1 n s 2 , K 1 = ∥x∥ 1 √ nH(X,0)
, and

K 2 = ∥x∥ 2 √ H(X,0)
, we have The results in the table are based on:

Textbook-Variance Algorithm

We present a numerical application of the textbook-variance algorithm for floating-point numbers chosen uniformly at random between 0 and 1. As expected, for a fixed n, figure 5.2 shows that the three bounds are close for a probability around 0.9. Nevertheless, AH and DM bounds are more accurate for higher probabilities than BC bound. The result is unsurprising because, generally, Azuma-Hoeffding inequality provides a bound for the deviation of the sum of a sequence of independent and bounded random variables, martingales in this instance, which gives tighter bounds for higher probabilities. In contrast, Bienaymé-Chebyshev inequality is a less restrictive result that provides an upper bound for the probability of deviation between the mean of a distribution and a particular value. The two-pass-variance algorithm exhibits analogous boundary behavior.

Textbook-Variance Against Two-pass-Variance

We now compare the forward errors of both algorithms under SR through two experiences.

0.008 0.016 0.031 0.062 0.125 0.250 0.500 1.000 2.000 4.000 8.000 n In figure 5.3, when the floating-point numbers are randomly chosen with zero mean distribution, the absorption errors cancel each other out because both positive and negative errors are uniformly distributed. Therefore, the computed mean is close to zero with low absolute error, and the two-passvariance algorithm degenerates into the textbook-variance algorithm. Interestingly, this effect is captured by the theoretical bounds because the condition term K Many computations are non-linear in various fields, such as numerical analysis. In this chapter, we have proposed a general framework that allows us to analyze problems under SR. We have demonstrated how an algorithm can be expressed as a martingale plus a drift. On the one hand, our generalization shows that the drift is zero for algorithms with multi-linear error, and the error can be described by a martingale for which we use AH or BC methods to bound it. On the other hand, for algorithms with non-linear error, the drift is non-zero. Regardless of the problem under study, we have proposed an expression for the bias and demonstrated that it is negligible at the first order over u.

In 1983, Chan, Golub, and LeVeque investigated the forward error of variance computation algorithms using RN. To the best of our knowledge, this is the first theoretical study of this problem using SR as well as of any algorithm with non-linear errors. In this chapter, we have investigated two variance computation algorithms that exhibit non-linear errors under SR. The study demonstrates that they are biased and using SR results in probabilistic

-Industrial Example

In the previous two chapters, we have demonstrated that SR has several advantages and positive effects on computations. For algorithms with multilinear error, SR is unbiased. Moreover, Chapter 5 demonstrates that when using SR, the forward error has a probabilistic bound in O( √ nu) instead of O(nu). To the best of our knowledge, previous SR theoretical studies have only considered algorithms with elementary operations, such as additions/subtractions or multiplications, while no study has explored the division.

In this chapter, we examine the division case with SR through two examples. First, an industrial application proposed by Anthony Scemama ("Laboratoire de Chimie et Physique Quantiques (LCPQ)") from the Quantum Package [START_REF] Garniron | Quantum package 2.0: An open-source determinant-driven suite of programs[END_REF], an open-source program designed for quantum chemistry. This example, MP2, involves calculating the Møller-Plesset Perturbation Theory, a method used to estimate the correlation energy of molecules.

The MP2 example is of great interest since it computes a summation of divisions such that each division involves a denominator with three operations. We have shown that their error impact is negligible compared to the summation errors under SR, and the probabilistic error bound is proportional to √ nu.

Table 6.1 illustrates the previous results and shows the SR benefit through an industrial code with two applications. Second, in the MP2 example, there are only three elementary operations involved in the denominator in each division. To analyze the SR effect in the division computation, we have examined a more complex problem, the inverse of a summation of n numbers. We have introduced a computational model for computing the error of the division in the context of SR. Using BC or AH methods, we have demonstrated probabilistic error bounds in O( √ nu). where

. Møller-Plesset Perturbation Theory

ψ ijab = k∈K ijab (1 + α k )(1 + β k )(1 + δ k )(1 + η k )(1 + µ k )
and K ijab is a subset of N that depends on n virt and n occ . For all k ∈ K ijab , α k ,β k , δ k , η k , and µ k represent the rounding errors from the subtractions a -b, the multiplications by a, the multiplications by 2, the divisions by εijab , and the summations, respectively. Using the BC method, the following theorem computes a probabilistic bound of the forward error of the computation of E.

Theorem 28. Denote N = n virt n occ . For all 0 < λ < 1, the computed Ê satisfies under SR-nearness = ϵ ijab + (ϵ i + ϵ j )θ 1 + ϵ i + ϵ j -ϵ a θ 2 + (ϵ i + ϵ j )θ 1 θ 2 (1 + θ 3 ) =ϵ ijab + (ϵ i + ϵ j )θ 1 + ϵ i + ϵ j -ϵ a θ 2 + (ϵ i + ϵ j )θ 1 θ 2 + ϵ ijab θ 3 + (ϵ i + ϵ j )θ 1 θ 3

Ê -E E ≤ K (1 + u) N f (u) 1 -f (u) + γ N (u 2 ) λ , (6.3
+ ϵ i + ϵ j -ϵ a θ 2 θ 3 + (ϵ i + ϵ j )θ 1 θ 2 θ 3 =ϵ ijab 1 + ϵ i + ϵ j ϵ ijab θ 1 + ϵ i + ϵ j -ϵ a ϵ ijab θ 2 + θ 3 + ϵ i + ϵ j ϵ ijab θ 1 θ 2 + ϵ i + ϵ j ϵ ijab θ 1 θ 3 + ϵ i + ϵ j -ϵ a ϵ ijab θ 2 θ 3 + ϵ i + ϵ j ϵ ijab θ 1 θ 2 θ 3 =ϵ ijab (1 + χ ijab ),
where

χ ijab = ϵ i + ϵ j ϵ ijab θ 1 + ϵ i + ϵ j -ϵ a ϵ ijab θ 2 + θ 3 + ϵ i + ϵ j ϵ ijab θ 1 θ 2 + ϵ i + ϵ j ϵ ijab θ 1 θ 3 + ϵ i + ϵ j -ϵ a ϵ ijab θ 2 θ 3 + ϵ i + ϵ j ϵ ijab θ 1 θ 2 θ 3 .
Since the rounding errors satisfy |θ k | ≤ u for k = 1, 2 and 3, and the quotients ϵ i +ϵ j ϵ ijab and ϵ i +ϵ j -ϵa ϵ ijab are between -1 and 1, we have

-f (u) ≤ χ ijab ≤ f (u), (6.4)
with f (u) = 3u + 3u 2 + u 3 = 3u + O(u 2 ). It follows that

1 1 + f (u) ≤ 1 1 + χ ijab ≤ 1 1 -f (u)
.

Ê -E E ≤ K(1 + u) N f (u) 1 -f (u) + K γ N (u 2 ) λ = K (1 + u) N f (u) 1 -f (u) + γ N (u 2 ) λ .
Remark 12. When N u ≪ 1,

• (1 + u) N = 1 + N u + O(u 2 ). • γ N (u 2 ) = √ N u + O(u 2 ).
• f (u) = 3u + O(u 2 ).

• In conclusion, the forward error of E has a probabilistic bound proportional to √ N u. Note that the deterministic bound of this forward error is given by

A = K (1 + u) N 1 -f (u) -1 .
Remark 12 implies

(1 + u) N 1 -f (u) -1 = (1 + N u)(1 + 3u) -1 + O(u 2 ) = (3 + N )u + O(u 2 ).
Denote the probabilistic bound of the forward error of E by

B = K (1 + u) N f (u) 1 -f (u) + γ N (u 2 ) λ .
In the following, we compare the number of significant digits of RN and SR for the computation of E in single and double precision through two input datasets: "benzene-dz" and "benzene-tz".

Note that

• The red numbers correspond to the digits lost compared to the reference.

121

• -log(A) represents the minimum of significant digits possible with RN.

• -log(B) represents the minimum of significant digits possible with SR.

In binary-32, SR is more favorable as it guarantees at least three significant digits even in the worst-case scenario, whereas RN can potentially lose all significance. These findings align with previous research on SR, emphasizing its advantage in low-precision. However, In binary-64, SR and RN are comparable.

Interestingly, we have demonstrated in this section that SR ensures a probabilistic error bound of O( √ nu) in the division case. Since each term in the denominator involves three elementary operations, we have shown that its effect is negligible compared to the overall summation. In the next section, we will examine the computation of the inverse of a summation, where the denominator involves a large number of operations.

. Inverse of Summation

Let s = n i=1 x i and y = 1 The Taylor series of 1 1+Φ around 0 yields

Ψ = 1 + θ 1 + Φ -1 = θ -Φ 1 + Φ = (θ -Φ)(1 -Φ + Φ 2 ) + O(Φ 3 ) = θ -θΦ + θΦ 2 -Φ + Φ 2 + O(Φ 3 )
= θ -Φ(θ + 1) + Φ 2 (θ + 1) + O(Φ 3 ) = θ + Φ(θ + 1)(Φ -1) + O(Φ 3 ).

Note that this development fails to hold when s is near 0 and then the error Φ will be away from 0. The following theorem computes a probabilistic bound of the relative error Ψ using the BC method and when nu 2 ≪ 1. The following theorem compute a probabilistic bound of the relative error Ψ using the AH method and when nu ≪ 1. for the error of any algorithm, but calculations can sometimes prove complex, as demonstrated for the two-pass-variance algorithm in Section 5.3.

We have shown that this generalization is highly flexible: the drift is zero for problems exhibiting multi-linear error, a characteristic applicable to all algorithms analyzed in Chapter 4. Subsequently, we use the AH or BC method to establish probabilistic bounds for the martingale term. For the textbookvariance algorithm in Section 5.2, a part of the error is similar to an inner product, which directly forms a martingale. Consequently, we apply the generalization to address the remaining part of the error. For the two-pass-variance algorithm in Section 5.3, we have considered the errors in the computation of m as inputs and formed the stochastic process from the remaining errors.

In general, once the stochastic process is determined, the applicability of this generalization is direct. We have demonstrated that the drift is negligible at the first order over u for the two variance computation algorithms. While for the martingale term, we have used the AH method to obtain probabilistic bounds proportional to √ nu. We also have shown that the previous results under SR remain valid in the pairwise case.

We have extended the properties of SR to the division problem in Chapter 6. We have demonstrated probabilistic error bounds in O( √ nu) for the MP2 computation and the summation inversion. This is of great interest because it illustrates the advantage of SR in all elementary operations.

. Perspectives

Artificial intelligence has motivated the use of SR in various applications, and the theoretical analysis in this thesis has confirmed the interest of using SR in computations. For an algorithm with a large problem size, using SR instead of RN improves the computation accuracy of this algorithm. However, it is known that SR is relatively expensive [START_REF] Croci | Stochastic rounding: implementation, error analysis and applications[END_REF]. This thesis did not address the implementation challenges, particularly in assessing the trade-off between the accuracy gained and the potential energy/time overhead when opting for SR over RN.

Moreover, the numerical experiments in Chapter 4 raise concerns regarding the use of SR as a model to estimate RN rounding errors [START_REF] Denis | Verificarlo: Checking floating point accuracy through Monte Carlo arithmetic[END_REF][START_REF] Sohier | Confidence intervals for stochastic arithmetic[END_REF], in particular for the inner product and a large number of operations. Further studies are required to assess precisely the limits of this model and possibly give criteria to detect them. Furthermore, in the context of RN, the necessary assumptions on errors to validate Wilkinson's intuition require exploration.

In this thesis, we have demonstrated probabilistic error bounds in O( √ nu)

through two methods based on a common tool: concentration inequalities.

The AH method uses Azuma-Hoeffding inequality, which is particularly effective when nu ≪ 1 and for higher probabilities. The BC method uses Bien-126 aymé-Chebyshev inequality and a bound on the variance, which is useful for large problem sizes n, particularly when nu 2 ≪ 1. Further research is needed to establish a connection between AH and BC methods, possibly through a unifying tool such as the Central Limit Theorem for instance. The idea is to ensure tight probabilistic bounds for higher probabilities and large problem sizes n.

Furthermore, using SR in low-precision formats, especially bfloat-16 is becoming increasingly attractive due to its higher speed and lower energy consumption. The BC method ensures tight probabilistic error bound when nu 2 ≪ 1. However, in low-precision, the first n values greater than u 2 can be relatively small, which can make the bound less significant in these formats when n is large and the condition nu 2 ≪ 1 does not hold. Therefore, further studies are required to examine other properties, such as higher-order moments, to ensure more general probabilistic bounds.

The generalization in Section 5.1 decomposes the error of an algorithm into a martingale plus a drift. For the martingale term, the concentration inequalities allow us to obtain probabilistic bounds in O( √ nu). Nevertheless, we are not quite there yet regarding the drift part. For the variance computation, we have demonstrated that at the first order over u, it is negligible compared to the martingale term. Furthermore, by construction, the drift corresponds to the summation biases accumulated in the computation. Since the martingale term captures the compensated errors (the errors of degree one) under SR, and the remaining errors (the errors of degree higher than one) are in the drift term, we believe that the drift is almost negligible compared to the martingale term at the first order. The results of this thesis have considered algorithms based on elementary operations and a fixed number of iterations. One numerical scheme is to investigate other situations, such as conditional statements, algorithms with while loop, and iterative methods such as Newton's method, which is a substantial concern, necessitating a thorough analysis under SR. We expect that the stopping condition can be modeled using a stopping time, which allows us to recover all the probabilistic properties used in this thesis. The challenge is to verify whether previous properties proved for SR hold in these situations. 127

Titre:

  Modèles stochastiques pour l'évaluation des erreurs numériques Mots clés: Arrondi Stochastique, Arithmétique en virgule flottante, Erreurs d'arrondi, Théorie des martingales Résumé: L'idée de considérer les erreurs d'arrondi comme des variables aléatoires n'est pas nouvelle. Basées sur des outils tels que l'indépendance des variables aléatoires ou le théorème central limite, plusieurs propositions ont démontré des bornes d'erreur en O( √ n).

tions 5 . 2 ,

 52 and 5.3 also show probabilistic error bounds in O( √ nu) for both algorithms. Moreover, Section 5.4 demonstrates that SR results extend to the pairwise case, with probabilistic error bounds in O( log(n)u).

For a given

  basis β and a working precision p, a floating-point number (FP) is a real number x characterized by (m, e) such that x = m × β e-p , (2.1) where • m is an integer (the significand) such that β p-1 ≤ |m| < β p .

Figure 2 . 2 :

 22 Figure 2.2: SR-up-or-down.

Figure 2 . 3 :

 23 Figure 2.3: SR-nearness.

Proof.

  It suffices to consider quantities a and b resulting from the computation of k -1 scalar operations that have produced rounding errors δ 1 , • • • , δ k-1 . Consider now the computation of c ← a op b for any scalar operation op ∈ {+, -, * , /}, resulting in ĉ = fl(a op b) = (a op b)(1 + δ k ). The rounding error δ k = (ĉ -c)/c is a random variable that depends on δ 1 , • • • , δ k-1 and is given by δ k = (⌈c⌉ -c)/c with probability p = (c -⌊c⌋)/(⌈c⌉ -⌊c⌋), (⌊c⌋ -c)/c with probability 1 -p.

  0.150... 0.7500 -3.725290e-09 1.490116e-08 3 0.200... 0.2500 3.725290e-09 1.490116e-08 4 0.250... 0.6250 -3.725290e-09 2.980232e-08 5 0.300... 0.6250 -3.725290e-09 2.980232e-08 6 0.350... 0.6250 -3.725290e-09 2.980232e-08 7 0.400... 0.6250 -3.725290e-09 2.980232e-08 8 0.450... 0.6250 -3.725290e-09 2.980232e-08 9 0.500... 0.3125 1.117587e-08 5.960464e-08 10 0.550... 0.8125 -1.862645e-08 5.960464e-08 11 0.600... 0.8125 -1.862645e-08 5.960464e-08 12 0.650... 0.8125 -1.862645e-08 5.960464e-08 13 0.700... 0.8125 -1.862645e-08 5.960464e-08 14 0.749... 0.8125 -1.862645e-08 5.960464e-08 15 0.799... 0.8125 -1.862645e-08 5.960464e-08 16 0.849... 0.8125 -1.862645e-08 5.960464e-08 17 0.899... 0.8125 -1.862645e-08 5.960464e-08 18 0.949... 0.8125 -1.862645e-08 5.960464e-08 19 0.999... 0.8125 -1.862645e-08 5.960464e-08 Table 3.1: s k , p(s k ), bias and ϵ for N = 20 in single precision.

Figure 3 . 1 :

 31 Figure 3.1: Round to nearest (RN-binary32) vs stochastic rounding SR-up-ordown and SR-nearness for Listing 1.

Figure 3 .

 3 Figure 3.1 illustrates that SR-nearness mode is unbiased regardless of the number N of rectangles. The unbiased nature is unsurprising since the SRnearness mode is a sub-case of Monte Carlo Arithmetic (MCA). Stott Parker proves [60, p. 46] that the expectation of a sum of terms with MCA is the exact mathematical result.

π 2 0

 2 cos(t) dt = 1. We run one time the C program in Listing 2 with each of the two previously defined stochastic rounding modes.

Figure 3 . 2 :

 32 Figure 3.2: Round to nearest (RN-binary32) vs stochastic rounding SR-up-ordown and SR-nearness for Listing 2.

Figure 3 .

 3 Figure 3.2 shows that SR-nearness is unbiased, while SR-up-or-down is biased, and the bias magnitude increases with the number of rectangles N .

Figure 3 . 3 :

 33 Figure 3.3: Absolute value of the error over the number of rectangles used to integrate the cosine function between 0 and π/2 using the rectangular rule.

Figure 4 . 1 :

 41 Figure 4.1: (4.9) and (4.10) bounds with λ = 0.1, η = 0.05, δ = 0.05, κ = 1 and u = 2 -23 .

Figure 4 . 2 :

 42 Figure 4.2: Last operation in A.

Theorem 7 .

 7 For all 0 < λ < 1, the computed ŝ in the Equation (4.1) satisfies under SR-nearness E(ŝ) = s and|ŝ -s| |s| ≤ K γ n-1 (u 2 )/λ,(4.13)with probability at least 1 -λ, where K = n i=1 |a i | | n i=1 a i|is the condition number.

y 1 .

 1 = a ⊤ b, where a, b ∈ R n |ŷ -y| |y| ≤ Kγ n (u), (Det-IP) |ŷ -y| |y| ≤ Kγ n ( 2 ln (2n/λ)) with probability at least 1 -λ, (AH1-IP) |ŷ -y| |y| ≤ K uγ 2n (u) ln(2/λ) with probability at least 1 -λ, (AH2-IP) |ŷ -y| |y| ≤ K γ n (u 2 ) 1/λ with probability at least 1 -λ, (BC-IP) where K = n i=1 |a i b i | | n i=1 a i b i | is the condition number and γn ( 2 ln (2n/λ)) = exp 2n ln (2n/λ)u + nu 2 1 -u -Note that (AH1-IP) is the bound obtained by Connelly et al. [12, thm 4.8], (AH2-IP) is the bound obtained by Ipsen, and Zhou [46, cor 4.7], and (BC-IP) is the bound obtained in Theorem 8.

then γn ( 2

 2 ln (2n/λ)) ≈ exp 2n ln (2n/λ)u .

γ n (u 2 )

 2 ≈ exp (nu 2 ) -1 = √ nu + O(u 2 ).

(4. 22 )

 22 Therefore, the previous analysis and in particular (4.20), (4.21), and(4.22)

(4. 23 )

 23 Now, let us compare 1/λ and 2 ln(2/λ) for λ ∈]0; 1[.

Figure 4 . 3 :

 43 Figure 4.3: Illustration of 1/λ and 2 ln(2/λ) behaviour for all λ ∈]0; 1[.

Figure 4 .

 4 Figure 4.3 and the inequality (4.23) show that whatever the problem size n and for a probability at most ≈ 0.7678, the BC method gives a tighter probabilistic bound than the AH2 method.

  73

Figure 4 . 4 :

 44 Figure 4.4: AH1, AH2 and BC bounds with probability 0.9 and u = 2 -23 for the inner product.

Figure 4 .

 4 Figure 4.4 confirms the discussion of Sub-section 4.3.1 and shows that with a probability 0.9, when nu ≫ 1, AH2 bound grows rapidly compared to AH1

Figure 4 . 5 : 5 .

 455 Figure 4.5: Probabilistic error bounds with probability 1 -λ = 0.5 (left) and 1 -λ = 0.9 (right) vs deterministic bound for the Horner's evaluation of T 20 (x) and u = 2 -23 . Triangles mark 30 instances of the SR-nearness relative errors evaluation in binary32 precision, a circle marks the relative errors of the 30 instances average, and a star represents the IEEE RN-binary32 value.

Figure 4 . 6 :

 46 Figure 4.6: Normalized forward error (error/K) with probability 1 -λ = 0.5 (left) and 1-λ = 0.9 (right) for Horner's evaluation of T N (24/26) and u = 2 -23 .

Figure 4 . 7 :

 47 Figure 4.7: Probabilistic bounds with probability 1 -λ = 0.9 vs deterministic bound of the computed forward errors of the inner product with u = 2 -23 .

Figure 4 . 8 :

 48 Figure 4.8: AH1 bound vs BC bound with probability 1 -λ = 0.9 and u = 2 -23for the inner product.

Figure 4 .

 4 Figure 4.8 illustrates the advantage of using (BC-IP) and shows that for a large n ≥ 10 13 and u = 2 -23 , the AH1 bound increases faster than the BC bound. This unsurprising result confirms the previous theoretical analysis in Section 4.3.

x 2 i

 2 (ψ i -1) and | s -s| represent the absolute errors of the inner product n i=1 x 2 i and the summation s = n i=1 x i , respectively. Then [23, sec 5.1] proves that n i=1

( 5 . 8 )Lemma 16 .

 5816 Note that |H(X, δ) -H(X, 0)| = |M + A| ≤ |M | + |A|. To bound the martingale M by Azuma-Hoeffding inequality, we need to bound the martingale steps. The martingale M 0 , . . . , M 2n-1 satisfies |M

9 Figure 5 . 1 :

 951 Figure 5.1: Probabilistic error bounds over n with probability 1 -λ = 0.9 vs deterministic bound for the textbook-variance algorithm and u = 2 -23 .

Figure 5 . 2 :

 52 Figure 5.2: Probabilistic error bounds over n with n = 10 6 and over λ vs deterministic bound for the textbook-variance algorithm and u = 2 -23 .

Figure 5 . 3 :

 53 Figure 5.3: The forward errors of textbook-variance and two-pass-variance algorithms in binary32 precision for floating-points chosen uniformly at random in [-1; 1].

Figure 5 . 4 : 2 i and 1 n s 2

 54212 Figure 5.4: The forward errors of textbook-variance and two-pass-variance algorithms in binary32 precision for floating-points chosen uniformly at random in [1024; 1025].

For 1 .

 1 a, b, i, j ∈ N and ϵ i , ϵ j , ϵ a , ϵ b ∈ R, the MP2 method involves the compuϵ ijab = ϵ i + ϵ j -ϵ a -ϵ b with ϵ i , ϵ j ≤ 0, and ϵ a , ϵ b ≥ 0. Then, the condition number cond(ϵ ijab ) =|ϵ i |+|ϵ j |+|ϵa|+|ϵ b | |ϵijab| =For a fixed i, j, a, b, we haveεijab = (ϵ i + ϵ j )(1 + θ ijab 1 ) -ϵ a (1 + θ ijab 2 ) -ϵ b (1 + θ ijab 3 ).

  -b) εijab k∈K ijab (1 + α k )(1 + β k )(1 + δ k )(1 + η k )(1 + µ k )

)

  with probability at least 1 -λ. With K = number of E using the 1-norm, and f (u) = 3u + 3u 2 + u 3 .Proof. For all 0 ≤ a, b ≤ n virt and 0 ≤ i, j ≤ n occ , denote θ l = θ ijab l where l = 1, 2, and 3. Then, Equation (6.2) yieldsεijab = (ϵ i + ϵ j )(1 + θ 1 ) -ϵ a (1 + θ 2 ) -ϵ b (1 + θ 3 ) = ϵ i + ϵ j -ϵ a + (ϵ i + ϵ j )θ 1 (1 + θ 2 ) -ϵ b (1 + θ 3 )

(6. 5 )

 5 Using SR-nearness and the mean independence, we have E(θ k ) = 0 and E(θ k θ l ) = 0 for k, l = 1, 2 or 3. Then, E(χ ijab ) = 0, andE(ε ijab ) = ϵ ijab . 2a(a -b) ϵ ijab (1 + χ ijab ) ψ ijab -2a(a -11 implies that E(ψ ijab ) = 1 and V (ψ ijab ) ≤ γ N (u 2 ), where N = n virt n occ . Therefore, Bienaymé-Chebyshev inequality yields|ψ ijab -1| ≤ V (ψ ijab ) λ ≤ γ N (u 2 ) λ ,with probability at least 1 -λ. Thus least 1 -λ. Finally

1 1 -

 1 f (u) = 1 + f (u) + O(u 2 ) = 1 + 3u + O(u 2 ).It follows that(1 + u) N f (u) 1 -f (u) + γ N (u 2 ) λ = (1 + N u)(1 + 3u)3u + N/λu + O(u 2 ) = (3 + N/λ)u + O(u 2 ).

Figure 6 . 1 :

 61 Figure 6.1: SR versus RN in single and double precision. The probabilistic bound B is given with probability at least 0.9.

x

  i ψ , where ψ i = n k=max (i,2) (1 + δ k ). Denote Φ = ŝ-s s and Ψ = ŷ-y y the relative errors of ŝ and ŷ, respectively. We have ŷ = 1 ŝ (1 + θ), where θ represents the rounding error of the division. Our goal is to compute a probabilistic bound for the relative error Ψ. Let us first present a model of computation for this problem. Note that ŝ = s(1 + Φ) and ŷ = y(1 + Ψ).

  Finally, |Ψ| = ŷ -y y = θ + Φ(θ + 1)(Φ -1) + O(Φ 3 ) .

1 . 29 .

 129 computes a deterministic bound of the relative error Ψ when nu ≪ Theorem The computed Ψ satisfies |Ψ| ≤ (K(n -1) + 1) u + O(u 2 ), We have shown in Equation (2.5) that |Φ| ≤ Kγ n-1 (u). Because |θ| ≤ u, using Equation (6.6) yields|Ψ| ≤ |θ| + |Φ| |θ + 1| (|Φ| + 1) + O(Φ 3 ) ≤ u + Kγ n-1 (u)(1 + u)(Kγ n-1 (u) + 1) + O (Kγ n-1 (u)) 3= (K(n -1) + 1) u + O(u 2 ).

Theorem 30 . 3 =

 303 For all 0 < λ < 1, the computed Ψ satisfies under SR-nearness|Ψ| ≤ K 2(n -1) 1/λ + 1 u + O(u 2 ),with probability at least 1 -λ, where K = We have shown in Theorem 7 that |Φ| ≤ K γ n-1 (u 2 )/λ with probability at least 1 -λ. Because |θ| ≤ u, using Equation (6.6) yields|Ψ| ≤ |θ| + |Φ| |θ + 1| (|Φ| + 1) + O(Φ 3 ) ≤ u + K γ n-1 (u 2 )/λ(1 + u)(K γ n-1 (u 2 )/λ + 1) + O (K 2(n -1) ln(2/λ)) K 2(n -1) 1/λ + 1 u + O(u 2 ),with probability at least 1 -λ.

Theorem 31 .

 31 For all 0 < λ < 1, the computed Ψ satisfies under SR-nearness|Ψ| ≤ K 2(n -1) ln(2/λ) + 1 u + O(u 2 ), with probability at least 1 -λ, where K = n i=1 |x i | | n i=1 x i|

  The following table summarizes various binary floating point formats. The precisions of FP numbers commonly used in numerical computations defined in the latest revision of the IEEE-754 standard are given by binary-k, where k = 32, 64 or 128.

	Precision Bits Sign Mantissa Exponent	e min	e max
	binary-128 128	1	112	15	-16382 16383
	binary-64	64	1	52	11	-1022	1023
	binary-32	32	1	23	8	-126	127
	binary-16	16	1	11	4	-14	15
	bfloat-16	16	1	8	7	-126	127
		Table 2.1: Binary floating-point formats

Table 2 .

 2 1 presents a variety of binary floating-point formats.

becomes smaller for zero-mean distributions. This is confirmed by the experiment of this figure, which shows a similar forward error for the two algorithms, whether for SR or RN.

Acknowledgments

Remark 9. The inequality (4.14) gives V (ŝ) ≤ ∥x∥ 2 1 γ n-1 (u 2 ). Then, the bias satisfies

where

. Thus E (H(X, δ)) ≥ H(X, 0) 1 -K 2 1 γ n-1 (u 2 ) .

We now turn to bound the forward error of this algorithm. We employ three methods to establish probabilistic bounds on the error: the DM method, which uses Equation (5.1), and the BC and AH methods, which use similar techniques to the previous general framework.

. DM Method

This sub-section uses a method based on the Doob-Meyer decomposition to provide a probabilistic bound on the forward error of the textbook-variance algorithm under SR-nearness. We need before to demonstrate the following lemma: Proof.

The proof of the second item uses the first point and the following property log(ab) = log(a) + log(b).

We have demonstrated in Equation (5.2) that the error of the textbookvariance algorithm under SR-nearness can be written as:

For the martingales n i=1 x 2 i (ψ i -1) and 1 n s 2 (ψ n+1 -1), we will apply the previous result on the inner product and summation. Let us examine the martingale X 1 , . . . , X n .

Finally,

We now have all the tools to state and demonstrate the main result of this sub-section:

with probability at least 1 -λ, where

.

Proof. From the Inequalities (5.9) and (5.11) we have

with probability at least 1 -λ.

Discussion: Interestingly, this application illustrates the power and richness of this method, showing that it can be applied to any algorithm. For the two-pass-variance algorithm, we did not consider the errors made in the computation of m but instead constructed our stochastic process from the random errors generated during the remaining elementary operations. This choice sufficed to establish a probabilistic bound of ( √ nu). Note that the generalization 5.1 can be applied to the entire algorithm, allowing us to construct a stochastic process that depends on all random errors in the computation.

In some applications, studying this stochastic process can be challenging to investigate.

. BC Method

We present a computational scheme for the proofs of the two-pass-variance algorithm errors in this sub-section. This computational scheme allows to use both BC and AH methods. Recall that H(X, δ) satisfy

where

Therefore, one needs first to separate the errors of order two.

The following equation will be used to compute a probabilistic bound on the two-pass-variance forward error by the BC and AH method:

(5.12)

The following theorem presents a probabilistic bound of the forward error of this algorithm through the BC method. • γ n (u) ≈ nu + O(u 2 ) when nu ≪ 1.

•

This table displays the advantage of the probabilistic bounds of the textbookvariance forward error in terms of O( √ nu) compared to the deterministic bounds in O(nu), when nu ≪ 1. Additionally, the BC method is far better when nu ≫ 1 and nu 2 ≪ 1. The previous discussion also holds for the twopass-variance forward error bounds.

. Numerical Experiments

We performed a series of numerical experiments comparing these new probabilistic bounds to the deterministic ones. We show that probabilistic bounds are tighter and accurately reflect the behavior of SR-nearness forward errors. Two types of plots are presented. Firstly, the plots are displayed over n, and show that for large values of n, BC bounds provide significant benefits compared to AH or DM bounds for the textbook-variance. Secondly, the plots are shown over λ, and show that AH bound holds a significant advantage for higher probabilities. All SR computations are repeated 30 times with verificarlo [START_REF] Denis | Verificarlo: Checking floating point accuracy through Monte Carlo arithmetic[END_REF]. All samples and the forward error of the average of the 30 SR instances are plotted. bounds on the forward error proportional to √ nu, which is better than the deterministic bounds in O(nu).

While introducing pairwise algorithm in textbook-variance and two-passvariance algorithms, SR leads to probabilistic bounds proportional to log(n)u, instead of O(log(n)u) for deterministic bounds. We also demonstrate that the two-pass-variance algorithm performs better than the textbook-variance algorithm under SR, as it does under RN.

The generalization in Section 5.1 completes the scope of the probabilistic analysis of problems under SR. With Equation (5.1), and under SR, we can decompose any algorithm into two parts: a martingale to which we apply AH or BC methods and a drift for which we provide the exact expression and bound it deterministically.

The two examples treated in this chapter demonstrate the flexibility of this generalization in obtaining the desired decomposition. When a part of the error directly provides a martingale, the focus can be on the remaining part, as illustrated in the textbook-variance algorithm. When the algorithm is more complex, we have shown, as demonstrated in the two-pass-variance algorithm, that we can construct a stochastic process without considering the previous errors in the average computation and ensure tight probabilistic bounds.

The generalization can be applied to the entire algorithm, and the stochastic process should be chosen carefully. Equation (5.1) ensures the exact expression of the drift and, consequently, the martingale term. However, studying the martingale to obtain a probabilistic bound can be challenging in some complex situations.

The scripts for reproducing the numerical experiments of this chapter are published in the repository: https://github.com/verificarlo/sr-non-linear-bounds.

probability at least 1 -λ. Because |θ| ≤ u, using Equation (6.6) yields

with probability at least 1 -λ.

Note that the probabilistic analysis in Section 4.3 remains valid for the previous three bounds.

In this chapter, our contribution is twofold: We have investigated an industrial application with SR that involves divisions, and we have studied the error of the inverse of a summation with SR. For the MP2 computation, we have shown an interest in using SR with low-precision and demonstrated a probabilistic error bound in O( √ nu). We have proposed a model for the inverse summation that enables the computation of probabilistic error bounds in O( √ nu).

-Conclusion

Stochastic rounding has drawn a lot of attention in various domains [START_REF] Croci | Stochastic rounding: implementation, error analysis and applications[END_REF] due to its efficiency compared to the deterministic IEEE-754 [START_REF]IEEE standard for floating-point arithmetic, IEEE Std[END_REF] default rounding mode. The numerical applications in domains such as machine learning [START_REF] Gupta | Deep learning with limited numerical precision[END_REF], or climate modeling [START_REF] Paxton | Climate modeling in low precision: Effects of both deterministic and stochastic rounding[END_REF] have shown that SR provides positive effects in computation, especially in low-precision formats, for example, by avoiding stagnation effects.

. Findings from Our Research

It is well known that the worst-case bound of a computation involving n elementary operations is given by O(nu). However, Wilkinson had the intuition that the error bound for this computation is typically proportional to O( √ nu). In this thesis, we proved the validity of this intuition in the context of SR through probabilistic error bounds. But, the intuition is not yet fully verified for RN.

SR errors satisfy the mean independence property, allowing us to derive tight probabilistic error bounds using two methods: AH or BC. The AH method is based on the martingale theory and Azuma-Hoeffding inequality, which is preferable for higher probabilities (1 -λ near 1). We have demonstrated in Section 4.1 the applicability of this method to several algorithms that exhibit multi-linear error. The BC method is based on Bienaymé-Chebyshev inequality and a bound of the variance, which is suitable for large problem sizes n.

We have demonstrated Lemma 11, a general framework applicable to a wide class of algorithms that allows to compute a variance bound. We have shown in Section 4.2 the strength of this new approach to the same previous algorithms studied by the AH method.

In Sub-sections 4.1.4 and 4.2.5, we have shown that the two previous methods can be generalized to any complete computation tree with a multi-linear sequence of elementary operations {+, -, * }. In the case of only additions/subtractions, the probabilistic error bound obtained is proportional to log(n) (Pairwise summation). At the same time, the presence of multiplications leads to a probabilistic error bound in O( √ n).

In this thesis, we have proposed a general framework for the probabilistic analysis of algorithms under SR. Using the Doob-Meyer decomposition, we have demonstrated in Chapter 5 that the error of any algorithm can be decomposed into a martingale plus a drift. It is essential to note that the choice of the stochastic process is crucial in the use of this generalization, as well as the quality of the final bound. The generalization suggests an analysis scheme