1 Introduction↩︎

At first glance, the Jacobi eigenvalue algorithm for computing the eigendecomposition of Hermitian matrices and the Gram-Schmidt procedure for computing the QR decomposition of tall matrices are quite distinct: they take different inputs, compute different factorizations, have different asymptotic complexities, etc. However, there is an important commonality: both methods are iterative and consider only pairs of columns of their inputs at a time. In this work, this connection is taken further, and both algorithms are realized as special cases of a more general matrix factorization algorithm.

Our primary result, Corollary 1, is a unified analysis of all special cases of this general procedure when using a particular randomized pivoting rule; we show that these algorithms are all linearly convergent in expectation, and further that the rate of convergence is unaffected by which particular factorization is being computed.

In finite arithmetic, we provide a framework for showing that instances of this general procedure are numerically stable. This framework is applied to several orthogonalization methods (Theorem 3), as well as the Jacobi eigenvalue algorithm (Theorems 4, 5, and Corollary 4), building on the work of [1] in the latter case.

We first give a brief overview of the history and context of these algorithms. As these algorithms have been extensively studied, we cannot give a comprehensive overview. We instead highlight their origins, practical relevance today, and interactions between them.

Eigendecomposition / Singular value decomposition: The classical Jacobi eigenvalue algorithm (which we henceforth refer to as the two-sided Jacobi eigenvalue algorithm²) is an iterative algorithm for computing an approximate eigendecomposition of a given Hermitian³ matrix that dates back to 1846 [3]. It was rediscovered and popularized in the 1940s and 50s when the advent of computers made large numerical computations feasible (we refer the reader to [4] for a detailed history). Today, it is included in popular textbooks [5], [6] and software packages [7]. Although it is often less preferred due to a slower runtime [5], it has several strengths when compared with other popular methods for diagonalization. First, if the given matrix \(B\) is already close to diagonal, the arithmetic cost is greatly reduced [8]. Second, the computation is highly parallelizable, and each thread only needs a small amount of memory (see, for example, [9]). Third, numerical experiments reveal classes of matrices for which all the eigenvalues appear to be computed with small relative error [1].

The corresponding algorithm for computing the singular value decomposition of a given matrix \(A\in{\mathbb{C}}^{d\times n}\), today called the one-sided Jacobi eigenvalue algorithm, was introduced as an independent algorithm by Hestenes [10], then rediscovered and connected to the classical two-sided Jacobi eigenvalue algorithm by Chartres [11]. It can be viewed as implicitly running the two-sided Jacobi eigenvalue algorithm on the matrix \(B=A^*A\), without ever forming \(B\). It enjoys many of the same advantages as the two-sided Jacobi eigenvalue algorithm [1], [12].

The numerical stability of both the one- and two-sided variants was investigated by Demmel and Veselić in a highly influential paper [1]. They primarily focus on the accuracy with which the small eigenvalues of positive definite matrices are computed. They find “modulo an assumption based on extensive numerical tests” that the Jacobi eigenvalue algorithm essentially produces the best approximation of the small eigenvalues one could hope for in the floating point arithmetic model. A family of matrices with more problematic behavior that those tested by [1] was published by [13], though they do not result in instability of the method overall.

Cholesky decomposition / QR decomposition: Gaussian elimination is the oldest method for factoring a given matrix \(B\) as \(B=LU\) where \(L\) is lower triangular and \(U\) is upper triangular [14]. When \(B\) is positive definite, one is guaranteed that \(U=L^*\) resulting in the Cholesky decomposition \(B=LL^*\). Gaussian elimination with various pivoting rules is the primary algorithm used in practice for the Cholesky decomposition [7], [15], [16].

There are many methods for the QR decomposition. Both classical and modified Gram-Schmidt applied to \(A\in{\mathbb{C}}^{d\times n}\) are mathematically equivalent to Gaussian elimination applied to the matrix \(B=A^*A\) [17]. Some other algorithms actively exploit the connection between the QR and Cholesky decompositions: [18] proposes computing the Cholesky decomposition \(A^*A=R^*R\) and then computing \(Q=AR^{-1}\), and the algorithm(s) of [19] and [20] explore this idea further, with modifications to improve stability.

Several works consider using these algorithms in conjunction to improve accuracy. For instance, [21] observes that one can accelerate the one-sided Jacobi eigenvalue algorithm for computing the SVD by first computing \(A=QR\) and then using one-sided Jacobi algorithm to find the SVD of the upper triangular factor \(R\). Similarly, [22] shows the stability of two-sided Jacobi can be improved by first running Cholesky on a positive definite matrix \(B\), then applying one-sided Jacobi to the Cholesky factor of \(B\).

