Sampling the eigenvalues of random orthogonal and unitary matrices

We develop an efficient algorithm for sampling the eigenvalues of random matrices distributed according to the Haar measure over the orthogonal or unitary group. Our technique samples directly a factorization of the Hessenberg form of such matrices, and then computes their eigenvalues with a tailored core-chasing algorithm. This approach requires a number of floating-point operations that is quadratic in the order of the matrix being sampled, and can be adapted to other matrix groups. In particular, we explain how it can be used to sample the Haar measure over the special orthogonal and unitary groups and the conditional probability distribution obtained by requiring the determinant of the sampled matrix be a given complex number on the complex unit circle.

In applications, one is often interested in random matrices with a given structure. In quantum mechanics, for example, the energy levels of a system are described by the eigenvalue of its Hamiltonian, a Hermitian operator on an infinite-dimensional Hilbert space. By approximating this space by a Hilbert space of finite dimension, one can reduce the problem of finding the energy levels to that of solving a Hermitian eigenvalue problem. The true Hamiltonian, however, is typically not known, thus it is customary to make statistical assumptions on the distribution of its entries, enforcing only the symmetry of the operator. The distribution of the eigenvalues of random symmetric and Hermitian matrices has been extensively studied [16,24,32] and an algorithm for sampling the eigenvalues of uniformly distributed Hermitian matrices has been developed by Edelman, Sutton, and Wang [12].
Here we are interested in matrices sampled from the orthogonal group Opnq " tQ P R nˆn | Q T Q " I n u, and the unitary group U pnq " tU P C nˆn | U˚U " I n u, where I n denotes the identity matrix of order n and A T and A˚denote the transpose and the conjugate transpose of A, respectively. It is easy to prove that the determinant of an orthogonal or unitary matrix lies on the unit circle, and that the special orthogonal group SOpnq " tQ P Opnq | det Q " 1u and the special unitary group SU pnq " tU P U pnq | det U " 1u are subgroups of Opnq and U pnq, respectively. Opnq is made of two connected components, the already mentioned SOpnq and one in which all matrices have determinant´1, which we denote by O´pnq. Clearly, the latter is not a group.
Random unitary matrices find application in quantum physics where they are employed, for example, to model scattering matrices and Floquet operators [15,Section 2.1]. Random orthogonal matrices, on the other hand, are used in statistical mechanics to characterize the behavior of certain log-gas systems [15,Section 2.9].
For a group G, the measure µ such that µpGq " 1 is a normalized left or right Haar measure if for any Q P G and any measurable G Ă G it satisfies µpQGq " µpGq or µpGQq " µpGq, respectively. For compact Lie groups, it can be shown that the left and right measures are unique and coincide. Hence, since Opnq, SOpnq, U pnq, and SU pnq are all compact Lie groups [22, Chapter 1], they have a unique normalized (left and right) Haar measure [23, § 58, § 60].
We consider the problem of sampling efficiently the joint eigenvalue distribution of unitary (or orthogonal) matrices distributed according to the Haar measure. Numerically, this may be obtained by sampling matrices from the desired group uniformly, and then computing their eigenvalues by relying, for instance, on the QR iteration. The latter step requires Opn 3 q floatingpoint operations (flops) to sample the n eigenvalues of a matrix of order n. The key observation is that for this task it is not necessary to explicitly sample matrices from the corresponding group, but it suffices to understand the distribution of their Hessenberg forms, which we analyze in detail in Section 4. The main advantage of this approach is that unitary or orthogonal matrices in Hessenberg form can be diagonalized in Opn 2 q flops. We will exploit this to derive the algorithm discussed in Section 5, which has quadratic complexity and linear storage requirements.
The algorithm we propose can efficiently sample the joint distribution of the eigenvalues of Haar-distributed matrices from any of the Lie groups Opnq, SOpnq, U pnq, and SU pnq. In Section 6 we show that the empirical phase and spacing of eigenvalues sampled by our algorithm follow the corresponding theoretical distributions for U pnq, and then we explore empirically the distribution of the eigenvalues of matrices from the Haar distribution of SU pnq, Opnq, SOpnq, and O´pnq, for which fewer theoretical results are available.
Our starting point is an algorithm proposed by Stewart [30] for sampling random matrices from the Haar distribution of Opnq. We recall this approach and the subsequent generalization to U pnq, due to Mezzadri [25], in Section 3. This technique exploits an algorithm for the QR factorization based on Householder transformations, which we revise in Section 2.
These techniques require the sampling of Opn 2 q random variables, and need Opn 2 q memory for storing the result. We provide an alternative and more compact formulation for the Hessenberg form obtained by the algorithms above, which requires the sampling of Opnq random variables and Opnq storage. We show that this formulation can be used to compute the eigenvalues in Opn 2 q floating point operations by leveraging the unitary QR algorithm in [5].
The use of a condensed factorization for storing random matrices has been already explored, for instance, by Edelman and Ure [13], who sample unitary matrices by taking random Schur parameters [20]. Methods similar to the technique presented in this work might be obtained by representing the Hessenberg forms of unitary Haar-distributed matrices using Schur parameters, or similarly condensed forms such as CMV matrices [9], and then using a quadratic method to compute their eigenvalues [2,3,4,6,7,18,20].
The approach discussed here is based on the unitary QR algorithm in [5,Section 5], The latter can be seen as a special case of the rootfinding algorithm of Aurentz et al. [2], which has been proven to be backward stable [4], and compares favorably with the methods above in terms of performance [2].
Finally, we introduce some notation. Throughout the manuscript, we use capital letters (A, B, . . . ) to denote matrices, lower case letters (u, v, . . . ) to denote vectors, and lower case Greek letters (α, β, . . . ) to denote scalars. We indicate the entries of matrices and vectors using a subscript notation, so that a ij denotes the entry in position pi, jq of the matrix A and v k refers to the kth element of the vector v. We use the same notation for random variables.
We denote by N R pµ, σ 2 q the Gaussian distribution centered at µ P R with variance σ 2 , and by N pm,nq R the distribution of mˆn random matrices with independently distributed Gaussian entries, that is, The chi-squared distribution with k degrees of freedom, denoted by χ 2 pkq, is the distribution of the sum of the squares of k independent Gaussian random variables, and is formally defined by These are real-valued distributions. The complex counterpart of N R pµ, σ 2 q is denoted by N C pµ, σ 2 q and defined by The distribution N pm,nq C is defined analogously.

