Oleg Balabanov
email: oleg.balabanov@inria.fr.

Laura Grigori
email: laura.grigori@inria.fr.

Randomized block Gram-Schmidt process for solution of linear systems and eigenvalue problems

Keywords: Gram-Schmidt, QR factorization, randomization, sketching, numerical stability, rounding errors, loss of orthogonality, multi-precision arithmetic, block Krylov subspace methods, Arnoldi iteration, Petrov-Galerkin, Rayleigh-Ritz, generalized minimal residual method

, our algorithm has an advantage of performing standard high-dimensional operations, that define the overall computational cost, with a unit roundoff independent of the dominant dimension of the matrix. This unique feature makes the methodology especially useful for large-scale problems computed on low-precision arithmetic architectures.

Block algorithms are advantageous in terms of performance as they are mainly based on cache-friendly matrix-wise operations, and can reduce communication cost in high-performance computing. The block Gram-Schmidt orthogonalization is the key element in the block Arnoldi procedure for the construction of Krylov basis, which in its turn is used in GMRES and Rayleigh-Ritz methods for the solution of linear systems and clustered eigenvalue problems. In this article, we develop randomized versions of these methods, based on the proposed randomized Gram-Schmidt algorithm, and validate them on nontrivial numerical examples.

Introduction

Let W ∈ R n×m be a matrix with a moderately large number of columns, so that m n. We consider column-oriented block Gram-Schmidt (BGS) algorithms for QR factorization of W:

W = QR
where Q ∈ R n×m has 2 -orthonormal (or very well-conditioned) columns such that range(Q) = range(W), and R ∈ R m×m is upper triangular with positive diagonal entries. Block algorithms are usually based on matrix-wise BLAS-3 operations allowing proper exploitation of modern cache-based and high-performance computational architectures. The BGS orthogonalization forms a skeleton for block Krylov subspace methods for solving clustered eigenvalue problems as well as linear systems with multiple right-hand sides. It is also used in s-step, enlarged and other communication-avoiding Krylov subspace methods [START_REF] Grigori | Enlarged Krylov subspace conjugate gradient methods for reducing communication[END_REF][START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF]. Please see [START_REF] Carson | An overview of block Gram-Schmidt methods and their stability properties[END_REF] and the references therein for an extensive overview of BGS variants, and [START_REF] Baker | On improving linear solver performance: A block variant of GMRES[END_REF][START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF][START_REF] Stewart | A Krylov-Schur algorithm for large eigenproblems[END_REF][START_REF] Zou | GMRES algorithms over 35 years[END_REF] for the underlying block Krylov methods.

A novel randomized Gram-Schmidt (RGS) process based on the random sketching technique (see [START_REF] Martinsson | Randomized numerical linear algebra: Foundations and algorithms[END_REF][START_REF] Vershynin | High-dimensional probability: An introduction with applications in data science[END_REF][START_REF] Woodruff | Sketching as a tool for numerical linear algebra[END_REF] and the references therein) was recently introduced in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF]. It requires nearly half as many flops and passes over the data than the classical Gram-Schmidt (CGS) process, and is very stable, in the sense that it yields a very well-conditioned Q factor if cond(W) = O(u -1 ), where u is the unit roundoff. While the RGS algorithm can be already very beneficial under unique precision computations, there are even more benefits that can be gained by working in two precisions: using

Preliminaries

Throughout this work, we use the following notations, which is an adaptation of the notations from [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF] to block linear algebra. As in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF], we denote Euclidean vectors x by bold lowercase letters. A matrix composed of several column vectors x 1 , . . . , x k is denoted with the associated bold capital letter and a subscript specifying the number of columns, i.e., X k . If k is constant, this notation can be simplified to X. The (i, j)-th block of a block matrix X is denoted by X (i,j) i.e, we have

X =      X (1,1) X (1,2) • • • X (1,p) X (2,1) X (2,2) • • • X (2,p)
. . . . . . . . .

X (l,1) X (l,2) X (l,p)      ,
for some p and l. When l = 1 (i.e., matrix X is partitioned column-wise), the notation X (1,j) can be simplified to X (j) . The sub-block of X composed of blocks X (i,j) with N 1 ≤ i ≤ N 2 and M 1 ≤ j ≤ M 2 , is denoted by X (N1:N2,M1:M2) . We denote the matrix X (N1:N2,M1:M1) by simply X (N1:N2,M1) . Moreover, if l = 1, then matrix X (1:1,M1:M2) is denoted by X (M1:M2) , and matrix X (M1:M1) by X (M1) . The minimal and the maximal singular values of X are denoted by σ min (X) and σ max (X), and the condition number by cond(X). We let •, • and • = σ max (•) be the 2 -inner product and 2 -norm, respectively. • F denotes the Frobenius norm. For two matrices (or vectors) X and Y, the relation X ≤ Y indicates that the entries of X satisfy x i,j ≤ y i,j . Furthermore, for a matrix (or a vector) X, we denote by |X| the matrix Y with entries y i,j = |x i,j |. We also let X T and X † respectively indicate the transpose and the Moore-Penrose inverse of X. Finally, we let I be the identity matrix of size appropriate to the expression where this notation is used. A quantity or an arithmetic expression X computed with finite precision arithmetic is denoted by fl(X) or X.

Block Gram-Schmidt process

Let matrix W ∈ R n×m be partitioned into p blocks W (i) ∈ R n×mp , with 1 ≤ i ≤ p, m = m p p:

W = W (1:p) = W (1) W (2) . . . W (p) .
BGS process proceeds recursively, at iteration i, selecting a new block matrix W (i) and orthogonalizing it with the previously orthogonalized blocks yielding a matrix Q (i) , followed by orthogonalization of Q (i) itself. This procedure is summarized in Algorithm 1.

Algorithm 1 Block Gram-Schmidt process

Given: n × m block matrix W = W (1:p) , m ≤ n Output: n × m factor Q = Q (1:p) and m × m upper triangular factor R = R (1:p,1:p) . for i = 1 : p do 1. Compute a projection Q (i) = Π (i-1) W (i) (also yielding R (1:i-1,i) ).

2. Compute QR factorization Q (i) R (i,i) = Q (i) with suitable efficient routine. end for For standard methods, the projector Π (i-1) in Algorithm 1 is taken as an approximation to 2 -orthogonal projector onto range(Q (1:i-1) ) ⊥ . For the classical BGS process (BCGS), one chooses Π (i-1) as

Π (i-1) = I -Q (1:i-1) Q T (1:i-1)
. Whereas, for the modified BGS (BMGS) process, we have

Π (i-1) = (I -Q (i-1) (Q (i-1) ) T )(I -Q (i-2) (Q (i-2) ) T ) . . . (I -Q (1) (Q (1) ) T ). (1)
In infinite precision arithmetic, these two projectors are equivalent, and are exactly equal to the 2 -orthogonal projector onto range(Q (1:i-1) ) ⊥ , since, by construction, Q (1:i-1) is an orthogonal matrix. However, in the presence of rounding errors these projectors may cause instability. In the first case, the condition number of Q can grow as cond(W) 2 or even worse and thus requires special treatment [START_REF] Carson | The Stability of Block Variants of Classical Gram-Schmidt[END_REF]. While in the second case, cond(Q) can grow as cond(W) max 1≤j≤p cond(W (j) ) unless Step 2 is unconditionally stable [START_REF] Barlow | Block Modified Gram-Schmidt Algorithms and Their Analysis[END_REF][START_REF] Carson | An overview of block Gram-Schmidt methods and their stability properties[END_REF]. The stability of these processes can be improved by re-orthogonalization, i.e., by running the inner loop twice. In particular, it can be shown that the BCGS process with re-orthogonalization, here denoted by BCGS2, yields an almost orthogonal Q factor as long as the matrix W is numerically non-singular [START_REF] Barlow | Reorthogonalized block classical Gram-Schmidt[END_REF].

Inter-block orthogonalization is an important step in BGS algorithm. For standard versions, this step can be performed with any suitable efficient and stable routine for ( 2 -)QR factorization of tall and skinny matrices, such as TSQR from [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF], a Gram-Schmidt QR applied l times, Cholesky QR applied l times, and others. Typically, such QR factorization takes only a fraction of the overall computational cost.

Random sketching

In this subsection we recall the basic notions of the random sketching technique, and then provide some useful results from [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF].

Let Θ ∈ K k×n with k ≤ n be a sketching matrix such that the associated sketched product Θ•, Θ• approximates well the 2 -inner product between any two vectors in the subspace (or subspaces) of interest V ⊂ R n . Such matrix Θ is referred to as a ε-embedding for V , as defined below.

Definition 1.1. For ε < 1, we say that the sketching matrix Θ is a ε-embedding for subspace V (or matrix V spanning V ), if it satisfies the following relation:

∀x, y ∈ V, | x, y -Θx, Θy | ≤ ε x y . (2)
Furthermore, we assume that Θ is chosen at random from a certain distribution, such that it satisfies (2) for any fixed d-dimensional subspace with high probability (see Definition 1.2).

Definition 1.2. A random sketching matrix

Θ is called a (ε, δ, d) oblivious 2 -subspace embedding, if for any fixed V ⊂ R n of dimension d, it satisfies P(Θ is a ε-embedding for V ) ≥ 1 -δ Corollary 1.3 (Corollary 2.4 in [3]). If Θ ∈ R k×n is a (ε, δ/n, 1
) oblivious 2 -subspace embedding, then with probability at least 1 -δ, we have

Θ F ≤ (1 + ε)n.
Corollary 1.4 (Corollary 2.2 in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF]). If Θ ∈ R k×n is an ε-embedding for V, then the singular values of V are bounded by

(1 + ε) -1/2 σ min (ΘV) ≤ σ min (V) ≤ σ max (V) ≤ (1 -ε) -1/2 σ max (ΘV).
In recent years, several distributions of Θ have been proposed that satisfy Definition 1.2 and have a small first dimension k that depends at most logarithmically on the dimension n and probability of failure δ. Among them, the most suitable distribution should be selected depending on the problem and computational architecture. The potential of random sketching is here realized on (rescaled) Rademacher matrices and Subsampled Randomized Hadamard Transform (SRHT). The entries of a Rademacher matrix are i.i.d. random variables satisfying

