Numerical Stability of Linear Algebra Operations over Relational Databases

Numerical Stability of Linear Algebra Operations over Relational Databases Andrei Draghici University of Zurich

Yuchen He University of Zurich

Dan Olteanu University of Zurich

andrei.draghici@uzh.ch

yuchen.he@uzh.ch

dan.olteanu@uzh.ch

arXiv:2609.13802v1 [cs.DB] 12 Sep 2026

Abstract A large body of work in the database literature develops efficient algorithms for linear algebra and machine learning over matrices defined by relational joins, yet the numerical stability of such computations has so far received no attention. This is a practical concern: a join matrix can be much larger than the input database, and the repeated copies of input values it contains compound the floating-point errors incurred by the numerical operations performed over it. This paper initiates a formal investigation of numerical stability for linear algebra over database joins. We first show that backward stability, the standard yardstick of numerical stability, loses its effectiveness in this setting: join matrices form a structured subspace of the ambient matrix space, so a perturbation explaining a computed result need not correspond to any perturbed input database. This failure already occurs for operations as simple as matrix-vector multiplication. To overcome this limitation, we introduce projected backward stability, a generalization of backward stability from the computation of one function to that of a composition of two functions, and establish its connection to classical backward stability. In our database setting, the two functions are the join query and the numerical operation. We further introduce the database condition number as the square root of the ratio of maximal to minimal number of copies of input data values into the join matrix, and show that it quantifies how a perturbation of the join matrix is amplified into a perturbation of the input database, independently of the computation used. The database condition number coincides with the classical condition number of the expansion matrix that replicates input values into the join matrix. Finally, we show that the QR decomposition and the Singular Value Decomposition are projected backward stable and that the FiGaRo algorithm, which pushes these decompositions past α-acyclic joins, admits an error bound that can be asymptotically smaller than the size of the join matrix. 2012 ACM Subject Classification Information systems Ñ Database management system engines; Mathematics of computing Ñ Numerical analysis Keywords and phrases backward stability, condition number, linear algebra over relational data

1

Introduction

In this paper we start an investigation into the numerical stability of linear algebra operations over matrices defined by database joins. A significant body of research in the database literature, e.g., [23, 28, 26, 20, 9, 27], put forward efficient algorithms for computing linear algebra operations and also for training a variety of machine learning models over such matrices, albeit no consideration was given to their numerical stability. The risk of numerical instability is exacerbated by the multiple copying of input data in the join result that compounds the effect of numerical errors that come with such operations.

1

Numerical stability is not merely a concern for traditional scientific computing; it directly affects the reliability of modern data systems. Database engines increasingly support numerical analytics, statistical estimation, and machine-learning workloads, where relational joins define large matrices over which downstream computations are performed. The need to integrate complex mathematical computations into data management systems has driven a rich line of research exploring the connection between relational algebra and linear algebra. To bridge the semantic gap between these domains, frameworks like LaRa have been proposed as minimalist unified algebras [4]. Foundational frameworks like MATLANG have been introduced to formalize the expressive power of matrix operations and their exact relationship to relational query languages [6, 7]. This theoretical bridge yields direct algorithmic benefits, such as mapping fragments of linear algebra to tractable conjunctive queries to enable constant-delay enumeration and efficient updates for matrix computations [25]. Building on this unified perspective, Galley applies modern relational query optimization techniques to compiling sparse tensor programs into highly efficient imperative kernels [11], LaraDB demonstrates that relational engines can efficiently execute linear algebra tasks [15], and the IFAQ [30] and SDQL [29] frameworks allow users to specify, optimize, and efficiently execute workloads that mix linear algebra and relational algebra operations. Fast matrix multiplication is integrated into query plans to lower the complexity of query evaluation [19]. However, as database systems take on complex numerical workloads, precision becomes a critical challenge. Prior work has shown that even the computation of simple aggregates can suffer from severe floating-point error: various implementations in modern database systems lose precision due to catastrophic cancellation, while floating-point aggregation can become non-reproducible because different execution orders produce different numerical results [16, 24]. Computing the marginal probability of a query in large probabilistic databases is known to incur large round-off errors [22]. Similar issues arise in modern machine-learning software, where mathematically equivalent implementations may have very different numerical-stability properties and generally inaccurate results [21]. This concern has become particularly acute in the modern era of artificial intelligence. Recent investigations into fundamental deep learning operators reveal significant vulnerabilities to finite-precision effects [10]. Notably, the transformer architecture has been demonstrated to suffer from numerical instability [8], driving the need for rigorous analyses of both forward and backward stability during model training [17]. These examples suggest that formal correctness of exact arithmetic over large data (such as query results) is insufficient in practice: we need formal guarantees controlling how finite-precision errors propagate through the computation. This motivates our study of numerical stability in computations with matrices defined by database queries. Example 1. The least squares problem argminx }Ax ´ b}22 is central to polynomial regression and a wide range of machine learning tasks. In our setting, the matrix A is defined by a database query: its rows represent tuples in the query result and its columns represent features. The vector b is also a column in the query result. The solution x gives the fitted model parameters [28]. Since practical algorithms operate in finite precision arithmetic, the computed result is generally not the exact solution of the mathematical problem [14, 12]. A key question is how far off the computed result is from the mathematical result. Backward stability is the main measure of numerical stability. A computation of a function f is backward stable if its output ŷ equals the exact value of f px ` ∆xq on a slightly perturbed input x ` ∆x [14, 12], as illustrated in Fig. 1. For instance, f can be the composition of the join query defining the matrix and the numerical operation computed over the join matrix, whereas x and ∆x are databases over the same schema and join keys for every join attribute and that may differ in the values for the data attributes. To our knowledge, this paper is the first to formally analyze the numerical stability of operations over matrices defined by join queries directly with respect to perturbations to the input 2

Input space x

backward perturbation x ` ∆x

Output space f

com put

atio n

f

y “ f pxq of f

forward error ŷ “ f px ` ∆xq

Figure 1: Illustration of backward stability for a computation of a function f . The computation output ŷ equals the exact value of f on a slightly perturbed input x ` ∆x. The smallest possible perturbation |∆x| is the backward error, while |ŷ ´ y| is the forward error. database. We consider matrices defined by natural join queries. Such join matrices can be much larger than the input database as they may have many copies of input values. These copies compound the effect of computational errors that come with numerical operations. We start with two simple operations, the Frobenius norm and matrix-vector multiplication, and then move to the QR decomposition and Singular Value Decomposition of the join matrix. Our setting calls for a refinement of Fig. 1, as we now have a composite function g ˝ f , where f denotes the join query and g denotes the subsequent matrix operation. The numerical stability of composed operations has been very recently recognized as an important challenge across disciplines: the numerical analysis community has investigated conditions under which forward stable algorithms compose stably [5], while the programming languages (PL) community has developed formal methods to check programs for numerical instability [31] and to statically derive backward error bounds for composed operations [18]. However, our database setting introduces a special type of composition that prior work does not address. Because the intermediate output (the join matrix) is structurally constrained to a specific subspace, we show that the classical notion of stability becomes insufficient. In our setting, we address two fundamental questions: 1. Does the backward stability of the join and numerical operation computations translate to the backward stability of a computation of their composition? 2. Can tighter error bounds be derived based on the size of the input database and the query structure rather than on the large size of the join matrix? Our answer to the first question is that, even though (i) any join computation is trivially backward stable as it does not involve arithmetic computations on the data and (ii) there are backward stable algorithms for matrix multiplication, QR, and SVD [14], there may be no backward stable computation of the composition of the join and the numerical operation. This is because a join matrix A belongs to a specific subspace of structured matrices, rather than to the full ambient matrix space. Consequently, a small perturbation ∆A to A that explains the computed result may not lie in the join-induced subspace. Thus, there may not exist a perturbed input database for which the join query precisely defines the matrix A ` ∆A, i.e., the corresponding perturbed input to g may not lie in the range of the join query f . To resolve the domain mismatch between the full matrix space and the join-induced subspace, we introduce the notion of projected backward stability, illustrated in Fig. 2. Projected backward stability generalizes classical backward stability to the composition g ˝ f of functions f and g, where f restricts its outputs to a structured subspace Y0 . As in general it is not the case that the perturbed input ỹ to g lies in Y0 , we decompose it into ỹ “ ỹ0 ` ∆y0K , where ỹ0 is the projection of ỹ onto Y0 , so that ỹ0 is the output of f for some perturbed input x ` ∆x, and ∆y0K is a small perturbation orthogonal to the subspace Y0 . That is, ỹ0 “ PY0 pỹq “ argmin }ỹ ´ ỹ0 }F , ỹ0 PY0

and ∆y0K “ ỹ ´ ỹ0 . This also means that ỹ0 is the closest to ỹ according to the Frobenius norm

3

Input space X x

projected backward perturbation x ` ∆x

f

Intermediate space Y y “ f pxq g comp

f ỹ0 “ f px ` ∆xq P Y0

utat

ion o

fg˝

Output space Z z “ gpyq “ gpf pxqq

f

g ỹ “ ỹ0 ` ∆y0K P Y

ẑ “ gpỹq “ gpf px ` ∆xq ` ∆y0K q

Figure 2: Illustration of projected backward stability for an algorithm that implements the composition g ˝ f of two functions, where the first function is f : X ÞÑ Y0 and the second function g : Y ÞÑ Z where Y Ě Y0 . The input perturbation ∆x is small relative to x and the projection error ∆y0K is also small relative to f px ` ∆xq. } ¨ }F . The projected backward perturbation ∆x is the difference between the true input x to f and a nearby input x ` ∆x whose image under f matches this projection. We argue that projected backward stability is the appropriate numerical stability concept in our setting: it accounts for the (join dependency) constraints while still permitting the reconstruction of a perturbed input database close to the original, in direct analogy to classical backward stability. In particular: f is the join query; X is the set of databases; Y0 is the set of join matrices; g is the numerical operation, e.g., matrix-vector multiplication or SVD; and Y is the set of all matrices. Then, the perturbation ∆y0K is the minimal change to a perturbed matrix to turn it into a join matrix for a given join query. Finally, the projected backward perturbation ∆x is the change to the input database that explains mathematically the result of g’s computation. The computation of g ˝ f is projected backward stable with respect to the input database if ∆x is small relative to x and ∆y0K is small relative to ỹ0 . We show that there is a backward stable computation for g with respect to its input matrix if and only if there is a projected backward stable computation for g ˝ f with respect to the input database for every join query f . We can also derive the backward error from the projected backward error and vice-versa. Towards an answer to the second question, we show how to bound the projected backward error using a new parameter that we call the database condition number. This is a natural analogue to the fundamental notion of condition number1 of matrices, which measures how much a problem’s output changes when the input is slightly perturbed [14]. The (database) condition number is a property of the problem and not of a particular computation. The database condition number of a relation is defined as the square root of the ratio of the maximal and minimal counts of an input data value in the join matrix. We show that the perturbation in the join matrix yields a perturbation in the input database that is amplified by the database condition number, completely independent from the choice of the computation. This is in addition to the computation-dependent numerical error. In particular, an error upper bound c in the join matrix translates to an error upper bound κjoin ¨ c in the input database, where κjoin denotes the database condition number. There is a natural relationship of the database condition number and the classical condition number: The former is precisely the condition number of the expansion matrix that transforms the input data columns into the corresponding columns in the join matrix. This transformation creates a number of copies of the input values 1

The condition number is defined as the value of the asymptotic worst-case relative change in output for a relative change in input. A problem with a low condition number (close to 1) is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned. A problem is ill-conditioned if a small change in the input leads to a large change in the output. A (projected) backward stable computation can be expected to accurately solve well-conditioned problems.

4

in the join matrix as dictated by the join query. The rows in this expansion matrix thus select and replicate entries from the input relation. For the QR and SVD decompositions of the join matrix, we answer the second question in the positive: We give an upper bound on the error of the FiGaRo algorithm [27] that depends on the size of the input database and not on the size of the join matrix. To sum up, the main contributions of this paper are as follows: 1. We initiate a formal investigation on numerical stability of linear algebra operations over database joins. We show that in this setting, the widely accepted yardstick of numerical stability, namely backward stability, loses its effectiveness already for simple operations like the matrix-vector multiplication (Sec. 3). 2. We introduce a natural generalization of backward stability called projected backward stability (Sec. 4) and we show its connection to the classical backward stability (Sec. 6). 3. We introduce the notion of database condition number that amplifies the error bound of the computation over the join matrix into the error in the input database. We also exhibit the connection between this number and the classical condition number (Sec. 5). 4. We show that the QR and SVD decompositions of the join matrix are projected backward stable in our database setting and their backward error bounds depend on the number of computational steps performed over the join matrix. We further show that the FiGaRo algorithm [27], which pushes the computation of the upper triangular matrix of the QR decomposition past α-acyclic joins, can have an asymptotically tighter error bound (Sec. 7).

2

Preliminaries

In this section, we recall standard relational algebra and linear algebra terminology and matrices defined by database joins, as used in prior work [27]. Matrices. We denote matrices by bold upper-case letters and column vectors by bold lower-case letters. For a real number x P R, |x| denotes its absolute value. By rns we denote the set t1, . . . , nu. For a matrix A P Rmˆn and indices i P rms, j P rns, we use Ari, js, Ari, ˚s, and Ar˚, js to denote its entry at row i andbcolumn j, its i-th row, and respectively its j-th column. řm řn def m 2 The Frobenius norm of A is ∥A∥F “ i“1 j“1 Ari, js . For a vector a P R , aris denotes def ařm mˆn , 2 its i-th component and its Euclidean norm is ∥a∥2 “ i“1 aris . For matrices A, B P R |A| ď |B| means that |Ari, js| ď |Bri, js| for all i P rms, j P rns. Databases and Queries A database consists of a set of relations S1 , . . . , Sr . Each relation Si has schema pXi , Yi q, where Xi is the tuple of its join attributes and Yi is the tuple of its data attributes. Each data attribute takes values in R. For relations Si and Sj , we write Xij for the tuple of join attributes shared by both relations. We consider full conjunctive (or join) queries expressed as Q “ S1 pX1 , Y1 q ’ ¨ ¨ ¨ ’ Sr pXr , Yr q, where two relations Si and Sj equi-join on Xij ‰ H. A query is α-acyclic if and only if it admits a join tree [1]. From Databases to Matrices A relation Si with schema pXi , Yi q defines a row-labeled, column-labeled real-valued matrix Si , where the rows of Si are indexed by the tuples of Si and the columns are indexed by the attributes in Yi . For every tuple si “ pxi , yi q P Si , the join keys xi are not entries of Si ; they appear only in the row label. Two tuples with the same join key but different data are distinct rows. We denote by Si r˚, Ys the matrix formed by the columns in Y of Si . A block-diagonal matrix Ablock with blocks A1 , . . . , Ab is compactly denoted by the tuple pA1 , . . . , Ab q. The matrix D “ pS1 , . . . , Sr q defined by the database D is block-diagonal, with the blocks S1 , . . . , Sr . The matrix defined by the result of a join query over a database is

5

called the join matrix, denoted by A. We fix an arbitrary order of the row labels of a matrix; row permutations do not affect the algebraic quantities considered in this paper. Given a relation S, we denote by ∆S a relation with the same schema and the same join keys, but with possibly different data values. The sum S ` ∆S is the entry-wise addition of the data values. For a database D “ pS1 , . . . , Sr q, we denote ∆D by p∆S1 , . . . , ∆Sr q and D ` ∆D by pS1 ` ∆S1 , . . . , Sr ` ∆Sr q. The join space of D is the set of all join matrices corresponding to a perturbed database D ` ∆D of D. Such perturbed databases induce a join space that is a subspace of the matrices in Rmˆn (Lemma 23). Example 2. Consider the join Q “ S1 pY1 , X12 q ’ S2 pX12 , Y2 , X23 q ’ S3 pX23 , Y3 q and the input relations and the join matrix shown below: S1

S2

X12 Y1 k11 a1 k12 a2

X12 X23 Y2 k11 k21 b1 k12 k21 b2

A

S3

row label in query result Y1 Y2 Y3 pk11 , k21 , a1 , b1 , c1 q a1 b1 c1 pk12 , k21 , a2 , b2 , c1 q a2 b2 c1

X23 Y3 k21 c1