Householder transformations and QR factorization
In this section we briefly recall some basic facts about Householder transformations, and discuss how they can be employed to compute the QR factorization of a square matrix.
Let v P C n be a nonzero vector. The matrix is a Householder transformation. It is easy to verify that P pvq is unitary and Hermitian, and in particular is orthogonal and symmetric if the entries of v are real. This implies that P pvq 2 " I n . Moreover, computing the action of P pvq on a vector requires only Opnq flops, instead of the Opn 2 q flops that would be needed for a generic nˆn matrix-vector product. Householder transformations are a convenient tool to zero out the trailing entries of a nonzero vector u P C n . For instance, let θ 1 " Arg u 1 , where Arg : C Ñ p´π, πs denotes the principal value of the argument function. Then the matrix P pvq for v " u`e iθ 1 u 2 e 1 is such that P pvqu "´e iθ 1 u 2 e 1 , where e i denotes the ith column of the identity matrix. In order to zero out only the last n´k´1 components of u P C nˆn , it suffices to consider the block matrix p P k puq :" where u pkq P C n´k collects the last n´k entries of the vector u.
Any matrix A P R nˆn has the QR factorization A ": QR, where Q P Opnq and R P R nˆn is upper triangular [19,Theorem 5.2.1]. This factorization is unique up to the sign of the diagonal entries of R. This result can be extended to complex matrices: any matrix A P C nˆn has a unique QR factorization A ": QR, where Q P U pnq and R P C nˆn is upper triangular with real positive entries along the diagonal [31,Theorem 7.2]. More generally, the QR factorization of a full-rank matrix is unique as long as the phases of the diagonal entries of R are fixed.
In many of the following proofs, it will be useful to assume that the matrix A under consideration has full rank. This is typically not restrictive, since rank-deficient matrices are a zero-measure set in N pn,nq R and N pn,nq C ; we will comment on this fact in further detail when needed.
We now explain how to compute efficiently the QR factorization of an nˆn complex matrix A by means of Householder reflections [19,Section 5.2.2]. The corresponding procedure for real matrices can be obtained by employing real Householder reflectors. Let the matrix A p0q :" A be partitioned by columns as and let r P 1 :" p P 1`A p0q 1˘. We obtain that . If we apply this procedure recursively to the trailing submatrix A p1q , after n´1 steps we obtain where the Q :" r P 1¨¨¨r P n´1 is unitary and R is upper triangular. This algorithm produces the triangular matrix R and the factors r P i for i " 1, . . . , n´1 in 4n 3 {3 flops [19, Section 5.2.2]. Computing the matrix Q explicitly requires an additional 4n 3 {3 flops [19, Section 5.1.6], but by exploiting the structure one can compute the action of Q on a vector or matrix with only n 2 and n 3 flops, respectively.
In order to normalize the factorization, note that if D is the diagonal matrix such that d ii "´e´i θ i where θ i is defined as in (2.2), then the matrix DR has positive real entries along the diagonal. Therefore, the normalized factorization with positive diagonal entries in the upper triangular factor is: A " r Q r R, r R :" D˚R, r Q :" r P 1 . . . r P n´1 D, D :"´diagpe iθ 1 , . . . , e iθ n q, (2.4) where r P 1 , . . . , r P n´1 are as in (2.3). A matrix H P C nˆn is in upper Hessenberg form if h ij " 0 when i ą j`1. Any square matrix is unitarily similar to an upper Hessenberg matrix, that is, for any A P C nˆn there exists a matrix U A P U pnq such that U A AUÅ is an upper Hessenberg matrix. The matrix U A can be constructed as the product of n´2 Householder transformations using a technique analogous to that discussed in this section to compute the QR factorization.
3 The Haar measure and Stewart's algorithm Birkhoff and Gulati [8,Theorem 4] note that if the QR factorization A ": QR of an nˆn matrix A ∼ N pn,nq R is normalized so that the entries along the diagonal of R are all positive, then Q is distributed according to the Haar measure over Opnq.
This observation suggests a straightforward method for sampling Haar distributed matrices from Opnq. One can simply generate an nˆn real matrix A ∼ N pn,nq R , compute its QR decomposition A ": QR, and normalize it as discussed in Section 2. This procedure is easy to implement, since efficient and numerically stable routines for computing the QR factorization are available in most programming languages.
The computational cost of this technique can be approximately halved by computing the QR factorization implicitly. Stewart [30] proposes an algorithm that does not explicitly generate the random matrix A, but produces the transpose of the matrix Q in factored form as D r P 1 . . . r P n´1 , where D is an nˆn diagonal sign matrix and r R . This algorithm readily generalizes to the complex case, as discussed by Mezzadri [25]. In order to sample Haar-distributed random matrices from U pnq, it suffices to generate vectors with entries drawn from the standard complex normal distribution N C p0, 1q, and replace the real sign function signpαq by its complex generalization e i Argpαq .
We outline this approach in Algorithm 3.1. The function Umult(X, F) computes the action of a Haar-distributed matrix from Opnq (if F " R) or U pnq (if F " C) on the rectangular matrix X P C nˆm . In order to determine the computational cost of the algorithm, note that asymptotically only the two matrix-vector products and the matrix sum on line 3.8 are significant. Therefore, each iteration of the for loop starting on line 3.2 requires 4km flops, and the computation cost of Algorithm 3.1 is approximately 2n 2 m flops.
In order to sample Haar-distributed matrices, it suffices to set X to I n . In this case, the computational cost of the algorithm can be reduced by taking into account the special structure of A: the cost of line 3.8 drops to 4k 2 flops per step, which yields an overall computational cost of 4n 3 {3 flops.

Sampling from the special groups
The ideas presented so far can be modified in order to sample matrices with prescribed determinant. Imposing that the determinant be 1 is of particular interest, as it implies sampling from the compact Lie groups SU pnq and SOpnq. As discussed in the previous section, the QR factorization of a random matrix A can be used to sample matrices distributed according to the Algorithm 3.1: Action of a matrix from the Haar distribution.
Haar measure over U pnq and Opnq. An analogous result holds for the special groups, if the last diagonal entry of R is chosen so that det Q " 1.
) and let A ": QR be its QR factorization, where Q and R are chosen so that , whenever A is nonsingular. Then, Q is distributed according to the Haar measure over SU pnq (resp. SOpnq).
Proof. We consider the complex case first. Note that the set of rank deficient matrices has measure zero in N pm,nq C ; therefore, the distribution of the unitary QR factor of such matrices is irrelevant for the distribution under consideration.
When A is nonsingular, fixing the phases of the diagonal entries of R makes the QR factorization unique. Hence, the random variables q ij , for i, j " 1, . . . , n are well-defined. In addition, the choice of the diagonal entries of R ensures that det R " det A and thus that det Q " 1.
In order to prove that Q is distributed according to the Haar measure over U pnq, we need to show that it has the same distribution as P Q for any constant matrix P P SU pnq. For any such P , the matrix P A has the QR factorization P A ": pP QqR.
Being independent Gaussian random variables, the entries of A are invariant under unitary transformations, thus P A has the same distribution as A. The triangular QR factor of P A is R, which necessarily satisfies the normalization constraints on the diagonal entries. Therefore, P Q has the same distribution as Q.
The proof for the real case is analogous and therefore omitted.
The above result yields a method for sampling the Haar distribution of the special unitary and orthogonal groups. In the next sections, we discuss how to make this method efficient for sampling the corresponding eigenvalue distribution. The approach we propose can be used for both U pnq and Opnq, and SU pnq and SOpnq.
Remark 3.2. Note that the Haar distribution of SU pnq coincides with the Haar distribution of U pnq conditioned to the event det Q " 1. This is easily verified by checking that the latter measure is invariant under the action of elements in SU pnq. This suggests that the above procedure can be further generalized and used to sample the probability µ ξ obtained by conditioning µ with det Q " ξ for some ξ P S 1 , where S 1 " tξ P C : |ξ| " 1u denotes the complex unit circle. If ξ ‰ 1, these matrices do not form a group, but we can write where P is any constant matrix such that det P " ξ. Sampling the matrices in SU pnq and then multiplying them by any fixed P yields the correct conditional probability distribution. More specifically, by choosing the diagonal matrix P " diagp1, . . . , 1, ξq we can readily adapt the algorithm discussed in the next section to sample unitary or orthogonal matrices with determinant ξ.