P(θ i,j = 1/ √ k) = P(θ i,j = -1/ √ k) = 1/2.
Rademacher matrices can be efficiently multiplied by vectors and matrices through the proper exploitation of computational resources such as cache or distributed machines. For n, which is a power of 2, SRHT is defined as a product of diagonal matrix of random signs with Walsh-Hadamard matrix, followed by uniform sub-sampling matrix and scaling factor 1/ √ k. For a general n, the SRHT has to be combined with zero padding to make the dimension a power of 2. Random sketching with SRHT can reduce the complexity of an algorithm. Products of SRHT matrices with vectors require only n log 2 n flops using the fast Walsh-Hadamard transform or 2n log 2 (k + 1) flops following the methods in [START_REF] Ailon | Fast dimension reduction using Rademacher series on dual BCH codes[END_REF]. Furthermore, for both distributions, the usage of a seeded random number generator can allow efficient storage and application of Θ. It follows from [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF] (based on [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF][START_REF] Woodruff | Sketching as a tool for numerical linear algebra[END_REF]) that the rescaled Rademacher distribution, and SRHT (possibly with zero padding) respectively are (ε, δ, d) oblivious 2 -subspace embeddings, if they have a sufficiently large first dimension: if

k ≥ 7.87ε -2 (6.9d + log(1/δ)) (3a)
for Rademacher matrices, or if

k ≥ 2(ε 2 -ε 3 /3) -1 √ d + 8 log(6n/δ) 2 log(3d/δ) (3b)
for SRHT.

Effect of random sketching on rounding errors

Let us now discuss a very useful result from [3, Section 2.2] characterizing the rounding errors in a sketched matrixvector product, which can be easily extended to matrix-matrix products. In short, this result states that multiplying Θ by a matrix-vector product x = fl(Yz) does not increase the rounding error bound of x by more than a small factor. This property is nothing else than a sketched version of the standard "rule of thumb" for randomized algorithms. It was rigorously proven in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF] for the probabilistic rounding model, [START_REF] Connolly | Stochastic Rounding and its Probabilistic Backward Error Analysis[END_REF]Model 4.7]. The extension of this result to matrix-matrix products is provided below. Let us fix a realization of an oblivious 2 -subspace embedding Θ ∈ R k×n of sufficiently large size, and consider a matrix-matrix product

X = YZ, with Y ∈ R n×m , Z ∈ R m×l ,
computed in finite precision arithmetic with unit roundoff u < 0.01/m. The following results can be derived directly from [3, Section 2.2] with the observation that each column of X is a matrix-vector product:

x i = Yz i , 1 ≤ i ≤ l. We have, |X -X| ≤ U, (4)
for some matrix U describing the "worse-case scenario" rounding error. In general, the standard analysis (e.g., see [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]) gives the bound

|X -X| ≤ mu 1 -mu |Y||Z| ≤ 1.02mu|Y||Z|, (5)
which means that U = 1.02mu|Y||Z| satisfies (4). Furthermore, as in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF], if YZ represents a sum of a matrix-matrix product with a matrix, i.e., YZ = Y Z + H, then one can take

U = 1.02u(|H| + m|Y ||Z |).
Let us now address bounding the norm of the rounding error after sketching X. We have the following "worse-case scenario" bound

Θ(X -X) F ≤ Θ X -X F ≤ Θ U F , (6)
which, combined with Corollary

1.3, implies that if Θ is (ε, δ/n, 1) oblivious 2 -subspace embedding, then Θ(X -X) F ≤ √ 1 + ε √ n U F holds with probability at least 1 -δ.
The next important point is that, as was argued in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF], this bound is pessimistic and can be improved by a factor of O( √ n) by exploiting statistical properties of rounding errors. We summarize the results from [3, Theorem 2.5 and Corollary 2.6] below.

Corollary 1.5. Consider probabilistic rounding model, where the rounding errors due to each elementary arithmetic operation are random variables possibly depending on each other, but are independently centered. Furthermore assume that the errors are bounded so that, it holds

|X -X| ≤ U. If Θ is a (ε/4, l -1 n d -1 δ, d) oblivious 2 -subspace embedding, with d = 4.2c -1 log(4/δ), where c ≤ 1 is some universal constant, then Θ(X -X) F ≤ √ 1 + ε U F (7)
holds with probability at least 1 -2δ.

By using (3) we deduce that for l ≤ n the relation ( 7) is satisfied with probability at least 1 -2δ if Θ is a Rademacher matrix with O(log(n) log(1/δ)) rows or P-SRHT matrix with O(log 2 (n) log 2 (1/δ)) rows. This result implies that in practice, the matrix-vector products in randomized algorithms can be performed with a unit roundoff independent of the (high) dimension n.

Randomized block Gram-Schmidt process

Consider BGS algorithms with projectors Π (i-1) that respect the relation

Π (i-1) W (i) = W (i) -Q (1:i-1) X, where X = R (1:i-1,i) is computed from W (i) and Q (1:i-1)
. Standard algorithms take X as an approximate solution to the following minimization problem: min

Y Q (1:i-1) Y -W (i) F . (8)
For instance, BCGS algorithm approximates the exact solution to (8) by Q T (1:i-1) W (i) , whereas BMGS projector (1) improves this approximation under finite precision arithmetic, though it may cause a computational overhead in terms of inter-processor communication and operation with cache/RAM.

As proposed in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF], a reduction of the computational cost and/or improvement of stability of Gram-Schmidt process can be obtained with a projector that gives a Q factor orthonormal with respect to Θ•, Θ• instead of the 2 -inner product as in standard methods. In our case, this corresponds to taking X as an approximate solution to

min Y ΘQ (1:i-1) Y -ΘW (i) F , (9)
which is a k-dimensional block least-squares problem. Since the sketching dimension satisfies k n, hence a very accurate solution to [START_REF] Carson | An overview of block Gram-Schmidt methods and their stability properties[END_REF] should be more efficient to compute than any, even the BCGS (i.e, the cheapest) solution to [START_REF] Carson | The Stability of Block Variants of Classical Gram-Schmidt[END_REF]. Furthermore, with the right choice of random sketching matrices, the precomputation of the sketches of Q (1:i-1) and W (i) should also have only a minor cost compared to that of standard high-dimensional operations. If [START_REF] Carson | An overview of block Gram-Schmidt methods and their stability properties[END_REF] is very small, its solution can be obtained with a direct solver based, for example, on Householder transformations or Givens rotations, requiring a cubic complexity. If the problem has a moderate size, so that the direct solvers have a considerable computational cost, one can use the fact that matrix ΘQ (1:i-1) is almost orthonormal, which implies applicability of iterative solvers requiring a quadratic complexity. Several such solvers are provided in Section 2.1. For the QR orthogonalization with respect to Θ•, Θ• we have to also use the corresponding inter-block orthogonalization routine in Step 2 of Algorithm 1, providing a Q factor orthonormal with respect to the sketched inner product. An efficient procedure for this task can be based on a standard 2 -orthogonalization followed by a generalized Cholesky QR factorization (see Section 2.2). On the other hand, this task can also be performed with the single-vector randomized Gram-Schmidt algorithm from [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF]. The RBGS process using the random sketching projector is depicted in Algorithm 2.

Algorithm 2 Randomized block Gram-Schmidt algorithm (RBGS)

Given: n × m block matrix W = W (1:p) , and k × n matrix Θ, m ≤ k n. Output: n × m factor Q = Q (1:p) and m × m upper triangular factor R = R (1:p,1:p) . for i = 1 : p do 1. Sketch W (i) : P (i) = ΘW (i) .
# macheps: u f ine 2. Solve small block least-squares problem:

R

(1:i-1,i) = arg min Y S (1:i-1) Y -P (i) F . # macheps: u f ine 3. Compute projection of W (i) : Q (i) = W (i) -Q (1:i-1) R (1:i-1,i) .
# macheps: ucrs 4-5. Compute QR factorization

Q (i) R (i,i) = Q (i) with Q (i) orthogonal with respect to Θ•, Θ• . # macheps: u f ine Compute sketch of Q (i) : S (i) = ΘQ (i) .
# macheps: u f ine 6. (Optional) compute S (i) = ΘQ (i) , I -S (i) T S (i) F and

S (i) -S (i) R (i,i) F S (i) F
. At the iteration i = 1 in Algorithm 2 we used the notation that R (1:i-1,i) is a 0 by m p matrix and Q (1:i-1) is a n by 0 matrix, so that Q (i) = W (i) . Algorithm 2 is presented under multi-precision arithmetic using two unit roundoffs, like its single-vector counterpart from [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF]. The working precision is represented by a coarse roundoff u crs . It is used for standard high-dimensional operations in Step 3, which determine the overall computational cost. All other (inexpensive) operations in Algorithm 2 are computed with a fine unit roundoff u f ine , u f ine ≤ u crs . It is shown in Section 3 that Algorithm 2 is stable if u crs ≤ O(cond(W)m -2 ), which is a very mild condition on u crs . The fact that this bound is independent of the high dimension n explains the potential of our methodology for large-scale problems computed on low-precision arithmetic architectures. Furthermore, according to our numerical experiments, the RBGS algorithm can be sufficiently stable even when u crs is larger than O(cond(W)m -2 ). In such cases, the stability of the algorithm can be certified by a posteriori bounds given by the quantities ∆ (p) , ∆(p) , etc., computed in (optional) Steps 6 and 7. These bounds provide stability certification if u crs = O(m -3/2 ), which is a milder condition than the one for a priori guarantees, and is independent not only of n but also cond(W).

The stability guarantees of Algorithm 2 executed in unique precision can be obtained directly from the analysis of the multi-precision algorithm by taking u f ine = u and u crs = F (m, n)u, where F (m, n) is some polynomial of low degree.

In terms of performance, the SRHT-based RBGS algorithm requires about half the flops and data passes of the cheapest standard block Gram-Schmidt algorithm, which is BCGS. Moreover, its computational cost is defined by p wellparallelizable BLAS-3 operations. As in BCGS, in RBGS only one global synchronization between distributed processors is required per iteration, not counting the cost of inter-block QR.

Remark 2.1. If necessary, the output of the RBGS algorithm can be post-processed with a Cholesky QR to provide a QR factorization with Q factor, which is not only well-conditioned but is 2 -orthonormal up to machine precision. More specifically, we can compute an upper-triangular matrix R such that R T R = Q T Q with Cholesky decomposition, and consider

Q ← QR -1 and R ← R R.
The numerical stability of this procedure follows directly from the fact that the Q factor produced by the RBGS algorithm (and therefore R ) is very well-conditioned. Its computational cost is dominated by the computation of Q T Q and, possibly, QR -1 . In terms of flops, the computational cost of computing Q T Q is similar to the cost of the entire RBGS algorithm.

Though as a single BLAS-3 operation, it is better suited for cache-based and parallel computing. The computation of QR -1 can be omitted when the application allows operating with the Q factor in an implicit form, for instance, in Arnoldi iteration or Rayleigh-Ritz algorithm. Otherwise, this product can be computed by well-parallelizable forward substitution, or even a direct inversion allowing a second BLAS-3 multiplication. In this case, the computation and application of R -1 should not cause any stability issues, since R is a very well-conditioned matrix.