The keys of the join attributes X12 and X23 are only retained in the row labels of A; the actual matrix entries are values for the data attributes Y1 , Y2 , Y3 . Definition 3 (Our DB Setting). • The input database D “ pS1 , . . . , Sr q defines the rowlabeled, column-labeled matrices D “ pS1 , . . . , Sr q. We denote by d the total number of rows in D and by ℓ (n) the total number of values (columns) in S1 , . . . , Sr . • A perturbation to D “ pS1 , . . . , Sr q is denoted by ∆D “ p∆S1 , . . . , ∆Sr q. • A bound |∆D| ď c ¨ |D| is a shorthand notation for |∆Si | ď c ¨ |Si | for i P rrs. • The join query is Q “ S1 ’ ¨ ¨ ¨ ’ Sr and has no self-joins. • The query result QpDq defines the join matrix A P Rmˆn that has m rows and n columns. Numerical Stability. In numerical computations, it is essential to quantify the effects of rounding errors that arise due to finite-precision arithmetic. Let flpxq denote the floatingpoint number obtained by rounding the real number x to the nearest floating-point number representable by the computer. The unit round-off, denoted by u, measures the relative error incurred in representing a real number in floating-point arithmetic. Then, flpx op yq “ px op yqp1 ` δq, where |δ| ď u, for some δ P R, op P t`, ´, ˆ, {u. A convenient and common def nu notation for bounding the accumulation of rounding errors is the quantity γn “ 1´nu . When nu ă 1, it follows from standard rounding models that a product of n floating-point factors satisfies flpx1 ¨ ¨ ¨ xn q “ px1 ¨ ¨ ¨ xn qp1 ` θn´1 q, for |θn´1 | ď γn´1 . Thus γn gives a compact way to express worst-case relative error growth in sequences of operations, especially products. It is often useful to allow in such bounds a small positive integer c constant independent of the def cnu problem size as a multiple of nu: γ̃n “ 1´cnu . This notation γ̃n enables concise statements of error bounds when several rounding terms combine. Both γn and γ̃n are treated as quantities of order nu whenever nu ! 1 and cnu ! 1 [14]. A numerical algorithm is backward stable if the computed solution can be interpreted as the exact solution to a slightly perturbed input. Definition 4 (Backward Stability [14]). Let C be a numerical computation of a function f : if y “ f pxq for some x then ŷ “ Cpxq. We say that C is backward stable if and only if there exists ∆x such that ŷ “ f px ` ∆xq and ∆x is small relative to x. The backward error is defined as min∆x:f px`∆xq“ŷ |∆x|. The computation C is ϵ-backward stable if for any input x, its backward error satisfies |∆x| ď ϵ ¨ |x|. 6

We use the convention that computed quantities wear a hat: ŷ denotes the computed approximation to y. A numerical computation is forward stable if the computed solution ŷ is close to the exact output: f pxq “ ŷ ` ∆y, where ∆y is small relative to y. A computation of f is mixed forward-backward stable if its output is close to the exact output corresponding to a nearby admissible input: f px ` ∆xq “ ŷ ` ∆y. Perturbation Bounds. For a matrix A and a perturbation ∆A, we distinguish three types of relative perturbation bounds: (1) Entry-wise: |∆Ari, js| ď ϵ ¨ |Ari, js|; (2) Column-wise: }∆Ar˚, js}2 ď ϵ ¨ }Ar˚, js}2 ; and (3) Frobenius-norm-wise: }∆A}F ď ϵ ¨ }A}F . These bounds are progressively weaker: an entry-wise ϵ-perturbation is also a column-wise ϵ-perturbation, which is also a Frobenius-norm-wise ϵ-perturbation. Entry-wise backward stability implies column-wise backward stability, which implies Frobenius-norm-wise backward stability. Condition Number of a Matrix [14]. Consider a rank-k matrix with its k non-zero singular values σ1 ě ¨ ¨ ¨ ě σk ą 0. Its spectral-norm condition number, denoted hereafter as the condition def number, is κpAq “ σ1 {σk .

3

Backward Stability with Respect to the Input Database

We initiate the investigation of backward stability of simple numerical operations over join matrices with respect to the input database. Whereas the stability of many numerical operations with respect to their input matrices is well understood [14], our setting, which investigates the stability of the composition of the numerical operation with the join, has not been addressed so far. In this section, we exemplify our setting for two simple numerical operations: the Frobenius norm and matrix-vector multiplication. Computing the Frobenius norm over an arbitrary matrix is backward stable [14] and remains backward stable also with respect to the input database. Proposition 5. Consider the setting in Def. 3. The computation C of an operation f , which first materializes the join matrix and then computes its Frobenius norm, is backward stable with respect to the input database. That is, there exists a perturbation ∆D such that, if ẑ “ CpDq, then ẑ “ f pD ` ∆Dq. Moreover, |∆D| ď γ̃mn ¨ |D|. Proof. Let z “ f pDq “ }A}F , where the join matrix A is fully materialized. If z “ 0, then every entry of A is zero, and we choose ∆D “ 0. Suppose now that z ą 0. The computation C of f needs one pass over A, so it does mn square operations followed by the square root of their sum. Thus the computation requires Opmnq arithmetic operations. Since all terms in the sum are nonnegative and the computation of the square root is forward stable, the forward error satisfies |ẑ ´ z| ď γ̃mn z. We write ẑ ´ z “ θz for some θ such that |θ| ď γ̃mn . We construct ∆D by scaling every numerical attribute included in the matrix A by the factor 1 ` θ, denoted compactly by ∆D “ θD. Then, D ` ∆D “ p1 ` θqD. Now, consider the result of f pD ` ∆Dq. The key attributes are unchanged and every entry of the materialized join matrix is scaled by the same factor 1 ` θ. Then, A ` ∆A “ p1 ` θqA and }A ` ∆A}F “ }p1 ` θqA}F “ |1 ` θ| ¨ }A}F “ |1 ` θ|z “ p1 ` θqz “ ẑ. In the previous chain of equalities, |1 ` θ| “ 1 ` θ since |θ| ď γ̃mn . Finally, ∆D “ θD ùñ |∆D| ď γ̃mn ¨ |D|. Prop. 5 establishes backward stability by analyzing the standard algorithm that first materializes the join matrix and then computes its Frobenius norm. However, materializing the join is unnecessary, as observed for aggregates over joins [28]: we can push the computation of the norm past the join. Besides avoiding the materialization, this yields a smaller error bound. Whereas the bound in Prop. 5 depends on the number of values in the join matrix (mn), the improved bound only depends on the number of values in the input database (ℓ). Since the materialized join can be much larger than the input database (ℓ ! mn), this means a substantially smaller backward error. 7

Proposition 6. Consider the setting in Def. 3. The computation C of an operation f , which pushes the computation of the Frobenius norm of the join matrix past the join, is backward stable with respect to the input database. That is, there exists a perturbation ∆D such that, if ẑ “ CpDq, then ẑ “ f pD ` ∆Dq. Moreover, |∆D| ď γ̃ℓ ¨ |D|. Even though standard matrix-vector multiplication is backward stable with respect to the join matrix [14], it is not with respect to the input database: there exist problem instances for which no perturbation of the input database can explain the computed result. Proposition 7. Consider the setting in Def. 3 and the operation f that defines the matrix-vector multiplication fx pDq “ Ax for a vector x. Then, there exist a join matrix A and vector x for which the standard computation of matrix-vector multiplication is not backward stable with respect to the input database. Proof. Consider the join matrix defined by the natural join of the relations S1 pX, Aq “ tpk, tq, pk, 0qu and S2 pX, Bq “ tpk, `1q, pk, ´1qu, where X is the join attribute and A and B are data attributes and t “ 2p`2 for p defined below. For any floating point number system [14] with base 2, each nonzero normalized floating point number can be represented by µ ˆ 2e´p , where e and µ are both integers, p is the positive integer precision, and 2p´1 ď |µ| ď 2p ´ 1.2 Let t “ 2p`2 ; t can be represented using µ “ 2p´1 and e “ p ` 3. The next floating point number has µ “ 2p´1 ` 1 and e “ p ` 3 and its value is p2p´1 ` 1q ˆ 2p`3´p “ t ` 8. The previous number has µ “ 2p ´ 1 and e “ p ` 2 and its value is p2p ´ 1q ˆ 2p`2´p “ t ´ 4. No numbers between t and t ` 8, as well as between t ´ 4 and t, exist in this floating point number system. Therefore, the integers t ` 1 and t ´ 1 will both be rounded to the closest floating point number, which is t. We therefore have: flpt ` 1q “ flpt ´ 1q “ t for t “ 2p`2 . Let x “ p1, 1qJ . With data columns ordered as pA, Bq, the join matrix A, exact multiplication Ax, and the row-wise floating-point computation ŷ “ flpAxq are: ¨ ˛ ¨ ˛ ¨ ˛ ¨ ˛ ¨ ˛ t `1 t`1 t a1 b1 a1 ` b1 ˚ t ´1‹ ˚t ´ 1‹ ˚ t ‹ ˚a1 b2‹ ˚a1 ` b2‹ 1 1 ‹ ˚ ‹ ˚ ‹ ˚ ‹ ˚ ‹ A“˚ ˝0 `1‚, Ax “ ˝ `1 ‚, ŷ “ ˝`1‚, A “ ˝a2 b1‚, A x “ ˝a2 ` b1‚. 0 ´1 ´1 ´1 a2 b2 a2 ` b2 We now show by contradiction that the computed vector ŷ cannot be obtained as the exact matrix-vector multiplication over any perturbation of the numerical attributes of the input database. Suppose there exist relations S11 and S21 that (i) are perturbations of S1 and S2 , respectively, and (ii) define the input matrices D1 “ pS11 , S12 q and the join matrix A1 , such that fx pD1 q “ ŷ. This means that A1 x “ ŷ. Since S11 and S21 are S1 and S2 , where the numerical values may be perturbed, they have the general form S11 “ tpk, a1 q, pk, a2 qu and S21 “ tpk, b1 q, pk, b2 qu for the same join key k and numerical values a1 , a2 , b1 , b2 P R. Then, the new join matrix A1 and the multiplication A1 x are as shown above. Since A1 x “ ŷ, we obtain the equalities a1 ` b1 “ t, a1 ` b2 “ t, a2 ` b1 “ `1, and a2 ` b2 “ ´1 by inspecting the first four rows of A1 x and ŷ. Subtracting the first two equations gives b1 ´ b2 “ 0, while subtracting the fourth equation from the third gives b1 ´ b2 “ 2. We reached a contradiction. Therefore, the numerical computation of Ax is not backward stable with respect to the input database. Although exemplified for two simple operations, the above observations apply to all linear algebra operations over the join matrix. The takeaway is that even though the join is trivially backward stable, since it does not touch the data values in the database, numerical operations that are backward stable with respect to the matrix may trivially lose this desirable property when analyzed with respect to the input database! So the composition of two backward stable 2

For IEEE double: p “ 53 and ´1021 ď e ď 1024. Then, t “ 255 and f lp255 ` 1q “ f lp255 ´ 1q “ 255 .

8

operations may not be backward stable in our setting. The main reason is that the join forces copies of the same input value in the join matrix, yet the numerical operations over the join matrix may introduce different numerical errors for these copies. In particular, if these operations produce an output whose size is more than the size of the input database, then backward stability is often lost since we may not be able to propagate the many independent errors to a smaller number of input data values. This is the case for workhorse linear algebra operations, such as QR. This calls for a new notion of stability for the composition of a query with a numerical operation and that allows us to quantify the numerical error introduced by the computation. We discuss this new notion in the next section.

4

Projected Backward Stability

Prop. 7 shows that the standard multiplication of the join matrix with a vector is not backward stable with respect to the input database. This failure, however, should not be interpreted as numerical instability, but it is due to the rigid structure of join matrices: their numerical perturbations may not necessarily be join matrices. Yet although perturbed join matrices may not be defined by the query over a perturbed input database, they remain close to valid join matrices. This observation motivates a generalized notion of backward stability that takes the query into account. We first define this notion for arbitrary functions f and g and then specialize it to our setting, where f is the join and g the numerical operation. Definition 8 (Projected backward stability). Let X be a normed space, Y be a finite-dimensional Hilbert space (so that projections are well-defined), and let Y0 Ď Y be a subspace. Let f : X Ñ Y0 , g : Y Ñ Z, and h “ g ˝ f . Let Y0K denote the orthogonal complement of the subspace Y0 in Y. Let C be a numerical computation of g ˝ f : if z “ gpf pxqq for some x, then ẑ “ Cpxq. We say that C is projected backward stable ` if there exists a˘ perturbation ∆x P X and a normal perturbation ∆y0K P Y0K such that ẑ “ g f px ` ∆xq ` ∆y0K , ∆x is small relative to x, and ∆y0K is small relative to f px ` ∆xq. The computation C is pα, βq-projected backward stable if for any input x, there exists ∆x and ∆y0K described above satisfying |∆x| ď α ¨ |x| and |∆y0K | ď β ¨ |f px ` ∆xq|. The standard notions of normed spaces and Hilbert spaces are recalled in Appendix A. Projected backward stability is depicted in Fig. 2. In our database setting (Def. 3), f is the join query, X is the set of matrices D “ pS1 , . . . , Sr q defined by relational databases, Y0 is the set of join matrices A such that f pDq “ A, Y Ě Y0 is the unconstrained set of all mˆn matrices, and g is the numerical operation applied to the join matrix A. The orthogonal complement Y0K captures the perturbations that break the join dependencies from the query. The normal perturbation ∆y0K records the component that cannot be represented by perturbing the database, and its magnitude is controlled separately. Notice that if the intermediate space lacks constraints (i.e., Y0 “ Y), then ∆y0K “ 0, recovering the standard notion of backward stability. C ` A computation ˘ K of g ˝ f is projected backward stable with respect to D if CpDq “ ẑ “ g à ` ∆A for some join matrix à P Y0 and perturbation ∆AK P Y0K , where ∆AK is relatively small compared to Ã. Furthermore, à must correspond to a small relative perturbation of the input database, i.e., à “ f pD ` ∆Dq, where ∆D is relatively small compared to D. We are now ready to apply this new stability lens to the matrix-vector multiplication Ax. The power of projected backward stability lies in its ability to reconcile floating-point errors with the structure of the join matrix. By projecting the inevitable rounding errors of the matrix-vector multiplication back onto the subspace of valid join matrices, the computed output maps to a small, structurally valid perturbation of the input database, as stated next. Proposition 9. Consider the setting in Def. 3, the join query f such that f pDq “ A, gx pAq “ Ax for a vector x, and the join space Y0 of A. Let C be the computation of gx ˝ f that forms A 9

2γ̃n q-projected backward stable if γ̃n ă 1. and computes Ax using dot product. Then, C is pγ̃n , 1´γ̃ n K K That is, there exist à P Y0 , ∆A P Y0 and a perturbation ∆D of D such that f pD ` ∆Dq “ Ã, 2γ̃n ¨ |Ã|, and |∆D| ď γ̃n ¨ |D|. gx pà ` ∆AK q “ CpDq, |∆AK | ď 1´γ̃ n