The eigenvalue distribution
Given a random matrix sampled according to one of the measures described so far, we are interested in describing the distribution of a generic eigenvalue. This can be computed as a marginal probability by integrating the joint eigenvalue distribution with respect to n´1 variables. For U pnq, the latter is known explicitly [24,Chapter 11], and can be used to prove that a generic eigenvalue is uniformly distributed over S 1 .
We are unaware of an analogous result for SU pnq, and we could not find any references stating the expected distribution for a generic eigenvalue. Nevertheless, using the fact that the eigenvalue distribution arises from Haar-distributed matrices, we can obtain the partial characterization in the following lemma. The well-known distribution for U pnq is for ease of comparison. Lemma 3.3. Let µ and µ 1 be the Haar distributions over U pnq and SU pnq, respectively, and let Λ µ and Λ µ 1 be the corresponding distributions for a generic eigenvalue. Then, Λ µ is the uniform distribution over S 1 , and Λ µ 1 has a 2π n -periodic phase, that is, Proof. We start considering the distribution of U pnq. Recall that µ is invariant under left multiplication in U pnq, and since ξI P U pnq for any ξ P S 1 , we have that Λ µ is invariant under multiplication by ξ P S 1 . Being S 1 a compact Lie group, Λ µ must be its Haar measure, which is the uniform distribution. We can use a similar argument for SU pnq. Since ξI P SU pnq for any ξ such that ξ n " 1, we have that Λ µ 1 is invariant under multiplication by a root of the unity, thus must be 2π n -periodic as in (3.1).
In Section 6 we will verify this claim experimentally, to test the correctness of our implementation. In particular, we will find that Λ µ is the uniform distribution, as expected, and that Λ µ 1 has the periodicity predicted by Lemma 3.3.