Solution of block least-squares problem in Step 2

As already pointed out, the stability of Algorithm 2 strongly depends on the stability of the least-squares solver used in Step 2. In particular, in our analysis in Section 3 we require the solution X = R (1:i-1,i) to satisfy the following backward-stability condition.

Assumptions 2.2 (Backward-stability of Step 2 of RBGS).

X = arg min Y ( S (1:i-1) + ∆S (i-1) )Y -( P (i) + ∆P (i) )
F for some matrices ∆S (i-1) and ∆P (i) such that

∆S (i-1) F ≤ 0.01u crs S (1:i-1) , ∆P (i) F ≤ 0.01u crs P (i) .
This condition can be met by standard direct solvers based on Householder transformation or Givens rotations and a sufficiently large gap between u crs and u f ine , which follows from [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]Theorems 8.5,19.10 and 20.3] and their proofs. However, direct solvers require a cubic complexity and can become too expensive even for relatively small m and k. In such cases, one has to appeal to iterative methods that can exploit the approximate orthogonality of S (1:i-1) to speedup the computations, as is discussed below.

Richardson iterations

Algorithm 3 describes, perhaps, the simplest iterative method for obtaining a solution in Step 2. When the least-squares problem is seen as orthogonalization of columns of P (i) to S (1:i-1) , this algorithm can be interpreted as nothing more than the classical Gram-Schmidt (CGS) algorithm with l re-orthogonalizations. In more general terms, Algorithm 3 is Richardson method applied to the normal system of equations

S T (1:i-1) S (1:i-1) X = S T (1:i-1) P (i) . (10)
Assuming that l = O(1), it requires O(m p mk) flops for computing R (1:i-1,i) = X, which makes the cost of Step 2 negligible in comparison to other costs of the RBGS algorithm and in particular to the Step 3 of RBGS, which has complexity O(m p mn). Finally, we note that the RBGS algorithm based on l = 1 Richardson iteration in Step 2 is equivalent to the BCGS process for the orthogonalization with respect to the sketched inner product. In general, there is a strong connection between the RBGS algorithm using l Richardson iterations and BCGS using l re-orthogonalizations. Let Y = [y 1 , . . . , y mp ] be the exact solution to [START_REF] Connolly | Stochastic Rounding and its Probabilistic Backward Error Analysis[END_REF] and

R (1:i-1,i) = X = [x 1 , . . . ,
x mp ] be the solution computed with Algorithm 3. In exact arithmetic, we have:

x j -y j ≤ I -S T (1:i-1) S (1:i-1) l y j = (∆ (i-1) ) l y j , 1 ≤ j ≤ m p ,
where ∆ (i-1) = I -S T (1:i-1) S (1:i-1) F measures the orthogonality of the sketch of Q factor. This result can be easily extended to finite precision arithmetic.

Algorithm 3

Step 2 of RBGS, using CGS with l re-orthogonalizations

Given: P (i) , S (1:i-1) Output: R (1:i-1,i) = X 1. Set X = 0 (1:i-1,1) . for 1 : l do 2. X ← X + S T (1:i-1) P (i) -S (1:i-1) X . end for

Iterations of the modified GS process

Another way to compute the solution in Step 2 of Algorithm 2 is to use the modified GS (MGS) process, applied l times. This algorithm can provide an accurate solution in fewer iterations than Algorithm 3, though it can have a higher cost per iteration from a performance standpoint. This drawback can be remedied by appealing to the block version of MGS depicted in Algorithm 4. The RBGS algorithm based on Algorithm 4 can be linked to BMGS algorithm, in which the Algorithm 4 Step 2 of RBGS, using BMGS with l re-orthogonalizations

Given: P (i) , S (1:i-1) Output: R (1:i-1,i) = X = [X (1) , . . . , X (i-1) ] 1. Set X = 0 (1:i-1,1) . for 1 : l do for j = 1 : i -1 do 2. P (i) ← P (i) -S (j) X (j) .
3. X (j) ← X (j) + S T (j) P (i) . end for end for

Other methods

In principle, we can compute X from the normal equation ( 10) using any suitable iterative method such as Conjugate Gradient or GMRES. In this case, the normal matrix can be operated with as an implicit map outputting products with vectors using O(m p mk) flops, making the cost of the iterative method similar to that of a Richardson iteration.

Inter-block QR factorization in Steps 4-5

Let us now provide three ways of efficient and stable inter-block orthogonalization in Steps 4-5 of the RBGS algorithm.

Single-vector RGS algorithm

First, the inter-block orthogonalization can be readily performed with the single-vector RGS algorithm from [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF] computed with unique roundoff u = u f ine = F (m, n)u crs , where F (m, n) is a low-degree polynomial. This algorithm provides a QR factorization that satisfies (11a), which follows directly by construction. Furthermore, it should also satisfy (11b), since Θ F ≤ √ 2n holds with high probability (according to Corollary 1.3).

| Q (i) -Q (i) R (i,i) | ≤ 0.1u crs | Q (i) || R (i,i) | (11a) Θ( Q (i) -Q (i) R (i,i) ) F ≤ 0.1u crs Q (i) R (i,i) (11b) As shown in [3], if Θ is ε-embedding for Q (i) , then RGS algorithm satisfies 1 -0.1u crs cond( Q (i) ) ≤ σ min (Θ Q (i) ) ≤ σ max (Θ Q (i) ) ≤ 1 + 0.1u crs cond( Q (i) ). (12)
Moreover, in this case Θ is guaranteed to be an ε -embedding for Q (i) , with ε = 2ε + u crs . Though, in practice, this property holds for smaller values ε , say, ε + u crs .

According to [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF], the single-vector RGS algorithm applied to an inter-block has the cost of 2m 2 p n flops and m p n bytes of memory, which can be considered negligible compared to the complexity and memory consumption of other computations such as the standard high-dimensional operations in Algorithm 2. The RGS algorithm, however, requires m p global synchronizations between distributed processors, which can dominate the computational costs in cache-based or massively parallel architectures. In such cases, it is necessary to appeal to different approaches for the inter-block QR factorization, some of which are described below.

Sketched CholQR factorization

Another way to perform inter-block QR factorization of matrix Q (i) with respect to the sketched inner product is to appeal to the following sketched version of the Cholesky QR (CholQR). We can first compute the R factor R (i,i) by performing a ( 2 -)QR factorization of a small matrix S (i) = ΘQ (i) with any suitable routine, possibly relying on Cholesky factorization of (S (i) ) T S (i) . Then the Q factor is retrieved by computing

Q (i) = Q (i) (R (i,i) ) -1
with forward substitution (see Algorithm 5). While sketched Cholesky QR is efficient in terms of communication, it can entail stability issues. They can be alleviated by re-orthogonalization (i.e, by running the algorithm l times), or by combination with modern efficient routines for 2 -orthogonalization.

Algorithm 5 Sketched CholQR

Given: Q (i) , Θ. Output: QR fact. Q (i) = Q (i) R (i,i) with Q (i) orthogonal with respect to Θ•, Θ• . Sketch S (i) = ΘQ (i) . 1. Compute S (i) = ΘQ (i) . 2. Compute R (i,i) as the R factor of 2-QR factorization of S (i) . 3. Compute Q (i) = Q (i) (R (i,i) ) -1 with forward substitution (or direct inversion, if R (i,i) is well conditioned). 4. Calculate the sketch as S (i) = S (i) (R (i,i) ) -1 or as S (i) = ΘQ (i) .

2.2.3

2 -QR+sketched CholQR factorization The idea here is to improve the stability of the sketched CholQR by pre-processing the matrix Q (i) with classical routines for 2 -orthogonalization, as shown in Algorithm 6. In principle, the pre-processing step can be done with any suitable routine for 2 -orthogonalization of tall-and-skinny matrices such as the classical or modified Gram-Schmidt with l reorthogonalizations, Householder algorithm, the tall-and-skinny QR from [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF] (TSQR) or any others. In each particular situation, the most suitable routine should be selected depending on the computational architecture and programming environment. For instance, because of their popularity and reliability, the Gram-Schmidt and Householder algorithms are often available in scientific libraries as greatly optimized high-level routines. On the other hand, TSQR is favorable in massively parallel environments, since it reduces the amount of messages/synchronizations between processors. It has to be noted that Algorithm 6 permits the usage of algorithms that provide the Q factor in an implicit form, e.g., equal to the first m p columns of I -YLY T , which is the output of the Householder method. The only requirement is that such implicit Q factor is as efficient to multiply by the sketching matrix Θ, and vectors/matrices, as an explicit Q factor would be. In the rest of the paper, however, the matrix Q (i) is assumed to be computed explicitly.

Algorithm 6 RBGS, Steps 4-5 ( 2 -QR+CholQR factorization of Q (i) ) Given: Q (i) , Θ. Output: QR fact. Q (i) = Q (i) R (i,i) with Q (i) orthogonal with respect to Θ•, Θ• . Sketch S (i) = ΘQ (i) . 1. Compute 2-QR fact. Q (i) = Q * R with any suitable efficient/stable routine. 2. Use Algorithm 5, taking Q * as Q (i) , to compute sketched CholQR fact. Q * = Q (i) R and sketch S (i) = ΘQ (i) . 3. Calculate R (i,i) = R R .
It can be shown that Algorithm 6 produces a stable QR factorization satisfying [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF] and [START_REF] Grigori | Enlarged Krylov subspace conjugate gradient methods for reducing communication[END_REF], if it is computed with roundoff u = u f ine = F (m, n)u crs , where F (m, n) is some low-degree polynomial. Let us briefly outline the main arguments that support this fact. First, since we use a sufficiently stable routine in Step 1, we have cond(

Q * ) = O(1) and | Q (i) -Q * R | = O(u crs | Q * || R |).
This property, for instance, is satisfied by the classical GS algorithm with l ≥ 2 re-orthogonalizations, Householder QR, TSQR, and others. Furthermore, since Step 2 is performed independently of Θ we can assume that Θ satisfies the ε-embedding property for Q * , and therefore we also have cond(Θ Q * ) ≈ cond( Q * ) = O(1). These properties imply that cond( R ) = O(1), which in its turn implies the stability of CholQR step. Notice that the fact that cond( R ) = O(1), gives an ability to perform Step 3 in the sketched CholQR algorithm with efficient direct inversion instead of the usual forward substitution. It is concluded that Θ Q (i) is approximately orthogonal (i.e., it satisfies [START_REF] Grigori | Enlarged Krylov subspace conjugate gradient methods for reducing communication[END_REF]) [START_REF] Demmel | Communication-optimal parallel and sequential QR and LU factorizations[END_REF].

