{'fraction_non_alphanumeric': 0.04932700536839648, 'fraction_numerical': 0.028825955029954095, 'mean_word_length': 4.119896434973113, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Projection method for eigenvalue problems of linear nonsquare matrix pencils

Eigensolvers involving complex moments can determine all the eigenvalues in a given region in the complex plane and the corresponding eigenvectors of a regular linear matrix pencil. The complex moment acts as a filter for extracting eigencomponents of interest from random vectors or matrices. This study extends a projection method for regular eigenproblems to the singular nonsquare case, thus replacing the standard matrix inverse in the resolvent with the pseudoinverse. The extended method involves complex moments given by the contour integrals of generalized resolvents associated with nonsquare matrices. We establish conditions such that the method gives all finite eigenvalues in a prescribed region in the complex plane. In numerical computations, the contour integrals are approximated using numerical quadratures. The primary cost lies in the solutions of linear least squares problems that arise from quadrature points, and they can be readily parallelized in practice. Numerical experiments on large matrix pencils illustrate this method. The new method is more robust and efficient than previous methods, and based on experimental results, it is conjectured to be more efficient in parallelized settings. Notably, the proposed method does not fail in cases involving pairs of extremely close eigenvalues, and it overcomes the issue of problem size.

A new wave of eigensolver's development emerged, initiated by the scaling method [18,19] for electronic structure analysis, the Sakurai-Sugiura method [54], and the FEAST algorithm [50]. In this class of projection methods, a complex moment consisting of a resolvent filters out undesired eigencomponents and extracts the desired ones in a pseudo-random matrix whose columns are supposed to have eigencomponents corresponding to the eigenvalues of interest. Thus, methods of this kind project a regular matrix pencil onto the eigenspace associated with eigenvalues in a prescribed region and give the eigenvalues and the corresponding eigenvectors of a regular matrix pencil. The complex moment results from a contour integral of a resolvent matrix. In numerical computations, the integral is approximated by a quadrature rule. Each quadrature point produces a linear system of equations to solve. Each linear system can be solved independently, and hence, this kind of methods can be efficiently implemented in parallel.
Another perspective of this class of projection methods is in the filter diagonalization, leading to the design of rational filters [47,4]. See [28] for mutual relationships among projection methods of this kind. Their advantage is in its hierarchal parallelization. This feature is suitable for the heterogeneous architectures of modern computers. The typical bottleneck in parallelization is the data transfer among memories in different layers and among processors. The following hierarchical procedures can be efficiently executed in parallel in the projection methods: • internal eigencomputations for partitioned regions in the top layer [43,40], • solutions of linear systems arising from quadrature points in the middle layer [20], and • computations of each linear solve in the bottom layer [53,60].
This work considers an extension of such a projection method to singular cases. To this end, motivated by a consideration of spectra of a single nonsquare matrix [6, Section 6.7], the resolvent matrix in the projection method is replaced with the pseudoinverse of the linear matrix pencil. The integral is also approximated by a numerical quadrature, similarly to the regular case. Each quadrature point gives a (block) least squares problem to solve, and it can be solved independently. This extension carries out those features of the above projection methods in parallel by nature. Note that the proposed method does not attempt to determine the sizes of Kronecker blocks, and thus, it is not naively comparable to GUPTRI. Also, GUPTRI is not designed to compute the eigenvalues in a prescribed region and the corresponding eigenvectors and is not designed to be efficient for sparse matrices A and B. Note also that the proposed method will be considered in theory for linear matrix pencils whose singular blocks are of size zero, as will be formally assumed later.

Preliminary
The spectral properties of the linear matrix pencil zB − A are given by the Kronecker canonical form [ where P ∈ C m×m and Q ∈ C n×n are nonsingular matrices,⊕ denotes the direct sum of matrices is a Jordan block corresponding to finite eigenvalues λ i , and N ρ i ∈ C ρ i ×ρ i is the shift matrix whose superdiagonal entries are 1 and the remaining entries are zero. Here, we denote the identity matrix of size η by I η ∈ R η×η and the ith standard basis vector by e i . Note that G(z) and K(z) have the finite and infinite eigenvalues of zB − A, respectively. We call the block matrices G i (z) and K i (z) regular blocks. Hence, matrix pencil zB − A has η finite eigenvalues λ η i , i = 1, 2, . . ., , with multiplicity η i , respectively, and ρ infinite eigenvalues. Meanwhile, we call L i (z) and U i (z) the right and left singular blocks, respectively. Hence, q and r are the numbers of right and left singular blocks, respectively. Let t be the number of finite eigenvalues of matrix pencil zB − A counting multiplicity in the prescribed region Ω. Then, without loss of generality, the eigenvalues in the region Ω are denoted by λ 1 , λ 2 , . . ., λ s , and the associated regular blocks are denoted by Throughout, matrix pencil zB − A is assumed to have left and right singular blocks of size µ = ν = 0 if not specified otherwise. Hence, the proposed method aims at determining the eigenvalues of J i for all i such that λ i ∈ Ω. Note that the proposed method does not aim to determine the Kronecker structures such as the minimal indices µ i and ν i . Further, the matrix pencil zB − A is assumed to be exact. In different contexts, perturbations in matrices A and B are assumed [9,10,30].

Organization
Section 2 gives the formulation of the proposed method and the condition such that the method determines the desired eigenpairs, and presents the implementation issues. Section 3 presents the results of numerical experiments on test matrix pencils. Finally, Section 4 concludes the paper.

Proposed method
This work considers extending the scope of problems that can be solved by a projection method for regular matrix pencils (CIRR [55,25]) to singular nonsquare matrix pencils. The extended method reduces a given nonsquare matrix eigenproblem (1.1) to a generalized square eigenproblem. The main principle of this approach lies in the kth order complex moment matrix, given by the contour integral where π is the circular constant, i is the imaginary unit, Γ is a positively-oriented closed Jordan curve of which Ω is the interior, † denotes the Moore-Penrose generalized inverse (pseudoinverse) of a matrix, and M is the order of moments such as M < n. (See [6, Section 6.7] for the spectral properties of nonsquare matrices.) Throughout, it is assumed that no eigenvalues of J i exist in Γ for all i = 1, 2, . . ., . We reduce (1.1) to a smaller eigenproblem having the eigenvalues λ ∈ Ω of interest. Let L ∈ Z >0 such as L < min(m, n), T ∈ C LM ×m and V ∈ C m×L be random matrices, and Then, the reduction of (1.1) to the square generalized eigenproblem is yielded by a Rayleigh-Ritz-like procedure [52, §4.3] over the subspace R(S) to obtain the desired eigenvalues, where R(·) denotes the range of a matrix.
We deal with both cases m ≥ n and m < n in a unified manner. To derive conditions such that (2.3) has eigenpairs of interest, we prepare a form of matrix pencil zB − A. From the Kronecker canonical form (1.2), the pseudoinverse of matrix pencil zB − A has the form Here, we used the fact [7,Fact 8.4.23] that if E ∈ C ×m has full-column rank and F ∈ C m×n has full-row rank, then E † E = F F † = I m and Hence, (zB − A) † is decomposed into three full-rank matrices (2.4): the first is related to the right singular blocks, the second is related to the regular blocks, and the third is related to the left singular blocks. Assume that the singular Kronecker blocks L(z), U (z) have size µ = ν = 0. Then, the decomposition (2.4) is reduced into where O r,q is the zero matrix of size r × q. The next lemma gives useful expressions of the first factor ([I η+ρ , O]Q −1 ) † and third factor (P −1 [ I η+ρ O ]) † . Here, Π S denotes the orthogonal projector onto a subspace S and S ⊥ denotes the orthogonal complement of S. Lemma 2.1. Let C ∈ C n×n be a nonsingular matrix. Partition columnwise into submatrices C 1 ∈ C n×k and C 2 ∈ C n×(n−k) . Then, we have Noting that the Schur complement of C H 2 C 2 in the 2 × 2 block matrix Partition Q = [Q r , Q s ] columnwise into submatrices Q r ∈ C n×(η+ρ) and Q s ∈ C n×r corresponding to the regular and right singular Kronecker blocks, respectively. Applying Lemma 2.1 to Q, we obtain Additionally, partition P = P r P s rowwise into submatrices P r ∈ C (η+ρ)×m corresponding to the regular Kronecker blocks and P s ∈ C r×m corresponding to the left singular Kronecker blocks. Transposing the statements in Lemma 2.1, we obtain For convenience, as per the Kronecker block structure, we partition rowwise into submatrices corresponding to the regular Kronecker blocks regarding the finite eigenvalues, the regular Kronecker blocks regarding the infinite eigenvalues, and the left singular blocks, respectively. Here, the first column of P i is the left eigenvector corresponding to λ i and the remaining columns of P i are the generalized right eigenvectors for i = 1, 2, . . ., . And, partition Furthermore, according to the Kronecker block structure, we partition columnwise into submatrices corresponding to the regular Kronecker blocks regarding the finite eigenvalues, the regular Kronecker blocks regarding the infinite eigenvalues, and the right singular blocks, respectively. Here, the first column of Q i is the right eigenvector corresponding to λ i and the remaining columns of Q i are the generalized right eigenvectors for i = 1, 2, . . ., . Further, partition Now, we express the moment matrix M k of order k using the Kronecker canonical form via the expansion The next lemma gives the localized moment matrix in a simple form.

Lemma 2.2.
LetP Proof. The expansion (2.9) gives in which two contour integrals of G i (z) and K i (z) are considered. The inverse of a regular Kronecker block corresponding to a finite eigenvalue is of the form From the assumption, G i (z) is regular for all z ∈ Γ. By using Cauchy's integral formula, we have Therefore, we obtain (2.12).
The decomposition (2.12) reduces to that in [25, equation (6)], when m = n and no singular Kronecker blocks L(z) and U (z) exist. Hence, the decomposition (2.12) is a generalization to the singular case. Also, (2.12) extends the formula [21, equation (5.8)] to singular linear matrix pencils.
In the following lemma, the localized moment matrix (2.12) is used to extract the Jordan blocks corresponding to the eigenvalues located in the region Ω. [26,Theorem 5]). Let L, M ∈ Z >0 such that min(m, n) > L ≥ t, W ∈ C L×n , V ∈ C m×L be arbitrary matrices, t be the number of eigenvalues of matrix pencil zB − A in a region Ω, andP Ω ,Q Ω be defined by (2.10), (2.11), respectively. If rank(WQ Ω ) = rank(P Ω V ) = t, then the nonsingular part of the size-reduced matrix pencil Proof. Because rank(WQ Ω ) = rank(P Ω V ) = t ≤ L, there exist nonsingular matrices P ∈ C L×L and Q ∈ C L×L such thatP Hence, we have Lemma 2.3 shows that solving the size-reduced eigenproblem W M 1 V y = λW M 0 V y, y = 0 gives the desired eigenvalues. Based on Lemma 2.3, the following theorem demonstrates an approach to project the matrix pencil zB − A onto the desired Ritz space to form a size-reduced matrix pencil having the desired eigenvalues.
Proof. The coefficient matrices are expressed as Then, we have Here, we used Then, we have Applying [I η+ρ , O] Q −1 † to both sides of the equation (2.13), we have (2.14) The first t columns of both sides of (2.14) form Hence, from the assumption, we have Therefore, because rank(WQ Ω ) = rank(P Ω V ) = t, Lemma 2.3 gives the assertion.
The following lemma prepares the proof of the main theorem.
Lemma 2.5. Let L, M ∈ Z >0 , J Ω ,P Ω ,Q Ω be defined as in Lemma 2.2, V ∈ C m×L , S be defined as in (2.2), and Then, the equalities S =Q Ω Y and rank(S) = rank(Y ) hold.
Proof. From the definition Eq. (2.2) of S and Lemma 2.2, it follows that i.e.,Q Ω = Π R(Qr) ⊥ Q Ω has full-column rank. Hence, rank(S) = rank(Y ) holds. Now, the main theorem can be proved. Remark 2.7. The desired eigenvector x can be obtained from y of (2.3) by x = Sy, where y is the eigenvector corresponding to an eigenvalue λ of (2.3).
Remark 2.8. The principle behind this formulation is the construction of a filter for eigencomponents. Combining Lemma 2.5 with Theorem 2.6 shows that the orthogonal projector Π R(Qs) ⊥ in S filters out the R(Q s ) component in the vector Q Ω Y y, i.e., stop undesired eigencomponents. Remark 2.9. Theorem 2.6 shows that the numbers of L and M must be chosen to satisfy LM ≥ t to obtain the eigenvalues in Ω. The number t of eigenvalues in Ω would be estimated using analogous techniques given in [16,39,48,57]. High-order complex moments with M > 1 enables one to efficiently enlarge the space R(S).
Remark 2.10. The condition rank(Y ) = t in Theorem 2.6 is not necessarily satisfied even ifP Ω V has full-column rank.

Implementation
In numerical computations, the contour integral (2.1) is approximated using the N -point trapezoidal quadrature ruleM where z j is a quadrature point and ω j is its corresponding weight. Thus, we obtain the approx-imationsS Moreover, to reduce computational costs and improve numerical stability, a low-rank approximation of the reduction (2.3) is applied using principal basis vectors of R(S) associated with the truncated singular value decompositions (TSVD) ofS where the columns of U S = [U S,1 , U S,2 ] are the left singular vectors ofS, the columns of V S = [V S,1 , V S,2 ] are the right singular vectors ofS, and Σ S = Σ S,1 ⊕ Σ S,2 ∈ R n×m is a diagonal matrix whose diagonal entries are the singular values ofS, arranged in decreasing order. Here, the diagonal entries of Σ S,1 are the dominating singular values ofS, and the columns of U S,1 and V S,1 are the corresponding left and right singular vectors, respectively. A practical way to determine the size τ of Σ S,1 ∈ R τ ×τ will be given in Section 3. Hence, U S,1 Σ S,1 V S,1 H is a low-rank approximation ofS.
Thus, by solving the reduced square generalized eigenproblem whereT ∈ C m×τ is a random matrix, the approximate eigenpair (λ,x) = (λ, U S,1ỹ ) can be obtained. These procedures are summarized in Algorithm 2.1.   [22,49,15,45,46], where v i is the ith column of the random matrix V and · 2 denotes the Euclidean norm, or the minimum-norm least squares problem with multiple right-hand side

Algorithm 2.1 Proposed method.
by using iterative solvers such as the block CGLS and LSQR methods [31,34] and the global CGLS method [38], where · F denotes the Frobenius norm. The L psuedoinverse solutions (z j B − A) † v i , i = 1, 2 . . ., L, can be computed in parallel by solving the problems (2.16).

Numerical experiments
Numerical experiments show that the proposed method (Algorithm 2.1) is superior to previous methods in terms of efficiency and robustness. The efficiency and robustness were evaluated in terms of CPU time and the relative residual norm (RRN) All computations were performed on a computer with an Intel Core i7-8565U 1.80 GHz central processing unit (CPU), 16 GB of random-access memory (RAM), and the Microsoft Windows 10 Pro 64 bit Version 20H2 operating system. All programs for implementing the proposed method were coded and run in MATLAB R2020b for double precision floating-point arithmetic with unit roundoff u = 2 −53 1.1·10 −16 . The compared implementation of GUPTRI was in the Matrix Canonical Structure (MCS) Toolbox [33].
For the proposed method, the quadrature points z i were set to z i = γ + R exp (iθ i ), θ i = (2i − 1)π/N , i = 1, 2, . . ., N , on the circle with center γ = 1 + i and radius R. The MATLAB function pinv was used to compute (z j B − A) † V for small sizes min(m, n) < 1000. On the other hand, the global CGLS method [31] was used to compute it for large sizes min(m, n) > 1000 for efficiency. Here, the initial guess was set zero and the stopping criteria for the global CGLS iteration were given by the threshold 10 −14 in terms of the relative residual norm and the maximum number of iterations min(m, n). In the row-rank approximation ofS (2.15), the rank was truncated such that the size of Σ S,1 was maximized subject to the constraint on its condition number not exceeding 1/u. Here, the condition number of a matrix is the ratio of the largest singular value to the smallest value of the matrix. The eigenpairs (λ, y) in line 4 of Algorithm 2.1 were computed by using the MATLAB function eig.
For comparison, the following folklore methods were tested. Let V ∈ C n×m be a pseudorandom matrix. On the one hand, suppose m < n. Then, the eigenpairs of square matrix pencil zBV − AV were computed by using the MATLAB function eig. This method is referred to as folklore method 1 herein. Additionally, we computed the eigenpairs of the square matrix pencil using the MATLAB function eig. This method is referred to as folklore method 2. On the other hand, suppose m > n. Then, the eigenpairs of square matrix pencil zV B − V A were computed using the MATLAB function eig. This method is referred to as folklore method 3. Further, the eigenpairs of square matrix pencil were computed using the MATLAB function eig. This method is referred to as folklore method 4. For further comparisons, we used the Rayleigh-Ritz-type contour integral method [25] for the above matrix pencils zBV − AV and zV B − V A, and call them folklore methods 1' and 3', respectively. The matrix resolvent involved in this methods may not exist and was replaced with the pseudoinverse. This pseudoinverse solution was computed similarly to the proposed method. Applying the Rayleigh-Ritz-type contour integral method [25] to the above matrix pencils (3.1) and (3.2) is essentially the same as applying the proposed method to zB − A, when replacing the matrix resolvent involved in the former with the pseudoinverse. Hence, we do not report on these approaches.

Case m < n
We illustrate experimental results for the case m < n. Test matrix pencils zB−A were generated as follows. First, a matrix pencil G(z) = zI η −Λ ∈ C η×η having finite eigenvalues was generated.
Here, Λ ∈ C η×η is a diagonal matrix whose diagonal entries are complex numbers having real and imaginary parts drawn from the standardized normal distribution. Second, a matrix pencil K(z) = z ⊕ p i=1 N ρ i − I ρ ∈ C ρ×ρ having infinite eigenvalues was generated. Here, N ρ i is a shift matrix whose size was randomly chosen such that ρ/2 superdiagonal entries of ⊕ p i=1 N ρ i are one and the remaining entries are zero. Finally, the test matrices were formed, where R 1 and R 2 were generated using the MATLAB function randn for (m, n) = (30, 100), (300, 1000); the products of randomly generated Givens rotation matrices were used to set the nonzero density of sparse matrices A and B to 0.001 for (m, n) = (3000, 10000), (30000, 100000). The block sizes η = ρ = r, q = n − η − ρ, and µ = ν = 0 were set. Table 3.1 gives information on the test matrix pencils. This includes the size of each matrix pencil and its Kronecker blocks, and the values of parameters for the proposed method, including the center and radius of the circle Γ of the prescribed region Ω, number of eigenvalues t in the region, number of columns L of matrix V , order of complex moments M , and number of quadrature points N . Table 3.2 gives the elapsed CPU time of each method in seconds, as well as the maximum relative error (max RERR) of eigenvalues in the prescribed region and the maximum RRN (max RRN) in terms of eigenvalues in the prescribed region and the corresponding eigenvectors. The symbol * indicates the least CPU time among the compared methods for each test problem. Table 3.2 shows that the proposed method was faster than GUPTRI in terms of CPU time and was more robust than other methods. Folklore methods 1, 1', and 2 were faster than the proposed method for the cases (m, n) = (30, 100), (300, 1000); however they did not terminate within 24 hours for the case (m, n) = (30000, 100000). Folklore method 2 computed only one of the three eigenvalues in the prescribed region for the case (m, n) = (3000, 10000). GUPTRI fails to give finite eigenvalues for the case (m, n) = (30000, 100000). This may be due to rounding errors. Note that GUPTRI does not compute the eigenvectors and hence the maximum RNN for the method is not given.
To examine the effect of the number of quadrature points N on the accuracy of the proposed method, we tested the method on different numbers of N = 8, 9, . . ., 48. Here, the values of the other parameters were the same as those given in Table 3.1. Figure 1 shows the maximum RERR and RRN and CPU time in seconds versus the number of quadrature points N on the above test matrix pencil with (m, n) = (3000, 10,000). This figure shows that the maximum RERR and RRN seem to exponentially converge to the machine epsilon regarding N , and the CPU time increased proportionally to N . The Supplementary Material describes the relationship between the quadrature error and the number of quadrature points in a diagonalizable case. We see that the number of quadrature points N = 48 chosen for the above tests is sufficiently large and taken as a safeguard. Similar trends were observed for the other problems in this study. However, it is not clear how to find the least number of quadrature points N required to achieve a specified accuracy in general. Table 3.3 gives the proportion of CPU time in seconds for each step of Algorithm 2.1 with quadrature points N = 48 on the above test matrix pencils. The CPU time for step 2 occupies  30  100  10  10  10  1 2  4  2 48  300  1,000  100  100  100  0.3 3  4  2 48  3,000 10,000 1,000 1,000 1,000 0.

Case m > n
We show experimental results for the case m > n. Test matrix pencils zB − A were generated as in (3.3) and (3.4) with block sizes η = ρ = q, r = m − η − ρ, and µ = ν = 0. As in Table 3.1,  Table 3.4 gives information on the test matrix pencils. Table 3.5 gives the elapsed CPU time of each method, as well as the maximum RERR and RRN, similarly to Table 3.2. The table shows that the proposed method was faster than GUP-TRI in terms of CPU time and was more robust than other methods. Folklore methods 3 and 4 were faster than the proposed method for the cases (m, n) = (100, 30), (1000, 300); however they took more than 24 hours for the case (m, n) = (100000, 30000). Folklore method 3 did not give eigenvalues in the prescribed region for the case (m, n) = (10000, 3000). Folklore method 4 computed only two of the three eigenvalues in the prescribed region for the case (m, n) = (10000, 3000). GUPTRI failed to give finite eigenvalues for (m, n) = (100000, 30000). This may be due to rounding errors.

Case ν > 0
We show experimental results for the case ν > 0, although the proposed method is not theoretically guaranteed to work in this case. Test matrix pencils zB − A ∈ C m×n were generated as

follows:
where Λ, ⊕ p i=1 N ρ i , R 1 , and R 2 were generated as in Section 3.1. We fixed m = 10000, η = ρ = 1000, r = 1, and R = 0.1, and varied ν = 1, 10, 100, and 1000 (n = 2001, 2010, 2100, and 3000, respectively). Then, the number of eigenvalues in the prescribed region was t = 3. The numbers of parameters were set to L = 8, M = 4, and N = 48. Table 3.6 gives the elapsed CPU time of each method and the maximum RERR and RRN, similarly to Table 3.2. The table shows that the proposed method was faster than other methods in terms of the CPU time. Folklore method 4 computed only two of the three eigenvalues in the prescribed region for the case ν = 1000. GUPTRI did not give the three eigenvalues in the prescribed region for ν = 10 and 100.

Conclusions
In this study, a projection method for computing interior eigenvalues and corresponding eigenvectors of square linear matrix pencils was extended to nonsquare cases. The proposed method was successfully demonstrated on matrix pencils with specific Kronecker structures. Numerical experiments showed that this approach is superior to previous ones in terms of efficiency and robustness for large problems. Experimental results conjecture that the proposed method works  100  30  10  10  10  1 2  4  2 48  1,000  300  100  100  100  0.3 3  4  2 48  10,000 3,000 1,000 1,000 1,000 0.1 3 8 4 48 100,000 30,000 10,000 10,000 10,000 0.05 2 8 4 48 m: number of rows of a matrix, n: number of columns of a matrix, η: size of the Kronecker block corresponding to finite eigenvalues, ρ: size of the Kronecker block corresponding to infinite eigenvalues, R: radius, t: number of eigenvalues in the curve Γ, L: number of columns of the pseudo-random matrix V , M : order of complex moments, and N : number of quadrature points. on matrix pencils with left singular blocks and theoretical support for this case was left open. This extension took the Rayleigh-Ritz-type approach and may be applied to other types of projection methods [54,50,8,27]. This study provides direction for tackling further extensions of the projection method to nonlinear nonsquare matrix eigenproblems, cf. [2,3,8,61].