Example 10. We use the counterexample in the proof of Prop. 7 to show how a small normal perturbation ∆AK maps the computed output ŷ to a valid input database. Consider the query and input database from the proof of Prop. 7: the natural join of the relations S1 pX, Aq “ tpk, tq, pk, 0qu and S2 pX, Bq “ tpk, `1q, pk, ´1qu, for t “ 2p`2 in a floating point number system with base 2 and precision p. A, x “ p1, 1qJ , Ax, and ŷ “ flpAxq are as before, whereas the new quantities à and its normal perturbation ∆AK are as shown below: ˛ ¨ ˛ ¨ ˛ ¨ ˛ ¨ ˛ ¨ t `1 t`1 t t `1 ´1 0 ˚ t ´1‹ ˚ ‹ ˚ ‹ ˚ ‹ ˚ ‹ ‹ , Ax “ ˚t ´ 1‹ , ŷ “ ˚ t ‹ , à “ ˚ t ´1‹ , ∆AK “ ˚`1 0‹. A“˚ ˝0 `1‚ ˝ `1 ‚ ˝`1‚ ˝0 `1‚ ˝ 0 0‚ 0 ´1 ´1 ´1 0 ´1 0 0 In this example, A “ Ã: The entries t ` 1 and t ´ 1 in à ` ∆AK originate from the same input data value t and are replaced by their average t in Ã. This defines the perturbation ∆AK as shown above. This replacement gives the desired projection onto the space of valid join matrices (Lemma 24 in Appendix B.2). It can be checked that pà ` ∆AK qx “ ŷ. Moreover, the entry-wise perturbations of ∆AK with respect to à at indices p1, 1q and p2, 1q are 1{t=2´p´2 , well below the unit round off 2´p . Since A “ Ã, ∆D “ 0 in this example.

5

The Database Condition Number

The perturbation of the database that explains the error in the output of a numerical operation over the join matrix is a function of two fundamental phenomena: (1) the inherent numerical error due to the numerical operation; and (2) the amplification of this error due to the join that may produce many copies of the same input data value. In this section, we capture the second by the new notion of database condition number and then use it to bound the projected backward error. The amplification of an error is controlled by how many times an entry from an input matrix is copied in the join matrix. The database condition number is a property of the database and of the join query. We show it to be a natural analogue to the classical condition number of matrices and exhibit the connection between the two. We first introduce a simple notation. For a matrix Si , we fix an order on its rows so we can refer uniquely to the j-th row in Si . We denote by ϕi pjq the number of occurrences of the j-th row of Si in the join matrix. This count holds for all entries in that row. Definition 11 (Relation/Database Condition Number). Consider the setting in Def. 3. For any matrix Si with di rows in the database D “ pS1 , . . . , Sr q, the relation condition number κjoin pSi q and the database condition number κjoin pDq with respect to the join query are: d d maxjPrdi s ϕi pjq maxSi PD maxjPrdi s ϕi pjq κjoin pSi q “ κjoin pDq “ . minjPrdi s: ϕi pjqą0 ϕi pjq minSi PD minjPrdi s: ϕi pjqą0 ϕi pjq The above holds in case the join matrix is not empty. If the join matrix is empty, then the relation/database condition number is 1. The relation (database) condition number is thus the square root of the ratio of the maximum and the minimum numbers of occurrences of a row from a relation (the database) in the query result. It ranges from 1 to the square root of the size of the query result. The relation and database condition numbers are equivalent to the classical condition number of so-called expansion matrices that transform the input matrices represented by the input 10

relations and the entire database into the join matrix. To show this connection, we first introduce the expansion matrix for a relation and the database. Definition 12 (Relation/Database Expansion Matrix). Consider the setting in Def. 3. Let D “ pS1 , . . . , Sr q, where Si P Rdi ˆni , and write the join matrix ¨ Ar s, where ř as A “ rA1 ¨ ¨ ř mˆn i Ai P R consists of the columns contributed by Si , with d “ iPrrs di and n “ iPrrs ni . For each i P rrs, the relation expansion matrix Ei P t0, 1umˆdi is (for k P rms and j P rdi s): # 1, if the row k in Ai is a copy of the row j in Si , Ei rk, js “ 0, otherwise. The database expansion matrix ED is the block-diagonal matrix pE1 , . . . , Er q. By construction, Ai “ Ei Si for every i P rrs and therefore ED D is a block-diagonal representation pA1 , . . . , Ar q of A, which we call Ablock . We show that the classical condition number of the expansion matrix of each relation (database) is equal to our newly introduced relation (database) condition number. Proposition 13. Consider the setting in Def. 3. If the join matrix is non-empty, the following hold: • κpEi q “ κjoin pSi q, where Ei is the expansion matrix for input matrix Si , and • κpED q “ κjoin pDq, where ED is the expansion matrix of the database matrix D. We next illustrate the connection between the database(relation) condition number and the condition number of the expansion matrix. Example 14. Consider the join Q “ S1 pX, Y1 q ’ S2 pX, Y2 q and the join matrix A: S1 X Y1 , x1 a1 x2 a2

S2 » X Y2 1 —0 x1 b1 , E1 “ — –0 x2 b2 0 x2 b3 x2 b4

fi » fi » 0 b1 a1 „ ȷ — ffi — 1ffi a b ffi , S1 “ 1 , E2 “ I4 , S2 “ — 2 ffi , A “ —a2 –b3 fl –a2 1fl a2 1 b4 a2

fi b1 b2 ffi ffi . b3 fl b4

It can be checked that: A1 “ Ar˚, Y1 s “ E1 S1 and A2 “ Ar˚, Y2 s “ E2 S2 . Also, Ablock “ ED D, where the block-diagonal matrices are: Ablock “ pA1 , A2 q, ED “ pE1 , E2 q, and D “ pS1 , S2 q. The two rows of S1 occur once and in A,˘whereas every row of S2 occurs ` ˘ thrice, respectively, ` once in A. ?Hence ϕ1 p1q, ϕ1 p2q “ p1, 3q, ϕ p1q, . . . , ϕ 2 p4q “ p1, 1, 1, 1q. It ?2 ? follows that κjoin pS1 q “ 3, κjoin pS2 q “ 1, and κjoin pDq “ 3. We next show that κpE1 q “ 3 “ κjoin pS1 q and κpE2 q “ 1 “ κjoin pS2 q, so the relation condition numbers are the condition numbers of the expansion matrices E1 and E2 . We first show how to compute the condition number of E1 . 2 J J Using the SVD decomposition E1 “ U1 Σ1 V1J , we have EJ 1 E1 “ V1 Σ1 V1 . Since E1 E1 is the 2 diagonal matrix diagp1, 3q in our example, we have V1 “ I2 and Σ1 “ diagp1, 3q. So 3 (1) is the max (min) singular value of EJ 1 E1 and its square root is the max (min) singular value of ? E1 . So κpE1 q “ 3. Since E2 “ I4 , the singular values of E2 are all 1,?so κpE2 q “ 1. Moreover, ED “ pE1 , E2 q and EJ 3 “ κjoin pDq. D ED “ diagp1, 3, 1, 1, 1, 1q. Therefore, κpED q “ The relation (database) condition number translates the column-wise (Frobenius-norm-wise) perturbation bound from the join space to the input database, as shown next. Proposition 15. Consider our setting in Def. 3 and let ∆A be a perturbation of the join matrix A on the join space. Then, there exists a perturbation ∆D “ p∆S1 , . . . , ∆Sr q of the input D such that A ` ∆A is the join matrix of the perturbed input D ` ∆D. Moreover, for every γ ě 0, the following statements hold: 11

A

Cg

}A ´ Ã}F ď ρ ¨ }A}F

}Ã ` ∆AK ´ A}F

ẑ “ Cg pAq “ gpà ` ∆AK q g

Ã

à ` ∆AK

}∆AK }F ď β ¨ }Ã}F ď β ¨ p1 ` ρq ¨ }A}F

Figure 3: Geometric interpretation of the Frobenius-norm-wise backward-error bound of Cg in Theorem 16. We let ρ “ κjoin pDq ¨ α. The quantities }A ´ Ã}F and }∆AK }F are the lengths of the perpendicular legs of a right triangle, whose hypotenuse length is the Frobenius norm of the total perturbation E “ Ã ` ∆AK ´ A. • For every j P rrs and Y P Yj , if ∥∆Ar˚, Y s∥2 ď γ ¨ ∥Ar˚, Y s∥2 then ∥∆Sj r˚, Y s∥2 ď κjoin pSj q ¨ γ ¨ ∥Sj r˚, Y s∥2 . • If ∥∆A∥F ď γ ¨ ∥A∥F , then ∥∆D∥F ď κjoin pDq ¨ γ ¨ ∥D∥F . The database condition number can be seen as a measure of the data skew in the query result (and the join matrix). A value of one means a well-conditioned database, where the input entries that contribute to the output are represented uniformly, and no amplification of relative perturbations.

Relationship of Projected Backward Stability and Backward Stability

6

There is a strong relationship between the projected backward stability of the computation C of g ˝ f with respect to the input database, where f is the join query and g is the numerical computation, and the backward stability of the computation Cg of g with respect to the join matrix. In particular, the projected backward stability of C implies the backward stability of Cg (Theorem 16), while the backward stability of Cg implies the projected backward stability of C (Theorem 17). Theorem 16. Consider the setting in Def. 3 and let C be the computation of the join query followed by the computation Cg of a numerical operation over the join matrix A. For α, β ě 0, the following hold: ` ˘ • If C is entry-wise pα, βq-projected backward stable, then Cg is entry-wise α ` β ¨ p1 ` αq backward stable on the subspace Y0 ; • If is column-wise pα, βq-projected backward stable, then Cg is column-wise bC ` ˘2 ` ˘2 κjoin pSi q ¨ α ` β 2 ¨ 1 ` κjoin pSi q ¨ α -backward stable on the subspace Y0 for each input data column in Si and i P rrs; • If C b is Frobenius-norm-wise pα, βq-projected backward stable, then Cg is Frobenius-norm` ˘2 ` ˘2 κjoin pDq ¨ α ` β 2 ¨ 1 ` κjoin pDq ¨ α -backward stable on the subspace Y0 . wise The Frobenius-norm-wise bound in Theorem 16 admits a geometric explanation (Fig. 3). Let ρ “ κjoin pDq ¨ α, where α satisfies }∆D}F ď α ¨ }D}F . Lemma 27 and the normal-error bound β give3 }à ´ A}F ď ρ ¨ }A}F and }∆AK }F ď β ¨ }Ã}F ď β ¨ p1 ` ρq ¨ }A}F . These two perturbations are orthogonal to each other because à ´ A P Y0 and ∆AK P Y0K . Their sum 3

We assume here that the input database is fully reduced (Lemma 27 applies to fully reduced databases). Theorem 16 holds for uncalibrated databases, too (cf. its proof in Appendix B.4).

12

E “ Ã ` ∆AK ´ A therefore represents the hypotenuse of`the right triangle˘ in Fig. 3. By the Pythagorean theorem, }E}2F “ }A ´ Ã}2F ` }∆AK }2F ď ρ2 ` β 2 ¨ p1 ` ρq2 ¨ }A}2F . Since a Cg pAq “a gpA ` Eq and }E}F ď ρ2 ` β 2 ¨ p1 ` ρq2 ¨ }A}F , the relative backward-error bound of Cg is ρ2 ` β 2 ¨ p1 ` ρq2 . The column-wise error bound follows from the same geometric argument. Theorem 17. Consider the setting in Def. 3 and let C be the computation of the join query followed by the computation Cg of a numerical operation over the join matrix A. For 0 ď ϵ ă 1, the following hold: 2ϵ • If Cg is entry-wise ϵ-backward stable, then C is entry-wise pϵ, 1´ϵ q-projected backward stable. ϵ • If Cg is column-wise ϵ-backward stable, then C is column-wise pκjoin pSi q ¨ ϵ, ?1´ϵ 2 q-projected backward stable for each input data column in Si and i P rrs;

• If Cg is Frobenius-norm-wise ϵ-backward stable, then C is Frobenius-norm-wise ϵ pκjoin pDq ¨ ϵ, ?1´ϵ 2 q-projected backward stable. The condition ϵ ă 1 in Theorem 17 is necessary. We explain this next for the Frobeniusnorm-wise bound by showing that when ϵ ě 1, there can be no finite bound β for the normal perturbation. Assume A P Y0 zt0u and Y0K ‰ t0u. Further, assume Cg is ϵ-backward stable for ? K K K 2 ϵ ą 1. Consider the scenario where à “ 0 and ∆A P Y0 with }∆A }F “ ϵ ´ 1 ¨ }A}F . By the Pythagorean Theorem, }∆AK ´ A}2F “ }A}2F ` }∆AK }2F “ ϵ2 ¨ }A}2F (the geometry is the same as shown in Figure 3, with the extra condition that à “ 0). Denote E “ à ` ∆AK ´ A and we have Cg pAq “ gpA ` Eq, with }E}F “ ϵ ¨ }A}F , satisfying the backward error bound of ϵ for Cg . But the projected backward stability of C requires a finite β satisfying }∆AK }F ď β ¨ }Ã}F , which cannot be fulfilled in our scenario, since the left-hand side is positive and the righthand side is? zero. If ϵ “ 1, consider the scenario where à “ tA and ∆AK P Y0K with }∆AK }F “ 2t ´ t2 }A}F`, where 0 ă t ă 1.˘By the Pythagorean Theorem, }à ` ∆AK ´ A}2F “ }∆AK }2F ` }A ´ Ã}2F “ p1 ´ tq2 ` 2t ´ t2 ¨ }A}2F “ }A}2F . Denote E “ à ` ∆AK ´ A and we have Cg pAq “ gpA b ` Eq, with }E}F “ }A}F , satisfying the backward error bound of 1 for K

}F “ 2t ´ 1 Ñ 8 as t Ñ 0. Thus, there is no finite bound β for the normal Cg . But }∆A }Ã}F perturbation for ϵ “ 1. Therefore, the condition ϵ ă 1 is necessary in Theorem 17. Prop. 9 is a direct corollary of Theorem 17. As the matrix vector product Ax is entry-wise γ̃n -backward stable with respect to A [14, Sec. 3.5, Eq. (3.11)], Theorem 17 states that the entire 2γ̃n computation of the join matrix A followed by Ax is entry-wise pγ̃n , 1´γ̃ q-projected backward n stable with respect to the input database.

7

Numerical Stability of QR and SVD

We now investigate the stability of the QR and SVD decompositions. In the standard setting, so for arbitrary input matrices, there are backward stable algorithms for these decompositions. In our database setting with structurally restricted join matrices, backward stability fails while projected backward stability holds. We first discuss the QR and then the SVD. Assume we want to compute the upper triangular matrix R̂ from the QR decomposition of a matrix A P Rmˆn using Givens Rotations [13]. Then, there exists an orthogonal matrix Q1 and a perturbation matrix ∆A, such that A ` ∆A “ Q1 R̂, where ∆A satisfies the column-wise bounds }∆Ar˚, js}2 ď γ̃m`n´2 }Ar˚, js}2 , for j P rns [14]. Backward stability with respect to the join matrix does not generally imply backward stability with respect to the input database. Prop. 28 in Appendix B.5 establishes this failure for a common algorithm for QR decomposition. Yet Theorem 17 together with the backward stability result for QR decomposition [14] implies projected backward stability. 13

Corollary 18. [Theorem 17 and [14]] Consider the setting in Def. 3 and let C be the computation that materializes the join matrix A and then applies Givens Rotations to compute the uppertriangular matrix R from the QR decomposition of A. If γ̃m`n´2 ă 1, then C is column-wise pκjoin pSj q ¨ γ̃m`n´2 , b γ̃m`n´2 q-projected backward stable for each input data column in Sj and 2 1´γ̃m`n´2

j P rrs. A similar guarantee holds for computing the singular values. LAPACK’s xGESVD returns the exact singular values of A ` ∆A, where }∆A}2 ď ppm, nq ¨ u ¨ }A}2 , ppm, nq is a modestly growing dimension-dependent factor, and u is the unit round-off [2, Sec. 4.9]. Since }∆A}F ď a minpm, nq ¨ }∆A}2 and }A}2 ď }A}F , xGESVD is Frobenius-norm-wise ϵsvd pm, nq-backward def a stable, where ϵsvd pm, nq “ minpm, nq ¨ ppm, nq ¨ u. Combining this classical result with Theorem 17, we obtain the following corollary. Corollary 19. [Theorem 17 and [2, Sec. 4.9]] Consider the setting in Def. 3 and let C be the computation that materializes the join matrix A and then applies xGESVD to compute Σ, the diagonal matrix containing ˙ singular values of A. If ϵsvd ă 1, then C is Frobenius-norm-wise ˆ κjoin pDq ¨ ϵsvd , ? ϵsvd2

1´ϵsvd

-projected backward stable.

Recall that pushing the computation of the Frobenius norm past the join down to the input relations improves the numerical stability by avoiding repeated arithmetic on copied values (Prop. 6). We establish a similar benefit for QR using the FiGaRo algorithm [27], which pushes the computation of the upper-triangular matrix R past the α-acyclic joins and thereby avoids the materialization of the join matrix A. In particular, FiGaRo computes an almost upper-triangular matrix R0 such that A “ Q0 R0 , where Q0 is some orthonormal matrix defined by a sequence of Givens rotations. The matrix R0 can then be transformed into the upper-triangular matrix R using, e.g., the Givens Rotations algorithm for QR decomposition. Note also that the singular values of R0 are precisely those of A. By pushing the linear operation past the join we obtain a much stronger numerical stability, where the error depends on the number d of rows in the input database instead of the number ˚ m of rows in the query result (and the join matrix). Note that m “ Opdρ pQq q for any input database and there are databases for which this is asymptotically tight. Here, ρ˚ pQq is the fractional edge cover number of Q [3]. It can be as large as the number r of relations in Q. Since d can easily be in the order of hundreds of millions and r in the order of 5-10, the error reduction from m to d brought by FiGaRo can be significant. Definition 20 (FiGaRo setting). In addition to the database setting from Def. 3, in the FiGaRo setting we require calibrated databases, α-acyclic queries, and given join trees as input to the FiGaRo algorithm. Theorem 21. Consider the FiGaRo setting in Def. 20 and let C be the computation that uses FiGaRo on the input database followed by Givens Rotations to compute ` ˘the upper-triangular matrix R from the QR decomposition of A. Let η “ γ̃ ` γ̃ ¨ 1 ` γ̃ d d`n´2 d . If η ă 1, then C is ´ ¯ η column-wise κjoin pSj q ¨ η, ? -projected backward stable for each input data column in Sj 1´η 2