and that | Q (i) -Q (i) R (i,i) | = O(| Q (i) -Q * R |) = O(u crs | Q * || R |) = O(u crs |Q (i) ||R (i,i) |), which yields
If needed, the output of Algorithm 6 can be certified a posteriori by computing I -S T (i) S (i) F and

Θ Q (i) -S (i) R (i,i) F Θ Q (i) F
, similarly to the RGS algorithm in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF].

Stability analysis

In this section we provide a rigorous stability analysis of the RBGS algorithm. In particular, we show that the RBGS algorithm has similar a priori as well as a posteriori stability guarantees as its single-vector counterpart (see [3, Section 3]).

Assumptions

Our analysis will be based on the following assumptions. We first assume that Θ has a bounded norm as in [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]. This condition is satisfied with probability at least 1 -δ, if Θ is (1/2, δ/n, 1) oblivious subspace embedding (see Corollary 1.3). Let the matrix ∆Q (i) define the rounding error in Step 3:

∆Q (i) = Q (i) -W (i) -Q (1:i-1) R (1:i-1,i) .
Then, the standard worse-case scenario rounding analysis gives [START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF]. Furthermore, we assume that Θ also satisfies [START_REF] Martinsson | Randomized numerical linear algebra: Foundations and algorithms[END_REF] which, according to Corollary 1.5, holds under the probabilistic rounding model with probability at least 1 Finally, it is assumed that in Steps 4-5, the QR factorization and the sketch of the Q factor satisfies ( 16) and ( 17) that can be achieved by using algorithms from Section 2.2.

-4δ, if Θ is (1/8, m -1 n d -1 δ, d)
Assumptions 3.1. It is assumed that for some ε ≤ 1/2, Θ F ≤ √ 1 + ε √ n. ( 13
)
Furthermore, we assume that

|∆Q (i) | ≤ 1.02u crs (| W (i) | + im p | Q (1:i-1) || R (1:i-1,i) |) (14)
and

Θ∆Q (i) F ≤ 1.02u crs √ 1 + ε | W (i) | + im p | Q (1:i-1) || R (1:i-1,i) | F ( 15
)
with 1 ≤ i ≤ p and 1 ≤ t ≤ m p . Finally, it is assumed that in in Steps 4-5 of Algorithm 2, we have

| Q (i) -Q (i) R (i,i) | ≤ u f ine m(| Q (i) | + | Q (i) || R (i,i) |) ≤ 0.1u crs | Q (i) || R (i,i) | (16a) Θ( Q (i) -Q (i) R (i,i) ) F ≤ u f ine m Θ F | Q (i) | + | Q (i) || R (i,i) | ≤ 0.1u crs Q (i) R (i,i) , (16b)
and, if Θ is ε -embedding for

Q (i) , 1 -0.1u crs cond( Q (i) ) ≤ σ min (Θ Q (i) ) ≤ σ max (Θ Q (i) ) ≤ 1 + 0.1u crs cond( Q (i) ), (17)
and

Q (i) ≤ (1 -ε ) -1/2 Θ Q (i) .

Stability guarantees of RBGS algorithm

Our stability analysis will rely on the condition that Θ satisfies the ε-embedding property for W and Q. See Section 3.3 for a characterization of this property.

A posteriori analysis of RBGS algorithm

Let us first give an a posteriori characterization of the RBGS algorithm. Such characterization can be performed by measuring or bounding coefficients ∆ (p) = I -S T S F and ∆(p) = P-S R F P F

, as done in [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF]. We have the following result, which essentially extrapolates [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF]Theorem 3.1] to the block version of the RGS algorithm.

Theorem 3.2. Consider Algorithm 2. Assume that

100m 1/2 n 3/2 u f ine ≤ u crs ≤ 0.01,
along with Assumptions 3.1, possibly excluding [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF].

If Θ is an ε-embedding for Q and W, with ε ≤ 1/2, and if ∆ (p) , ∆(p) ≤ 0.1, then the following inequalities hold:

(1 + ε) -1/2 (1 -∆ (p) -0.1u crs ) ≤ σ min ( Q) ≤ σ max ( Q) ≤ (1 -ε) -1/2 (1 + ∆ (p) + 0.1u crs ) and W -Q R F ≤ 4u crs m 3/2 W F .
Proof. See appendix.

Remark 3.3. In Theorem 3.2, we also have

Θ( W -Q R) F ≤ 5u crs m 3/2 W F .
Proof. See appendix.

Theorem 3.2 implies the numerical stability of the RBGS algorithm if ∆ (p) and ∆(p) are ≤ 0.1. These coefficients can be efficiently computed a posteriori from the sketches S and P and the R factor R, thus providing a way for the certification of the solution. Such certification, in particular, does not involve any assumptions on cond( W) and the accuracy of the least-squares solution in Step 2. Furthermore, we would like to highlight the very mild condition u crs = O(m -3/2 ) on the working (coarse) unit roundoff to guarantee accuracy of the algorithm, which is in particular independent of the high-dimension n.

A priori analysis of RBGS algorithm

Clearly, to get a priori bounds for ∆ (p) and ∆(p) we need more assumptions than in Theorem 3.2. In particular, it is necessary to impose a stability condition on the least squares solver used in Step 2. We also need W to be numerically non-singular, i.e., to satisfy u crs ≤ O(cond( W) -1 ).

The following a priori guarantee of stability of RBGS algorithm is, in fact, a restating of [3, Theorem 3.2] proposed for the single-vector algorithm.

Theorem 3.4. Consider Algorithm 2 with a backward-stable solver (e.g., based on Richardson iterations) satisfying Assumptions 2.2. Under Assumptions 3.1, if Θ is an ε-embedding for Q (1:p-1) and W, with ε ≤ 1/2, and if

u crs ≤ 10 -3 cond( W) -1 m -2 , u f ine ≤ 10 -2 m -1/2 n -3/2 u crs ,
then ∆ (p) and ∆(p) are bounded by

∆(p) ≤ 4.2u crs m 3/2 W F / P F ≤ 6u crs m 3/2 , (18)
∆ (p) ≤ 20u crs m 2 cond( W). (19)
Proof. See appendix.

Theorem 3.4 states that the RBGS algorithm is stable unless the input matrix W is numerically singular. This stability guarantee is optimal and is seen in other stable algorithms such as MGS, CGS2, and others. Furthermore, the stability is proven for working unit roundoff independent of the high dimension n. This property of RBGS algorithm is unique and can be especially interesting for large-scale problems solved on low-precision arithmetic architectures.

Epsilon embedding property

The stability analysis in Section 3 holds if Θ satisfies the ε-embedding property for Q and W. In this section we analyze this property.

We consider the case when W and Θ are independent of each other. Then, if Θ is (ε, δ, m) oblivious 2 -subspace embedding, it satisfies the ε-embedding property for W with high probability. Below, we show that, in this case Θ will also satisfy an ε-embedding property for Q with moderately increased value of ε. This result is basically a restatement of [3, Proposition 3.5] proven for the single-vector RGS algorithm. Proposition 3.5. Consider Algorithm 2 with a backward-stable solver satisfying Assumptions 2.2 and

u crs ≤ 10 -3 cond( W) -1 m -2 , u f ine ≤ 10 -2 m -1/2 n -3/2 u crs .
Under Assumptions 3.1, if Θ is an ε-embedding for W, with ε ≤ 1/4, then it satisfies the ε -embedding property for Q, with ε = 2ε + 180u crs m 2 cond( W).

Proof. See appendix.

The ε-embedding property will likely hold even when matrix W depends on Θ. In such a case the quality of Θ can be certified a posteriori by computing additional sketches Φ Q and Φ W, associated with a new sketching matrix Φ, in addition to the sketches Θ Q and Θ W. Then one may characterize the quality of Θ by measuring the orthogonality of (Θ W)X and (Θ Q)Y, where X and Y are inverses of Q factors of Φ W and Φ Q, as is described in [3, Propositions 3.6-3.8] extrapolated from [START_REF] Balabanov | Randomized linear algebra for model reduction. Part II: minimal residual methods and dictionary-based approximation[END_REF]. For efficiency in terms of cache or communication, at each iteration, the products Φ W (i) and Φ Q (i) can be computed together with Θ W (i) in Step 2, and Θ Q (i) in Steps 4-5, respectively.

Randomized block Krylov methods

In this section we discuss practical applications of the methodology. We particularly focus on the improvement of efficiency of popular block Krylov methods, such as the block GMRES method and the Rayleigh-Ritz method, for solving block linear systems of the form

AX = B, (20)
and eigenvalue problems of the form

AX = XΛ, (21)
where A is large, possibly non-symmetric n × n matrix, B is n × m p matrix, and Λ is m p × m p diagonal matrix of the extreme eigenvalues of A.

The block Krylov methods proceed with approximation of X by a projection X (j) onto Krylov space K (j) (A, B), defined as

K (j) (A, B) := span{B, AB, • • • , A j-1 B},
with j being the order of the subspace.

The GMRES method computes X (j) that minimizes the Frobenius norm of the residual, while the Rayleigh-Ritz method seeks X (j) that is optimal in the Galerkin sense. Both methods first construct an orthonormal basis for K (p) (A, B) with GS orthogonalization, called Arnoldi iteration, and then determine the coordinates of the columns of X (j) in the computed basis.

Krylov basis computation: randomized block Arnoldi iteration

The Arnoldi algorithm produces orthonormal matrix Q = Q (1:p) satisfying the Arnoldi identity

AQ (1:p-1) = Q (1:p) H,
where H = H (1:p,1:p-1) is block upper Hessenberg matrix. The Arnoldi algorithm can be viewed as a block-wise QR factorization of matrix [B, AQ (1:p-1) ]. In this context, the Hessenberg matrix H can be seen as the R factor R = R (1:p,1:p) without the first column of subblocks, that is, as R (1:p,2:p) .

Below, we propose a randomized Arnoldi process based on RBGS algorithm. Note that unlike standard methods, Algorithm 7 produces a Krylov basis orthonormal with respect to the sketched product Θ•, Θ• .

Algorithm 7 RBGS-Arnoldi algorithm

Given: n × n matrix A, n × mp matrix B, k × n matrix Θ with k n, parameter p. Output: n × m factor Q = Q (1:p) and m × m upper triangular factor R = R (1:p,1:p) . 1. Set W (1) = B and perform 1-st iteration of RBGS (Algorithm 2). for i = 2 : p do 2. Compute W (i) = AQ (i-1) .

3. Perform i-th iteration of RBGS (Algorithm 2). end for 4. (Optional) compute ∆ (p) and ∆(p) . Use Theorem 3.2 to certify the output.