The Hessenberg form of Haar-distributed matrices
As mentioned in the introduction, the unitary QR algorithm of [5] cannot be applied directly to the representation of the upper Hessenberg form of a random matrix proposed by Stewart and Mezzadri. In this section, first we show how to sample a factorization of the upper Hessenberg form of Haar-distributed matrices using only Opnq random variables, then we explain how to rewrite this factorization in a form that is suitable for computing the eigenvalues with corechasing algorithms, which we briefly review in Section 5. The main result of this section is the following.
Theorem 4.1. Let w 1 , . . . , w n´1 P C 2 be independent random vectors such that and let H " P 1 . . . P n´1 D P C nˆn (4.1) be the unitary Hessenberg matrix such that where θ j " Arg α j for j " 1, . . . , n´1 and θ n ∼ U p´π, πs. Then, the joint eigenvalue distribution of H is that of Haar-distributed unitary matrices of U pnq.
Proof. In view of the discussion in Section 3, we can sample the joint eigenvalue distribution of Haar-distributed matrices by computing the eigenvalues of the unitary factor of the normalized QR factorization in (2.4), which has the form Q " r P 1 . . . r P n´1 D.
We now proceed to reduce the matrix Q to upper Hessenberg form using similarity transformations. Let q :" Qe 1 be the first column of Q and let U 1 " p P 2 pqq, where the matrix p P k is defined in (2.2). The first step of the reduction gives where r D " diagp1, d 22 , . . . , d nn q. In the formula above we used the fact that D r D´1 commutes with U 1 and r P 2 , . . . , r P n´1 . Now note that for some u P C n , which implies that the first column of Q depends only on r P 1 and on the first entry of D, and that by construction the trailing n´2 entries of U 1 u are zeros. Therefore, the matrix P 1 :" U 1 r P 1 U1 D r D´1 is of the required form (4.2). The matrix where r Q P C n´1ˆn´1 , is unitary by construction and is also Haar distributed, since r U 1 P U pn´1q is independent of r Q. The process can be repeated recursively on the smaller matrix r U 1 r Q r U 1 , and the result follows by induction.
The analogue of Theorem 4.1 for real matrices is obtained by replacing the first element of w j by α j ∼ N R p0, 1q and by sampling θ n uniformly from the set t0, πu. The result can be easily modified in order to sample the joint eigenvalues distribution of matrices from the Haar measure over SU pnq and SOpnq: setting guarantees that det H " 1, while Lemma 3.1 ensures that the matrices are sampled according to the Haar measure.
Remark 4.2. In a similar way, we may set the last diagonal entry of D to obtain det Q " ξ, for any ξ P S 1 . According to Remark 3.2, this procedure samples the Haar distribution conditioned on the event det Q " ξ.  [20], and later investigated by numerous authors, see for instance [1,17,21]. Here, in particular, we are interested in the approach proposed by Aurentz, Mach, Vandebril, and Watkins [5]. This algorithm, briefly described in Section 5.2, is implemented in eiscor [3], a Fortran 90 package for the solution of eigenvalue problems by core-chasing methods available on GitHub. 1 The software computes the eigenvalues of the Hessenberg matrix where the unitary matrices G 1 , . . . , G n´1 are plane rotations of the form Because of its special structure, the matrix G j in (5.2) is said to be "essentially 2ˆ2", and the 2ˆ2 matrix p G j is called a core block. In principle, the core chasing algorithm in [5] could be applied to any factorization of H involving only "essentially 2ˆ2" unitary matrices, even though the particular implementation described in [3] involves only the special family of plane rotations. In practice, however, the key operation in the QR algorithm-the so-called turnover -has to be implemented with care in order to ensure backward stability. In order to leverage the analysis done for rotations of the form (5.2), it is thus convenient to refactorize H given in the form (4.1) as a product of the the form (5.1).
This section is structured as follows. First, we show how to refactorize a representation in terms of 2ˆ2 Householder transformations into one consisting only of plane rotations with real sines. Then we briefly recall the main ideas underlying the unitary QR algorithm implemented in eiscor.