Other works take inspiration from algorithms for one task for the other. For instance, [23] proposes a Jacobi-like method for the QR decomposition, applying Givens rotations to transform a matrix \(A\) into an upper triangular one \(R\). The accumulation of all the Givens rotations then forms \(Q\) (which is akin to the method for the QR decomposition using Householder reflections [15]). This method for QR now commonly appears in numerical linear algebra textbooks, e.g. [5].

The central focus of this paper is the striking similarity between the operations performed by the once-sided Jacobi eigenvalue algorithm and the modified Gram-Schmidt procedure. This similarity has been previously noted, for instance [24] observes that the operations are actually the same under a certain kind of infinite matrix limit. We show that this similarity can be exploited to obtain a unified analysis of these algorithms and more.

1.1 Paper overview↩︎

In Section 2, we present the Jacobi eigenvalue algorithm and modified Gram-Schmidt procedure (MGS) in a way that highlights their similarity. In Section 3 we formally describe an abstract algorithm which encompasses the Jacobi eigenvalue algorithm and MGS as special cases.

Our results about this generalized algorithm are presented in the final two sections. Section 4 focuses on exact-arithmetic implementations of the algorithm, in the real RAM model. Modulo a mild condition number assumption (which is satisfied by small perturbations of fixed matrices), we show that \(O(n^2\log(n/\delta))\) iterations of this procedure suffice to produce \(\delta\)-approximate factorizations (stated precisely in Corollary 1). For many reasonable special cases of this generalized algorithm, each iteration costs only \(O(n)\) arithmetic operations, resulting in many \(\widetilde{O}(n^3)\)-time factorization algorithms. Additionally, this section confirms a recent conjecture of Steinerberger about the convergence of particular method for orthogonalization [25] (see Remark 6).

Section 5 studies the general algorithm in finite arithmetic. Here, we have two primary results: the first is that a large class of orthogonalization algorithms are numerically stable (including an algorithm for computing the QR decomposition specifically). The second is that with randomized pivoting, the Jacobi eigenvalue algorithm is stable. Specifically, the “assumption based on extensive numerical tests” needed by Demmel and Veselić in [1] is proved to hold with high probability over the randomness of the pivoting (see Corollary 2).

2 Connecting SVD and QR↩︎

2.1 The core connection↩︎

Here we describe a simple connection between the SVD and QR decomposition (and the corresponding connection between the eigendecomposition and Cholesky decomposition) and algorithms for computing them. To see the connection, we first express the QR decomposition and a factorization that is almost the SVD in a common form: given \(A\) with full column rank, write \[{\label{a11}A=QT^{-1}}\tag{1}\] where \(Q\) has orthogonal columns and \(T\) is nonsingular. When \(T\) is unitary, by further decomposing \(Q=U\sqrt{Q^*Q}\), one obtains the singular value decomposition. When \(T\) is upper triangular and \(Q^*Q=I\), this is exactly the QR decomposition. Analogously, the eigendecomposition and Cholesky decomposition of \(B=A^*A\) can both be expressed as \(B=(T^{-1})^*DT\) where \(D=Q^*Q\). The four cases are summarized in the following table.