In Step 2 of Algorithm 7, the computation of the matrix-vector product can be executed either with roundoff u f ine or u crs depending on the situation. The stability guarantees of Algorithm 7 can be obtained directly from Theorems 3.2 and 3.4 with standard stability analysis similar to that from [3, Section 4.1]. In particular, it can be shown that, under the stability conditions of RBGS, the computed Q and H satisfy

(A + ∆A) Q (1:p-1) = Q (1:p) H, (22)
with ∆A close to machine precision, and cond( Q) = O(1). We leave the precise analysis of this fact outside of the scope of this manuscript.

Linear systems: randomized block GMRES method

Let us now discuss solution of block linear systems [START_REF] Vershynin | High-dimensional probability: An introduction with applications in data science[END_REF] with GMRES.

The GMRES method computes the approximate solution X (p-1) = U = [u 1 , . . . , u mp ] as

u i = Q (1:p-1) arg min z Hz -R (1:p,1) e i , (23)
performed in sufficient precision, where e i is the i-th column of the identity matrix. Let us now characterize the quasioptimality of such a projection when Q and H were obtained with the RBGS-Arnoldi algorithm. Under the stability conditions of RBGS, the Arnoldi identity [START_REF] Zou | GMRES algorithms over 35 years[END_REF] implies that

(A + ∆A) Q (1:p-1) z -b i = Q (1:p) ( Hz -R (1:p,1) e i ) ≤ Q Hz -R (1:p,1) e i ,
and, similarly,

(A + ∆A) Q (1:p-1) z -b i ≥ σ min ( Q) Hz -R (1:p,1) e i .
These two relations imply that

(A + ∆A)u i -b i ≤ cond( Q) min v∈Qp-1 (A + ∆A)v -b i , with Q p-1 = range( Q (1:p-1) ) = K p-1 (A + ∆A, B).
Consequently, the randomized version of the GMRES provides a solution which minimizes the norm of the residual associated with a slightly perturbed matrix, up to a factor of order 1.

Remark 4.1. Randomized Full Orthogonalization Method. In contrast to GMRES, the Full Orthogonalization Method (FOM) obtains solution X (p-1) = U = [u 1 , . . . , u mp ] by imposing a Galerkin orthogonality condition on the residuals (in exact arithmetic):

v, r(u i , b i ) = 0, ∀v ∈ Q p-1 , 1 ≤ i ≤ m p , (24)
where Q p-1 = range(Q (1:p-1) ) and r(u i , b i ) = Au i -b i is the residual associated with the i-th right-hand-side in [START_REF] Vershynin | High-dimensional probability: An introduction with applications in data science[END_REF].

The Galerkin projection may be a more appropriate choice than the minres projection, when the quality of the solution is measured with an energy error rather than the residual error. When Q (1:p-1) is obtained with the randomized Arnoldi method, the solution to (24) is given by

U = Q (1:p-1) (H (1:p-1,1:p-1) ) -1 R (1:p-1,1) . (25)
Suppose now that Q (1:p-1) was obtained by the randomized Arnoldi method. Then, by noticing that

H (1:p-1,1:p-1) = [I 0 (1:p-1,p) ]H = (ΘQ (1:p-1) ) T ΘQ (1:p-1) H = (ΘQ (1:p-1) ) T ΘAQ (1:p-1) ,
where 0 (1:p-1,p) is a null matrix of size of H T (p,1:p-1) , we deduce that the solution (25) satisfies

Θv, Θr(u i , b i ) = 0, ∀v ∈ Q p-1 , 1 ≤ i ≤ m p . ( 26
)
The relation (26) can be viewed as the sketched version of the Galerkin orthogonality condition. To our knowledge, it was first used in model order reduction community to obtain a reduced-basis solution of parametric linear systems [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF]. By using similar considerations as in [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF], one can show that (26) preserves the quality of the classical Galerkin projection when Θ is a ε-embedding for Q p with εcond(A) < 1, though, in practice, this condition can be too pessimistic [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF]. We leave the development of the randomized FOM method for future research.

Remark 4.2. It should be noted that the reasoning of this section could be reversed. We could first derive the sketched minres projection u i = arg min v∈Qp-1 Θr(v, b i ) and sketched Galerkin orthogonality condition (26) by replacing the

Eigenvalue problems: randomized Rayleigh-Ritz method

Next, we consider the computation of the extreme eigenpairs of A. Note that the methodology proposed below can be readily used for finding eigenpairs in the desired region by introducing a shift to the matrix. In addition, our methodology can be extended to inverse methods by replacing A with A -1 . In this case, the application of the inverse can be readily performed with (possibly preconditioned) randomized block GMRES method from Section 4.2 or other efficient methods.

To simplify the presentation, the analysis of the Rayleigh-Ritz method and its randomized version will be provided only in infinite precision arithmetic. The stability of these approaches follows directly from the stability of the Arnoldi iteration.

Let Q = Q (1:p) be the Krylov basis generated with the standard or randomized block-Arnoldi algorithm. The standard Rayleigh-Ritz method, based on the Arnoldi iteration, approximates the extreme eigenpairs (λ, x) of A by

(λ, x) ≈ (µ, u) = (µ, Q (1:p-1) y), (27)
where (µ, y) are the corresponding extreme eigenpairs of H (1:p-1,1:p-1) . Furthermore the residual error of an approximate eigenpair (µ, u) is estimated by H (p,1:p-1) y .

For the standard block-Arnoldi algorithm, we have

H (1:p-1,1:p-1) = [I 0 (1:p-1,p) ]H = (Q (1:p-1) ) T QH = (Q (1:p-1) ) T AQ (1:p-1)
,

and therefore v, r(u, µu) = 0, ∀v ∈ Q p-1 , (28)
where Q p-1 := range(Q (1:p-1) ) and r(u, µu) = Au -µu is the residual associated with (µ, u). The relation ( 28) is known as the Galerkin orthogonality condition. At the same time, we have

r(µ, u) ≥ σ min (Q) Hy -µ y 0 = σ min (Q) H (p,1:p-1) y , (29a)
and, similarly,

r(µ, u) ≤ σ max (Q) H (p,1:p-1) y . (29b)
Since, the classical Arnoldi iteration produces an 2 -orthogonal Q factor, hence in this case, the quantity H (p,1:p-1) y represents exactly the residual error of (µ, u).

On the other hand, the randomized block-Arnoldi algorithm produces a matrix Q orthogonal with respect to Θ•, Θ• . This implies that

H (1:p-1,1:p-1) = [I 0 (1:p-1,p) ]H = (ΘQ (1:p-1) ) T (ΘQ (1:p-1) )H = (ΘQ (1:p-1) ) T ΘAQ (1:p-1) ,
or equivalently, that (µ, u) satisfies the following sketched version of the Galerkin orthogonality condition:

Θv, Θr(u, µu) = 0, ∀v ∈ Q p-1 , (30)
similar to the sketched Galerkin condition for linear systems in Remark 4.1. Unlike for the sketched minres projection in GMRES, the optimality for the sketched Galerkin projection does not trivially follow from the ε-embedding property of Θ (or equivalently from the fact that cond(Q) = 1 + O(ε)). Consequently, for completeness of our paper, in Section 4.3.1 we establish a characterization of the accuracy of this projection. Note that, independently of this work, the sketched Galerkin orthogonality condition was also used by Nakatsukasa and Tropp in their recent paper [START_REF] Nakatsukasa | Fast & Accurate Randomized Algorithms for Linear Systems and Eigenvalue Problems[END_REF].

Since the ε-embedding property of Θ implies that cond(Q) = 1 + O(ε), hence according to (29), the quantity H (p,1:p-1) y estimates well the residual error.

The Rayleigh-Ritz algorithm based on RBGS-Arnoldi with restarting is depicted in Algorithm 8.

Remark 4.3. When the classical Galerkin projection is preferable to its sketched variant, say, because of its strong optimality properties for Hermitian operators, the RBGS-Arnoldi algorithm can be applied directly to the classical Galerkin methods. In this case, the 2 -orthogonal Krylov basis matrix Q and the Hessenberg matrix H in (27) can be readily obtained with the RBGS-Arnoldi algorithm followed by an additional Cholesky QR step, as is depicted in Remark 2.1. In particular, this situation can be accounted for in Algorithm 8 by adding the following line between Step 1 and Step 2: "Compute Q T Q and obtain its Cholesky factor R . Set

Q ← QR -1 , R ← R R."
The cost of this additional step is dominated by the computation of product Q T Q (since we can omit computing QR -1 ), which requires similar number of flops as the entire RBGS-Arnoldi algorithm. Nevertheless, as a single BLAS-3 operation, it is very well suited to modern computational architectures.

Algorithm

Brief analysis of the sketched Galerkin projection

We proceed with reformulation of the methodology in terms of projection operators similarly to [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF]Section 4.3] for classical methods. Our analysis will be based on the fact that Θ satisfies an ε-embedding property for some finite collection V of fixed low-dimensional subspaces V . This property can be satisfied with probability at least 1 -δ#V, if Θ is an (ε, δ, d) oblivious 2 -subspace embedding, with d = max V ∈V dim(V ).

Let K be an approximation space, which does not necessarily have to be a Krylov space. Let Π K denote the 2orthogonal projector onto K:

∀w ∈ C n , Π K w = arg min v∈K w -v . ( 31
)
Then the Galerkin orthogonality condition (28) (with K = Q p-1 ) can be expressed as

Π K r(u, µu) = 0, or, equivalently, Π K AΠ K u = µu.
In other words, the classical Rayleigh-Ritz method can be interpreted as approximation of eigenpairs of A by the eigenpairs of approximate operator Π K AΠ K . Similarly, the sketched Galerkin orthogonality condition (30) (with K = Q p-1 ) can be expressed as

Π Θ K r(u, µu) = 0, (32)
or, equivalently,

Π Θ K AΠ Θ K u = µu, where Π Θ
K is an orthogonal projector onto K with respect to the sketched inner product Θ•, Θ• , i.e.,

∀w ∈ C n , Π Θ K w = arg min v∈K Θ(w -v) , (33)
Proof. See appendix.

Corollary 4.7. Assume that Θ is an ε-embedding for K. If K is invariant under A then each pair (µ, u) that satisfies the sketched Galerkin orthogonality condition (32) is an eigenpair of A.

Theorem 4.6 implies two important results. The first one is that the sketched Galerkin projection provides an accurate result in terms of the residual error when the approximation space K approximates well the range of A or when K is close to an invariant space. In particular, it follows that if we have (I -Π K )AΠ K ≤ τ |µ m |, then the dominant m computed eigenpairs (µ, u) have relative error less than τ :

(A -µI)u µu ≤ (1 + O(ε))τ. (35)
At the same time, Theorem 4.6 also proves that if K represents a Krylov subspace Q p-1 , then the convergence of eigenvectors implies convergence of their residuals to zero. More rigorously, if

