{'fraction_non_alphanumeric': 0.045888257680216205, 'fraction_numerical': 0.01392642435607523, 'mean_word_length': 4.113947257086448, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Randomized subspace iteration: Analysis of canonical angles and unitarily invariant norms

This paper is concerned with the analysis of the randomized subspace iteration for the computation of low-rank approximations. We present three different kinds of bounds. First, we derive both bounds for the canonical angles between the exact and the approximate singular subspaces. Second, we derive bounds for the low-rank approximation in any unitarily invariant norm (including the Schatten-p norm). This generalizes the bounds for Spectral and Frobenius norms found in the literature. Third, we present bounds for the accuracy of the singular values. The bounds are structural in that they are applicable to any starting guess, be it random or deterministic, that satisfies some minimal assumptions. Specialized bounds are provided when a Gaussian random matrix is used as the starting guess. Numerical experiments demonstrate the effectiveness of the proposed bounds.

1. Introduction. The computation of low-rank approximations of large-scale matrices is a vital step in many applications in data analysis and scientific computing, such as principal component analysis, facial recognition, spectral clustering, model reduction techniques such as proper orthogonal decomposition (POD) and discrete empirical interpolation method (DEIM), approximation algorithms for PDEs and integral equations. The celebrated Eckart-Young theorem [10] says that the optimal low-rank approximation can be obtained by means of the Singular Value Decomposition (SVD); however, computing the full or truncated SVD can be computationally challenging, or even prohibitively expensive for many applications of interest.
Randomized algorithms for computing low-rank approximations have become increasingly popular in the last two decades, see the survey papers [12,17]. Randomized methods have gained in popularity since they are easy to implement, computationally efficient, and numerically robust. While randomized algorithms tend to have the same asymptotic cost comparable to classical methods, randomized algorithms are able to exploit parallel computing efficiently. Additionally, for large datasets that are too large to fit in memory, the randomized algorithms are efficient in the number of times they access the data, compared to classical algorithms. We focus on a particular method known as randomized subspace iteration. The main idea of this method is to identify a subspace of the matrix by randomly mixing its columns, i.e., by multiplying the matrix by a random matrix of appropriate dimensions. A low-rank approximation is obtained by projecting the matrix onto this subspace. A post-processing step here involves compressing the matrix by a deterministic method, and a conversion step to obtain an equivalent representation in the desired format (typically, a truncated SVD representation). In practice, to obtain a higher quality approximation, the mixing step requires repeated applications of the action of a matrix.
Randomized algorithms also have excellent numerical robustness and the resulting approximation come with guarantees of very high probability. Many advances have been made in the analysis of randomized algorithms for low-rank approximations. The analysis typically has two stages: the structural, or deterministic stage, in which minimal assumption about the distribution of the random matrix is made, and probabilistic stage, in which the distribution of the random matrix is taken into account to derive bounds for expected and tail bounds for the error distribution. As mentioned earlier, existing literature only targets the error in the low-rank representation [11,12]. When the low-rank representation is in the SVD format, it is desirable to understand the quality of the approximate subspaces and the individual singular triplets. This paper aims to fill in some missing gaps in the literature by a rigorous analysis of the accuracy of approximate singular values, vectors and subspaces obtained using randomized subspace iteration. This analysis will be beneficial in applications where an analysis beyond the low-rank approximation is desired. Examples include Model Reduction techniques [9,2], Leverage Score computation [14], Spectral Clustering [6], FEAST eigensolvers [22], Canonical Correlation Analysis [1].
1.1. Contributions and overview of paper. We survey the contents and the main contributions of this paper.
Canonical angles We have developed bounds for all the canonical angles between the spaces spanned by the exact and the approximate singular vectors. Several different flavors of bounds are provided: 1. The bounds in § 3.1 relate the canonical angles between the exact and the approximate singular subspaces. Analysis is also provided for unitarily invariant norms of the canonical angles. 2. In applications where lower dimensional subspaces are extracted from the approximate singular subspaces, the bounds in § 3.2 quantifies the accuracy in the extraction process. 3. Also in § 3.2 are bounds for the angles between the individual exact and approximate singular vectors, extracted from the appropriate subspaces. Our bounds suggest that the accuracy of the singular values and vectors, in addition to the low-rank approximations, is high provided (1) singular values decay rapidly beyond the target rank k, and (2) the larger the singular value gaps, the higher is the accuracy to be expected. Furthermore, the truncation step to extract the k dimensional subspaces does not significantly increase the canonical angles.
Low-rank approximation This paper provides the first known analysis of the randomized subspace iteration for an arbitrary unitarily invariant norm, with stronger, specialized results for Schatten-p norms. Bounds for the special cases of the Schattenp norm, namely the spectral and Frobenius norms, have already appeared in the literature-our result for the Schatten-p norm recovers these results as special cases.
Singular values We derive upper and lower bounds on the approximate singular values obtained by the randomized subspace iteration. Similar bounds also appear in [11]; however, our proof technique is different. We also present Hoffman-Weilandt type bounds for the accuracy of the singular values.
The conclusion of the bounds for the low-rank approximations and the singular values are similar to those of the conclusions for the canonical angles. We also make the following contribution, maybe of independent interest beyond the study of randomized algorithms.
Generalization of sin theta theorem The sin theta theorem [24] is a well known result in numerical analysis and relates the canonical angles between the true and approximate singular subspaces in the unitarily invariant norms. We derive a generalization of the sin theta theorem that derives bounds for the individual canonical angles between the two subspaces. The sin theta theorem is recovered as a special case.

Background and preliminaries.
2.1. Notation. Denote the target rank by k. Let the matrix A ∈ C m×n , with m ≥ n, have the singular value decomposition Here Σ 1 ∈ C k×k and Σ 2 ∈ C (n−k)×(n−k) ; the columns of U 1 and U 2 are the corresponding left singular vectors, and columns of V 1 and V 2 are the corresponding right singular vectors. We also denote by [[A]] k = U 1 Σ 1 V * 1 as the best rank−k approximation to the matrix A, in any unitarily invariant norm (for a definition, see below).
Singular values and ratios. Let · 2 denote the spectral norm, so that Σ 2 2 = σ k+1 and Σ −1 The singular values of A can be arranged in decreasing order as For later use, we define the singular value ratios Since the singular values are monotonically decreasing, the singular value ratios are monotonically increasing, i.e., γ 1 ≤ · · · ≤ γ k . Norms. We have already defined the spectral norm. The Frobenius norm of a matrix is A F = trace (A * A). We use the symbol |||·||| to denote any unitarily invariant norm. An example of the unitarily invariant norms is Schatten-p class of norms, defined as the vector ℓ p norm of the singular values of A, i.e., With this definition it can be readily seen that A 2 = |||A||| ∞ and A F = |||A||| 2 . Another example is the Ky-Fan-k class of norms defined A (k) = k j=1 σ j for every k = 1, . . . , min{m, n}. Associated with every unitarily invariant norm is a symmetric gauge function acting on the singular values of the matrix that it acts on.
Projection matrices. Suppose the matrix Z has full column rank Z † is a left multiplicative inverse and where † represents the Moore-Penrose inverse. We define the projection matrix P Z = ZZ † . For a matrix Q with orthonormal columns, the formula simplifies and P Q = QQ * .
Canonical angles. The separation between subspaces can be measured by the principal or canonical angles. Let M and N be two subspaces of C n , such that dim M = k, dim N = ℓ and ℓ ≥ k. Then the principal angles between the subspaces M and N are recursively defined to be the numbers 0 ≤ θ i ≤ π/2 such that subject to the constraints u i 2 = v i 2 = 1, and The canonical angles are arranged in increasing order as It can also be shown that sin θ i are also the singular values of P M − P N . We denote ∠(M, N ) to be the canonical angles between subspaces M and N . Let M and N be matrices with orthonormal columns, which form bases for subspaces M and N respectively. Then, sin ∠(M, N ) are obtained as singular values of (I − M M * )N and cos ∠(M, N ) are the singular value of M * N [5]. For ease of notation, in the rest of this paper, we write ∠(M, N ) instead of ∠(M, N ).

2.2.
Randomized subspace iteration. The basic version of the randomized subspace iteration is summarized in Algorithm 1. Given a starting guess, denoted by Ω ∈ C n×(k+ρ) , the algorithm performs q steps of the randomized subspace iteration to obtain the matrix Y , also known as the "sketch." The algorithm to compute Algorithm 1 Idealized version of Subspace iteration for Singular Value Decomposition Require: Matrix A, Starting guess Ω ∈ C n×(k+ρ) , an integer q ≥ 0.
an approximate singular value decomposition, given starting guess Ω ∈ C n×(k+ρ) is summarized in Algorithm 1. We say this is an idealized version, since the algorithm can behave poorly in the presence of round-off errors. A practical implementation of this algorithm alternates the QR factorization with matvecs involving A; the reader is referred to [20,12]. In Algorithm 1 the output has rank k + ρ. If a rank-k approximation to A is desired, then it can be obtained by discarding the p smallest singular values ofÂ. We denote this low-rank representation by 1 . This is summarized in Algorithm 2.
Partition the matrix V * Ω, which captures the influence of the starting guess on the right singular matrix, as where Ω 1 ∈ C k×(k+ρ) and Ω 2 ∈ C (n−k)×(k+ρ) . The following assumptions will be required for our analysis.
Let Ω 1 be defined as above. We assume that (2) rank (Ω 1 ) = k, and the singular value gap is inversely proportional to The first assumption guarantees that the starting guess Ω has a significant influence over the right singular vectors, whereas the second assumption ensures that the k dimensional subspace R (U 1 ) is well defined.
3. Accuracy of singular vectors. We want to understand how well R ( U ) approximates the subspace spanned by the first k singular vectors U 1 , measured in terms of the canonical angles between the subspaces.
3.1. Bounds for canonical angles. Our first result derives bounds for the canonical angles ∠(U 1 , U ). The analysis is based on the perturbation of projectors and the tools used here are similar to [12].
Theorem 1. Let U and V be obtained from Algorithm 1. With Assumption 1, the canonical angles θ j and ν j satisfy for j = 1, . . . , k This theorem has several interesting features worth pointing out. First, if the matrix has exact rank k, then all of the canonical angles are uniformly zero; that is, the randomized subspace iteration identifies the subspace exactly. Suppose the rank ρ ≥ k. The bounds for the canonical angles show explicit dependence on the singular value ratios γ j . In particular, the smaller the singular value ratio, smaller the canonical angles. Furthermore, the canonical angles converge to zero at different rates.
Second, the term Ω 2 Ω † 1 2 can be written in terms of the singular vector matrix V and the starting guess Ω as When the columns of Ω is linearly independent, this quantity is nothing but the tangent of the largest canonical angle between R (V 1 ) and R (Ω). This term appears frequently in randomized linear algebra and can be interpreted as a measure of the subspace overlap between the starting guess and the right singular vectors. In the ideal case, Ω contains the singular vectors in V 1 . A discussion of the meaning and interpretation of this term, is provided in [8, Section 2.5]. In particular, when Ω is a Gaussian random matrix Ω 2 Ω † 1 2 is roughly on the order of (n − k)k.
Third, the influence of Ω 2 Ω † 1 2 is subdued by the singular value ratios γ 2q+1 j . the bounds are informative in the sense that the sine of the canonical angles are at most 1, which is a desirable feature. With sufficiently large number of iterations, the denominator is O(1), the numerator is a sufficiently good approximation.
Lastly, it is worth mentioning here that the bounds for the canonical angles θ j are smaller than ν j because of the presence of an additional power of γ j . The reason for this higher accuracy is because the columns of V are the right singular vectors of Q * A. The multiplication step with Q amounts to an additional step of subspace iteration and gives the extra factor. Theorem 1 gives the sine of the canonical angles; these bounds can also be used to obtain upper bounds for the tangents and lower bounds for the cosines.
Remark 1. For j = 1, . . . , k, the relationship between the tangent and sine implies Lower bounds for cosine of the canonical angles follow similarly.
Unitarily invariant norms. The following result derives bounds for the canonical angles in any unitarily invariant norm, in contrast to Theorem 1 which bounds the individual canonical angles.
Let Ω ∈ R n×(k+ρ) . Let the approximate singular vectors U and V for a matrix A be computed according to Algorithm 1. Under Assumption 1 , for every unitarily invariant norm, The interpretation of this theorem is similar to that of Theorem 1. The connection between the two theorems follows from the identity sin θ k = sin ∠(U 1 , U ) 2 . It is immediately seen that the result in the spectral norm is weaker.

Extraction of k-dimensional subspaces.
In the previous subsection, the columns of U and V spanned ℓ = k + ρ dimensional subspaces. Many applications, however, require the extraction of k dimensional singular subspaces from the low-rank approximation A ≡ QQ * A. One way to extract the appropriate subspaces is to first compute the optimal rank-k truncation of A, denoted by [ We would like to understand how much additional error is incurred during this extraction process, for which we present several results. The important conclusion is that the accuracy of the extracted subspaces of dimension k is comparable to the accuracy of the k + ρ dimensional subspace up to a modest constant. The approach we take is different from that of the previous section. The starting point of our analysis is the well-known sin-theta theorem for singular subspaces [24]. Let A, A be two matrices of conformal dimensions. Assuming where the two matrices E 12 and E 21 are However, this version of the sin theta theorem does not provide us with a way to obtain bounds for the individual canonical angles. To this end, we first present a new generalization of the sin theta theorem.
Theorem 3. Let A ∈ C m×n with rank (A) ≥ k and let A be the perturbed matrix with same dimensions. Suppose the singular value gap satisfies (4) Let U 1 and V 1 be the left and right singular vectors of A corresponding to the top k singular values.
This theorem states that the sine of the canonical angles sin θ ′ j are bounded by sin θ ′ k up to a multiplicative factor, which is at the most 1.
Our main result provides the following bounds canonical angle when both the subspaces have the same dimension. The proof involves simplifying every term in equations (6).
Theorem 4. Let U and V be obtained from Algorithm 1, and matrices U 1 and V 1 from Algorithm 2. Under assumption 1, • for every unitarily invariant norm The factor φ = √ 2 for an arbitrary unitarily invariant norm, and φ = 1 for the spectral and Frobenius norms.
The interpretation of this theorem is: as the number of iterations q increase, the largest canonical angle converges to 0 quadratically and a larger singular value gap means that the subspace is computed more accurately. Compared to Theorem 2, for any unitarily invariant norm, we have gained the factor max{1, √ 2γ k }/(1 − γ k ) in the denominator. For the spectral and Frobenius norms, this bound is only 1/(1 − γ k ).
Both factors are greater than 1, and independent of the number of iterations q. This additional factor is the (one-time) price to be paid by restricting the subspace R ( U ) which is ℓ = k + ρ dimensional, to R ( U 1 ) which is k dimensional. The bound is devastating when γ k ≈ 1; however, in this case, we are in deep trouble to begin with, since the subspaces may not be well defined.
Individual singular vectors. The previous results give insight into the accuracy measured using the canonical angles between the exact and approximate singular subspaces. When individual singular vectors need to be extracted, does the extraction process introduce additional error? The following result quantifies the accuracy of the extraction process.
Theorem 5. Let the approximate singular vectors U and V be computed according to Algorithm 1. With assumption 1, we have the following inequalities for j = 1, . . . k Denote the approximate singular tripets (σ j ,û j ,v j ) for j = 1, . . . , k. Under assumption 1 The first result bounds the angles between the exact singular vector and the corresponding approximate singular subspaces. The second result compares the angles of the exact and the approximate singular vectors. This result also says that the extraction process does not adversely increase the error in the singular subspaces, provided the singular values are well-separated. The convergence of the individual singular vectors tell a similar story to that of Theorem 1. The singular vectors corresponding to the largest singular values converge earlier than the singular vectors corresponding to the smaller singular vectors. This is a consequence of the fact that the singular value ratios are non-decreasing.

3.3.
Comparison with other bounds. The subspace iteration dates to a 1957 paper by Bauer [3] for eigenvalue problems. The analysis of the subspace iteration has also been well-established, for we refer to [19,Chapter 14]. The effect of randomization has been discussed in [6] which has similar results as ours. One major exception in previous works, however, is that they assumed that the oversampling parameter ρ = 0.
Since Ω 1 ∈ C k×k , the assumption rank (Ω 1 ) = k, requires that Ω 1 is invertible, which is a rather strong assumption. They were able to show (in our notation) Returning to this assumption that Ω 1 has full row rank, when Ω is standard Gaussian matrix, Ω −1 with probability at least 1 − δ. By contrast, when Ω 1 ∈ C k×(k+ρ) and p > 0, then when Ω is a Gaussian random matrix, with probability at least 1 − δ .
It is clear that oversampling has an impact on the accuracy of the randomized subspace iteration. When Ω is a random matrix that arises from the SRHT, or Rademacher distributions, a more aggressive form of oversampling ℓ ∼ k log k is necessary to ensure that rank (Ω 1 ) = k. Therefore, by allowing for oversampling, our bounds are applicable to starting guesses that are not just Gaussian random matrices. Not only that, our bounds are also informative for matrices with decaying singular values and significant singular value gap. A recent paper by Nakatsukasa [18] considered the issue of accuracy of extracting singular subspaces for general projection-based approximation methods. In our notation, these refer to relating bounds for ∠(U 1 , U ) to ∠(U 1 , U 1 ). Our bounds for the canonical angles appear to be tighter than the result [18, Corollary 1]; here we do not go into a detailed comparison. This may be because the analysis was applicable to arbitrary subspace projections, whereas ours is specialized to randomized subspace iteration. Furthermore, our analysis is able to bound the individual canonical angles which is missing in [18].

Probabilistic bounds.
Thus far, we have not made specific assumptions on the matrix Ω, as long as it satisfies rank (Ω 1 ) = k. In particular, Ω need not be even be random, and may be deterministic. However, more can be said about the bounds when Ω is random is drawn from a specific distribution.
In many applications, the matrix Ω ∈ R n×(k+ρ) is taken to be the standard Gaussian random matrix. That is, the entries of Ω are i.i.d. N (0, 1) random variables. Here we give a taste of the probabilistic results that can be derived. Let ρ ≥ 2 and define the constant (9) C e = k ρ − 1 + e (k + ρ)(n − r) ρ and for 0 < δ < 1 define the constant Theorem 6 (Probabilistic bounds). Let Ω ∈ R n×(k+ρ) be a standard Gaussian random matrix with ρ ≥ 2. Assume that the singular value ratio γ k < 1. For j = 1, . . . , k, the expected value of the canonical angles satisfy Let 0 < δ < 1 be a user defined failure tolerance. With probability, at least 1 − δ, the following inequalities hold independently for j = 1, . . . , k Probabilistic results of this form can be derived for other distributions, following the strategy in [12]. We will not pursue them here.

Low-rank approximation and Singular values.
4.1. Low-rank approximation. Several results are available for estimating the error in the low-rank approximation A ≈ QQ * A in the spectral and Frobenius norms, when the matrix Q is obtained from the randomized subspace iteration [12,11,26]. As was mentioned earlier, the spectral and Frobenius norms are special cases of the Schatten-p norm, which are examples of unitarily invariant norms.
Here we present the first known analysis of randomized subspace iteration in a unitarily invariant norm.
Theorem 7. Let A ∈ C m×n be computed using Algorithm 1. Under assumption 1, the following inequalities hold in every unitarily invariant norm If a rank−k approximation is desired, then the error in the low-rank approximation is As in Theorem 4, φ = 1 for spectral and Frobenius norms, and √ 2 for an arbitrary unitarily invariant norm.
In this theorem, as the number of iterations q → ∞, the error in the low-rank approximation goes to zero.
We present a variant of the error in the low-rank approximation for the special case that a Schatten-p norm is used. The proof for the special case of the Frobenius norm was provided in [26].
Theorem 8. Let A be computed using Algorithm 1. Under assumption 1, we have The error bound in Theorem 7 is weaker than Theorem 8 for the Schatten-p norm since for α, β ≥ 0, we have α 2 + β 2 ≤ α + β. More generally, Theorem 8 is applicable to any unitarily invariant norm that is also a Q-norm [4, Definition IV.2.9]. A unitarily invariant norm |||·||| Q is a Q-norm, if there exists another unitarily invariant norm |||·||| a such that |||A||| 2 Q = |||A * A||| a . Note that the Schatten-p norms satisfy this property for p ≥ 2, since |||A||| It can be readily seen that the large singular values are computed more accurately since the singular value ratio corresponding to larger singular values is smaller. This result also shows that in the worst case scenario the singular values are a factor 1/ √ 2 smaller.
Rather than quantify the accuracy of the individual singular values, the next results are of the Hoffman-Weilandt type and account for all the singular values together. Define the two matrices of conformal sizes Under Assumption 1, the error in the singular values satisfies The proof combines [4, III.6.13] with Theorem 7. For the Schatten-p norm, with p ≥ 2, we can derive the bound The proof is similar, and is therefore omitted.

Proofs.
We recall some results here that will be useful in our analysis, see [

Proofs of Section 3.1 Theorems.
Proof (Theorem 1). We tackle each case separately. Bounds for sin θ j The proof is lengthy and proceeds in four steps. We give a great level of detail here, since the proof technique will be applicable to the subsequent proofs.
1. Converting an SVD to an EVD. We compute the thin SVD of (I − P U )U 1 = KS U G * . The matrix S U = diag (sin θ k , . . . , sin θ 1 ) ∈ R k×k contains the sine of the canonical angles between the subspaces spanned by the columns of U and U 1 [5, Equation (13)]. It is readily seen that 2. Shrinking space. In Algorithm 1, we had defined Y = (AA * ) q AΩ. It follows that where Ω 1 ≡V * 1 Ω and Ω 2 ≡V * 2 Ω. Next, by Assumption 1, Ω 1 has full row rank and therefore it has a right multiplicative inverse. Define Recall that Y = QR is the thin-QR factorization. Let Q 1 R 1 be the thin-QR factorization of RΩ 2 Ω † 1 ; here Q 1 ∈ C ℓ×k , R 1 ∈ C k×k . From Q 1 Q * 1 I, the conjugation rule implies P Z = U * QQ 1 Q * 1 Q * U U * QQ * U = P U * Q . Since R (U * Y ) = R (U * Q) = R (U * U ), they have the same projectors, so (20) P Z P U * U I − P U * U I − P Z .
Plug in U U * = I into (19), and use (20) to obtain 3. Simplifying P Z . Since P Z = ZZ † , and Z = I F , we have from which, it can be readily seen that To summarize the story so far, GS 2 U G * H.

Applying singular value inequalities.
A straightforward SVD argument shows that the j-th singular value of H satisfies The singular value inequalities (18) imply Plugging this inequality into σ j (H) Since GS 2 U G * H, Weyl's theorem implies sin 2 θ k−j+1 ≤ σ 2 j (H). Take square roots on both sides and rename j ← k − j + 1 to get the desired result.
Bounds for sin ν j Let GS 2 V G * be the eigenvalue decomposition of V * 1 (I − P V )V 1 . Note that the diagonals of S V are the sine of the canonical angles ∠(V 1 , V ). Since V is obtained from the thin SVD of A * Q, R (A * Q) = R( V ) and P V = P A * Q , since an orthogonal projection matrix is uniquely determined by the range. Next, consider Z defined as , from (AV ) * Q = Σ * U * Q, it can be verified that Using an argument similar to (20), we obtain The right hand side simplifies to I − (I + F * F ) −1 . The rest of the proof is similar to that of the proof for sin θ j .
A similar calculation shows that A similar argument gives Combining these relations Recognize that B 2 = σ k+1 ( A). Applying (5) in the spectral norm simplifies the expression since Therefore, Now σ j (X) = sin θ ′ k−j+1 and σ j (Y ) = sin ν ′ k−j+1 . Rename j ← k − j + 1 to finish. Proof (Theorem 4). We tackle each case independently. Unitarily invariant norms Our proof involves simplifying each term in (5), and (6) and has several steps.  Then applying (12) along with submultiplicativity gives 3. Simplifying |||E 21 |||. We can simplify E 21 in (6). First, E 21 = P U1 (I − QQ * )AP V1 , and since P U1 because of strong sub-multiplicativity. Applying Theorem 1 and (12) . Then for β ≥ 0, since γ k < 1 Putting everything together. Plugging in the intermediate quantities into (5), we have Dividing the numerator and denominator by σ k proves the stated result for unitarily invariant norms. Spectral/Frobenius norms Let · ξ denote the spectral and Frobenius norms. The first two steps are identical to the case for unitarily invariant norms. For the third step, using (12) With β defined as before, since γ k < 1, γ 2 k + γ 2 k β 2 / 1 + γ 2 k β 2 ≤ 1. Therefore, . The rest of the proof is the same.
Canonical angles Use Theorem 3. The term on the right hand side is nothing but max sin ∠(U 1 , U 1 ) The rest of the proof involves some simple manipulations.
Proof (Theorem 5). We first address (7). Following the steps of Theorem 1 we have where e j is the j-th column of the k × k identity matrix. Therefore, we have sin 2 ∠(u j , U ) ≤ e * j He j , where H was defined in (21). The inequality H F * F implies Taking square-roots on both sides gives the desired results. The strategy for bounding the canonical angles sin ∠(v j , V ) is very similar and will be omitted. We now address (8), which is a straightforward application of [13, Theorem 2.5]. Let P U = QQ * and P V = I. Then, in our notation, this result takes the form whereγ ′ = max{0, (I − QQ * )A 2 } andδ is as defined in the statement of the theorem. Theorem 8 for the spectral norm impliesγ ′ ≤γ, whereas Theorem 5 implies

Plug in the intermediate steps to obtain the desired bound.
Proof (Theorem 6). With the constant C e defined in [11,Theorem 5.8] says Let α > 0 be a constant. The map x → x/ √ 1 + αx 2 is convex. Therefore, by Jensen's inequality the results in expectation follow.

Proofs of Section 4 Theorems.
It is useful to recall some properties of the unitarily invariant norms. Every unitarily invariant norm |||·||| on C n is associated with a symmetric gauge function on R n . The |||·||| satisfies |||M ||| = (M * M ) 1/2 , since both matrices have the same nonzero singular values. The following inequality for unitarily invariant norms, also known as strong sub-multiplicativity, will be useful [4, (IV.40)] We will also need Proof. Combining the properties of the partial ordering, the eigenvalues of the scaled matrices satisfy λ j (D * AD) 1/2 ≤ λ j (D * BD) 1/2 for all j = 1, . . . , n. Since the matrices are positive semidefinite, the eigenvalues are the singular values and (D * AD) 1/2 (k) ≤ (D * BD) 1/2 (k) for every Ky-Fan-k norm k = 1, . . . , n. By the Fan dominance theorem [4, Theorem IV.2.2], the advertised inequality is true for every unitarily invariant norm.
We are finally ready to prove our main result for error analysis in unitarily invariant norms.
Proof (Theorem 7). Using the unitary invariance of the norms We use (20) combined with Lemma 10 to obtain The square root function is concave on [0, ∞) and Σ ⊤ (I−P Z )Σ is positive semidefinite. Therefore, an extension to Rotfel'd's theorem says [16, Theorem 2.1] Use the inequalities M 1 F * F and M 2 I, along with Lemma 10 gives and the submultiplicativity to obtain the advertised bounds.
Proof (Theorem 8). The proof is similar to Theorem 7. Consider the term of interest |||(I − QQ * )A||| 2 p , which can be simplified to The first equality holds only for p ≥ 2, whereas the last equality follows because of the unitary invariance. As in Theorem 1, we have The use of (20) and Lemma 10 ensures We apply [16, Theorem 2.1] to (23) with f (t) = t to obtain We have used The rest of the proof is similar to that of Theorem 7.
Proof (Theorem 9). From the inequality I QQ * , the conjugation rule gives Then Weyl's theorem implies λ i (A * A) ≥ λ i (A * QQ * A) for i = 1, . . . , k. Relating the eigenvalues to the singular values proves the first inequality.

Test matrices.
To demonstrate the performance of the bounds, we use the following test matrices 1. Controlled gap The first set of test matrices A ∈ R 3000×300 are constructed using the formula where x j ∈ R 3000 and y j ∈ R 300 are sparse random vectors with non-negative entries generated using the MATLAB commands sprand(3000,1,0.025) and sprand(300,1,0.025) respectively. While the formula above is not an SVD, since the vectors are not orthogonal. Nonetheless, the singular values decay like 1/j and the gap between the singular values between 15 and 16 is controlled by the parameter gap. We consider three cases: 2. Low-rank plus noise The matrices are of the form where G ∈ R n×n is a random Gaussian matrix. We consider three cases: (a) Small noise (NoiseSmall) γ n = 10 −2 , (b) Medium noise (NoiseMedium) γ n = 10 −1 ,

19
(c) Large noise (NoiseLarge) γ n = 1. 3. Low-rank plus decay The matrices take the form The unitary matrices U, V are obtained by drawing a random Gaussian matrix, and taking its QR factorization. We distinguish between the following cases (a) Slow decay (DecaySlow): d = 0.5, (b) Medium decay (DecayMedium): d = 1.0, (c) Fast decay (DecayFast): d = 2.0. The first example is adapted from [21], whereas the second and third examples are drawn from [23]. In all the examples, the random matrices were fixed by setting the random seed and we set r = 15. The singular values of all the test matrices are plotted in Figure 1.

Canonical angles.
For the first numerical example, we use the 9 test matrices in Section 6.1. For each matrix, we chose an oversampling parameter ρ = 20 and the target rank k was chosen to be 25. The starting guess Ω was taken to be a random Gaussian matrix. We plot the canonical angles sin ∠(U 1 , U ) in solid lines, the corresponding bounds from Theorem 1 are also plotted in dashed lines. The results are displayed in Figure 2. We make the following observations • The influence of the subspace iterations on the canonical angles is clear: the angles become smaller, as the number of iterations q increases. This implies that the subspace is becoming more accurate. • If there is a large singular value gap in the spectrum, this means that all the canonical angles below that index are captured accurately. This is prominently seen in Figure 2(e), in which there is a large gap between singular values 15 and 16. Similar observations can be made in the other figures. • Similarly, as the decay in the singular values becomes more rapid, the corresponding canonical angles become smaller. The canonical angles for Decay-Fast is much smaller than DecaySlow. • In all the figures the bounds are qualitatively informative, but in some figures, the bounds are also quantitatively accurate. • Similar results were observed for sin ∠(V 1 , V ) and, therefore, omitted. Likewise, the bounds for Theorems 5 and 2 yield similar looking plots and, therefore, we did include them.
Our next experiment is to test the performance of the extraction step. We now compute sin θ ′ j and sin ν ′ j for the test matrices described above. We plot the quantities max{sin θ ′ j , sin ν ′ j } for j = 1, . . . , k in solid lines. The corresponding bounds from Theorem 4 are plotted in dashed lines. Here the target rank was chosen to be k = 15, to exploit the singular value gap here. We make the following observations • The extraction step did not significantly affect the canonical angles. Conclusions similar to Theorem 1 can, therefore, be drawn. The subspaces are more accurate as the number of iterations increase, and if there is a large singular value gap at index j, then the canonical angles with index j ′ < j are captured accurately. • Compared to Theorem 1, the bounds in Theorem 4 are only qualitatively good. One reason is that the error bounds in Theorem 1 are strictly less than 1, but this is not true for the bounds in Theorem 4. Furthermore, the bound in Theorem 4 has the factor 1/(1 − γ k ) in the denominator, which can be quite large. Better bounds may be possible, but we could not immediately see how to derive them.
• We also compared the accuracy of the individual singular vectors (not shown here). The results and the conclusions are similar.

Singular Values.
We now consider the accuracy of the singular values. We use the same test matrices and the remaining parameters are kept fixed. The computed singular values are plotted against the upper and lower bounds. We make the following observations.
• were the worst for GapSmall.
• As the number of iterations increase, the singular values are increasingly accurately estimated and are practically indistinguishable from the upper bounds (the exact singular values). However, the lower bounds are not great approximations for the singular values close to the target rank. However, the bounds improve as the number of iterations increase. • The bounds for the singular values captured the behaviour of the singular values far better than the bounds for the canonical angles.
7. Acknowledgements. The author would like to thank Ilse Ipsen and Andreas Stathopoulos for helpful conversations. He would also like to acknowledge Ivy Huang for her help with the figures.