A Globally Convergent Iterative Method for Matrix Sign Function and Its Application for Determining the Eigenvalues of a Matrix Pencil

: In this research article, we propose a new matrix iterative method with a convergence order of five for computing the sign of a complex matrix by examining the different patterns and symmetry of existing methods. Analysis of the convergence of the method was explored on a global scale, and attraction basins were demonstrated. In addition to this, the asymptotic stability of the scheme was explored.Then, an algorithm for determing thegeneralized eigenvalues for the case of regular matrix pencils was investigated using the matrix sign computation. We performed a series of numerical experiments using numerous matrices to confirm the usefulness and superiority of the proposed method.


Introduction
Over the past several years, the efficient computing of matrix functions has emerged as a prominent area of research.Symmetry plays a significant role in the study of matrix functions, particularly in understanding their properties, behavior, and computational methods.Symmetry principles are widely used in physics and engineering to analyze and solve problems involving matrix functions.For instance, in quantum mechanics, symmetric matrices often arise in the context of Hamiltonians, and efficient computation of matrix functions is crucial for simulating quantum systems and predicting their behavior.In order to compute matrix functions, various iterative schemes have been suggested, which include methods for the evaluation of matrix sign function.The scalar sign function defined for z ∈ C such that Re(z) ̸ = 0 which is given as, is extended to form the matrix sign function.Roberts [1] introduced the Cauchy integral representation of the matrix sign function in 1980, given by which can be applied in model reduction and in solving algebraic Riccati and Lyapunov equations.Recall the definition of matrix sign function sign(A) via the Jordan canonical form [2], in which if A = PJP −1 is the Jordan decomposition of A such that J = diag(J 1 , J 2 ), then where A ∈ C n×n has no eigenvalue on the imaginary axis, P is a nonsingular matrix, and p and n − p are sizes of Jordan blocks J 1 , J 2 corresponding to eigenvalues lying in left half plane (C − ) and right half plane (C + ), respectively.
In the pursuit of evaluating sign(A) as a solution of Roberts applied Newton's method (N M) resulting into the following iterative formula: with the initial matrix W 0 = A and convergence of order two.Another alternative is accelerated Newton's method (AN M) that can be achieved via norm scaling parameter in (5) given by, In order to improve the convergence, acceleration and scaling of iteration (5), Kenney and Laub [3,4] suggested a family of matrix iterations by employing the Padé approximant to the function where x = 1 − z 2 and taking into the consideration Let (r, s)-Padé approximant to H(x) with r + s ≥ 1 is denoted as P r,s (x) Q r,s (x) , then converges to ±1 if r ≥ s − 1 and it has convergence rate equal to r + s + 1.So, the (r, s)-Padé approximant iteration (PA r,s ) for computing sign(A) is These iterations were later investigated by Gomilko et al. [5].Other popular iterative approaches for sign(A), such as the Newton-Schulz method (NSM) [6] and Halley's method (HM) [2] are either members of the Padé family (10) or its reciprocal family.
A recent fifth order iterative method [7] (represented as Lotfi's Method (LM)) for sign(A) can be investigated and is given by For more recent literature on matrix sign function, one can go through references [8][9][10].Motivated by constructing efficient iterative methods like the methods discussed above and to overcome on some of the drawbacks of the existing iterative methods, here we aim to propose a globally convergent iterative method for approximation of sign(A).Matrix sign function has several application in scientific computing (see e.g., [11,12]).Another applications of matrix sign function include solving various matrix equations [13], solving stochastic differential equations [14] and many more.
In this article, we present a fifth order iterative scheme for evaluation of sign(A).The outline for the rest of the article is defined here: In Section 2, the construction of new iterative scheme for evaluating matrix sign function has been concluded.In Section 3, basins of attraction guarantees the global convergence of developed method.Also, convergence analysis for the fifth rate of convergence is discussed along with stability analysis.Then the proposed method is extended to determine the generalized eigenvalues of a regular matrix pencil in Section 4. Computational complexity is examined in Section 5. Numerical examples are provided to evaluate the performance and to demonstrate the numerical precision of suggested method in Section 6.The final conclusion is presented in Section 7.

A Novel Iterative Method
In this section, we have developed an iterative method for evaluation of sign(A).The primary goal is to develop the scheme by examining the factorization patterns and symmetry of existing schemes rather than using a non-linear equation solver [13].We aim to design an iterative formula that does not conflict with the Padé family members (9) or its reciprocals.
If we assume the eigenvalues to converge to 1 or −1, and we aim to get a higher convergence order five, then the iterative formula will satisfy the equation If we solve Equation ( 14) for w l+1 , then we get Nevertheless, the formula mentioned above is not a novel one, as it is a member of Padé family (9).So, we make some modifications in (14) so that the formula becomes better in terms of converegnce as follows: provided that −4+w l −4−w l < 1. Solving (16) for w l+1 , we get which contains a polynomial of degree 5 in the numerator and a polynomial of degree 6 in the denominator.From (9), it is clear that Padé family has a unique member with a polynomial of degree 2r + 1 in numerator and a polynomial of degree 2s in denominator for each r, s ≥ 0. So 2r + 1 = 5 and 2s = 6 will result into r = 2 and s = 3, i.e., PA 2,3 and will have order 2 + 3 + 1 = 6 which is different from the order 5 that we are achieving.So, the above formula (17) is not a part of the family (9) or reciprocals of (9) and therefore our proposed method (PM) is as follows: where W 0 is the appropriate initial guess which we discuss in Section 3.

Remark 1.
Convergence can be slow if ∥W l ∥ >> 1 for some l.In that situation, it is advisable to introduce the scaling parameter µ l [15] given as, (determinant scaling). ( to speed up the convergence [16] of the scheme. Hence, we also present the accelerated proposed method (APM) which is accelerated version of (18) given by where lim l→∞ µ l = 1 and lim l→∞ W l = S = sign(A).

Convergence and Stability Analysis
In this section, we discuss the convergence analysis to ensure that the sequence generated from (18) converges to sign(A) for a matrix A ∈ C n×n .We draw the basins of attraction to confirm the global convergence of the proposed scheme and compared with the existing schemes for solving w 2 − 1 = 0 [17]. We The region has grid size of 256 × 256 and has two simple roots −1 and 1.The convergence was determined by the stopping criteria ∥ f (w l )∥ ≤ 10 −2 .The exact location of roots is shown by two white dots within the basins.The dark blue color (left side of basins) represent the convergence region corresponding to the root w = −1 and the light blue color (right side of basins) represent the convergence region corresponding to the root w = 1.
Figure 1 illustrates the basins for N M (5), NSM (11), HM (12), PA r,s (10), and the proposed method PM (18).It is evident from Figure 1b,d that methods NSM and PA 3,1 exhibit local convergence, while the other methods exhibit the global convergence.Moreover, Figure 1e,f illustarte broader and lighter basins of attraction because of higher convergence order of PA 2,2 and PM as compared to the basins of N M and HM as shown in Figure 1a,c.One can notice the symmetry in convergence region of each method.The convergence regions are symmetric about imaginary axis.
After achieving global convergence of the proposed method, it is important to examine certain critical factors related to matrix iterations, including convergence rate and numerical stability.

Convergence Analysis
Theorem 1. Suppose that the initial matrix W 0 = A ∈ C n×n has no eigenvalue on the imaginary axis, then the sequence of matrices {W l } ∞ l=0 generated from ( 18) converges to sign(A) = S.
Proof.Let W = PJP −1 be Jordan decomposition of matrix W ∈ C n×n and consider where R is the rational function linked to (18).Consequently, this will lead to the mapping of an eigenvalue λ of W l to an eigenvalue R(λ) of W l+1 .Also, R will satisfy the properties given as follows: 1.
To achieve this objective, consider A = PΛP −1 as the Jordan decomposition of A that can be modified as [2] where P is a non singular matrix while Jordan blocks C and N are associated to the eigenvalues in C − and C + respectively.Let the principal diagonals of blocks C and N have values λ 1 , λ 2 , . . ., λ p and λ p+1 , λ p+2 , . . ., λ n respectively.Then Therefore, Now, define a sequence D l = P −1 W l P with D 0 = P −1 AP.Then from (18), we get Mathematical induction states that all D l 's for l = 1, 2, 3... will be diagonal if D 0 is diagonal.When D 0 is not diagonal, the proof may be considered in the similar methodology as stated at the end of the proof.Now, (25) can be written for scalar case f (w) = w 2 − 1 = 0 as follows: where d i l = (D l ) i,i and 1 ≤ i ≤ n.Furthermore, sign(λ i ) = s i , with s i = ±1 and therefore from (26), The factor 4s i −d i l 4s i +d i l is bounded for 1 ≤ i ≤ n and does not affect the convergence.
Furthermore, by selecting the right initial matrix W 0 = A and implying that lim l→∞ (d i l ) = s i = sign(λ i ) and hence lim l→∞ (D l ) = sign(Λ).In the case when D 0 is not a diagonal matrix, it is essential that we investigate the relation between the eigenvalues of W i 's which was briefly explained at the starting of the proof.
In a comparable approach, it is clear from ( 29) that λ i 's will converge to s i = ±1, i.e., lim As a result, it is easy to conclude that which finishes the proof.
Theorem 2. Considering the assumptions of Theorem 1, the matrix sequence {W l } ∞ l=0 generated from (18) has convergence rate five to sign(A) = S.
Proof.Since, W l depends on matrix A ∀ l ≥ 0 and the matrix A commutes with S, so W l also commutes with S ∀ l ≥ 0. Taking into consideration the substitution we get, Apply any type of matrix norm to the both sides of (33), This completes the proof.

Stability Issues
Theorem 3. Considering the assumptions of Theorem 1, the sequence of matrix iterations {W l } ∞ l=0 generated from ( 18) is stable asymptotically.
Proof.The iterate obtained from ( 18) is a function of A because of W 0 = A, and hence commutes with A. Let ∆ l be the perturbation generated at l th iterate in (18).Then Based on the results of the error analysis of order one, we may consider that (∆ l ) j ≈ 0, j ≥ 2. As long as ∆ l is small enough, this usual adjustment will work.Also suppose that W l ≈ sign(A) = S for larger l, then we get Here, the identity [18] given by has been applied, where B is any nonsingular matrix and C is any matrix.Also, considering ∆ l+1 = W l+1 − W l+1 , we get At this point, it is possible to deduce that ∆ l+1 is bounded, i.e., which completes the proof.

Extension to Determine Generalized Eigenvalues of a Matrix Pencil
Let A, B ∈ C n×n .Then λ ∈ C is said to be a generalized eigenvalue of a matrix pencil A − λB if there exists a non-zero vector x ∈ C n such that Ax = λBx. (40) Here x will be a generalized eigenvector corresponding to λ and (40) is referred as a generalized eigenvalue problem.This can be determined by finding the zeros of characteristic polynomial det(A − λB).
Complex problem solving in physics was the first scenario for the appearance of the characteristic polynomial and the eigenproblem.For a recent coverage of the topic one can see [19].When more generality is needed, generalized eigenvalue problem is taking the place of the classical approach [20].
Problem (40) arises in several fields of numerical linear algebra.However, the dimensions of matrices A and B are typically quite enormous, hence further complicating the situation.There are a lot of different approaches that claim to be able to solve (40) including the necessary eigenvectors in a smaller subspace or series of matrix pencils [21][22][23] or using matrix decomposition [24-26].Our interest is in employing the matrix sign function which involves evaluating certain matrix sign computations based on the given matrix pencil to compute the required decomposition and hence to solve the problem (40).
To solve the problem (40) for the case of regular matrix pencils using matrix sign function, we consider rank revealing QR (RRQR) decomposition [27] after computing sign of two matrices depending on the given matrix pencil.
Theorem 4. Let A, B ∈ C n×n such that A − λB be a regular matrix pencil and let U 1 and U 2 are defined as where S 1 and S 2 are sign functions given by where λ i 's are eigenvalues of A − λB, α ii and β ii are the diagonal entries of Q * 2 AQ 1 and Q * 2 BQ 1 respectively for i = 1, 2, . . ., n.Moreover, if B is non-invertible, then β ii = 0 for some i which implies that some eigenvalues of A − λB will be infinity.
Let τ denotes tolerance value and max be the maximum iterations.Then, in accordance with the proposed method (18) and Theorem 4, we propose Algorithm 1 to find out the eigenvalues of a regular matrix pencil A − λB.

Computational Complexity
The computational efficiency index (CEI) [28] of an iterative method can be defined as follows: Here p denotes convergence order and C denotes computational cost per iteration.The approximate cost of matrix multiplications in current programming packages is 2.373 [29].For a matrix of size n, we therefore consider that the cost of one matrixmatrix multiplication is n 2373 1000 and cost of one matrix inverse is 2n 3  3 .Additionally, we do not take into account the cost per iteration of norm computation (in AN M (6)), addition and subtraction.Then, CEI's for the methods N M (5), AN M (6), HM (12) and PM (18) to compute sign(A) are given in Table 1.

Method
Convergence Order CEI N M 2 2 1 n 2373/1000 + 2n 3 3 Figure 2 presents a comparison of the efficiencies of several approaches, demonstrating the greater efficiency index of the suggested method.

Numerical Aspects
In this section, we have examined a variety of examples to evaluate our proposed method by comparing it with various existing techniques.For this study, we have assumed the following stopping criteria: where τ stands for the tolerance value and ∥.∥ * will be suitable matrix norm.Also, l ∞ and l 2 norms are considered for real and complex input matrices, respectively [14].
For the sake of comparison, we only consider the methods that exhibit global convergence and we do not include any methods that have local convergence.The methods that are being compared will be N M (5), AN M (6), HM (12), PA 2,2 ((10) for r = s = 2), LM (13), PM (18) and APM (20) with spectral scaling (19).All simulations are done in Mathematica using system specifications "Windows 11 Home Intel(R) Core(TM) i7-1255U CPU running at 1.70 GHz processor having 16.0 RAM (64-bit operating system)".

Example 1.
Here, we perform the eigenvalue clustering of a random real 500 × 500 matrix given as The outcomes are shown in Figure 3 using the stopping criteria (45) with l ∞ norm and τ = 10 −4 .We can see that eigenvalues of W 0 are scattered at the initial stage as shown in Figure 3a, while proceeding for the next iterations W i in Figure 3b-g for i = 1, 2, . . ., 6 respectively using PM (18), the eigenvalues are converging towards 1 or −1.Then in Figure 3h, all eigenvalues converge to 1 or −1 (based on the stopping criteria) implying that the sequence of matrices {W l } ∞ l=0 converges to sign(A).The result of Theorem 1 is therefore confirmed.Example 2. In this example, convergence of many different methods for computing sign(A) has been compared, where A is a random 1000 × 1000 complex square matrix given as follows: The comparison is shown in Figure 4 with τ = 10 −4 and l 2 norm which illustrates the results with regard to the number of iterations and the absolute error and demonstrates the consistent behavior of both PM and APM.The comparison is presented in Figure 5 with τ = 10 −4 and l ∞ norm which illustrates the results with respect to the number of iterations and the absolute error and demonstrates the consistent behavior of both PM and APM for evaluating sign(A).The comparative analysis based on iteration number and CPU timings for matrices of various sizes, are presented in Table 2 and Table 3 respectively using τ = 10 −5 and l 2 -norm.The implementation of new methods (PM) and (APM) result in an improvement in the average iteration number and the amount of time spent by CPU.The outcomes demonstrate that the suggested method (PM) is significantly better as compared to other fifth order methods like PA 2,2 and LM.Now, in Examples 5 and 6, we aim to evaluate generalized eigenvalues of matrix pencils using matrix sign function for which we use Algorithm 1 for different methods by replacing (18) with existing methods in Step 3.
where I 20 is the identity matrix of size 20.
The eigenvalues of given pencil A − λB are λ i = i−1 100 for i = 1, 2, . . .20 and remaining eigenvalues are infinite.We verify these eigenvalues using proposed method and Algorithm 1.We compute the matrices Q Example 6 (Bounded finline dielectric waveguide problem).The generalized eigenvalue problem (40) arises in the finite element analysis for finding the propagating modes of a rectangular waveguide.We consider three cases for the pairs of matrices A and B as BFW62A, BFW62B; BFW398A, BFW398B and BFW782A, BFW782B from the Matrix Market [31].
Table 5 shows the outcomes depending on the iteration numbers l 1 and l 2 for computing S 1 = sign((A − B) −1 (A + B)) and S 2 = sign((A + B)(A − B) −1 ) with corresponding absolute errors ϵ l 1 = ∥W 2 l 1 − I∥ ∞ and ϵ l 2 = ∥W 2 l 2 − I∥ ∞ using stopping criterion (45) with τ = 10 −10 .Also, total computational time in seconds taken for computing S 1 and S 2 has been compared in Figure 6.The outcomes demonstrate that the suggested method is better with respect to errors and timings.

Conclusions
This paper has discussed an iterative approach for evaluation of sign of a matrix.Our proposed iterative method attains a convergence order of five, showcasing its efficiency in evaluation of the sign of complex matrices.Through comprehensive convergence analysis discussed globally and the demonstration of asymptotic stability, we have established the robustness and reliability of the proposed approach.Then an algorithm has been designed for determining the generalized eigenvalues of a regular matrix pencil utilizing matrix sign computation.Furthermore, extensive numerical experiments performed on matrices of different sizes have validated the effectiveness and superiority of our method over existing techniques.

Figure 2 .
Figure 2. Comparative analysis of efficiency indices for different methods.

Figure 5 .
Figure 5. Convergence history of different methods for Example 3.

Example 5 .
Let A − λB be a matrix pencil, where matrices A and B are given by[30]

Figure 6 .
Figure 6.Comparison of CPU time for various methods to compute S 1 and S 2 in Example 6.
2 are RRQR decompositions of U 1 and U 2 such that Q * 2 AQ 1 and Q * 2 BQ 1 are upper triangular, then

Table 2 .
Comparative analysis for the iteration number in Example 4.

Table 3 .
Comparative analysis for the CPU-time (in seconds) in Example 4.

Table 5 .
Comparison of results obtained from various methods in Example 6.