u (p) -u (p-1) ≤ τ |µ m | u (p) ,
where u (p) denotes the sketched Galerkin solution in Q p-1 , then (35) holds.

Unfortunately, the fact that the residual norm is small does not guarantee that the eigenpair is accurate due to possibility of bad conditioning of the eigenvalue. If A is Hermitian, this issue can be circumvented for classical Galerkin approximation, by appealing to the Min-Max principle and the Courant characterization [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF]. Perhaps, similar ideas can be also used for the sketched Galerkin projection. The sketched Rayleigh quotient R(w) can be naturally defined as a number that minimizes Θ(Aw -R(w)w) . It is easy to see that R(w) = Θw,ΘAw Θw,Θw . Furthermore, it can be shown that the solution (µ, u) to (32) and the exact solution (λ, x) satisfy µ = R(u) and λ = R(x). We can also show that |R(Π K x) -λ| is small, if I -Π K x is small. The main difficulty now becomes how to use the Min-Max principle, since the sketched operator Π Θ K AΠ Θ K can be non-Hermitian, even when A is Hermitian.

Further applications

s-step Krylov methods

The proposed RBGS algorithm can be used to improve the s-step Krylov methods. These methods use efficient matrix power kernels that can output a matrix F (v) of basis vectors for K (s+1) (A, v) , of the form

F (v) = [p 0 (A)v, p 1 (A)v, . . . , p s (A)v],
where v is some input vector, p 0 (A) is usually taken as 1, and p 1 (A), . . . p s (A) are some suitable polynomials that aim to make F (v) not too badly conditioned. Popular options for p 0 (A), . . . p s (A) are the monomials and Newton, or Chebyshev polynomials. A relatively large s (say s = 10 or 30 [START_REF] Hoemmen | Communication-avoiding Krylov subspace methods[END_REF]) may lead to stability problems even when a polynomial basis is used. Consequently, to obtain a Krylov basis suitable for practical applications, it becomes necessary to proceed with a block-wise generation, at each iteration, computing a small block of a Krylov basis by using the matrix power kernel and subsequently orthogonalizing it against the previously computed vectors with a BGS approach. In more concrete terms, we have to perform a block-wise QR factorization of the matrix W = W (1:p) generated recursively as

W (i) = b if i = 1 F * (q s(i-2)+1 ) if i = 2, . . . p
where F * (v) corresponds to the matrix F (v) without the first column, and q s(i-2)+1 is the last Krylov basis vector computed at iteration i -1. Such block-wise orthogonalization can readily be done with the RBGS algorithm, as is depicted in Algorithm 9, with all the computational benefits of this algorithm.

An important point is that the random sketching allows not only to effectively construct a well-conditioned Krylov basis, but also to determine the Hessenberg matrix H, without any additional high-dimensional operations nor the need to appeal to cumbersome recurrence relations. With the output of Algorithm 9, an approximate solution to the linear system can be obtained in a classical way presented in Section 4.2 (by setting B = b). We leave the numerical analysis of the proposed s-step RBGS-Arnoldi algorithm outside the scope of this manuscript.

Algorithm 9 s-step RBGS-Arnoldi algorithm

Given: n × n matrix A, n × 1 vector b, k × n matrix Θ with k n, parameters p and s. Output: n × m Krylov basis matrix Q = Q (1:p) and Hessenberg matrix H. 1. Set W (1) = b and perform 1-st iteration of RBGS (Algorithm 2). for i = 2 : p do 2. Compute W (i) = F * (q s(i-2)+1 ) and perform i-th iteration of RBGS (Algorithm 2). end for 3. Obtain H = S † P (1:p-1) (possibly with setting the below-subdiagonal elements to zeros) by QR or SVD of S, or solving a small least-squares problem. 4. (Optional) compute ∆ (p) and ∆(p) . Use Theorem 3.2 to certify the output.

Numerical experiments

We test the methodology on three numerical examples: QR factorization of a synthetically generated matrix from [3, Section 5.1], solution of a linear system with block GMRES, and solution of an eigenvalue problem with Rayleigh-Ritz method. The RBGS process is validated by comparison with standard methods such as BCGS, BMGS and BCGS2. Depending on the example, the finite precision arithmetic is performed in float32 or float64 format. To ensure good comparison, we decided to perform the inter-block orthogonalization routines in float64, even in the unique float32 precision algorithms. Furthermore, for standard methods we perform the inter-block orthogonalization with an efficient explicit Householder QR routine while for the RBGS, this routine is combined with the sketched Cholesky QR according to Algorithm 6. Several solvers are considered for the sketched least-squares problems depending on the experiment, but in all of them this step has a negligible complexity and memory requirements. Finally, in all the numerical examples the SRHT and Rademacher matrices give very similar results, so we present the results for SRHT only.

Orthogonalization of a numerically singular matrix

Take [W] i,j = f µj (x i ), 1 ≤ i ≤ n, 1 ≤ j ≤ m, where f µ (x) =
sin (10(µ + x)) cos (100(µ -x)) + 1.1 , and x j and µ j are chosen as, respectively, n = 10 6 and m = 300 equally distanced points with x 0 = µ 0 = 0 and x 10 6 = µ 300 = 1. The matrix W is partitioned into blocks of columns of size m p = 10. Then a block-wise QR factorization of W is performed with the RBGS process (Algorithm 2) executed either in unique float32 precision or in multi-precision, using Θ with k = 3000 rows. The least-squares solver in Step 2 of Algorithm 2 is chosen either as the Householder solver, or as 20 iterations of CG applied to the normal equation. Thereafter the stability behavior of RBGS is compared to the standard BCGS, BMGS, and BCGS2 algorithms executed in float32 arithmetic. We observe a similar picture as in [3, Section 5.1] comparing single-vector Gram-Schmidt algorithms. According to Figure 1a, the BCGS and BCGS2 methods become dramatically unstable respectively at iterations i ≥ 8, and i ≥ 17, with the latter iterations corresponding to W (1:i) being numerically singular. The BMGS algorithm exhibits a better stability than the other two standard approaches. However, it still yields a Q factor with a condition number of two orders of magnitude. The unique precision RBGS algorithm using Householder least-squares solver in Step 2, has similar stability profile as BMGS algorithm. On the other hand, the usage of 20 iterations of CG does not provide a sufficient accuracy in Step 2 and as a consequence, reduces stability of the unique precision RBGS. In contrast to all tested methods, the multi-precision RBGS algorithm remains perfectly stable, even at iterations where W (1:i) is numerically singular, and outputs a Q factor with the condition number close to O(1). Unlike in the unique precision algorithm, here the CG solver in Step 2 turns out to be sufficiently accurate. Moreover, Figure 1b examines the behavior of the approximation error W (1:i) -Q (1:i) R (1:i,1:i) / W (1:i) . We see that for all tested algorithms, besides BCGS2, this error remains close to machine precision at all iterations. 5.2 Solution of a linear system with block GMRES

Consider linear system, (A Ga + αI)X Ga = B, (36)
where the matrix A Ga is taken as the "Ga41As41H72" matrix of dimension n = 268096 from the SuiteSparse matrix collection, I is the n × n identity matrix, and α = 0.2 introduced to improve conditioning of A Ga . The right-hand-side matrix B is taken as n × m p random matrix with entries being i.i.d normal random variables and m p = 100. Solving such system could be a part of the inverse subspace iteration for computing eigenvalues of A Ga of the smallest magnitude. The shifted "Ga41As41H72" matrix has many clustered, possibly negative, eigenvalues, which can bring stability issues to the solvers. Furthermore, the system (36) is preconditioned from the right by the incomplete LU factorization P Ga of A Ga + αI with zero level of fill-in and symmetric reverse Cuthill-McKee reordering. With this preconditioner the final system of equations has the following form AX = B, where A = (A Ga +αI)P Ga and X Ga = P Ga X. Here matrix A is not computed explicitly but considered as an impicit map outputting product with vectors. This system is approximately solved with GMRES method based on different versions of BGS process. We restart the GMRES every 30 iterations, i.e., when the dimension of the Krylov space becomes m = 3100.

Here we examine the behavior of BCGS, BCGS2, BMGS and the unique precision RBGS under float32 arithmetic. There is no need to test the multi-precision RBGS as the unique precision algorithms already are able to provide nearly optimal solution. The products with matrix A and solutions of reduced GMRES least-squares problems (23) are computed in float64 format. The BGS iterations and other operations are performed in float32 format. In RBGS, the sketched leastsquares problems are solved with 5 Richardson iterations as is explained in Section 2.1.

Figure 3a provides convergence of the residual error max j=1,...,100 Au j -b j / b j . The condition number of the computed Krylov basis Q (1:i-1) at each iteration i is given in Figure 3b. We see that already starting from the first iterations, the BCGS algorithm entails dramatic instability and stagnation of the residual error. At the same time, the BCGS2 algorithm remains stable at all iterations providing a Q factor that is orthonormal up to machine precision. The BMGS algorithm has partial instability resulting in a deteriorated convergence of the error. Finally, the RBGS algorithm is as stable as the BCGS2 algorithm. It provides a well-conditioned Q factor and an optimal convergence of the residual error for all tested sizes of the sketching matrix.

Solution of an eigenvalue problem with randomized Rayleigh-Ritz

In this test we seek the smallest (negative) eigenpairs of the "Ga41As41H72" matrix A Ga from Section 5.2. For this we first transform the eigenproblem to the following positive-definite one:

(αI -A Ga )x = λx, (37)
where α = 1500 ≈ A Ga , and seek the dominant eigenpairs of A = αI -A Ga . Then one can use the fact that A Ga and A have same eigenvectors and the associated eigenvalues related as λ Ga = α -λ. In this case the first 50 eigenvalues of A are clustered inside [1500.9 1501.3]. To compute the dominant eigenpairs of A we either apply a subspace iteration or the Rayleigh-Ritz method based on different BGS algorithms. The initial guess matrix B is taken as an n × m p Gaussian matrix of dimension m p = 50. The Arnoldi algorithm is restarted every 50 iterations (with columns of B chosen as the dominant Ritz vectors) which corresponds to the Krylov space of dimension m = 2550. In this numerical example all arithmetic operations are performed in float64. The Step 2 of RBGS algorithm is again performed with 5 Richardson iterations as in the previous example.

We measure the approximation error by the maximum relative residual of the first 80% of computed eigenpairs (i.e., 40 eigenpairs out of 50). The convergence of the approximation error max j=1,...,40 Au j -µ j u j / µ j u j is depicted in Figure 3a. The condition number of the computed Krylov basis is depicted in Figure 3b. It is revealed that the subspace iteration method yields early stagnation of the error and is unsuitable for this numerical example. On the other hand, the Rayleigh-Ritz method, if stable, provides convergence of the error to machine precision.