Notably, both Modified Gram Schmidt (MGS) and one-sided Jacobi can be interpreted as methods for producing an orthogonal basis \(Q=\begin{bmatrix}q_1&\cdots&q_n\end{bmatrix}\) for a given collection of vectors, expressed as columns of \(A=\begin{bmatrix}a_1&\cdots&a_n\end{bmatrix}\) (along with the corresponding column operation sending \(A\) to \(Q\)). The algorithms themselves closely resemble each other too: they both iteratively pick pairs of indices \((i,j)\) and replace the columns \(a_i,a_j\) with a different orthogonal basis for their span. By repeating this many times with appropriate selection of \((i,j)\) each time, the collection of vectors will approach an orthogonal set. More specifically let \(A^{(t)}=\begin{bmatrix}a_1^{(t)}&\cdots&a_n^{(t)}\end{bmatrix}\) be the state after this has been done for \(t\) pairs \((i,j)\). For each \(t\), the Jacobi eigenvalue algorithm selects a pair \((i,j)\) and then updates \[\label{a12}\begin{align} a_i^{(t+1)}&=\cos(\theta)a_i^{(t)}+ \sin(\theta)a_j^{(t)} \\ a_j^{(t+1)}&=-\sin(\theta)a_i^{(t)}+ \cos(\theta)a_j^{(t)} \\ a_\ell^{(t+1)}&=a_\ell^{(t)}\quad\forall\ell\not\in{\left\{{i,j}\right\}}\end{align}\tag{2}\] where \(\theta\) is such that the resulting columns \(a_i^{(t+1)}\) and \(a_j^{(t+1)}\) are orthogonal. MGS selects \(i<j\) then updates \[\label{a13}\begin{align} a_{j}^{(t+1)}&=\frac{a^{(t)}_{j}-\left\langle a^{(t)}_{j},a^{(t)}_{i}\right\rangle a^{(t)}_{i}}{\left\|a^{(t)}_{j}-\left\langle a^{(t)}_{j},a^{(t)}_{i}\right\rangle a^{(t)}_{i}\right\|} \\ a_\ell^{(t+1)}&=a_\ell^{(t)}\quad\forall\ell\neq j.\end{align}\tag{3}\] The update equation 2 is a Givens rotation of \(A^{(t)}\) (also called a plane rotation); importantly it is unitary. The product of these rotations forms \(T\) in the decomposition 1 . The update equation 3 is an elementary matrix operation followed by a scaling; importantly, it is upper triangular. The product of these elementary operations forms \(T\) in the decomposition 1 .

In both cases, the result of this update is that the \(i\)th and \(j\)th columns of \(A^{(t+1)}\) will be orthogonal. Despite their similarity, the existing analysis of these algorithms is quite different, as is their empirical performance measured in terms of computational and communication costs and numerical stability.

2.2 Pivoting rules↩︎

In both cases, the strategy for picking \((i,j)\) at each time \(t\) is determined by a “pivoting rule”. Multiple pivoting rules have been proposed; we group them into three categories to discuss: fixed, adaptive, and randomized. Before doing so, we must clarify some terminology. There is some subtlety in what the literature refers to as “one iteration” of either procedure. For the Jacobi eigenvalue algorithm, one iteration is sometimes taken to mean one “sweep” (which we describe in a moment) of all \(n(n-1)/2\) pairs \((i,j)\). For MGS (more specifically, left-looking MGS), one iteration typically refers to the grouping of updates \[{\label{a14} (1,i),(2,i),\ldots,(i-1,i)}\tag{4}\] for each \(i\). For Gaussian elimination (and right-looking MGS), it typically refers to the grouping of updates \[{\label{a15} (i,i+1),(i,i+2),\ldots,(i,n)}\tag{5}\] for each \(i\). Because of this, the word “pivot” in those contexts typically refer to the single column index \(i\). In this paper, we refer to each selection of \((i,j)\) and update of the corresponding column(s) as one iteration, i.e. \(A^{(t)}\) defined by either 2 , 3 is the state after \(t\) iterations. One pivot then corresponds to a pair of column indices \((i,j)\).

Fixed pivot rules: These are rules that fix a sequence of \((i,j)\) ahead of time. The ones we highlight for the Jacobi eigenvalue algorithm are called cyclic or sweep rules, due to [26], see also [27]. The row-cyclic version uses the sequence \[(1,2),(1,3),(1,4),(1,5),\ldots,(1,n),(2,3),(2,4),(2,5),\ldots,(2,n),\ldots (n-1,n)\] and repeats it over and over again. The column-cyclic version uses the sequence \[(1,2),(1,3),(2,3),(1,4),(2,4),(3,4),(1,5),(2,5),(3,5),(4,5),\ldots,(n-1,n)\] and repeats it over and over again. Each run through the sequence of \(n(n-1)/2\) pairs is referred to as one “sweep.” In the context of Gaussian elimination and MGS, those sweeps correspond to applying the groupings 4 or 5 for \(i=1,\ldots,n-1\). Running MGS for just one or two sweeps are both popular proposals.

For appropriately chosen rotations, [27] show that row-cyclic and column-cyclic Jacobi converges, and [8] improve this to (eventual) quadratic convergence. However, there is no rigorous theory to predict the number of sweeps required to achieve a desired distance to diagonal [18]. In practice, it is observed to be a constant number of sweeps [5]. Better theoretical results are available for the next category of pivoting rules.