and j P rrs. FiGaRo can improve the numerical stability of computing the singular values, too. Theorem 22. Consider the FiGaRo setting in Def. 20 and let C be the computation that uses FiGaRo on the input database followed by xGESVD to compute the singular values of A. Let ´ ¯ ` ˘ η “ γ̃d ` ϵsvd pd, nq ¨ 1 ` γ̃d . If η ă 1, then C is Frobenius-norm-wise κjoin pDq ¨ η, ? η 2 1´η

projected backward stable. 14

8

Conclusion

This paper initiates the study of numerical stability for computations that combine relational algebra and linear algebra. It introduces the fundamental concepts of projected backward stability and database condition number to allow making formal statements about the stability of numerical operations over matrices defined by database join queries. The framework introduced in this paper is a fertile ground for future work. We mention below a few directions: • Analyze the numerical stability in the context of more complex queries, e.g., different join types and even self-joins. Every column-wise error bound in the formal statements of this paper fails in the presence of self-joins. The presence of self-joins does not make a difference, however, for our entry-wise and Frobenius-norm-wise error bounds. • Determine tight or optimal bounds for projected backward stability. Our current analysis establishes upper bounds on the database and orthogonal perturbations, but does not show whether these bounds are tight or characterize the best achievable pair of bounds. • Extend the framework to analytical pipelines, e.g., series-parallel graphs whose nodes are relational and linear algebra operators. How do the structural and numerical perturbations compose across such pipelines? Develop a syntax-directed, compositional method that derives global error bounds from the operators in the pipeline, while avoiding a numerical analysis for every combination of relational and linear-algebraic operations. • Our work examines numerical stability with respect to perturbations of data values. An alternative view is that it investigates the sensitivity of (relational and linear) algebraic data transformations: How much the input data must change to account for a change in the output of a given scale. This perspective can inform what-if analysis of computational pipelines: understanding whether a desired change in the output can be realized through an input perturbation, and how large this perturbation must be. This view relates to prior work [32] that studies how varying query parameters on a fixed database affects query results, with applications to computational fact-checking.

References [1] Serge Abiteboul, Richard Hull, and Victor Vianu. Foundations of Databases. Addison-Wesley, 1995. URL: http://webdam.inria.fr/Alice/. [2] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, 3 edition, 1999. URL: https://www. netlib.org/lapack/lug/. [3] Albert Atserias, Martin Grohe, and Dániel Marx. Size bounds and query plans for relational joins. In FOCS, pages 739–748, 2008. doi:10.1109/FOCS.2008.43. [4] Pablo Barceló, Nelson Higuera, Jorge Pérez, and Bernardo Subercaseaux. On the expressiveness of lara: A proposal for unifying linear and relational algebra. Theor. Comput. Sci., 935:105–127, 2022. doi:10.1016/J.TCS.2022.09.003. [5] Carlos Beltrán, Vanni Noferini, and Nick Vannieuwenhoven. When can forward stable algorithms be composed stably? IMA Journal of Numerical Analysis, 44(2):886–919, 03 2024. doi:10.1093/imanum/drad026.

15

[6] Robert Brijder, Floris Geerts, Jan Van den Bussche, and Timmy Weerwag. MATLANG: matrix operations and their expressive power. SIGMOD Rec., 48(1):60–67, 2019. doi: 10.1145/3371316.3371331. [7] Robert Brijder, Floris Geerts, Jan Van den Bussche, and Timmy Weerwag. On the expressive power of query languages for matrices. ACM Trans. Database Syst., 44(4):15:1–15:31, 2019. doi:10.1145/3331445. [8] Stanislav Budzinskiy, Wenyi Fang, Longbin Zeng, and Philipp Petersen. Numerical stability analysis of large language models, 2026. arXiv:2503.10251. [9] Lingjiao Chen, Arun Kumar, Jeffrey F. Naughton, and Jignesh M. Patel. Towards linear algebra over normalized data. Proc. VLDB Endow., 10(11), 2017. [10] Xinye Chen. Automated numerical stability analysis of deep learning operators, 2026. arXiv:2607.25494. [11] Kyle Deeds, Willow Ahrens, Magdalena Balazinska, and Dan Suciu. Galley: Modern query optimization for sparse tensor programs. Proc. ACM Manag. Data, 3(3):164:1–164:24, 2025. doi:10.1145/3725301. [12] James W. Demmel. Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics, 1997. doi:10.1137/1.9781611971446. [13] Wallace Givens. Computation of plain unitary rotations transforming a general matrix to triangular form. J. of the Soc. for Industrial and Applied Maths, 6(1):26–50, 1958. doi:10.1137/0106004. [14] Nicholas J. Higham. Accuracy and stability of numerical algorithms, Second Edition. SIAM, 2002. doi:10.1137/1.9780898718027. [15] Dylan Hutchison, Bill Howe, and Dan Suciu. Laradb: A minimalist kernel for linear and relational algebra computation. In BeyondMR@SIGMOD, pages 2:1–2:10, 2017. doi: 10.1145/3070607.3070608. [16] Niranjan Kamat and Arnab Nandi. A closer look at variance implementations in modern database systems. SIGMOD Rec., 45(4):28–33, 2016. doi:10.1145/3092931.3092936. [17] Kelvin Kan, Xingjian Li, Benjamin J. Zhang, Tuhin Sahai, Stanley Osher, Krishna Kumar, and Markos A. Katsoulakis. Stability of transformers under layer normalization, 2026. arXiv:2510.09904. [18] Ariel E. Kellison, Laura Zielinski, David Bindel, and Justin Hsu. Bean: A language for backward error analysis. Proc. ACM Program. Lang., 9(PLDI), June 2025. doi: 10.1145/3729324. [19] Mahmoud Abo Khamis, Xiao Hu, and Dan Suciu. Fast matrix multiplication meets the submodular width. Proc. ACM Manag. Data, 3(2):98:1–98:26, 2025. doi:10.1145/3725235. [20] Mahmoud Abo Khamis, Hung Q. Ngo, Xuanlong Nguyen, Dan Olteanu, and Maximilian Schleich. Learning models over relational data using sparse tensors and functional dependencies. ACM Trans. Database Syst., 45(2), June 2020. doi:10.1145/3375661. [21] Eliska Kloberdanz, Kyle G. Kloberdanz, and Wei Le. DeepStability: A study of unstable numerical methods and their solutions in deep learning. In ICSE, pages 586–597, 2022. doi:10.1145/3510003.3510095.

16

[22] Christoph Koch and Dan Olteanu. Conditioning probabilistic databases. Proc. VLDB Endow., 1(1):313–325, 2008. URL: http://www.vldb.org/pvldb/vol1/1453894.pdf, doi: 10.14778/1453856.1453894. [23] Arun Kumar, Jeffrey Naughton, and Jignesh M. Patel. Learning generalized linear models over normalized data. In SIGMOD, page 1969–1984, 2015. doi:10.1145/2723372.2723713. [24] Ingo Müller, Andrea Arteaga, Torsten Hoefler, and Gustavo Alonso. Reproducible floatingpoint aggregation in RDBMSs. In ICDE, pages 1049–1060, 2018. doi:10.1109/ICDE.2018. 00098. [25] Thomas Muñoz, Cristian Riveros, and Stijn Vansummeren. Enumeration and updates for conjunctive linear algebra queries through expressibility. Log. Methods Comput. Sci., 22(1), 2026. doi:10.46298/LMCS-22(1:1)2025. [26] Dan Olteanu. The relational data borg is learning. Proc. VLDB Endow., 13(12):3502–3515, 2020. doi:10.14778/3415478.3415572. [27] Dan Olteanu, Nils Vortmeier, and Djorđe Živanović. Givens rotations for QR decomposition, SVD and PCA over database joins. VLDB J., 33(4):1013–1037, 2024. doi:10.1007/ s00778-023-00818-9. [28] Maximilian Schleich, Dan Olteanu, and Radu Ciucanu. Learning linear regression models over factorized joins. In SIGMOD, pages 3–18, 2016. doi:10.1145/2882903.2882939. [29] Amir Shaikhha, Mathieu Huot, Jaclyn Smith, and Dan Olteanu. Functional collection programming with semi-ring dictionaries. Proc. ACM Program. Lang., 6(OOPSLA1):1–33, 2022. doi:10.1145/3527333. [30] Amir Shaikhha, Maximilian Schleich, Alexandru Ghita, and Dan Olteanu. Multi-layer optimizations for end-to-end data analytics. In CGO, pages 145–157, 2020. doi:10.1145/ 3368826.3377923. [31] Shaila Sharmin, Anwar Hossain Zahid, Subhankar Bhattacharjee, Chiamaka Igwilo, Miryung Kim, and Wei Le. Automatically detecting numerical instability in machine learning applications via soft assertions. Proc. ACM Softw. Eng., 2(FSE), June 2025. doi:10.1145/ 3729394. [32] Brett Walenz and Jun Yang. Perturbation analysis of database queries. Proc. VLDB Endow., 9(14):1635–1646, October 2016. doi:10.14778/3007328.3007330.

17

This appendix complements the main paper with further preliminaries and proofs of the formal statements.

A

Further Preliminaries

We review the notions of normed and Hilbert spaces used in our definition of projected backward stability (Definition 8). A normed space is a vector space X equipped with a norm } ¨ } : X Ñ Rě0 satisfying, for all x, y P X and α P R, }x} “ 0 ðñ x “ 0,

}x ` y} ď }x} ` }y}.

}αx} “ |α| }x},

The norm provides a notion of magnitude, and hence allows us to quantify the size of perturbations such as ∆x through quantities such as }∆x}{}x}. A finite-dimensional Hilbert space is a finite-dimensional vector space Y equipped with an inner product x¨, ¨y. The inner product induces a norm a }y} “ xy, yy. The inner product also defines orthogonality: two vectors y1 , y2 P Y are orthogonal if xy1 , y2 y “ 0. For a subspace Y0 Ď Y, its orthogonal complement is Y0K “ ty P Y : xy, y0 y “ 0 for all y0 P Y0 u. Because Y is finite-dimensional, every subspace Y0 is closed and every y P Y has a unique orthogonal decomposition y “ y0 ` y0K ,

y0 P Y0 ,

y0K P Y0K .

Equivalently, there is a unique orthogonal projection PY0 : Y Ñ Y0 such that pI ´ PY0 qy “ y0K .

PY0 y “ y0 ,

Thus, in our setting, the finite-dimensional Hilbert space structure provides both a well-defined notion of magnitude through the norm and a well-defined notion of orthogonal projection, which we use to separate the component in Y0 from the off-range component in Y0K .

B

Proofs of Statements from the Main Body of the Paper

B.1

Proofs of Statements from Section 3