It is revealed that the BCGS-based Arnoldi method becomes unstable already at early iterations and produces an ill-conditioned Krylov basis with condition number close to O(u -1 ). Therefore this method needs to be restarted, say every 5 iterations. Figure 3a plots the error also for this case. It is seen that the early restarting does not help as it causes a dramatic effect on the convergence of the error. In contrast to BCGS, the BCGS2 and RBGS algorithms show perfect stability and yield an approximation that converges to machine precision in 200 iterations. The BMGS method is unstable in the sense that it outputs an ill-conditioned Krylov basis. Despite this, it yields an approximation error that converges to machine precision, though, with somewhat deteriorated rates compared to BCGS2 and RBGS.

Concluding remarks

In this work we developed a block version of the randomized Gram-Schmidt process, called RBGS, to compute QR factorization of a large-scale matrix. It was shown that this algorithm inherits main properties of its single-vector analogue from [START_REF] Balabanov | Randomized Gram-Schmidt process with application to GMRES[END_REF], and in particular that it is at least as stable as the single-vector MGS process, and requires nearly half the cost of the classical BGS process in terms of flops and data passes. At the same time, our RBGS algorithm is well suited to cache-based and highly parallel computational architectures because it mainly relies on the matrix-matrix BLAS-3 operations, in contrast to its single vector counterpart relying on matrix-vector BLAS-2 operations. Like the single-vector RGS algorithm, it can be implemented using multi-precision arithmetic allowing to perform the dominant large-scale operations in precision independent of the dimension of the problem. This can become especially handy for simulations on low-precision arithmetic architectures. Different strategies for treating the inter-block orthogonalization step, as well as solution of the sketched least-squares problems have been proposed. Special care has been taken to ensure that the computational cost of these steps is negligible compared to other steps in the RBGS algorithm.

The stability of our algorithms was verified in numerical experiments. In particular, it was seen that the multiprecision RBGS can provide a stable QR factorization of a numerically singular matrix, where even the most stable standard algorithms, including BCGS2, do not work. Furthermore, the randomized Arnoldi algorithm based on RBGS showed excellent stability in the context of GMRES and Rayleigh-Ritz approximation in numerical examples where BMGS showed instability and convergence degradation. These factors indicate the robustness of the RBGS algorithm.

Our next goal is to combine RBGS with model order reduction and compression techniques for even more, possibly asymptotic, reduction of the cost of Gram-Schmidt orthogonalization and the associated Krylov methods. Another direction is the application of random sketching to increase the robustness and/or efficiency of Krylov methods that use short recurrences such as CG, BCG, Lanczos, and others. Furthermore, we want to improve not only the Krylov methods but also other methods that involve an orthogonalization of the approximation basis, such as block LOBPCG. As for the theoretical analysis, besides improving the characterization of the sketched Galerkin projection in the randomized Rayleigh-Ritz method, we also plan to investigate the ε-embedding property of the sketching matrix for the computed Krylov space in the presence of rounding errors. Despite the fact that this property has been thoroughly verified in numerical experiments, it still remains unproven for the randomized Arnoldi algorithm.

Appendix

Here we provide proofs of some propositions and theorems from the paper.

Proof of Theorem 3.2. We have

Θ Q -S F ≤ 1.02u f ine n Θ F Q F ≤ 1.02 √ 1 + εu f ine n 3/2 Q F ≤ 0.02u crs Q := F 1 , (38)
which implies that σ min ( S) -

F 1 ≤ σ min (Θ Q) ≤ σ max (Θ Q) ≤ σ max ( S) + F 1 . By noticing that 1 -∆ (p) ≤ 1 -∆ (p) ≤ σ min ( S) ≤ σ max ( S) ≤ 1 + ∆ (p) ≤ 1 + ∆ (p) ,
and Corollary 1.3 we deduce that Q ≤ 1.57 and that

(1 + ε) -1/2 (1 -∆ (p) -F 1 ) ≤ σ min ( Q) ≤ σ max ( Q) ≤ (1 -ε) -1/2 (1 + ∆ (p) + F 1 ).
Next is proven the second statement of the theorem, which is obvious for p = 1. Assume that p ≥ 2 and notice that R F ≤ σ min ( S) -1 S R F ≤ 1.12 S R F ≤ 1.12( P F + P -S R F ) ≤ 1.12(1 + ∆(p) )

P F ≤ 1.4(1 + ∆(p) ) W F ≤ 1.6 W F , (39)
We also have, for 1 ≤ i ≤ p,

W (i) = Q (1:i) R (1:i,i) + ∆W (i) ,
with

∆W (i) F ≤ Q (i) -Q (i) R (i,i) F + 1.02ucrs( W (i) F + imp Q (1:i-1) F R (1:i-1,i) F)
≤ 0.1ucrs Q (i) F R (i,i) F + 1.02ucrs( W (i) F + 1.57i Proof of Remark 3.3. We have,

W (i) = Q (1:i) R (1:i,i) + ∆W (i) ,
with

Θ∆W (i) F ≤ Θ( Q (i) -Q (i) R (i,i) ) F + 1.02 √ 1 + εucrs( W (i) F + imp Q (1:i-1) F R (1:i-1,i) F)
≤ 1.02ucrs(1.25 W (i) F + 0.16m Proof of Theorem 3.4. The proof is done by induction on p. Assume that Theorem 3.4 holds for p = i -1 ≥ 1. Below, we show that the statement of the theorem then also holds for p = i. Clearly, S (1:i-1) F ≤ (i -1)m p 1 + ∆ (i-1) ≤ 1.01 (i -1)m p and P (i) ≤ 1.02 Θ W (i) ≤ 1.25 W (i) .

(40)

Moreover, we have for 1 ≤ j ≤ i, Θ Q (j) -S (j) F ≤ 1.02u f ine √ 1 + εn 3/2 Q (j) F ≤ 0.02u crs Q (j)

From the fact that σ min ( S (1:i-1) ) ≥ √ 1 -∆ (i-1) ≥ 0.989, and Assumptions 2.2 and (40), notice that σ min ( S (1:i-1) + ∆S (i-1) ) ≥ σ min ( S (1:i-1) ) -∆S (i-1) ≥ 0.98. (41)

Consequently, we have R (1:i-1,i) F ≤ P (i) + ∆P (i) F /σ min ( S (1:i-1) + ∆S (i-1) ) ≤ 1.4

W (i) F (42)
This fact combined with the fact (which holds by the induction hypothesis: ∆ (i-1) ≤ 0.02 and Theorem 3.2) that Q (1:i-1) ≤ 1.5 and Q (1:i-1) F ≤ 1.5(i -1)

1/2 m 1/2 p , (43)
and ( 14), leads to the following result:

∆Q (i) F = Q (i) -( W (i) -Q (1:i-1) R (1:i-1,i) ) F ≤ 1.02ucrs | W (i) | + imp| Q (1:i-1) || R (1:i-1,i) | F ≤ 1.02ucrs( W (i) F + imp Q (1:i-1) F R (1:i-1,i) F) ≤ 3.2ucrsi 3/2 m 3/2 p W (i) F, (44)
and

Q (i) F ≤ W (i) F + Q (1:i-1) R (1:i-1,i) F + ∆Q (i) F ≤ 3.2 W (i) F. ( 45
)
Then by [START_REF] Martinsson | Randomized numerical linear algebra: Foundations and algorithms[END_REF],

Θ∆Q (i) F ≤ 3.2 √ 1 + εucrsi 3/2 m 3/2 p W (i) F ≤ 4ucrsi 3/2 m 3/2 p W (i) F,
and

Θ Q (i) F ≤ Θ W (i) F + Θ Q (1:i-1) R (1:i-1,i) F + Θ∆Q (i) F ≤ 3.2 W (i) F. ( 46
)
Denote residual matrix Q (1:i-1) R (1:i-1,1:i-1) -W (1:i-1) by B (i-1) . Then, we have B (i-1) ≤ 4ucrsi 3/2 m 3/2 p W (1:i-1) F ≤ 0.01σmin( W), and Θ B (i-1) ≤ 5ucrsi 3/2 m 3/2 p W (1:i-1) F ≤ 0.01σmin( W).

(47)

Now we are all set to derive the ε -embedding property of Θ for Q (i) , which is needed to characterize the quality of the inter-block orthogonalization step. Let us first notice that,

σmin( Q (i) ) = σmin( W (i) -Q (1:i-1) R (1:i-1,i) + ∆ Q (i) ) ≥ σmin( W (i) -Q (1:i-1) R (1:i-1,i) ) -∆ Q (i)
≥ σmin( W (i) -( W (1:i-1) + B (i-1) ) R -1 (1:i-1,1:i-1) R (1:i-1,i) ) -∆ Q (i)

= σmin [ W (1:i-1) + B (i-1) , W (i) ] -R -1 (1:i-1,1:i-1) R (1:i-1,i)

I -∆ Q (i)
≥ σmin([ W (1:i-1) + B (i-1) , W (i) ]) -∆ Q (i)

≥ σmin( W) -B (i-1) -∆ Q (i) ≥ 0.98σmin( W).

(48)

Consequently, cond( Q (i) ) ≤ 5m

1/2 p cond( W). For any a ∈ R mp , it holds Θ Q (i) a = Θ( W (i) -Q (1:i-1) R (1:i-1,i) )a ± Θ∆ Q (i) a = Θ( W (i) -W (1:i-1) R -1 (1:i-1,1:i-1) R (1:i-1,i) )a ± Θ B (i-1) R -1 (1:i-1,1:i-1) R (1:i-1,i) a ± Θ∆ Q (i) a = 1 ± σmin(Θ W) -1 ( Θ B (i-1) + Θ∆ Q (i) ) Θ( W (i) -W (1:i-1) R -1 (1:i-1,1:i-1) R (1:i-1,i) )a = (1 ± 0.015) Θ( W (i) -W (1:i-1) R -1 (1:i-1,1:i-1) R (1:i-1,i) )a = (1 ± 0.015) √ 1 ± ε ( W (i) -W (1:i-1) R -1 (1:i-1,1:i-1) R (1:i-1,i) )a = (1 ± 0.015)

√ 1 ± ε ( W (i) -Q (1:i-1) R (1:i-1,i) )a ± B (i-1) R -1 (1:i-1,1:i-1) R (1:i-1,i) a = (1 ± 0.015)(1 ± 0.01) √ 1 ± ε ( W (i) -Q (1:i-1) R (1:i-1,i) )a = (1 ± 0.015)(1 ± 0.01) √ 1 ± ε( Q (i) a + ∆ Q (i) a ) = (1 ± 0.015)(1 ± 0.01)(1 ± 0.01) √ 1 ± ε Q (i) a = (1 ± 0.036) √ 1 ± ε Q (i) a ,
which, combined with the parallelogram identity, implies that Θ is an ε -embedding for Q (i) with ε ≤ 1.1ε + 0.1 ≤ 0.65. By combining this fact with Assumptions 3.1, we obtain:

1 -u crs m 1/2 p cond( W) ≤ σ min (Θ Q (i) ) ≤ σ max (Θ Q (i) ) ≤ 1 + u crs m 1/2 p cond( W),
and We also have R (i,i) F ≤ 1.02 Θ Q (i) F ≤ 3.3 W (i) F , and σ min ( R (i,i) ) ≥ σ min (Θ Q (i) R (i,i) )/ Θ Q (i) ≥ 0.99(σ min (Θ Q (i) ) -0.1u crs Q (i) R (i,i) ) ≥ 0.58σ min ( W), and S T (1:i-1) S (i) F = S T (1:i-1) (Θ Q (i) + ∆ S (i) ) R -1 (i,i) + ∆ S (i) F ≤ ( S T (1:i-1) Θ Q (i) F + 1.01 ∆ S (i) ) F )/0.58σ min ( W) + ∆ S (i) F , where ∆ S (i) = Θ( Q (i) R (i,i) -Q (i) ) and ∆ S (i) = S (i) -Θ Q (i) . By noting that, ∆ S (i) F ≤ 0.1ucrs Q (i) R (i,i) ≤ 0.6ucrs W (i) F ∆ S (i) F ≤ 0.02ucrs Q (i) ≤ 0.04ucrs

Q (i) ≤ 1.8. Since Θ Q (i) -S (i) F ≤ 0.02u crs Q (i) ,

end for 7 .

 7 (Optional) compute ∆ (p) = I -S T (1:p) S (1:p) F and ∆(p) = P (1:p) -S (1:p) R (1:p.1:p) F P (1:p) F . # macheps: u f ine

  oblivious subspace embedding, with d = O(log(m/δ)). In its turn this condition is met by Rademacher matrices with k = O(log(n) log(m/δ)) rows or SRHT with k = O(log 2 (n) log 2 (m/δ)) rows.

Figure 1 :

 1 Figure 1: Block QR factorization of the matrix from [3, Section 5.1]. In the plots, the u.p. CG-RBGS and m.p. CG-RBGS respectively refer to the unique precision RBGS and the multi-precision RBGS algorithms, performing Step 2 with 20 iterations of CG. On the other hand, the u.p. HH-RBGS refers to the unique precision RBGS utilizing Householder solver.

20

 20

Figure 2 :

 2 Figure 2: Solution of a linear system with GMRES.

Figure 3 :

 3 Figure3: Solution of an eigenvalue problem with Rayleigh-Ritz method. In the plots, BCGS-restarted refers to BCGS-Arnoldi algorithm that restarts every 5 iterations.

we have 1 - 1 . 2 pI| x 2 - 2 ≤

 11222 05u crs m 1/2 p cond( W) ≤ σ min ( S (i) ) ≤ σ max ( S (i) ) ≤ 1 + 1.05u crs m 1/2 p cond( W).Consequently,I -S T (i) S (i) F ≤ m 1/-S T (i) S (i) ≤ m 1/2 p max x S (i) x 2 |x 3u crs m p cond( W).

2 pW

 2 i-1) Θ Q (i) F ≤ S T (1:i-1) Θ( W (i) -Q (1:i-1) R (1:i-1,i) ) F + 4.04ucrsi 3/2 m 3/(i) F ≤ S T (1:i-1) ( P (i) -S (1:i-1) R (1:i-1,i) ) F + 4.2ucrsi 3/2 m 3/2 p W (i) F ≤ 5ucrsi 3/2 m 3/2 p W (i) F,

  8 Randomized Rayleigh-Ritz algorithm for extreme eigenpairs with restarting Given: n × n matrix A, n × mp matrix B, k × n matrix Θ with k n, parameters p and Niter. Output: Λ and X.