Adaptive pivot rules: These are pivoting rules that use the current value of \(A^{(t)}\) to decide the next pivots. The original rule proposed by Jacobi [3], which we call the greedy pivoting rule, was to pick \[{\label{a16}(i,j)=\mathop{\mathrm{arg\,max}}_{i,j:i\neq j}\mleft|\left\langle a_i,a_j\right\rangle\mright|.}\tag{6}\] This is somewhat similar to column-pivoting rule of Gram Schmidt (implicitly equivalent to complete pivoting rule of Gaussian elimination). This rule for Gram Schmidt maintains a list of indices \(\mathcal{I}\) initialized to \([n]\). Then \[{\label{a17}i=\mathop{\mathrm{arg\,max}}_{i\in\mathcal{I}}\left\|a_i\right\|}\tag{7}\] is computed, \(i\) is removed from \(\mathcal{I}\), and the updates \({\left\{{(i,j):j\in\mathcal{I}}\right\}}\) are performed (in any order). This is repeated until \(\mathcal{I}\) is empty. Strictly speaking, these rules will not necessarily compute the QR decomposition of \(A\), since we no longer are assured that the update pairs \((i,j)\) satisfy \(i<j\). However, if one lets \(\sigma\) be the permutation which records the order in which indices are removed from \(\mathcal{I}\), then we are assured \(\sigma(i)<\sigma(j)\) and therefore the method computes QR decomposition of \(AP_\sigma\) where \(P_\sigma\) is the permutation matrix associated with \(\sigma\).

The same pivoting strategy has been considered for one-sided Jacobi [28]. Numerical experiments showed that this pivoting strategy can improve the rate of convergence in several (but not all) cases (see [28] for more details).

Randomized pivot rules: There seems to be much less attention given to randomized pivot rules. A recent work of Chen et al. [29] considers the partial Cholesky and QR decompositions for computing low-rank approximations. They study a randomized variant of the complete pivoting rule 7 where \(i\in\mathcal{I}\) is sampled with probability proportional to \(\left\|a_i\right\|^2\). This and similar ideas had been considered before as well [30]–[32].

One exceedingly simple randomized rule is just to sample \((i,j)\) uniformly at random among all pairs of indices. We were unable to find a reference that explicitly describes this rule (or any randomized rule) for the Jacobi eigenvalue algorithm; one reason that this is surprising is that the proof of global convergence of the greedy pivoting rule 6 essentially applies the probabilistic method to the hypothetical algorithm which uses this uniformly random pivot (see, for example, Section 5 of [33]).

In the context of orthogonalization, we note an important distinction between sampling \((i,j)\) uniformly at random among all pairs of indices and just pairs of indices with \(i<j\) (or \(\sigma(i)<\sigma(j)\) for some permutation \(\sigma\)). In the latter case, we are assured the computed decomposition (if the method converges) is the QR decomposition \(AP_\sigma=QR\). However, in the former case, we are given no assurances about \(R\) being upper triangular. That rule (without the \(i<j\) constraint) was proposed by Steinerberger [25], and both rules (both with and without the \(i<j\) constraint) were concurrently proposed by the authors in [34]. Both those works were motivated by a designing a preconditioner for the Kaczmarz linear system solver, and made analogies of the method to the Kaczmarz solver itself; neither work identified the connection to the Jacobi eigenvalue algorithm. Steinerberger conjectured that the rate of convergence of this procedure is \(\kappa(A^{(t)})\approx\exp(-O(t/n^2))\kappa(A)\) based on numerical experiments and a heuristic. The authors proved that \((-\log\det|A^{(t)}|)\le\exp(-O(t/n^2))(-\log\det\mleft|A\mright|)\).

2.3 Block variants↩︎

Both Gram Schmidt and Jacobi admit block variants, which can be abstractly described as follows: rather than pick a pair of indices \((i,j)\) at each step and updating columns \(i\) and \(j\) of \(A^{(t)}\), one can instead pick any subset of indices \(J\subset[n]\) and update all the corresponding columns. We call \(J\) the pivot set. In the non-block case, pivoting on \((i,j)\) ensures that \[\left\langle a^{(t+1)}_i,a^{(t+1)}_j\right\rangle=0.\] In the block variant, pivoting on \(J\) ensures that \[\left\langle a^{(t+1)}_i,a^{(t+1)}_j\right\rangle=0\quad\forall i,j\in J,i\neq j.\] These methods are often motivated by considering the practical performance of these algorithms on modern hardware. For example, one may be able to take advantage of large cache sizes to perform one iteration for \(\mleft|J\mright|>2\) very quickly.