Refactoring core transformations
The assumption that all core blocks in the factorization of the Hessenberg matrix H be plane rotations with real sines is not restrictive, as it is always possible to rewrite H as a product of the form (5.1). This refactorization can be performed by noting that any 2ˆ2 unitary matrix U can be written as where the diagonal entries are given by d 1 " e´i θ |u 11 | 2`eiθ u 2 21 , d 2 " e iθ pu 11 u 22´u21 u 12 q.
The procedure above can be performed in a backward stable manner, as it coincide with the computation of the QR decomposition of U [3].
In addition, the product of a plane rotation with real sines G and a unitary diagonal 2ˆ2 matrix D can be refactored so to swap the order of the two operations. In fact, there exist a unitary 2ˆ2 diagonal matrix r D and a plane rotation with real sines r G such that GD " r D r G. This property is easy to verify, and represents the foundation of most core-chasing algorithms [3]. Combining these two observations gives the following. Proof. The proof is by induction on n. For n " 2 we have that H " P 1 D, and the refactorization can be performed directly by relying on (5.3). If n ą 2, then there exist a plane rotation G 1 P C nˆn as in (5.2) and a unitary diagonal matrix such that P 1 " G 1 D 1 . Since the matrix " d 11 I n´1 ı commutes with D 1 , P 2 , . . . , P n´1 , we can write where r d ij " 1 if i " j " 1 and r d ij " d ij otherwise. We note that this refactorization has the form H " where r H has the same structure as H but order n´1. The inductive hypothesis yields This procedure provides an algorithm for refactoring H from the form (4.1) to the form (5.4).
Noting that each step requires Op1q flops, for a total of Opnq flops for the complete refactorization, concludes the proof.