for i = 1 : Niter do 1. Perform RBGS-Arnoldi algorithm (Algorithm 7) returning Q and R. Set H = R (1:p,2:p) . 2. Compute diagonal matrix Λ with mp extreme eigenvalues of H (1:p-1,1:p-1) on the diagonal, and the matrix of associated eigenvectors Y. 3. (Optional) Compute the residuals H (p,1:p-1) y to characterize the approximation error. 4. Compute B = Q (1:p-1) Y. end for 5. Normalize B with respect to 2-norm and return it as X.

  3/2 m 3/2 ≤ 1.02 2 u 2 crs 2(1.26 2 m W 2 F + 1.57 2 m 3 R 2 F ) ≤ 1.02 2 u 2 crs 2(1.26 2 m W 2 F + 2.52 2 m 3 W 2 F )

				p	R (1:i-1,i) F)
	≤ 1.02ucrs( W (i) F + 0.16m 1/2 p	R (i,i) F + 1.57i 3/2 m 3/2 p	R (1:i-1,i) F)
	≤ 1.02ucrs(1.26m 1/2 p	W F + 1.57i 3/2 m 3/2
	Consequently,
	∆W 2 F ≤ 1.02 2 u 2 crs	2(1.26 2 mp W 2 F + 1.57 2 i 3 m 3 p R (1:i-1,i)	2 F )
	1≤i≤p
	≤ 4ucrsm 3/2 W F	2	.

p R (1:i-1,i) F).

  Θ∆W 2 F ≤ 1.02 2 u 2 crs 2(1.51 2 m W 2 F + 2 2 m 3 R 2 F ) ≤ 1.02 2 u 2 crs 2(1.51 2 m W 2 F + 3.2 2 m 3 W 2 F ) ≤ 5ucrsm 3/2 W F

		1/2 p	R (i,i) F + 2i 3/2 m 3/2 p	R (1:i-1,i) F)
	≤ 1.02ucrs(1.51m 1/2 p	W F + 2i 3/2 m 3/2 p	R (1:i-1,i) F)
	Consequently,

2

.

-inner products are replaced by Θ•, Θ• . In particular, for the case l = 1 the two algorithms are essentially equivalent. The numerical analysis of Algorithm 4 is beyond the scope of this manuscript.

-inner products and 2 -norms in the standard minres equation u i = arg min v∈Qp-1 r(v, b i ) and the standard Galerkin orthogonality condition (24) by sketched ones, similarly as was done in[START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF][START_REF] Balabanov | Randomized linear algebra for model reduction. Part II: minimal residual methods and dictionary-based approximation[END_REF]. After that, it could be realized that the solution to the sketched minres and Galerkin equations can be obtained by orthogonalizing the Krylov basis with respect to Θ•, Θ• and using the classical identities (23) and (25).

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 810367).

or, in matrix form,

where Q is a matrix whose columns compose a basis for K. We see that the randomized Rayleigh-Ritz method corresponds to taking the approximate operator as Π Θ K AΠ Θ K instead of Π K AΠ K . Next we shall provide an upper bound for the residual norm of the exact eigenpair with respect to such approximate operator. For this, it is necessary to first establish a characterization of the sketched orthogonal projector, given in Lemma 4.4. lemma 4.4. Let w ∈ C n . If Θ is ε-embedding for V = K + span(w), then

Moreover, if w ∈ K, then

Proof. See appendix.

Theorem 4.5. (extension of [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF]Theorem 4.3]) Assume that Θ is an ε-embedding for V = K + AK + span(x). Let

Then the sketched residual norms of the pairs (λ, Π Θ K x) and (λ, x) for the linear operator

Proof. See appendix.

Theorem 4.5 together with Lemma 4.4 imply that

Furthermore, according to Lemma 4.4, we also have

Thus, Theorem 4.5 guarantees a good approximation if the distance (I -Π K )x is small and if the approximate eigenproblem is well conditioned. Interestingly, by taking Θ as an identity in Theorem 4.5, we exactly recover Theorem 4.3 in [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF].

On the contrary, the following result bounds the residual norm of the approximate eigenpair (µ, u) in (32) with respect to the exact operator. Theorem 4.6. Let (µ, u) be a solution to (32). If Θ is ε-embedding for V = K + AK, then we have

and,

Consequently,

The proof is finished by noting that

Proof of Proposition 3.5. We proceed with an induction on p. Clearly, the statement of the proposition holds for p = 1. Assume that Θ is an ε-embedding for Q (1:p) , p = i -1 ≥ 1. This condition is sufficient for the following results in Theorem 3.4 and its proof to hold for p = i:

In addition, we have

3, and

which can be proven similarly to (39) and (44) in the proofs of Theorems 3.2 and 3.4. Finally, the following holds

By using the fact that S (1:i) ≤ 1.02 and σ min ( P (1:i) ) ≥ σ min (Θ W (1:i) ) -0.02u crs W (1:i) , along with (49) and the ε-embedding property of Θ, we get

Properties (49), ( 52) and (53) imply that

Due to the fact that Θ is an ε-embedding for W (1:i) , we deduce that

We also have from ( 51) and (53),

By combining this property with (54), we get

We conclude that

where

, which in particular implies that Q (1:i) ≤ 1.6 and σ min ( Q (1:i) ) ≥ 0.86. Furthermore, from (52), we get

where 55) and (56) it is deduced that for any vector a ∈ R m ,

where ε = 2ε + 180u crs m 2 p i 2 cond( W (1:i) ). By using the parallelogram identity this relation can be brought to form (2) where V = range( Q (1:i) ) and ε = ε . It is deduced that Θ is a ε -embedding for Q (1:i) .

Proof of Lemma 4.4. By definitions of Θ, Π K and Π Θ K , we have

which gives the first inequality. Furthermore, we also have the following relation

we obtain (34b). Finally, if w ∈ K, then, by definition, Π Θ K w = arg min v∈K Θ(wv) = w, which finishes the proof.

Proof of Theorem 4.5. The proof directly follows that of [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF]Theorem 4.3] replacing 2 -inner products and norms by the sketched ones. We have,

Furthermore, by using the fact that the two vectors on the right hand side are orthogonal with respect to the sketched inner product, we obtain

Proof of Theorem 4.6. We have, Θr(u, µu) 2 = ΘΠ Θ K r(u, µu) 2 + Θ(I -Π Θ K )r(u, µu) 2 = Θ(I -Π Θ K )r(u, µu) 2 .

Consequently, by an ε-embedding property of Θ, we obtain

which gives theorem's first inequality. Moreover, notice that (I -Π Qp-1 )AΠ Qp-1 u = (I -Π Qp-1 )A(I -Π Qp-2 )u, which gives theorem's second inequality.