The block Jacobi eigenvalue algorithm of Drmač [35] is one such method. It fixes a partition \(I_1\sqcup\cdots\sqcup I_m\subset[n]\) and divides the given matrix \(B\) into blocks corresponding to those index sets. Then, the pivot set \(J\) is selected to be \(I_{\ell}\cup I_{\ell'}\), where the pair \((\ell,\ell')\) can be selected analogously to how \((i,j)\) were selected in the non-block version. Both [36] and [37] show convergence of block Jacobi under various pivoting strategies. There are many block versions of Gram Schmidt along similar lines, see [38] for a recent survey.

3 A generalized factorization↩︎

In describing the connection between Gram Schmidt and the one-sided Jacobi eigenvalue algorithm, we’ve already hinted at what the generalized algorithm is. Once a pivot set \(J\) has been selected, we may abstract away the process by which the corresponding columns are updated: all we assume is that it is an invertible column operation which results in orthogonal columns. We will use the letter \(S\) for the \(\mleft|J\mright|\times \mleft|J\mright|\) matrix performing this column operation. The Jacobi eigenvalue algorithm corresponds to the special case where \(S\) is unitary, and Gram Schmidt to the case where \(S\) is upper triangular. The general pseudo-code for this and the two-sided version is below. We first clarify some notation. \(\text{GL}_k({\mathbb{C}})\) is the set of all invertible \(k\times k\) matrices over \({\mathbb{C}}\). \(\binom{[n]}k\) is the set of all subsets of \([n]\) of size \(k\). For an index set \(J={\left\{{i_1,\ldots,i_k}\right\}}\in\binom{[n]}k,i_1<\cdots<i_k\), denote the submatrices \[A_J=\begin{bmatrix} a_{1,i_1} &\cdots& a_{1,i_k}\\ \vdots&&\vdots\\ a_{d,i_1} &\cdots& a_{d,i_k}\\ \end{bmatrix}\in{\mathbb{C}}^{d\times k} \quad\And\quad B_{JJ}=\begin{bmatrix}b_{i_1i_1} &\cdots& b_{i_1i_k} \\ \vdots && \vdots \\ b_{i_ki_1} &\cdots& b_{i_ki_k}\end{bmatrix}\in{\mathbb{C}}^{k\times k}.\] For a permutation \(\sigma\) on \([n]\) let \(P_\sigma\) be the associated permutation matrix, i.e. \(P_\sigma e_{j}=e_{\sigma(j)}\).

We present both versions of the algorithms side-by-side to highlight their mathematical equivalence: one can start with \(B=A^*A\) and maintain \(A^{(t)*}A^{(t)}=B^{(t)}\) by selecting the same \(T\) in line 7 of either algorithm. Because of this correspondence, we will refer to various lines of the algorithm without needing specify which version we are talking about. We also highlight that \(T\) need not be formed and inverted explicitly.

There are many pivot rules one can use for the selection in line 5. This paper uses a randomized pivoting rule where \(J\) is sampled uniformly at random from \(\binom{[n]}k\). We refer to this as “randomized size \(k\) pivoting”.

We obtain a unified analysis of all special cases of Algorithm 1 with randomized size \(k\) pivoting. The mechanism by which \(S\) is selected in line 6 and the additional structure one wishes to impose upon \(S\) do not affect our analysis in any way—in expectation, the rate of convergence depends only on \(k\) and \(n\). For concreteness, 2 lists several decompositions one can obtain by imposing different constraints on \(S\) and \((A_J^{(t)}S)^*(A_J^{(t)}S)=S^*B^{(t)}S\).

4 Exact arithmetic analysis↩︎

This section shows that one- and two-sided versions of Algorithm 1 quickly converge to the factorizations \(A=QT_{\textrm{acc}}\) and \(B=T_{\textrm{acc}}^*DT_{\textrm{acc}}\) where \(Q\) is nearly orthogonal and \(D\) is nearly diagonal. Concretely, we have the following results, stated in terms of the diagonal-normalized condition number \[\label{a20} \hat{\kappa}(B) = \kappa(D_B^{-1/2}BD_B^{-1/2})\tag{8}\] where \(\kappa(B)=\left\|B\right\|\cdot\|B^{-1}\|\) and \(D_B = B \odot I\). This condition number can be arbitrarily smaller than the usual condition number \(\kappa(\cdot)\) (consider, for example, any diagonal matrix), but cannot be much larger, \[\hat{\kappa}(B)\le n\kappa(B).\]

4.1 Standard Jacobi analysis↩︎