Here we provide the proofs of the formal statements in Section 3. Proposition 6. Consider the setting in Def. 3. The computation C of an operation f , which pushes the computation of the Frobenius norm of the join matrix past the join, is backward stable with respect to the input database. That is, there exists a perturbation ∆D such that, if ẑ “ CpDq, then ẑ “ f pD ` ∆Dq. Moreover, |∆D| ď γ̃ℓ ¨ |D|. Proof. Let z “ f pDq “ }A}F and ẑ “ CpDq. Compute the multiplicity ϕk piq of each row i of Sk in the join result. Then g fÿ f r ÿÿ z“e ϕk piqSk ri, js2 . k“1 i

j

18

Once the multiplicities have been computed exactly, evaluating this expression requires Opℓq arithmetic operations. Since all terms are nonnegative and the square root is forward stable, the standard rounding model gives |ẑ ´ z| ď γ̃ℓ z. If z “ 0, then ẑ “ 0, and we set ∆D “ 0. Otherwise, we write ẑ ´ z “ θz, where |θ| ď γ̃ℓ . Since ẑ ě 0, we have 1 ` θ “ ẑ{z ě 0. Define the database perturbation ∆D “ θD. This scales every numerical data entry by θ to form its perturbation, while preserving the join keys, row counts, and row order. Thus D ` ∆D “ p1 ` θqD on the numerical data entries. The join structure and all multiplicities ϕk piq remain unchanged. Consequently, g fÿ f r ÿÿ f pD ` ∆Dq “ e ϕk piqp1 ` θq2 Sk ri, js2 “ |1 ` θ|z “ p1 ` θqz “ ẑ. k“1 i

j

Finally, entry-wise, |∆D| “ |θ| |D| ď γ̃ℓ ¨ |D|. This proves backward stability with respect to the input database.

B.2

Proofs of Statements from Section 4

Here we provide the proofs of the formal statements in Section 4. Lemma 23. Consider our setting in Definition 3 and let A be the query matrix, and Y0 Ď Rmˆn be the query space. Then Y0 is a linear subspace of Rmˆn . Proof. With the row labels, key columns, and join structure fixed, every entry of a join matrix is a copy of a fixed numerical entry of one of the input relations in the database. Hence the map from the numerical entries of S1 , . . . , Sr to the materialized join matrix is linear. Therefore, if A, B P Y0 and α, β P R, then αA ` βB is the join matrix obtained by applying the same linear combination to the corresponding numerical entries of the input database. Thus αA ` βB P Y0 , and Y0 is a linear subspace. Lemma 24. Consider our setting in Definition 3 and let D be an input database and A be the query matrix, Y0 be the join space, and Y “ Rmˆn . For any à ` ∆AK P Y with à P Y0 and ∆AK P Y0K , its projection on Y0 , Ã, satisfies the following property: For any fixed entry s in a relation in the input database D that appears in the join, denote Ps “ tpi, jq|Ari, js is a copy of su as the set of indices of entries ř in A parameterized by the same input entry s, then, we have: for any pi1 , j 1 q P Ps , Ãri1 , j 1 s “

K pi,jqPPs pÃ`∆A qri,js |Ps |

Proof. Fix an entry s in a relation in the input database D, and define the matrix Ms P Rmˆn by # 1, pi, jq P Ps , Ms ri, js “ 0, pi, jq R Ps . Because the join keys are fixed and the data entry s is a free parameter, perturbing s by any scalar changes precisely the entries of the query matrix indexed by Ps by that scalar. Consequently, Ms P Y0 . Moreover, since all entries indexed by Ps are parameterized by the same input entry s, every matrix in Y0 has the same value at all indices in Ps . In particular, there exists a scalar cs such that Ãri, js “ cs for every pi, jq P Ps . Since ∆AK P Y0K and Ms P Y0 , ∆AK is orthogonal to Ms and their inner product is 0, so we have ÿ D @ 0 “ ∆AK , Ms F “ ∆AK ri, js. pi,jqPPs

19

It follows that ÿ ` ˘ 0“ p∆AK ` Ãqri, js ´ Ãri, js “ pi,jqPPs

ÿ

p∆AK ` Ãqri, js ´ |Ps |cs .

pi,jqPPs

Therefore, pi,jqPPs pà ` ∆A

ř cs “

K qri, js

|Ps |

.

Thus, for every pi1 , j 1 q P Ps , ř K pi,jqPPs pà ` ∆A qri, js 1 1 . Ãri , j s “ |Ps | Since s was arbitrary, the claimed property holds for every input entry s. Proposition 9. Consider the setting in Def. 3, the join query f such that f pDq “ A, gx pAq “ Ax for a vector x, and the join space Y0 of A. Let C be the computation of gx ˝ f that forms A 2γ̃n and computes Ax using dot product. Then, C is pγ̃n , 1´γ̃ q-projected backward stable if γ̃n ă 1. n K K That is, there exist à P Y0 , ∆A P Y0 and a perturbation ∆D of D such that f pD ` ∆Dq “ Ã, 2γ̃n gx pà ` ∆AK q “ CpDq, |∆AK | ď 1´γ̃ ¨ |Ã|, and |∆D| ď γ̃n ¨ |D|. n Proof. For each ř i P rms, the computed entry ŷi is obtained by evaluating the length-n dot product Ai,: x “ nj“1 Aij xj in floating-point arithmetic. Under the standard model flpa op bq “ pa op bqp1 ` δq,

|δ| ď u,

and assuming that no overflow or exceptional underflow occurs, the standard dot-product error analysis shows that there exist numbers θij satisfying |θij | ď γ̃n such that ŷi “

n ÿ

Aij p1 ` θij qxj .

j“1

Define the perturbation matrix E by Eij :“ θij Aij . Then ŷ “ pA ` Eqx, and |E| ď γ̃n |A|. This means that the computation of the matrix vector product Ax is γ̃n -backward stable with respect to A. Using Theorem 17, we immediately obtain: the computation for gx pf pDqq is 2γ̃n pγ̃n , 1´γ̃ q-projected backward stable. This completes the proof. n

B.3

Proofs of Statements from Section 5

Here we provide the proofs of the formal statements in Section 5. Proposition 13. Consider the setting in Def. 3. If the join matrix is non-empty, the following hold: • κpEi q “ κjoin pSi q, where Ei is the expansion matrix for input matrix Si , and • κpED q “ κjoin pDq, where ED is the expansion matrix of the database matrix D. Proof. Consider a relation expansion matrix Ei . Its j-th column contains exactly ϕi pjq ones, so the inner product of this column with itself is ϕi pjq. Two distinct columns have inner product zero: each join row uses exactly one tuple` from Si , so each˘ row of Ei contains exactly one entry J equal to 1. Consequently, EJ i Ei “ diag ϕi p1q, . . . , ϕi pdi q . Let Ei “ Ui Σi Vi be a full SVD. J Since Ui Ui “ Im , we obtain J J EJ i Ei “ Vi pΣi Σi qVi .

20

The matrix ΣJ i Σi is diagonal and nonnegative, with positive diagonal entries equal to the squares of the positive singular values of Ei . Since Vi is orthogonal, this is an SVD of EJ i Ei . On the other hand, EJ E is itself diagonal and nonnegative, so its singular values are its diagonal i i entries, upato ordering. Comparing these two descriptions, the positive singular values of Ei are precisely ϕi pjq for j P rdi s with ϕi pjq ą 0. Dangling tuples contribute zero diagonal entries, which do not enter the positive-singular-value ratio. Since the join is nonempty, at least one tuple of Si has positive multiplicity. Thus, c c ` σmax pEi q “ max ϕi pjq, σmin pEi q “ min ϕi pjq. jPrdi s

jPrdi s: ϕi pjqą0

` Here σmin pEi q denotes the smallest positive singular value of Ei . Taking their ratio gives d maxjPrdi s ϕi pjq σmax pEi q κpEi q “ ` “ “ κjoin pSi q, minjPrdi s: ϕi pjqą0 ϕi pjq σmin pEi q

where the last equality is the definition of the relation condition number. For the database expansion matrix, block-diagonal multiplication gives ` J ˘ J EJ D E D “ E1 E1 , . . . , E r Er . Each block is diagonal, so the entire matrix is diagonal, with entries ϕi pjq over all i P rrs and J , we also have j P rdi s. Writing a full SVD as ED “ UD ΣD VD J J EJ D ED “ VD pΣD ΣD qVD .

As above, this is an SVD of EJ D ED , whose positive singular values are therefore the squares of the positive singular values of ED . Since the same matrix is diagonal and nonnegative, these squared singular values are precisely the positive tuple multiplicities across the database. Consequently, d maxiPrrs maxjPrdi s ϕi pjq κpED q “ “ κjoin pDq, miniPrrs minjPrdi s: ϕi pjqą0 ϕi pjq as required. Lemma 25. Let s “ ps1 , . . . , sn q ‰ 0, and form the vector J a “ psloooomoooon s2 , . . . , s2 , . . . , sloooomoooon 1 , . . . , s1 , loooomoooon n , . . . , sn q , λ1 times

λ2 times

λn times

where the positive integers λi satisfy 1 ď λ1 ď λ2 ď ¨ ¨ ¨ ď λn . For a perturbation ∆s “ p∆s1 , . . . , ∆sn qJ , define the corresponding repeated-entry perturbation J ∆a “ p∆s ∆s2 , . . . , ∆s2 , . . . , ∆s 1 , . . . , ∆s1 , looooooomooooooon n , . . . , ∆sn q . looooooomooooooon looooooomooooooon λ1 times

λ2 times

λn times

Then c

λ1 }∆a}2 }∆s}2 ¨ ď ď λn }a}2 }s}2

c

λn }∆a}2 ¨ . λ1 }a}2

21

Proof. By the definitions of a and ∆a, }a}22 “

n ÿ

λi |si |2 ,

}∆a}22 “

i“1

n ÿ

λi |∆si |2 .

i“1

Since λ1 ď λi ď λn for every i, λ1 }∆s}22 ď }∆a}22 ď λn }∆s}22 and λ1 }s}22 ď }a}22 ď λn }s}22 . Therefore, λ1 }∆a}22 }∆s}22 λn }∆a}22 ¨ ¨ ď ď . 2 2 λn }a}2 λ1 }a}22 }s}2 Taking square roots gives the result. Now we provide two lemmas that translate the perturbation bounds between the input database and the join matrix in calibrated databases. We then generalize the results to any arbitrary database. Lemma 26. Consider the setting in Def. 3 with the condition that the input database is fully reduced and let A ` ∆A be the join matrix obtained from the perturbed input database D ` ∆D “ pS1 ` ∆S1 , . . . , Sr ` ∆Sr q. For every j P rrs, Y P Yj , and γ, γ 1 ě 0, the following statements hold: • If ∥∆Ar˚, Y s∥2 ď γ ¨ ∥Ar˚, Y s∥2 , then ∥∆Sj r˚, Y s∥2 ď κjoin pSj q ¨ γ ¨ ∥Sj r˚, Y s∥2 . • If ∥∆Sj r˚, Y s∥2 ď γ 1 ¨ ∥Sj r˚, Y s∥2 , then ∥∆Ar˚, Y s∥2 ď κjoin pSj q ¨ γ 1 ¨ ∥Ar˚, Y s∥2 . Proof. Let sp “ Sj rp, Y s. In the column Ar˚, Y s, the entry sp occurs exactly λp “ ϕj ppq times. Since the relations are fully reduced, λp ą 0 for every p P rdj s where dj denotes the number of rows in Sj . Moreover, because the perturbations preserve the key columns and the join structure, ∆Sj rp, Y s is repeated exactly λp times in ∆Ar˚, Y s. If Sj r˚, Y s “ 0, then Ar˚, Y s “ 0. In the first case, the assumption implies ∆Ar˚, Y s “ 0 and hence ∆Sj r˚, Y s “ 0. In the second case, the assumption implies ∆Sj r˚, Y s “ 0 and hence ∆Ar˚, Y s “ 0. Thus, both claims hold. We may therefore assume that Sj r˚, Y s ‰ 0. Applying Lemma 25 gives d d minpPrdj s λp }∆Ar˚, Y s}2 maxpPrdj s λp }∆Ar˚, Y s}2 }∆Sj r˚, Y s}2 ¨ ď ď ¨ . maxpPrdj s λp }Ar˚, Y s}2 }Sj r˚, Y s}2 minpPrdj s λp }Ar˚, Y s}2 By the definition of the relation condition number and using the fact that the database is fully reduced and hence the multiplicity ϕj ppq ą 0 for every row p in relation Sj , d d maxpPrdj s λp maxpPrdj s ϕj ppq “ “ κjoin pSj q. minpPrdj s λp minpPrdj s ϕj ppq Consequently, 1 }∆Ar˚, Y s}2 }∆Sj r˚, Y s}2 }∆Ar˚, Y s}2 ¨ ď ď κjoin pSj q ¨ . κjoin pSj q }Ar˚, Y s}2 }Sj r˚, Y s}2 }Ar˚, Y s}2 22

If }∆Ar˚, Y s}2 ď γ}Ar˚, Y s}2 , then the right-hand inequality gives }∆Sj r˚, Y s}2 ď κjoin pSj q ¨ γ ¨ }Sj r˚, Y s}2 . Similarly, if }∆Sj r˚, Y s}2 ď γ 1 ¨ }Sj r˚, Y s}2 , then the left-hand inequality gives }∆Ar˚, Y s}2 ď κjoin pSj q ¨ γ 1 ¨ }Ar˚, Y s}2 . This proves both claims. Lemma 27. Consider the setting in Def. 3 with the condition that the input database is fully reduced and let A ` ∆A be the join matrix obtained from the perturbed input database D ` ∆D “ pS1 ` ∆S1 , . . . , Sr ` ∆Sr q. For γ, γ 1 ě 0, • If }∆A}F ď γ ¨ }A}F , then }∆D}F ď κjoin pDq ¨ γ ¨ }D}F . • If }∆D}F ď γ 1 ¨ }D}F , then }∆A}F ď κjoin pDq ¨ γ 1 ¨ }A}F . Proof. If D “ 0, the join matrix A “ 0. Conversely, if the join matrix A “ 0, the input database D “ 0 since it is fully reduced. In both cases, the statements trivially hold. We consider cases where A ‰ 0 and D ‰ 0. For the p-th row in Sj , each of its data entries occurs in A exactly λj,p “ ϕj ppq times. Since the relations are fully reduced, all such multiplicities are positive. Moreover, because the perturbations preserve the key values and the join structure, the corresponding perturbed entries have the same multiplicities in ∆A. Therefore, the minimum and maximum parameter multiplicities are λmin “ min min ϕj ppq, jPrrs pPrdj s

λmax “ max max ϕj ppq. jPrrs pPrdj s

Applying Lemma 25 gives bř c c r 2 j“1 }∆Sj }F λmin }∆A}F λmax }∆A}F b ¨ ď ď ¨ . ř r λmax }A}F λmin }A}F }Sj }2 j“1

Since

b

F

λmax λmin “ κjoin pDq, we obtain

1 }∆A}F }∆D}F }∆A}F ¨ ď ď κjoin pDq ¨ . κjoin pDq }A}F }D}F }A}F If }∆A}F ď γ ¨ }A}F , then the right-hand inequality gives }∆D}F ď κjoin pDq ¨ γ ¨ }D}F . Similarly, if }∆D}F ď γ 1 ¨ }D}F , then the left-hand inequality gives }∆A}F ď κjoin pDq ¨ γ 1 ¨ }A}F . This proves both claims. We now show that Lemmas 26 and 27 imply that if we can bound the perturbation on the join matrix, we can prove the existence of a bounded perturbation in the input database corresponding to the perturbation in the join matrix, even if the database is not fully reduced. Proposition 15. Consider our setting in Def. 3 and let ∆A be a perturbation of the join matrix A on the join space. Then, there exists a perturbation ∆D “ p∆S1 , . . . , ∆Sr q of the input D such that A ` ∆A is the join matrix of the perturbed input D ` ∆D. Moreover, for every γ ě 0, the following statements hold: 23

• For every j P rrs and Y P Yj , if ∥∆Ar˚, Y s∥2 ď γ ¨ ∥Ar˚, Y s∥2 then ∥∆Sj r˚, Y s∥2 ď κjoin pSj q ¨ γ ¨ ∥Sj r˚, Y s∥2 . • If ∥∆A∥F ď γ ¨ ∥A∥F , then ∥∆D∥F ď κjoin pDq ¨ γ ¨ ∥D∥F . Proof. For every j P rrs, let Ij` “ tp : ϕj ppq ą 0u, and let S` j denote the restriction of Sj to the tuples indexed by Ij` . Thus, S` consists exactly of the non-dangling tuples of Sj . Let j ` D` “ pS` 1 , . . . , Sr q.

The database D` is fully reduced and has the same join matrix A as D. For each p P Ij` and Y P Yj , define ∆Sj rp, Y s to be the common value of the entries of ∆Ar˚, Y s corresponding to occurrences of Sj rp, Y s. This is well defined because ∆A is on the join space. For each dangling tuple, set ∆Sj rp, Y s “ 0

whenever

ϕj ppq “ 0.

Because the perturbations preserve the join-key columns and the join structure, the join matrix obtained from D ` ∆D is exactly A ` ∆A. ` ` ` ` Let ∆S` j be the restriction of ∆Sj to Ij , and let ∆D “ p∆S1 , . . . , ∆Sr q. Since removing ` dangling tuples does not change any positive tuple multiplicity, κjoin pSj q “ κjoin pSj q for every j P rrs, and κjoin pD` q “ κjoin pDq. Suppose first that ∥∆Ar˚, Y s∥2 ď γ¨∥Ar˚, Y s∥2 for every j P rrs and Y P Yj . Since D` is fully ` ď κjoin pS` reduced, Lemma 26 gives ∆S` j r˚, Y s j q ¨ γ ¨ Sj r˚, Y s . The perturbation vanishes 2

2

` on the dangling tuples, so ∥∆Sj r˚, Y s∥2 “ ∆S` j r˚, Y s 2 . Moreover, Sj r˚, Y s 2 ď ∥Sj r˚, Y s∥2 . Consequently, ∥∆Sj r˚, Y s∥2 ď κjoin pSj q ¨ γ ¨ ∥Sj r˚, Y s∥2 . Now suppose that ∥∆A∥F ď γ ∥A∥F . Since D` is fully reduced, Lemma 27 gives ∥∆D` ∥F ď κjoin pD` q ¨ γ ¨ ∥D` ∥F . Because the perturbation vanishes on all dangling tuples, ∥∆D∥F “ ∥∆D` ∥F , while ∥D` ∥F ď ∥D∥F . It follows that ∥∆D∥F ď κjoin pDq ¨ γ ¨ ∥D∥F . This proves both claims.

B.4

Proofs of Statements from Section 6

Here we provide the proofs of the formal statements in Section 6. Theorem 16. Consider the setting in Def. 3 and let C be the computation of the join query followed by the computation Cg of a numerical operation over the join matrix A. For α, β ě 0, the following hold: ` ˘ • If C is entry-wise pα, βq-projected backward stable, then Cg is entry-wise α ` β ¨ p1 ` αq backward stable on the subspace Y0 ; • If is column-wise pα, βq-projected backward stable, then Cg is column-wise bC ` ˘2 ` ˘2 κjoin pSi q ¨ α ` β 2 ¨ 1 ` κjoin pSi q ¨ α -backward stable on the subspace Y0 for each input data column in Si and i P rrs; • If C b is Frobenius-norm-wise pα, βq-projected backward stable, then Cg is Frobenius-norm` ˘2 ` ˘2 wise κjoin pDq ¨ α ` β 2 ¨ 1 ` κjoin pDq ¨ α -backward stable on the subspace Y0 . Proof. We first show that it suffices to consider calibrated input databases. Fix an arbitrary A P Y0 , and let D be any database with the fixed join keys such that f pDq “ A. Let D̄ be obtained from D by setting all data entries of every dangling tuple to zero, without changing any tuple identity or join key. Since dangling tuples do not contribute to the join, f pD̄q “ f pDq “ A, and hence CpD̄q “ Cg pAq. 24

Apply the assumed pα, βq-projected backward stability of C to D̄, and let ∆D̄ and ∆AK be the corresponding perturbations, with à “ f pD̄ ` ∆D̄q. Because perturbations do not change the join keys, every tuple that is dangling in D̄ remains dangling in D̄ ` ∆D̄. Therefore, setting the entries of ∆D̄ corresponding to dangling tuples to zero does not change à and cannot increase any of the entry-wise, column-wise, or Frobenius-norm-wise perturbation bounds. Now let D` be obtained from D̄ by removing all dangling tuples, and let ∆D` be the corresponding restriction of ∆D̄. Then D` is calibrated and f pD` q “ A, f pD` ` ∆D` q “ Ã. Moreover, since all removed data entries are zero, the relevant input norms are unchanged. For every relation Sj and column k, }S` j r˚, ks}2 “ }S̄j r˚, ks}2 ,

}D` }F “ }D̄}F ,

while restricting the perturbation can only decrease these norms: }∆S` j r˚, ks}2 ď }∆S̄j r˚, ks}2 ,

}∆D` }F ď }∆D̄}F .

The entry-wise perturbation bound is preserved directly by restriction. Finally, removing dangling tuples does not change any positive tuple-occurrence number in the join, and therefore κjoin pS` j q “ κjoin pS̄j q,