Computing the eigenvalues of unitary Hessenberg matrices
We have described how to generate unitary upper Hessenberg matrices whose joint eigenvalue distribution follows the Haar measure, and we have shown how to write it in the factored form (5.4).
In order to compute the spectrum of H in the form (5.4) in Opn 2 q flops, we rely on the method proposed in [5], which belongs to the family of core-chasing algorithms [3]. Here we provide a high-level overview of this technique, and refer the interested reader to the original paper [5] for a detailed discussion.
With the term core transformation we indicate an essentially 2ˆ2 unitary matrix such as, for example, a plane rotation. The factorization (5.4) is an example of a matrix expressed as product of n´1 core transformation and a diagonal matrix. In this particular case, the facotrization can also be used to give a compact representation of H that uses only Opnq parameters, as opposed to the Opn 2 q that would be necessary if all the entries of H were explicitly stored. Each core transformation acts on a pair of adjacent indices. The matrix G j in (5.2), for example, acts on the indices j and j`1.
The standard single-shift bulge chasing QR algorithm works as follows. Given an upper Hessenberg matrix H, we determine a first core transformation Q 1 acting on the indices 1 and 2 such that Q 1 pH´ρIqe 1 " αe 1 . The parameter ρ is the shift, and has to be carefully chosen in order to ensure fast and reliable convergence of the method [34]. The implementation considered here uses a Wilkinson shift, which is projected onto S 1 as the matrix H is unitary [5].
We use the core transformation Q 1 to compute Q 1 HQ1, which is not upper Hessenberg having a nonzero element in position p3, 1q. Another core transformation Q 2 acting on the indices 2 and 3 is used to restore the upper Hessenberg structure and obtain Q 2 Q 1 HQ1. The similarity Q 2 Q 1 HQ1Q2, however, yields a matrix that is not upper Hessenberg because of a nonzero element in position p4, 2q, and the process is repeated until the nonzero element, the so-called bulge, is eliminated from the last row of the matrix. The focus on the nonzero element that breaks the upper Hessenberg structure and is "chased to the bottom" until it disappears from the matrix, justifies the name bulge-chasing QR used for this algorithm.
Core-chasing algorithms have a similar formulation, and indeed are mathematically equivalent [3], but do not handle the entries of the matrix directly, as we now explain.
An upper Hessenberg matrix can always be written as H " QR, where R is upper triangular and Q " G 1 . . . G n´1 is the product of unitary plane rotations. The core-chasing step starts by computing a core transformation Q 1 such that Q 1 pH´ρIqe 1 is a scalar multiple of e 1 . Keeping H " QR in factored form, the similarity transformation with Q 1 gives We now make the following two key observations.
• Given an upper triangular matrix R and a core transformation Q 1 , it is always possible to find another upper triangular matrix r R and core transformation r Q 1 such that RQ 1 " r Q 1 r R. Using the terminology of core-chasing algorithms, the computation of r Q 1 and r R from Q 1 and R is a passthrough operation, and can be represented pictorially as .
• Given the matrices G i and K i acting on the ith and i`1st indices, and J i`1 acting on the i`1st and the i`2nd indices, the product G i J i`1 K i can be refactored as r G i`1 r J i r K i`1 , where r G i`1 and r K i`1 act on the i`1st and i`2nd indices, and r J i acts on the ith and i`1st indices. Similarly, this operation can be represented pictorially as " .
In the context of core-chasing algorithms, it is more natural to reinterpret the refactorization above as moving the rightmost core transformation to the left, which can be displayed as .
Clearly, all the rotations involved in the step above change, but from the point of view of the structure, it is as if only one rotation had moved. This operation is called turnover.
With this notation, we can rephrase the factorization obtained after introducing the core transformation Q 1 as In the rightmost factorization, we have first fused the top-left rotations, and then used the passthrough and turnover to take the rotation that was on the right to the left. If we call this new core transformation Q 2 , we can now perform the similarity transformation Q 2 Q 1 HQ1Q2 and obtain the matrix The structure of this matrix is similar to that of Q 1 HQ1, in (5.5), but the rightmost core transformation has moved down one step. Indeed, it now acts on indices 2 and 3 instead of 1 and 2. Carrying on this process will move it further down, until it is fused at the bottom. At the very end, the core transformation will hit the bottom rotation in the sequence G 1 . . . G n´1 , and they will be fused together. This completes the chasing, and is mathematically equivalent to chasing the bulge into the bottom-right corner.
There are a few more technical details to address in order to obtain a complete algorithm. One key point is how to detect deflations, that is, eigenvalues that have converged. In the usual bulge-chasing setting, we monitor subdiagonal elements, setting them to zero as soon as they become "small enough". For core-chasing algorithms, we observe that a subdiagonal element is small if and only if the corresponding rotation in the sequence G 1 . . . G n´1 is close to the identity. In fact, this technique is often more accurate than the customary criterion in practice [3].
The main computational step in this process is the refactorization RQ 1 " r Q 1 r R, which requires Opnq flops in general. Therefore, a single core-chasing step requires Opn 2 q flops, and Opn 3 q flops will be necessary to compute the eigenvalues of a generic Hessenberg matrix, if about Opnq steps are required for the QR iteration to converge. All the other operations require only Op1q flops, and contribute only a low-order term to the total cost. However, if H is unitary to begin with, then so is the upper triangular matrix R. As upper triangular unitary matrices must be diagonal, the passthrough operation can be performed in Op1q flops, making the cost of the QR algorithm quadratic instead of cubic in n. For a more detailed analysis, we refer the reader to the paper where the algorithm was first introduced [5].

Sampling the eigenvalues of random unitary and orthogonal matrices
The approach underlying the discussion in Sections 4 and 5 is summarized in Algorithm 5.1.
The function SampleEigs samples the joint distribution of orthogonal or unitary matrices from a specific distribution determined by the value of third parameter ξ. If ξ is 0, then the function samples the eigenvalues of Haar distributed matrices from the orthogonal group, if F " R, or from the unitary group, if F " C. If ξ ‰ 0, the algorithm samples the eigenvalue distribution of matrices whose determinant has the same phase as ξ. We recall that these matrices form a group only if ξ " 1, in which case the algorithm samples the eigenvalue distribution of Haar-distributed matrices from the special orthogonal group SOpnq, if F " R, or special unitary group SU pnq, if F " C.
In order to achieve this, we note that for H in the form (5.4), we have that det H " det D " e iθ for some θ P p´π, πs, since the determinant of plane rotations is 1. Therefore, once the first n´1 entries of D are chosen, it suffices to choose d n " e i Arg ξ´i Argp ś n´1 j"1 d j q , which ensures that det D " e i Arg ξ . In the pseudocode, the function UnitaryQR denotes a call to the eiscor routine, which computes the eigenvalues of a product of plane rotations of the form (5.2).
The computational cost of the algorithm can be determined by noticing that each step of the for loop on line 5.3 requires only a constant number of operations, which implies that the whole preprocessing taking place between line 5.2 and line 5.18 requires only Opnq flops.

Experimental results
In this section we first validate the new algorithm experimentally, and then compare its performance with that of the naïve method for sampling the joint eigenvalues distribution of Haardistributed unitary matrices. The experiments were run in MATLAB 9.8.0 (R2020a) Update 4 on a GNU/Linux machine equipped with an Intel Xeon E5-2640 v3 CPU running at 2.60GHz.
In our tests we compare the following implementations.
• samplemat, an algorithm that generates a Haar-distributed unitary matrix by calling the function Umult in Algorithm 3.1 on the identity matrix, and then computes its eigenvalues by using the built-in MATLAB function eig.
• sampleeig, an implementation of Algorithm 5.1 that exploits the eiscor package to run the QR algorithm on the unitary matrix in factored form.
Our implementations of samplemat and sampleeig are available on Github. 2 For reproducibility, the repository also includes the scripts we used to run the tests reported here.
Algorithm 5.1: Sample the joint distribution of random matrices matrices. 5.1 function SampleEigspn P N, F P tR, Cu, ξ P Fq Sample eigenvalues of orthogonal (if F " R) or unitary (if F " C) matrices. If ξ ‰ 0, the determinant of the sampled matrices has the same phase as ξ. If ξ " 0, the matrices are sampled from Opnq (if F " R) or U pnq (if F " C).
ϕ Ð e i Arg u 21

Unitary matrices
We start by considering the joint distribution of the eigenvalues of Haar-distributed matrices in U pnq and SU pnq. Figure 1 reports the phase distribution and the spacing of the eigenvalues of 1,000,000 unitary matrices of order 10 sampled from the unitary group (top row) and from the special unitary group (bottom row) using sampleeig. The histograms in the four plots are normalized so that the total area of the columns is one. In this way, the histograms can be interpreted as empirical probability densities and can be compared directly with the probability density functions they are expected to match.
Let e iθ 1 , . . . , e iθ n be the eigenvalue of the matrix A P U pnq normalized so that for i between 1 and n the phase θ i lies in the interval r0, 2πq. The histogram in Figure 1a shows the distribution of the phases of the 10,000,000 sampled eigenvalues, whereas the dashed lines indicates the probability density function of the uniform distribution over the interval r0, 2πq. As the eigenvalues of unitary matrices lie on the unit circle, our results indicates that the eigenvalues sampled by the procedure are uniformly distributed.
Next we investigate the statistical correlation among the eigenvalues sampled by sampleeig. In Figure 1b     of consecutive eigenvalues, defined by where the eigenvalues are ordered so that θ 1 ď¨¨¨ď θ n . The empirical distribution of the sampled eigenvalues matches closely the theoretical one, confirming that the matrices whose eigenvalues sampleeig samples are in fact Haar distributed.
To the best of our knowledge, the probability distribution for the phase and spacing of Haardistributed matrices in the special unitary group are not known in closed form, but we can use sampleeig to obtain a relative frequency distribution based on 10,000,000 samples. The results in Figure 1c and Figure 1d suggest that the phase of the eigenvalues of these matrices is not uniformly distributed, but the spacing appears to be the same as for Haar-distributed matrices in U pnq, as the empirical distribution matches the Wigner surmise in (6.1).
We note that the invariance under multiplication by elements in SU pnq implies that Q and diagpξ, . . . , ξqQ must have the same distribution for any ξ such that ξ n " 1, and the phase of the density of the eigenvalue distribution needs to be 2π{n-periodic, as proven in Lemma 3.3. This is clearly visible in Figure 1c.

Orthogonal matrices
The joint eigenvalue probability density functions for SOpnq and O´pnq are reported explicitly in [15,Section 2.6]. The corresponding joint eigenvalue distribution for the orthogonal group can be obtained easily, since a matrix in Opnq belongs with equal probability to SOpnq or O´pnq. However, we are not aware of a closed form expression for the distribution of the eigenvalues of such matrices, which may be obtained by integrating out n´1 variables from the expressions of the joint eigenvalue distribution.
We can use sampleeig to get the empirical distribution of the phase and spacing of the eigenvalues of these matrices. In Figure 2, we report the relative frequency distribution of the phase and spacing of 1,000,000 random matrices of order 10 sampled from the orthogonal group (top row), from the special orthogonal group (middle row), and from the set of orthogonal matrices with determinant´1 (bottom row). In Figure 3 we report the same data for matrices of order 9, as the behavior of these distributions changes dramatically depending on the parity of n.
The distribution of phase and spacing for the eigenvalues of matrices sampled from Opnq appears identical for both matrix dimensions we consider. In particular, we note that in Figure 2a and Figure 3a there is a mass of probability corresponding to the eigenvalues 1 and´1, which is a consequence of the fact that the eigenvalues of real matrices always appear in conjugate pairs. Therefore, if n is even then matrices with determinant´1 must always have both eigenvalues 1 and´1 (see Figure 2e), whereas if n is odd then all matrices with determinant 1 must have the eigenvalue 1 (see Figure 3c) and all those with determinant´1 must have the eigenvalue´1 (see Figure 3e).

Timings and computational complexity
Now we compare the performance of our MATLAB implementations of sampleeig and samplemat. Figure 4 shows the time, in seconds, required by the two algorithms to sample the eigenvalues of matrices of order n between 2 and 2 15 . For matrices of order up to 16, sampleeig is slightly slower than samplemat; this is due to the fact that normalizing the rotations amounts to a large portion of the overall execution time of the algorithm.
As the computational cost of this operation scales linearly, however, its contribution becomes negligible as n grows: for matrices of order 32 and above the execution time grows much faster for samplemat than for sampleeig. This is expected, since the two algorithms have cubic and quadratic computational cost, respectively.

Conclusions
We have presented a method for sampling the joint distribution of the eigenvalues of Haardistributed orthogonal and unitary matrices. The two ingredients of our approach are a tech-       Our experimental results show that the new technique is more efficient than the naïve method that first samples a matrix from the Haar distribution and then computes its eigenspectrum numerically. We used this algorithm to investigate experimentally the distribution of the phase and spacing of the eigenvalues of Haar-distributed matrices from SU pnq, Opnq, SOpnq, and O´pnq, groups for which these distributions are not known explicitly.