Our analysis of Algorithm 1 follows a similar strategy as the classical analysis of the Jacobi eigenvalue algorithm. To highlight the connection, we review the classical analysis here. Say that \(B=B^{(0)},B^{(1)},B^{(2)},\ldots\) is the sequence of matrices produced by the iterations of the Jacobi eigenvalue algorithm. The proof strategy is to identify a monotone invariant, i.e. some quantity which strictly decreases as iterations are performed. In the classical analysis, this monotone invariant is taken to be the distance of the current iterate \(B^{(t)}\) to the closest diagonal matrix in Frobenius norm, denoted \(\mathop{\mathrm{off}}(\cdot)\). This quantity can be compactly expressed as \[\mathop{\mathrm{off}}(B) :=\inf_{D\text{ diagonal}}\left\|B-D\right\|_F^2 =\sum_{i,j\in[n],i\neq j}\mleft|b_{ij}\mright|^2.\] A straightforward calculation reveals that that\[{\label{a21}\mathop{\mathrm{off}}(B^{(t+1)})=\mathop{\mathrm{off}}(B^{(t)})-2|b^{(t)}_{ij}|^2}\tag{9}\] where \((i,j)\) is the pivot. By appealing to the probabilistic method, if \((i,j)\) is picked according to the greedy pivot rule (i.e. \(|b^{(t)}_{ij}|^2=\mathop{\mathrm{arg\,max}}_{i,j:i\neq j}|b^{(t)}_{ij}|^2\)), then one has \[|b^{(t)}_{ij}|^2\ge\frac{1}{n(n-1)}\mathop{\mathrm{off}}(B^{(t)})\] which ensures global convergence of \(\mathop{\mathrm{off}}(B^{(t)})\to0\) at a linear rate. By using a randomized size 2 pivoting strategy, one has equality in expected value, \[{\label{a22}\mathop{\mathrm{\mathbb{E}}}\mathop{\mathrm{off}}(B^{(t+1)})=\mleft(1-\frac{2}{n(n-1)}\mright)\mathop{\mathrm{off}}(B^{(t)}),}\tag{10}\] which again ensures global convergence in expectation at a linear rate.

In the more general setting where \(B^{(t)}\) is the sequence produced by non-unitary transformations, we no longer have 9 (and therefore not 10 either). Instead, we obtain an expression similar to 10 but with a different function in place of \(\mathop{\mathrm{off}}(\cdot)\).

4.2 Our analysis↩︎

Our potential function is the following, \[{\label{a23} \Gamma(B):=\mathop{\mathrm{tr}}(B\odot B^{-1}-I) =\sum_{i\in[n]}(b_{ii}\cdot(B^{-1})_{ii}-1) =-\sum_{i,j\in[n],i\neq j}b_{ij}\cdot(B^{-1})_{ji}.}\tag{11}\] It is not immediately clear that convergence of \(\Gamma(\cdot)\to0\) implies convergence of positive definite \(B^{(t)}\) to diagonal, and in fact this is not true for indefinite matrices. For definite matrices, this implication does hold; see Section 6 for several proofs and quantitative statements of this implication.

The analysis to show \(\Gamma(B^{(t)}) \to 0\) consists of two lemmas: the first is an update formula analogous to 9 for \(\Gamma(\cdot)\), and the second is the computation of its expectation under the randomized pivoting rule analogous to 10 .

As promised, we now compute the expected value of the bound from Lemma 1 under the randomized size \(k\) pivot rule.

An immediate consequence is global linear convergence of Algorithm 1.

5 Finite arithmetic analysis↩︎

The previous section showed that every instantiation of Algorithm 1 (i.e. every process for selecting \(S\) in line 6) with randomized size \(k\) pivoting has the same expected behavior in terms of \(\Gamma(\cdot)\). To obtain a convergence result in the presence of round-off error, one needs to assume the updates line 8 and 9 are performed with small mixed forward-backward error (see the diagram 13 ).

In this setting, \(\Gamma(\cdot)\) will play two roles: the first role, as before, is that the convergence of \(\Gamma(\cdot)\to0\) tracks the convergence of the algorithm. The second role is in controlling the accumulation of the error acquired during each iteration.