κjoin pD` q “ κjoin pD̄q.

Thus the projected-backward-stability result on an uncalibrated database can always be represented using a calibrated database, given by the subset of tuples that contribute to the join matrix, without changing A, Ã, ∆AK , the parameters α, β, or the relevant join condition numbers. Hence, without loss of generality, we assume throughout the remainder of the proof that the input database D is calibrated. We write C “ pCf , Cg q, where Cf denotes forming the join, and Cg denotes the computation of g on the join matrix. Case 1: C “ pCf , Cg q is entry-wise pα, βq-projected backward stable. Let ẑ “ Cg pAq “ CpDq be the computed result of g ˝ f . Since C “ pCf , Cg q is entry-wise pα, βq-projected backward stable, there exists a perturbation on the input database ∆D satisfying |∆D| ď α|D|, and the join matrix obtained from D`∆D can be denoted as Ã. Further, there exists ∆AK P Y0K satisfying |∆AK | ď β|Ã|, such that gpà ` ∆AK q “ ẑ. Denote A to be the join matrix of the input database D. Since |∆D| ď α|D|, we have |A ´ Ã| ď α|A| since the entries in A and à are copies of the entries in D and D ` ∆D, respectively. We then have: |A ´ pà ` ∆AK q| ď |A ´ Ã| ` |∆AK | ď α|A| ` β|Ã| ď α|A| ` βp|à ´ A| ` |A|q ´ ¯ ď α|A| ` βpα|A| ` |A|q “ α ` βpα ` 1q |A|. K In other words, ´ for any input ¯ A P Y0 , write E “ à ` ∆A ´ A, if Cg pAq “ ẑ then gpA ` Eq “ ẑ

with |E| ď α ` βpα ` 1q |A|. Hence Cg is entry-wise α ` βpα ` 1q-backward stable when the input lies in Y0 . Case 2: C “ pCf , Cg q is column-wise pα, βq-projected backward stable. Let ẑ “ Cg pAq “ CpDq be the computed result of g ˝ f . Since C “ pCf , Cg q is column-wise pα, βq-projected backward stable, there exists a perturbation on the input database ∆D “ t∆S1 , ..., ∆Sr u satisfying }∆Sj r˚, ks}2 ď α}Sj r˚, ks}2 for every relation Sj and column k, and the join matrix obtained from D ` ∆D can be denoted as Ã. Further, there exists ∆AK P Y0K satisfying }∆AK r˚, is}2 ď β}Ãr˚, is}2 for every column i, such that gpà ` ∆AK q “ ẑ. Denote by A the join matrix for the input database D. }∆Sj r˚, ks}2 ď α}Sj r˚, ks}2 for every relation Sj and column k. Using Lemma 26, we have }pA ´ Ãqr˚, is}2 ď κjoin pSj qα}Ar˚, is}2 for column i in A that comes from the relation Sj , as the entries in A and à are copies of the 25

entries in D and D ` ∆D, respectively. Notice that the query Q on D does not contain self joins, projecting à ` ∆AK on Y0 is equivalent to projecting on Y0 with each column projected independently. The Pythagorean theorem therefore tells us that }pA´ Ãqr˚, is}22 `}∆AK r˚, is}22 “ ›` › ˘ › A ´ pà ` ∆AK q r˚, is›2 for every column i. 2 We then have: ›` › ˘ › A ´ pà ` ∆AK q r˚, is›2 “ }pA ´ Ãqr˚, is}22 ` }∆AK r˚, is}22 2 ď α2 κjoin pSj q2 }Ar˚, is}22 ` }∆AK r˚, is}22 ď α2 κjoin pSj q2 }Ar˚, is}22 ` β 2 }Ãr˚, is}22 . Further, using the triangle inequality, ` ˘ }Ãr˚, is}2 ď }Ar˚, is}2 ` }pA ´ Ãqr˚, is}2 ď 1 ` ακjoin pSj q }Ar˚, is}2 , we obtain ›` › ˘ › A ´ pà ` ∆AK q r˚, is›2 ď α2 κjoin pSj q2 }Ar˚, is}22 ` β 2 }Ãr˚, is}22 2 ` ˘2 ď α2 κjoin pSj q2 }Ar˚, is}22 ` β 2 1 ` ακjoin pSj q }Ar˚, is}22 , which gives us the bound for the backward error of g: b › ›` › ˘ ` ˘ › › A ´ pà ` ∆AK q r˚, is› ď κjoin pSj q2 α2 ` β 2 1 ` κjoin pSj qα 2 ››Ar˚, is›› 2

2

K In other words, for any input bA P Y0 , write E` “ à ` ∆A ˘´ ›A, if C›g pAq “ ẑ then 2 › › gpA ` Eq “ ẑ with }Er˚, is}2 ď κjoin pSj q2 α2 ` β 2 1 ` κjoin pSj qα ›Ar˚, is› . Hence Cg is 2 b ` ˘2 2 2 2 column-wise κjoin pSj q α ` β 1 ` κjoin pSj qα -backward stable when the input lies in Y0 . Case 3: C “ pCf , Cg q is Frobenius-norm-wise pα, βq-projected backward stable. Let ẑ “ Cg pAq “ CpDq be the computed result of g ˝ f . Since C “ pCf , Cg q is Frobeniusnorm-wise pα, βq-projected backward stable, there exists a perturbation on the input database ∆D “ t∆S1 , ..., ∆Sr u satisfying }∆D}F ď α}D}F , and the join matrix obtained from D ` ∆D can be denoted as Ã. Further, there exists ∆AK P Y0K satisfying }∆AK }F ď β}Ã}F such that gpà ` ∆AK q “ ẑ. Denote A to be the join matrix of the input database D. Since }∆D}F ď α}D}F , using Lemma 27, we have }pA ´ Ãq}F ď κjoin pDqα}A}F . Moreover, the Pythagorean theorem tells us ›` ˘›2 that }pA ´ Ãq}2F ` }∆AK }2F “ › A ´ pà ` ∆AK q ›F , we then have: ›` ˘› › A ´ pà ` ∆AK q ›2 “ }pA ´ Ãq}2F ` }∆AK }2F ď α2 κjoin pDq2 }A}2F ` }∆AK }2F F

ď α2 κjoin pDq2 }A}2F ` β 2 }Ã}2F . Further, using the triangle inequality, ` ˘ }Ã}F ď }A}F ` }pA ´ Ãq}F ď 1 ` ακjoin pDq }A}F , we obtain ›` ˘› › A ´ pà ` ∆AK q ›2 ď α2 κjoin pDq2 }A}2F ` β 2 }Ã}2F F ` ˘2 ď α2 κjoin pDq2 }A}2F ` β 2 1 ` ακjoin pDq }A}2F which gives us the bound for the backward error of g: b ›` ˘› ` ˘ › › › A ´ pà ` ∆AK q › ď κjoin pDq2 α2 ` β 2 1 ` κjoin pDqα 2 ››A›› F

F

K In other words, bfor any input A P`Y0 , write E “˘Ã›` ›∆A ´ A, if Cg pAq “ ẑ then gpA ` Eq “ ẑ 2› › with }E}F ď κjoin pDq2 α2 ` β 2 1 ` κjoin pDqα ›A› . Hence Cg is Frobenius-norm-wise F b ` ˘2 2 2 2 κjoin pDq α ` β 1 ` κjoin pDqα -backward stable when the input lies in Y0 .

26

Theorem 17. Consider the setting in Def. 3 and let C be the computation of the join query followed by the computation Cg of a numerical operation over the join matrix A. For 0 ď ϵ ă 1, the following hold: 2ϵ • If Cg is entry-wise ϵ-backward stable, then C is entry-wise pϵ, 1´ϵ q-projected backward stable. ϵ • If Cg is column-wise ϵ-backward stable, then C is column-wise pκjoin pSi q ¨ ϵ, ?1´ϵ 2 q-projected backward stable for each input data column in Si and i P rrs;

• If Cg is Frobenius-norm-wise ϵ-backward stable, then C is Frobenius-norm-wise ϵ pκjoin pDq ¨ ϵ, ?1´ϵ 2 q-projected backward stable. Proof. We write C “ pCf , Cg q, where Cf denotes forming the join, and Cg denotes the computation of g on the join matrix. Case 1: Cg is entry-wise ϵ-backward stable. Let ẑ “ Cg pAq “ CpDq be the computed result of g ˝ f . Since Cg is entry-wise ϵ-backward stable, there exists a perturbation matrix E satisfying ẑ “ gpA ` Eq, such that |E| ď ϵ|A|. We can write A ` E as à ` ∆AK , where à is the projection of A ` E on the join space Y0 , and ∆AK “ A ` E ´ à denotes the normal perturbation. Now, consider any entry in the input database s “ Sk rp, qs that occurs in the join matrix A. In this case, the set Ps “ tpi, jq | Ari, js is parameterized by su is nonempty. Lemma 24 tells us that, in the projection of A ` E on Y0 , namely Ã, for any pi1 , j 1 q P Ps , ř ř ř K pi,jqPPs pà ` ∆A qri, js pi,jqPPs pA ` Eqri, js pi,jqPPs Eri, js 1 1 Ãri , j s “ “ “s` , |Ps | |Ps | |Ps | where ř ˇ ř ˇř pi,jqPPs |Eri, js| pi,jqPPs ϵ|Ari, js| ˇ pi,jqPPs Eri, js ˇ ď “ ϵ|s|. ˇď ˇ |Ps | |Ps | |Ps | This means |Ãri, js´Ari, js| ď ϵ|s| for every pi, jq P Ps . We construct the perturbed database D`∆D “ pS1 `∆S1 , . . . , Sr `∆Sr q as follows. For every database entry s “ Sk rp, qs that occurs in the join matrix, set its perturbation so that its perturbed value equals the corresponding value in Ã. By the argument above, this perturbation satisfies |∆Sk rp, qs| ď ϵ|Sk rp, qs|. If s “ Sk rp, qs belongs to a dangling tuple and therefore does not occur in the join matrix, set ∆Sk rp, qs “ 0. Such an entry does not affect the join result, and its perturbation trivially satisfies |∆Sk rp, qs| “ 0 ď ϵ|Sk rp, qs|. Therefore, f pD ` ∆Dq “ à and |∆Sk rp, qs| ď ϵ|Sk rp, qs| for every relation k and every entry pp, qq. Hence, ϵ is a valid entry-wise error bound on the input database. 2ϵ Next, we show that 1´ϵ is a valid bound for the normal perturbation ∆AK . Notice that ř

Ãri, js “ s `

pi,jqPPs Eri,js |Ps |

ř

and therefore s “ Ãri, js ´

pi,jqPPs Eri,js |Ps |

. We further have:

ˇř ˇ ˇ pi,jqPPs Eri, js ˇ |s| ď |Ãri, js| ` ˇ ˇ ď |Ãri, js| ` ϵ|s|, |Ps | which means |Ãri, js| ě p1 ´ ϵq|s|. On the other hand, |∆AK ri, js| “ |pA ` E ´ Ãqri, js| where |pA ` E ´ Ãqri, js| ď |Eri, js| ` |Ari, js ´ Ãri, js| ď ϵ|s| ` ϵ|s| “ 2ϵ|s|. 27

2ϵ If |s| “ 0, we immediately have |∆AK ri, js| “ 0 and |∆AK ri, js| ď 1´ϵ |Ãri, js| holds true. If 2ϵ|s| 2ϵ |s| ą 0, we also immediately have |∆AK ri, js| ď p1´ϵq|s| |Ãri, js| “ 1´ϵ |Ãri, js|. Applying this 2ϵ argument entry-wise gives |∆AK | ď 1´ϵ |Ã|. 2ϵ Therefore, the computation C “ pCf , Cg q is pϵ, 1´ϵ q-projected backward stable on g ˝ f . Case 2: Cg is column-wise ϵ-backward stable. Let ẑ “ Cg pAq “ CpDq be the computed result of g ˝ f . Since Cg is column-wise ϵ-backward stable, and the query has no self-joins, there exists a perturbation matrix E satisfying ẑ “ gpA ` Eq, such that }Er˚, is}2 ď ϵ}Ar˚, is}2 for every column i. Since the query has no self-joins, any data value in the input database will not parameterize entries from different columns in A. In other words, projecting A ` E on Y0 is equivalent to projecting each column in A ` E onto the vector space defined by each column in Y0 individually. We can write A ` E as à ` ∆AK , where à is the projection of A ` E on the join space Y0 , and ∆AK “ A ` E ´ à denotes the normal perturbation. For each column i, we apply the Pythagorean theorem, which states: }pA ´ Ãqr˚, is}22 ` }∆AK r˚, is}22 “ }Er˚, is}22 . Therefore }pA ´ Ãqr˚, is}2 ď }Er˚, is}2 ď ϵ}Ar˚, is}2 for every column i. Applying Proposition 15, we obtain that there exists a perturbed database D ` ∆D “ tS1 ` ∆S1 , ..., Sr ` ∆Sr u satisfying à is the join matrix of D ` ∆D, and for every relation Sj and column k, }∆Sj r˚, ks}2 ď κjoin pSj q ¨ ϵ ¨ }Sj r˚, ks}2 . We now analyze the normal perturbation bound. Fix a column i. We first consider the case where }Ar˚, is}2 “ 0. Since Cg is column-wise ϵ-backward stable, we have

}Er˚, is}2 ď ϵ}Ar˚, is}2 “ 0. Thus, Er˚, is “ 0. By the Pythagorean theorem, }pA ´ Ãqr˚, is}22 ` }∆AK r˚, is}22 “ }Er˚, is}22 , we obtain pA ´ Ãqr˚, is “ 0 and ∆AK r˚, is “ 0. Consequently, Ãr˚, is “ Ar˚, is “ 0, and hence ϵ }∆AK r˚, is}2 ď ?1´ϵ 2 }Ãr˚, is}2 holds trivially. Now suppose that }Ar˚, is}2 ą 0, and define λ“

}pA ´ Ãqr˚, is}2 . }Ar˚, is}2

Since }pA ´ Ãqr˚, is}2 ď }Er˚, is}2 ď ϵ}Ar˚, is}2 , we have 0 ď λ ď ϵ ă 1. By the triangle inequality, }Ãr˚, is}2 “ }Ar˚, is ´ pA ´ Ãqr˚, is}2 ě }Ar˚, is}2 ´ }pA ´ Ãqr˚, is}2 “ p1 ´ λq}Ar˚, is}2 ą 0. Therefore, division by }Ãr˚, is}2 and by p1 ´ λq}Ar˚, is}2 is well defined. Using the Pythagorean theorem again, we obtain }∆AK r˚, is}22 “ }Er˚, is}22 ´ }pA ´ Ãqr˚, is}22 . Since }Er˚, is}2 ď ϵ}Ar˚, is}2 , it follows that }∆AK r˚, is}22 ď ϵ2 }Ar˚, is}22 ´ λ2 }Ar˚, is}22 “ pϵ2 ´ λ2 q}Ar˚, is}22 . Combining this bound with the lower bound on }Ãr˚, is}2 , we obtain ? ? }∆AK r˚, is}2 ϵ2 ´ λ2 }Ar˚, is}2 ϵ2 ´ λ2 ď “ . p1 ´ λq}Ar˚, is}2 1´λ }Ãr˚, is}2 28

Moreover, ϵ2 ϵ2 ´ λ2 pλ ´ ϵ2 q2 ´ “ ě 0. 1 ´ ϵ2 p1 ´ λq2 p1 ´ ϵ2 qp1 ´ λq2 Since both sides are nonnegative, this implies ? ϵ2 ´ λ2 ϵ . ď? 1´λ 1 ´ ϵ2 Therefore, }∆AK r˚, is}2 ď ?

ϵ }Ãr˚, is}2 . 1 ´ ϵ2

The inequality ´ holds for every column ¯ i, including the case }Ar˚, is}2 “ 0. Thus, C “ pCf , Cg q