We let \(\widetilde{A}^{(t)}\), \(\widetilde{B}^{(t)}\) denote the state of the one- and two-sided algorithm respectively after \(t\) iterations are performed in finite precision with randomized size \(k\) pivoting (to contrast with \(A^{(t)}\), \(B^{(t)}\) without tildes computed exactly). The stability of each iteration can be defined in terms of the following commutative diagram, corresponding to a mixed forward-backward guarantee. \[\label{a31} \begin{tikzcd} {\widetilde{A}^{(t)}} & {\widetilde{A}^{(t+1)}} \\ {\widetilde{A}^{(t+1/3)}} & {\widetilde{A}^{(t+2/3)}} \arrow["\textrm{floating}", from=1-1, to=1-2] \arrow["\approx"{description}, tail reversed, from=1-1, to=2-1] \arrow["\textrm{exact}", from=2-1, to=2-2] \arrow["\approx"{description}, tail reversed, from=1-2, to=2-2] \end{tikzcd} \quad\quad\quad\quad\quad \begin{tikzcd} {\widetilde{B}^{(t)}} & {\widetilde{B}^{(t+1)}} \\ {\widetilde{B}^{(t+1/3)}} & {\widetilde{B}^{(t+2/3)}} \arrow["\textrm{floating}", from=1-1, to=1-2] \arrow["\approx"{description}, tail reversed, from=1-1, to=2-1] \arrow["\textrm{exact}", from=2-1, to=2-2] \arrow["\approx"{description}, tail reversed, from=1-2, to=2-2] \end{tikzcd}\tag{13}\] To avoid ambiguity, we always take \(t\) to be an integral index and \(s\) an index in \(\frac{1}{3}{\mathbb{N}}\). The diagram commutes in the following sense: if \(\widetilde{A}^{(t+1)}\) is the result of performing one iteration on \(\widetilde{A}^{(t)}\) in floating point arithmetic, then we are guaranteed the existence of “intermediate states” \(\widetilde{A}^{(t+1/3)}\) which is close to \(\widetilde{A}^{(t)}\) and \(\widetilde{A}^{(t+2/3)}\) which is close to \(\widetilde{A}^{(t+1)}\) such that \(\widetilde{A}^{(t+2/3)}\) is the result of applying one iteration to \(\widetilde{A}^{(t+1/3)}\) in exact arithmetic. Analogous states exist for the two-sided variant. We define “closeness” by \[A\approx A'\iff \mathop{\mathrm{dist}}_{\textrm{o}}(A,A')\le\varepsilon\quad\And\quad B\approx B'\iff \mathop{\mathrm{dist}}_{\textrm{t}}(B,B')\le\varepsilon\] for pseudo-metrics \[\mathop{\mathrm{dist}}_{\textrm{o}}(A,A')= \sup_j \left\|\frac{a_j}{ \left\|a_j\right\| }-\frac{a_j'}{ \|a_j'\| }\right\| \quad\And\quad \mathop{\mathrm{dist}}_{\textrm{t}}(B,B')=\sup_{ij}\mleft| \frac{b_{ij}}{\sqrt{b_{ii}b_{jj}}} - \frac{b'_{ij}}{\sqrt{b'_{ii}b'_{jj}}}\mright|\] where \(a_j,a'_j\) is the \(j\)th column of \(A,A'\) respectively. Note that we are implicitly assuming here that the diagonal entries of \(B\) and \(B'\) are positive. Importantly, \(\mathop{\mathrm{dist}}_{\textrm{o}}(A,A')\le\varepsilon\) implies that \(\mathop{\mathrm{dist}}_{\textrm{t}}(A^T A,A'^TA')\le2\varepsilon+\varepsilon^2=O(\varepsilon)\) so it suffices to prove our theorem under the hypothesis of just the two-sided assumption of 13 .

5.1 Control on \(\Gamma(\widetilde{B}^{(t)})\)↩︎

We begin by proving a finite-arithmetic analog of Theorem 1, i.e., we wish to show \(\Gamma(\widetilde{B}^{(t)})\) exhibits linear convergence. In the last section, specifically, 12 , we showed that the stochastic process \(C_{n,k}^{-t}\Gamma(B^{(t)})\) was a martingale where \(C_{n,k}=\mleft(1-\frac{k(k-1)}{n(n-1)}\mright)\). Thus \(C^{t}\to0\) roughly implied \(\Gamma(B^{(t)})\to0\) at the same rate. In this section, we define a new martingale and show \(C_{n,k}^{-t}\Gamma(\widetilde{B}^{(t)})\) is (nearly) bounded by it with high probability. To define the martingale, we introduce the following random variables.\[X_t=\frac{\Gamma(\widetilde{B}^{(t+2/3)})}{\Gamma(\widetilde{B}^{(t+1/3)})} \quad\And\quad Y_{t_1t_2}=\prod_{t=t_1}^{t_2}X_t\] where we use the convention \(Y_{t_1t_2}=1\) if \(t_2<t_1\).

Our next lemma allows us to control the change in \(\Gamma(\cdot)\) between the states \(\widetilde{B}^{(t\pm 1/3)}\) and the state \(\widetilde{B}^{(t)}\).

To stitch these two lemmas together, we essentially argue that by selecting \(\varepsilon\) sufficiently small, one can ensure that the additive error one obtains from Lemma 4 between states \(\widetilde{B}^{(t-1/3)}\) to \(\widetilde{B}^{(t+1/3)}\) does not accumulate too much.

As we stated at the outset, \(\Gamma(\cdot)\) plays two roles: it must converge to \(0\) to show the algorithm halts, and it also must never see a large intermediate value to ensure the numerical error does not blow up. To these ends, we have two consequences of Theorem 2. In both, we set \[\hat{B}^{(t)}= \mathop{\mathrm{diag}}\mleft(\widetilde{b}_{11},\ldots,\widetilde{b}_{nn}\mright)^{-1/2} \widetilde{B}^{(t)} \mathop{\mathrm{diag}}\mleft(\widetilde{b}_{11},\ldots,\widetilde{b}_{nn}\mright)^{-1/2}\] to be the diagonal-normalization of a matrix.

5.2 Orthogonalization is stable↩︎

The most general result is that if each iteration of the one-sided algorithm satisfies 13 , then the overall algorithm stably computes an orthonormal basis of the column-space of \(A=\begin{bmatrix}a_1&\cdots&a_n\end{bmatrix}\in{\mathbb{C}}^{d\times n}\). This can be extended to a basis of the flag \(\mathop{\mathrm{span}}(a_1), \mathop{\mathrm{span}}(a_1,a_2), \mathop{\mathrm{span}}(a_1,a_2,a_3),\cdots.\) We give a handful of examples of implementations of the algorithm which satisfy 13 .

We quantify error in terms of the following metric on subspaces: \[\mathop{\mathrm{dist}}(V,W)=\max\mleft( \sup_{v\in V,\left\|v\right\|=1}\inf_{w\in W}\left\|v-w\right\|, \sup_{w\in W,\left\|w\right\|=1}\inf_{v\in V}\left\|v-w\right\|\mright).\]

5.3 Diagonalization and SVD are stable↩︎

Demmel and Veselić study the stability of the Jacobi eigenvalue algorithm with an arbitrary pivoting rule. Their results make two assumptions. The first is that the pivoting rule is such that the method converges after some specified number of steps. The second is that the diagonal-normalized condition number of any iterate remains bounded. Under our randomized pivoting rule, both properties are satisfied. We restate the results of [1] concerning eigenvalue and singular value error in terms of our notation.

Our results allow us to provide a concrete upper bound for randomized size 2 pivoting.

Similar bounds for the eigenvectors and singular vectors can also be found in [1] (see their Theorems 3.3 and 4.3). We omit repeating the precise statements, but note that the application of Corollary 4 provides concrete bounds analogous to the bounds for eigenvalues and singular values.

Acknowledgments↩︎

We thank Nikhil Srivastava who brought the Jacobi eigenvalue algorithm to our attention when discussing extensions of our earlier work [34]. We also thank Jim Demmel and Ryan Schneider for helpful conversations about the Jacobi eigenvalue algorithm.

6 Characterization of \(\Gamma(\cdot)\)↩︎

Note that \(\Gamma(\cdot)\) is invariant under diagonal scaling, i.e. \[\Gamma(B)=\Gamma(DBD)\] for any diagonal \(D\). In this section, we fix a positive definite \(B\in{\mathbb{C}}^{n\times n}\) and let \[B'= \mathop{\mathrm{diag}}\mleft(b_{11},\cdots,b_{nn}\mright)^{-1/2} B \mathop{\mathrm{diag}}\mleft(b_{11},\cdots,b_{nn}\mright)^{-1/2}\] be the “diagonal normalization” of \(B\), so that \(B'\) has all 1s along the diagonal and \[\Gamma(B)=\Gamma(B')=\mathop{\mathrm{tr}}(B'^{-1})-n.\] We let \(\lambda_j(B')\) be the \(j\)th largest eigenvalue of \(B'\).

A few quantifiable versions of this lemma are available. They each can be derived by considering an optimization program of the eigenvalues \(\lambda_j\) of \(B'\). The constraints are the following: \[\label{a43}\begin{align} \lambda_1\ge\cdots\ge\lambda_n>0\quad\And\quad\lambda_1+\cdots+\lambda_n=n\quad\And\quad\lambda_1^{-1}+\cdots+\lambda_n^{-1}=n+\Gamma(B').\end{align}\tag{18}\]

References↩︎