ϵ is column-wise κjoin pSj q ¨ ϵ, ?1´ϵ 2 -projected backward stable on g ˝ f for columns originating from relation Sj . Case 3: Cg is Frobenius-norm-wise ϵ-backward stable. Let ẑ “ Cg pAq “ CpDq be the computed result of g ˝ f . Since Cg is Frobenius-norm-wise ϵ-backward stable, there exists a perturbation matrix E satisfying ẑ “ gpA ` Eq, such that }E}F ď ϵ}A}F . We can write A ` E as à ` ∆AK , where à is the projection of A ` E on the join space Y0 , and ∆AK “ A ` E ´ à denotes the normal perturbation. We apply the Pythagorean theorem, which states: }pA ´ Ãq}2F ` }∆AK }2F “ }E}2F . Therefore, }pA ´ Ãq}F ď }E}F ď ϵ}A}F . Applying Proposition 15, we obtain that there exists a perturbed database D ` ∆D “ tS1 ` ∆S1 , ..., Sr ` ∆Sr u satisfying à is the join matrix of D ` ∆D, and }∆D}F ď κjoin pDqϵ}D}F . We now analyze the normal perturbation bound. We first consider the case where }A}F “ 0. Since Cg is Frobenius-norm-wise ϵ-backward stable, we have }E}F ď ϵ}A}F “ 0. Thus, E “ 0. By the Pythagorean theorem,

}A ´ Ã}2F ` }∆AK }2F “ }E}2F , we therefore obtain A ´ Ã “ 0 and ∆AK “ 0. Consequently, Ã “ A “ 0, and hence }∆AK }F ď ?

ϵ }Ã}F 1 ´ ϵ2

holds trivially. Now suppose that }A}F ą 0, and define λ“

}A ´ Ã}F . }A}F

Since }A ´ Ã}F ď }E}F ď ϵ}A}F , we have 0 ď λ ď ϵ ă 1. By the triangle inequality, }Ã}F “ }A ´ pA ´ Ãq}F ě }A}F ´ }A ´ Ã}F “ p1 ´ λq}A}F ą 0. Therefore, division by }Ã}F and by p1 ´ λq}A}F is well defined. Using the Pythagorean theorem again, we obtain }∆AK }2F “ }E}2F ´ }A ´ Ã}2F . Since }E}F ď ϵ}A}F , it follows that }∆AK }2F ď ϵ2 }A}2F ´ λ2 }A}2F “ pϵ2 ´ λ2 q}A}2F . 29

Combining this bound with the lower bound on }Ã}F , we obtain ? ? }∆AK }F ϵ2 ´ λ2 }A}F ϵ2 ´ λ2 “ ď . p1 ´ λq}A}F 1´λ }Ã}F Moreover, ϵ2 ´ λ2 pλ ´ ϵ2 q2 ϵ2 ´ “ ě 0. 1 ´ ϵ2 p1 ´ λq2 p1 ´ ϵ2 qp1 ´ λq2 Since both sides are nonnegative, this implies ? ϵ2 ´ λ2 ϵ . ď? 1´λ 1 ´ ϵ2 Therefore, }∆AK }F ď ?

ϵ }Ã}F . 1 ´ ϵ2

This inequality also holds in the case ¯ }A}F “ 0. Thus, we have proved that C “ pCf , Cg q is ´ ϵ Frobenius-norm-wise κjoin pDqϵ, ?1´ϵ 2 -projected backward stable on g ˝ f .

B.5

Proofs of Statements from Section 7

Here we provide the proofs of the formal statements in Section 7. Proposition 28. There exists an α-acyclic join query Q with no self-joins over fully reduced relations for which the standard floating-point computation of the upper triangular matrix R, obtained from QR decomposition via Givens Rotations, is not backward stable with respect to the input database. Proof. Consider the four input relations Sj “

key Yj k 1 , k ´1

j “ 1, 2, 3, 4.

Their join is the full Cartesian product over the common key. We order the resulting tuples lexicographically, with 1 preceding ´1. Thus the materialized join matrix A P R16ˆ4 contains all 16 rows py1 , y2 , y3 , y4 q,

yj P t1, ´1u.

We compute the QR factorization of A using Givens rotations in IEEE double, whose unit roundoff is u “ 2´53 . The algorithm processes columns from left to right, so j “ 1, 2, 3, 4. Denote the intermediate matrix in Givens Rotations as B. In column j, it eliminates the entries below the diagonal from bottom to top. Thus, for i “ 16, 15, . . . , j ` 1, the algorithm eliminates the current entry Bri, js using the entry immediately above it, Bri ´ 1, js. For one such elimination, set x “ Bri ´ 1, js, y “ Bri, js. If y “ 0, no rotation is needed and the algorithm proceeds to the next entry. Otherwise, it computes ´a ¯ r “ fl f lpx2 ` y 2 q , c “ f lpx{rq, s “ f lpy{rq. The corresponding Givens rotation is then applied to the two rows i ´ 1 and i. Only columns k “ j, j ` 1, . . . , 4 need to be updated, since columns 1, . . . , j ´ 1 have already been reduced and 30

have zero on those rows and the values will stay zero after any rotation. For each such k, the update is ˆ ˙ ˆˆ ˙ˆ ˙˙ Bri ´ 1, ks c s Bri ´ 1, ks Ð fl . Bri, ks ´s c Bri, ks After the update, the entry Bri, js, which we target to annihilate, is explicitly set to zero. In exact arithmetic, the columns of A are mutually orthogonal and each column has squared norm 16. Hence AJ A “ 16I4 , and the exact QR factor with positive diagonal is R “ 4I4 . Let R̂ be the upper-triangular factor computed by the Givens algorithm described above. A direct IEEE double evaluation of this specified algorithm gives the following leading 4 ˆ 4 triangular factor: ¨ ¨ 3 ˛ 5 3 ˛ 4 ´ 52 u 0 ´ 52 45 2 4u 2u ‹ ˚0 16 ´1 ´1‹ ˚0 4 ` 16u ´u ´u ‹“˚ ‹. R̂ “ 4I4 ` u ˚ ˝0 0 8 ´1‚ ˝0 0 4 ` 8u ´u‚ 0 0 0 0 0 0 0 4 We now show that this R̂ cannot be the exact R factor of the join of any perturbed input relations, no matter how large the perturbation is. A perturbation of the numerical attributes of the input relations replaces the two numerical values in S1 , S2 , S3 , S4 by pa11 , a12 q, pb11 , b12 q, pc11 , c12 q, pd11 , d12 q, respectively. Thus the perturbed relations are S11 pK, Y1 q “ tpk, a11 q, pk, a12 qu,

S21 pK, Y2 q “ tpk, b11 q, pk, b12 qu,

S31 pK, Y3 q “ tpk, c11 q, pk, c12 qu,

S41 pK, Y4 q “ tpk, d11 q, pk, d12 qu.

Let A1 be the join matrix obtained from the natural join of these perturbed relations. Suppose, for contradiction, that R̂ is the exact R factor of A1 . Then there exists a matrix 1 Q with orthonormal columns such that A1 “ Q1 R̂. Therefore A1J A1 “ R̂J R̂. Let G “ A1J A1 . Also define σa “ a11 ` a12 ,

σb “ b11 ` b12 ,

σc “ c11 ` c12 ,

σd “ d11 ` d12 .

Since the join is the full Cartesian product, the off-diagonal entries of G satisfy Gr1, 2s “ 4σa σb ,

Gr1, 3s “ 4σa σc ,

Gr1, 4s “ 4σa σd ,

Gr2, 3s “ 4σb σc ,

Gr2, 4s “ 4σb σd ,

Gr3, 4s “ 4σc σd .

Hence every Gram matrix obtained from any such perturbed input relations must satisfy the polynomial invariant Gr1, 2sGr3, 4s “ Gr1, 3sGr2, 4s “ Gr1, 4sGr2, 3s. Now define the Gram matrix implied by the computed triangular factor as Ĝ “ R̂J R̂. Using the exact real values of the displayed entries of R̂, we obtain Ĝr1, 2sĜr3, 4s ´ Ĝr1, 3sĜr2, 4s “ 30u2 p2 ` 5uq ‰ 0. Thus Ĝ violates the invariant that every Gram matrix arising from perturbed input relations must satisfy. This contradicts A1J A1 “ R̂J R̂. Therefore there are no real perturbed input values a11 , a12 , b11 , b12 , c11 , c12 , d11 , d12 whose join matrix has R̂ as an exact QR factor. Thus, in this example, the computed Givens QR factor is not even backward-engineerable from any perturbation of the input relations, let alone from a small backward perturbation. This proves that Givens QR, although backward stable with respect to an arbitrary materialized join matrix, need not be backward stable with respect to the input relations that generate the join matrix. 31

The following lemma recalls the backward stability result of the QR decomposition using Givens rotations by Higham [14, Theorem 19.10]. Lemma 29. Let R̂ be the computed upper triangular matrix obtained by performing the QR decomposition of A P Rmˆn using a sequence of Givens rotations. Then there exists an orthogonal matrix Q1 and a perturbation matrix ∆A, such that A ` ∆A “ Q1 R̂, where ∆A satisfies the column-wise bounds }∆Ar˚, js}2 ď γ̃m`n´2 }Ar˚, js}2 , j P rns, where the constant inside γ̃m`n´2 is a modest integer constant independent of m and n. Therefore, the computation of the upper triangular factor in the QR decomposition is backward stable. It produces the exact upper triangular matrix of a matrix that differs from A by a small relative perturbation. Corollary 18. [Theorem 17 and [14]] Consider the setting in Def. 3 and let C be the computation that materializes the join matrix A and then applies Givens Rotations to compute the uppertriangular matrix R from the QR decomposition of A. If γ̃m`n´2 ă 1, then C is column-wise q-projected backward stable for each input data column in Sj and pκjoin pSj q ¨ γ̃m`n´2 , b γ̃m`n´2 2 1´γ̃m`n´2

j P rrs. Proof. By Lemma 29, there exist a perturbation E P Rmˆn and an orthogonal matrix Q1 such that A`E “ Q1 R̂, with }Er˚, is} ď γ̃m`n´2 }Ar˚, is} for every column i. Theorem 17 immediately tells us that the process of computing gpf q is column-wise pκjoin pSj q ¨ γ̃m`n´2 , b γ̃m`n´2 q-projected 2 1´γ̃m`n´2

backward stable. This completes our proof.

Lemma 30 (Data independence of the FiGaRo rotations). Consider the FiGaRo setting in Def. 20. There exists an orthogonal matrix Q0 , depending only on the fixed join keys and the join tree, such that, for every choice of the data columns, „ ȷ R0 A “ Q0 , 0 where A is the corresponding materialised join matrix and R0 is the exact output of FiGaRo. Proof. FiGaRo is equivalent, in exact arithmetic, to applying a collection of ordinary and generalized Givens transformations to blocks of the materialised join matrix [27]. For an ordinary head-and-tail transformation applied to a block with m rows, the corresponding Givens rotations have coefficients c 1 i´1 sin θi “ ´ ? , cos θi “ , i “ 2, . . . , m. i i These coefficients depend only on the number m of rows in the block and not on any data value. For a generalized head-and-tail transformation with weight vector v “ pv1 , . . . , vm qT , the corresponding rotations have coefficients sin θi “ ´

vi , }v1:i }2

cos θi “

}v1:i´1 }2 , }v1:i }2

i “ 2, . . . , m.

These coefficients depend only on the weight vector v. Once the join keys and the join tree are fixed, every key group processed by FiGaRo is fixed, as is the number of rows in each group. Moreover, all entries of the vectors scales, and hence all weight vectors used by the generalized head-and-tail transformations, are determined by group sizes and by the count aggregates Φ˝i , ΦÒi , and ΦÓi . These quantities depend only on the fixed join keys and not on the data columns. 32

It follows that every local Givens rotation represented by FiGaRo, as well as the order in which these rotations act on the rows of the materialised join, is fixed once the join keys, join tree, and row orderings are fixed. Let G0 be the product of all these rotations, together with any fixed row permutations required by the chosen ordering. Then G0 is orthogonal, depends only on the join keys and the join tree, and satisfies „ ȷ R0 G0 A “ 0 for every choice of the data columns. Setting Q0 :“ GT 0 gives „ ȷ R0 A “ Q0 . 0 Therefore, Q0 depends only on the fixed join keys and the join tree, and not on the data values. We recall the following result, stated as Lemma 3.6 by Higham [14], adapted here to our notation: Lemma 31. If Gj ` ∆Gj P Rnˆn satisfies }∆Gj } ď δj }Gj } for all j, where } ¨ } is a consistent norm, then › › ˜ ¸ m m m m ›ź › ź ź ź › › }Gj }. Gj › ď p1 ` δj q ´ 1 › pGj ` ∆Gj q ´ ›j“0 › j“0 j“0 j“0 Before presenting the proof of Theorem 21, first of all, we show that the FiGaRo algorithm guarantees column-wise forward stability. Lemma 32. Consider the FiGaRo setting (Def. 20). Let R0 and R̂0 be the exact and computed outputs of FiGaRo applied to the input database D, respectively. Then for every column i, we have }pR̂0 ´ R0 qr˚, is}2 ď γ̃d }pR0 qr˚, is}2 Proof. If Ar˚, is “ 0, then all values computed for this column are zero and the result follows immediately. We therefore assume that Ar˚, is ‰ 0. Lemma 30 states that there is an orthogonal matrix Q0 such that „ ȷ R0 A “ Q0 . 0 Consequently, ›„ ȷ› › R0 r˚, is › › › “ }R0 r˚, is}2 . }Ar˚, is}2 “ › › 0 2 We adopt the following implementation choice of FiGaRo: we assume that all tuple-count aggregates and scale radicands are computed exactly, such that every scale is stored through its exact integer radicand. ? ?Every scale ? appearing in FiGaRo is the square root of a positive integer. In particular, λ1 ¨ ¨ ¨ aλk “ λ1 ¨ ¨ ¨ λk . Moreover, if the weights of a generalized head-tail computation are vj “ λj , then each squared prefix norm is the integer v12 ` ¨ ¨ ¨ ` vj2 “ λ1 ` ¨ ¨ ¨ ` λj . 33

Hence products of count-based scales and squared prefix norms can be maintained exactly as integer radicands. A square root is rounded only when the corresponding coefficient is formed and used. This consolidates each product of scales that occurs at one point in the algorithm; scale applications separated by other transformations are not moved or combined, and are counted separately below. We now use the stagewise interpretation of the exact FiGaRo computation. At any point, conceptually restore the repeated rows represented by the current scale factors, append the rows already written to Out, and pad the rows not stored by FiGaRo with zeros. This gives a vector ? with the same number of entries as Ar˚, is. For example, multiplying a stored value x by λ represents the ordinary head-tail transformation »? fi » fi λx x — 0 ffi — .. ffi — ffi – . fl ÞÝÑ — .. ffi – . fl x 0 on λ equal virtual copies of x. Generalized head-tail operations similarly combine blocks whose sizes are represented by their scale factors. Thus each exact local operation is an orthogonal transformation of the padded virtual vector. Order these transformations as G1 , . . . , Gs , absorbing exact permutations between stages. Each Gt is orthogonal and, after a suitable exact permutation of coordinates, is block diagonal with ordinary or generalized head-tail blocks and identity blocks. Define z0 “ Ar˚, is,

zt “ Gt zt´1 ,

t “ 1, . . . , s.

Every Gt is orthogonal, and therefore }zt }2 “ }zt´1 }2 “ }Ar˚, is}2 . After the final stage, „ ȷ R0 r˚, is zs “ . 0 Apply the same virtual expansion to the floating-point values stored in Data and Out. Denote the resulting state after stage t by ẑt . Then „ ȷ R̂0 r˚, is ẑ0 “ z0 , ẑs “ . 0 We next analyze one stage. Consider one active generalized head-tail block with input vector bř j 2 x and positive weight vector v. Let ρj “ k“1 vrks . The exact head and tails computed by this block are p

1 ÿ Hpx, vq “ vrksxrks, ρp k“1 and, for j “ 1, . . . , p ´ 1, ρj xrj ` 1s ´ Tj px, vq “

vrj ` 1s řj k“1 vrksxrks ρj . ρj`1

An ordinary head-tail block is the special case vr1s “ ¨ ¨ ¨ “ vrps “ 1. Tj px, vq is also the formula evaluated by the sequential running-sum implementation. Importantly, a fixed input xk contributes through only one algebraic term to each nonzero output coefficient. Thus, after expanding the rounded running sums, multiplications, divisions, square 34

roots, and subtractions, each computed coefficient equals the corresponding exact coefficient multiplied by a product of elementary rounding factors p1 ` δq˘1 . Let qt be an upper bound on the number of such elementary rounding factors occurring along any coefficient path in stage t. The standard θ-notation therefore gives ÿ ` ptq ˘ ptq ẑt rjs “ Gt rj, ks 1 ` θjk ẑt´1 rks, |θjk | ď γqt . k

Equivalently, there exists a perturbation matrix ∆Gt such that ẑt “ pGt ` ∆Gt qẑt´1 ,

|∆Gt | ď γqt |Gt |,

where the absolute values and inequality are entry-wise. We now bound the norm of the entrywise absolute value of an exact head-tail block. Let „ Tȷ h G“ T be one such block. If the size of the head-tail transformation is 1, then G “ r1s and }|G|}2 “ 1. Suppose therefore that the size is at least 2. Since G is orthogonal, its head has unit norm and its tail rows are orthonormal: }h}2 “ 1,

}T}2 “ 1.

Every generalized tail row contains one nonnegative coefficient corresponding to the newly introduced input and nonpositive coefficients corresponding to preceding inputs. Hence the tail matrix can be written as T “ P ´ L, where P and L are entrywise nonnegative. The matrix P contains the single positive coefficient in each tail row. It has at most one nonzero entry in each row and column, and every such entry has magnitude at most one. Therefore, }P}2 ď 1. Since L “ P ´ T, }L}2 ď }P}2 ` }T}2 ď 2. The nonzero entry positions of P and L are disjoint, so |T| “ P ` L. It follows that › › ›|T|› ď }P}2 ` }L}2 ď 3. 2 Combining the head and tail bounds gives b a ? › › › › ›|G|› ď }h}2 ` ›|T|›2 ď 1 ` 32 “ 10. 2 2 2 The same estimate holds for Gt , because Gt is block diagonal and the spectral norm of a block-diagonal matrix is the maximum spectral norm of its blocks. Identity and permutation › › ? › › blocks have norm one. Thus, |Gt | 2 ď 10. We now have the following stagewise bound: ? › › › › }∆Gt }2 ď ›|∆Gt |›2 ď γqt ›|Gt |›2 ď 10 γqt ď 4γqt ď γ4qt , provided 4qt u ă 1. Since Gt is orthogonal, }Gt }2 “ 1. Hence }∆Gt }2 ď γ4qt }Gt }2 . Since ẑ0 “ z0 “ Ar˚, is, repeated application of ẑt “ pGt ` ∆Gt qẑt´1 gives ẑs “ pGs ` ∆Gs q ¨ ¨ ¨ pG1 ` ∆G1 qz0 , whereas the exact computation satisfies zs “ Gs ¨ ¨ ¨ G1 z0 . 35

Applying Lemma 31 in the spectral norm therefore gives }ẑs ´ zs }2 ď }pGs ` ∆Gs q ¨ ¨ ¨ pG1 ` ∆G1 q ´ Gs ¨ ¨ ¨ G1 }2 }z0 }2 ˜ ¸˜ ¸ s s ź ź ď p1 ` γ4qt q ´ 1 }Gt }2 }z0 }2 ˜

t“1 s ź

t“1

¸

p1 ` γ4qt q ´ 1 }z0 }2 .

“ t“1

Let νi :“

řs śs 1 t“1 qt . Using 1 ` γx “ 1´xu and t“1 p1 ´ 4qt uq ě 1 ´ 4 t“1 qt u, we obtain

řs

s ź p1 ` γ4qt q ď 1 ` γ4νi . t“1

Consequently, }ẑs ´ zs }2 ď γ4νi }z0 }2 . A head-tail computation on a group with p stored rows has arithmetic depth at most a fixed constant times p: the running sums are updated once per row, and every head or tail uses only a fixed number of additional floating-point operations. For each relation Sj , the ordinary head-tail groups partition the rows of Sj , so the sum of their sizes is dj , where dj denotes the number of rows in Sj . Immediately before projecting Sj toward its parent, the current Data matrix has at most one row for each distinct join-key value occurring in Sj , and therefore at most dj rows. The generalized head-tail groups partition these rows, so the sum of their sizes is also at most dj . Finally, for a fixed column, every scale product required for one stored entry is first consolidated into one exact integer radicand. Its floating-point application therefore has constant arithmetic depth, independent of the number of factors in that product. Under the no-self-join assumption, a fixed data column is carried through a unique path toward the root, and the number of stored entries of that column handled at relation Sj is at most a fixed constant times dj . The tail-scaling operations are bounded in the same way. ř Consequently, there is an absolute implementation constant c such that νi ď c rj“1 dj “ cd. In particular, this bound has no additional dependence on the number of relations or on the size of the materialized join. Therefore, the total number of rounded operations along all stages affecting column i satisfies νi ď cd, where d is the total number of rows in the input database and c is independent of the materialized join size. Therefore, }ẑs ´ zs }2 ď γ4cd }Ar˚, is}2 “ γ̃d }Ar˚, is}2 , provided 4cdu ă 1. Finally, „ ȷ R0 r˚, is zs “ , 0

„ ẑs “

ȷ R̂0 r˚, is . 0

Hence }R̂0 r˚, is ´ R0 r˚, is}2 “ }ẑs ´ zs }2 ď γ̃d }Ar˚, is}2 . Since A “ Q0 rR0 ; 0s with Q0 orthogonal, }Ar˚, is}2 “ }R0 r˚, is}2 . Therefore, }R̂0 r˚, is ´ R0 r˚, is}2 ď γ̃d }R0 r˚, is}2 , which proves the result. 36

Theorem 21. Consider the FiGaRo setting in Def. 20 and let C be the computation that uses FiGaRo on the input database followed by Givens Rotations to compute ` ˘the upper-triangular matrix R from the QR decomposition of A. Let η “ γ̃ ` γ̃ ¨ 1 ` γ̃ . If η ă 1, then C is d d`n´2 d ´ ¯ η column-wise κjoin pSj q ¨ η, ? -projected backward stable for each input data column in Sj 1´η 2

and j P rrs. Proof. Let f pDq “ A be the join, and gpAq “ R be the function mapping the join matrix to the upper-triangular R, and we use R̂ to denote the computed result from FiGaRo followed by Givens Rotations. Notice that this computation computes g ˝ f without ever materializing the join matrix A. Let R0 and R̂0 denote the exact and computed outputs of FiGaRo on D, respectively. Whenever a FiGaRo output is multiplied by an m ˆ m matrix, we implicitly pad it with zero rows to obtain an m ˆ n matrix. By Lemma 30, there exists an orthogonal matrix Q0 , determined only by the fixed key columns and the join tree τ , such that „ ȷ R0 A “ Q0 . 0 Define the FiGaRo forward error by EF :“ R̂0 ´ R0 . By Lemma 32, for every column i, }EF r˚, is}2 ď γ̃d }R0 r˚, is}2 . The matrix R̂0 has at most d nonzero rows. Therefore, by Lemma 29, there exist a perturbation EQ and an orthogonal matrix Q1 such that R̂0 ` EQ “ Q1 R̂, with }EQ r˚, is}2 ď γ̃d`n´2 }R̂0 r˚, is}2 for every column i. Since }R̂0 r˚, is}2 ď }R0 r˚, is}2 ` }EF r˚, is}2 ď p1 ` γ̃d q}R0 r˚, is}2 , we obtain }EQ r˚, is}2 ď γ̃d`n´2 p1 ` γ̃d q}R0 r˚, is}2 . Define the total perturbation in the FiGaRo-output space by E :“ EF ` EQ . Since R̂0 “ R0 ` EF , we have R0 ` E “ R̂0 ` EQ “ Q1 R̂. Moreover, for every column i, }Er˚, is}2 ď }EF r˚, is}2 ` }EQ r˚, is}2 ď η}R0 r˚, is}2 , where η “ γ̃d ` γ̃d`n´2 p1 ` γ̃d q. Define „ ȷ E EA :“ Q0 . 0 Because Q0 is orthogonal, for every column i, }EA r˚, is}2 “ }Er˚, is}2 ď η}R0 r˚, is}2 “ η}Ar˚, is}2 . Furthermore, „ A ` EA “ Q0

ȷ „ 1 ȷ R0 ` E Q R̂ “ Q0 . 0 0

After padding Q1 with an identity matrix on the omitted zero rows, the matrix multiplying R̂ on the right-hand side is orthogonal. Hence R̂ is the exact upper triangular QR factor of A ` EA , and therefore R̂ “ gpA ` EA q. Let PY0 denote the orthogonal projector onto the join space Y0 , and define à :“ PY0 pA ` EA q,

∆AK :“ pI ´ PY0 qpA ` EA q. 37

` ˘ Then à P Y0 , ∆AK P Y0K , and à ` ∆AK “ A ` EA . Consequently, R̂ “ g à ` ∆AK . Since A P Y0 , we have à ´ A “ PY0 EA and ∆AK “ pI ´ PY0 qEA . Since the query has no self-joins, the columns of Y0 are parameterized independently and the projection acts independently on every column. Hence, by the Pythagorean theorem, }pà ´ Aqr˚, is}22 ` }∆AK r˚, is}22 “ }EA r˚, is}22 . It follows that }pà ´ Aqr˚, is}2 ď }EA r˚, is}2 ď η}Ar˚, is}2 . Since à P Y0 , there exists a perturbation ∆D “ p∆S1 , . . . , ∆Sr q such that f pD ` ∆Dq “ Ã. Applying Proposition 15 to à ´ A, we obtain, for every relation Sj and every data column k of Sj , }∆Sj r˚, ks}2 ď κjoin pSj qη}Sj r˚, ks}2 . It remains to bound the normal perturbation. Fix a column i with Ar˚, is ‰ 0, and let λi :“

}pà ´ Aqr˚, is}2 . }Ar˚, is}2

Then 0 ď λi ď η ă 1. By the triangle inequality, }Ãr˚, is}2 ě }Ar˚, is}2 ´ }pà ´ Aqr˚, is}2 “ p1 ´ λi q}Ar˚, is}2 . The Pythagorean theorem and the bound on EA give }∆AK r˚, is}22 ď pη 2 ´ λ2i q}Ar˚, is}22 . Notice that if }Ãr˚, is}2 “ 0, the two inequalities above directly lead to }Ar˚, is}2 “ }∆AK r˚, is}2 “ 0, and the bound on the normal perturbation trivially holds. Now if }∆AK r˚, is}2 ą 0, we have: a η 2 ´ λ2i η }∆AK r˚, is}2 ď ďa , 1 ´ λi }Ãr˚, is}2 1 ´ η2 ? η 2 ´λ2 where the last inequality follows because the function hpλq “ 1´λ reaches its maximum over 0 ď λ ď η at λ “ η 2 . Thus, R̂ “ gpà ` ∆AK q, where à is the join matrix of a columnwise-small perturbation of the input database and ∆AK P Y0K satisfies the claimed ´ normal-perturbation ¯ bound. Therefore, η FiGaRo followed by Givens rotations is column-wise κjoin pSj q ¨ η, ? 2 -projected backward 1´η

stable for each input data column in Sj and j P rrs. Theorem 22. Consider the FiGaRo setting in Def. 20 and let C be the computation that uses FiGaRo on the input database followed by xGESVD to compute the singular values of A. Let ´ ¯ ` ˘ η “ γ̃d ` ϵsvd pd, nq ¨ 1 ` γ̃d . If η ă 1, then C is Frobenius-norm-wise κjoin pDq ¨ η, ? η 2 1´η

projected backward stable. Proof. Let f pDq “ A be the join, and gpAq “ Σ be the function mapping the join matrix to the diagonal matrix containing its singular values, and we use Σ̂ to denote the computed result from FiGaRo followed by xGESVD. Notice that this computation computes g ˝ f without ever materializing the join matrix A. Let R0 and R̂0 denote the exact and computed outputs of FiGaRo on D, respectively. Whenever a FiGaRo output is multiplied by an m ˆ m matrix, we implicitly pad it with zero rows to obtain an m ˆ n matrix. By Lemma 30, there exists an orthogonal matrix Q0 , determined only by the fixed key columns and the join tree τ , such that „ ȷ R0 A “ Q0 . 0 38

Define the FiGaRo forward error by EF :“ R̂0 ´ R0 . Lemma 32 claims the forward error of FiGaRo in every column is upper-bounded by γ̃d , which directly gives us a Frobenius-norm-wise forward error bound γ̃d too. Namely, }EF }F ď γ̃d }R0 }F . The matrix R̂0 has at most d nonzero rows. Therefore, computing Σ̂ with xGESVD from R̂0 means there exist a perturbation EQ and orthogonal matrices U1 and V1 such that R̂0 ` EQ “ U1 Σ̂V1 J , with }EQ }F ď ϵsvd pd, nq ¨ }R̂0 }F . Since }R̂0 }F ď }R0 }F ` }EF }F ď p1 ` γ̃d q}R0 }F , we obtain }EQ }F ď ϵsvd pd, nq ¨ p1 ` γ̃d q}R0 }F . Define the total perturbation in the FiGaRo-output space by E :“ EF ` EQ . Since R̂0 “ R0 ` EF , we have R0 ` E “ R̂0 ` EQ “ U1 Σ̂V1J . Moreover, }E}F ď }EF }F ` }EQ }F ď η}R0 }F , where η “ γ̃d ` ϵsvd pd, nq ¨ p1 ` γ̃d q. Define „ ȷ E EA :“ Q0 . 0 Because Q0 is orthogonal, }EA }F “ }E}F ď η}R0 }F “ η}A}F . Furthermore, „ A ` EA “ Q0

ȷ „ 1 ȷ R0 ` E U Σ̂V1 J “ Q0 . 0 0

After padding U1 with an identity matrix on the omitted zero rows, the matrix multiplying Σ̂V1J on the right-hand side is orthogonal. Hence Σ̂ is the exact diagonal matrix containing the singular values of A ` EA , and therefore Σ̂ “ gpA ` EA q. Let PY0 denote the orthogonal projector onto the join space Y0 , and define à :“ PY0 pA ` EA q,

∆AK :“ pI ´ PY0 qpA ` EA q.

` ˘ Then à P Y0 , ∆AK P Y0K , and à ` ∆AK “ A ` EA . Consequently, Σ̂ “ g à ` ∆AK . Since A P Y0 , we have à ´ A “ PY0 EA and ∆AK “ pI ´ PY0 qEA . By the Pythagorean theorem, }pà ´ Aq}2F ` }∆AK }2F “ }EA }2F . It follows that }pà ´ Aq}F ď }EA }F ď η}A}F . Since à P Y0 , there exists a perturbation ∆D “ p∆S1 , . . . , ∆Sr q such that f pD ` ∆Dq “ Ã. Applying Proposition 15 to à ´ A, we obtain }∆D}F ď κjoin pDqη}D}F . It remains to bound the normal perturbation. Notice that if the join matrix is a zero matrix, then the computed singular values will be all zero and we trivially recover the input database with all zero entries with no perturbation. Now suppose A ‰ 0, and let λ :“

}pà ´ Aq}F . }A}F

Then 0 ď λ ď η ă 1. By the triangle inequality, }Ã}F ě }A}F ´ }pà ´ Aq}F “ p1 ´ λq}A}F . 39

The Pythagorean theorem and the bound on EA give }∆AK }2F ď pη 2 ´ λ2 q}A}2F . Notice that if }Ã}F “ 0, the two inequalities above directly lead to }A}F “ }∆AK }F “ 0, and the bound on the normal perturbation trivially holds. Now if }∆AK }F ą 0, we have: a }∆AK }F η 2 ´ λ2 η ď , ďa 1´λ }Ã}F 1 ´ η2 ?

η 2 ´λ2

where the last inequality follows because the function hpλq “ 1´λ reaches its maximum over 0 ď λ ď η at λ “ η 2 . Thus, Σ̂ “ gpà ` ∆AK q, where à is the join matrix of a Frobenius-norm-wise-small perturbation of the input database and ∆AK P Y0K satisfies the claimed normal-perturbation bound. Therefore, FiGaRo followed by xGESVD is Frobenius-norm-wise ´ ¯ ` ˘ η κjoin pDq ¨ η, ? 2 -projected backward stable, where η :“ γ̃d ` ϵsvd pd, nq ¨ 1 ` γ̃d and when η ă 1.

1´η

40

Related documents

• Numerical Stability of Linear Algebra Operations over Relational Databases#1083300
• Algorithmic stability via ensembling#997311
• DB-SpMSpV: Dual-View Blocked Sparse Matrix-Sparse Vector Multiplication for Dynamic GPU Workloads#620708
• DB-SpMSpV: Dual-View Blocked Sparse Matrix-Sparse Vector Multiplication for Dynamic GPU Workloads#619029
• DB-SpMSpV: Dual-View Blocked Sparse Matrix-Sparse Vector Multiplication for Dynamic GPU Workloads#459612