A fast and efficient Newton-type iterative scheme to find the sign of a matrix

1.   Introduction

The author in ^[10] first proposed a matrix function for the sign to be used and applied in algebraic Riccati equations. The sign function of a matrix can be defined using the following Cauchy integral definition:

Assume that $A\in \mathbb{C}^{n\times n}$ is given, and $g$ is a scalar function. Then, we define $g(A)$ as a function of a matrix with the same size as the input matrix. As long as $g$ is given on the spectrum of the input matrix [5,Chapter 5], we can find the following features for $g(A)$:

● If $Y$ commutes with $A$, then $Y$ commutes with $g(A)$,

● The eigenvalues of $g(A)$ are $g(\lambda_i)$, where $\lambda_i$ stand for the eigenvalues of $A$,

● $g(Y A Y^{-1}) = Y g(A) Y^{-1}$,

● $g(A^T) = g(A)^T$,

● $g(A)$ commutes with $A$.

Note that both Jordan decomposition and Schur decomposition with block Parlett recursion can be used under some conditions to compute matrix functions. The matrix sign function can give us the ability to decompose an appropriate matrix into two components whose spectra lie on opposite sides of the axis of imaginary. For more information, see ^[3].

Practically speaking, we are motivated to derive a new Newton-type fast iteration method to find (1.1), which is efficient. That is, due to its higher order and wider convergence radii, fewer iterations will be required. This results in lower computational cost of the whole algorithm since fewer number of stopping terminations should be computed per cycle for the iterative structure. This is useful especially when the entries of the input matrix are changing a bit in a loop and each time, it is needed to find (1.1). The goal of this research is to boost convergence speed, resulting in high computational efficiency.

Here, we concentrate on functions of matrices and especially the sign function. Toward this aim, we first derive an iterative scheme for finding the matrix sign function. We derive a novel numerical procedure to find (1.1) for which the convergence rate is high via the number of matrix-matrix products while only one matrix inversion is needed per full stage.

For this research study, the remaining sections of this article are unfolded as comes next. In Section 2, iteration methods for finding (1.1) are discussed. Section 3 investigates a higher order method, which is under the umbrella of iteration methods to compute the sign of an invertible matrix. The convergence of the scheme and its asymptotical stability are brought forward. To test the correctness, efficiency, and stability, we run numerical tests on various problems in Section 4. Finally, based on the numerical results obtained, the suggested technique is determined to be efficient. Some remarks are furnished in Section 5.

2.   Iteration methods in the literature

A typical procedure for finding iteration methods to calculate matrix functions is via employing a suitable nonlinear solver and selecting an initial value so as to guarantee the (local) convergence. It is here requisite to remind that the results originated from scalar nonlinear solvers do not generalize necessarily to the matrix environment. As an instance, standard convergence conditions given via the derivatives of $g$ at a fixed point in the scalar case do not translate explicitly to similar conditions on the Frechét derivatives in the matrix environment.

The second-order Newton's method (NM) has the following structure to calculate the sign of a square invertible matrix:

where

is initial value. Authors in ^[7] proposed an important family of iterative schemes for finding (1.1) employing the following Padé approximants:

Let us consider that the $(m_1, m_2)$-Padé approximate to $g(\zeta)$ is defined by:

where $P$, $Q$ are polynomials of the approximate degree in the Padé approximation and $m_1+ m_2 \geq 1$. Then, ^[7] showed that the iterative scheme

converges with convergence speed $m_1+m_2+1$ to $\pm 1$.

By considering (2.5), the well-known locally convergent inversion-free Newton-Schulz method

and Halley's iteration scheme

belong to this family. It is noted that Newton's method (2.1) is a member of the reciprocal of (2.5), ^[4].

Authors in the work ^[2] investigated a fourth-order family of methods to find (1.1) with the following structure:

wherein $v$ is a free real parameter.

Two other fourth-order methods from (2.5) having global convergence behavior are given by:

For a recent work on this topic, an interested reader may refer to ^[12].

3.   A novel fourth-order method

3.1.   A Newton-type nonlinear equation solver

Recalling that iteration methods can be constructed for matrix sign function if they can solve the matrix equation below (see, e.g., ^[11]):

wherein $I$ is the unit matrix. In the scalar case, we consider $g(x) = x^2-1 = 0$. To get a better background on nonlinear equation solvers and their recent findings, one may refer to ^[1].

Without loss of generality, consider the scalar case. Now, we propose the following multi-step Newton-type iteration scheme for solving (3.1):

Theorem 3.1. Let $\xi \in D$ be a simple solution (root) of $g:D\subseteq\mathbb{C}\rightarrow \mathbb{C}$, which is a sufficiently smooth nonlinear function. Also assume the initial value $w_{0}$ is close enough to the solution. Then, the sequence generated by (3.2) has at least fourth-order convergence and satisfies the following equation of error:

where $k_j = \frac{g^{(j)}(\xi)}{j!g'(\xi)}$ and $e_m = w_m-\xi$.

Proof. The proof is via Taylor expansions. See e.g., ^[13,14]. Hence, it is omitted. □

3.2.   Extension to matrix sign computation

By employing (3.2) to solve the scalar version of (3.1), we attain the following iterative method in a unique way:

having the initial value (2.2). Similarly, one may obtain the reciprocal version of (3.4) as follows:

3.3.   Convergence

Now, the convergence properties of (3.5) are investigated. Recall that for any $H$ such that $A+H$ is positive definite, we have

as a fixed point of (3.5).

Theorem 3.2. In computing the sign of the matrix $A$ without any eigenvalues on the imaginary axis, consider that $W_0$ is sufficiently close to $S$ or is chosen by (2.2) and commutes with $A$. Then, the scheme (3.5) (or equivalently (3.4)) converges to the sign matrix $S$, and the convergence rate is four.

Proof. Let us use the Jordan block matrix $J$ and decompose $A$ using an invertible matrix $T$ of the same size as follows:

Employing this along with the iterative method yields an iteration scheme similar to the original matrix iterative method but on the eigenvalues from the step $m$ to the step $m+1$ as follows:

where

In general and after some mathematical simplifications, the expression (3.8) shows that the eigenvalues tend to $s_i = \pm1$, that is:

This shows the convergence and states that the eigenvalues tends to $\pm1$ after each iteration, and there would be eigenvalue clustering using (3.5). After studying the convergence, it is requisite to compute the convergence speed. Toward this, it is considered that:

Using (3.11), we can write:

Using (3.12), it is possible to get that:

which shows the convergence rate of four. This completes the proof. The analysis of error for (3.4) might be similarly performed. The proof is complete now. □

We remark that $W_0 = A$ is the best choice as the initial approximation since no eigenvalues are located on the axis of imaginary and could lead to convergence. However, another sufficiently close initial approximation based on $\mu A$ where $\mu$ is an appropriate non-zero parameter can be chosen as long as no eigenvalues are located on the axis of imaginary.

The presented method needs computation of matrix multiplications, as do the other competitors. On the other hand, most of the convergence theory for the presented method is via the calculation of the eigenvalues ^[8] from an iterate to the next one. Nonetheless, we have showed the fourth order is obtained only if the initial approximation is chosen properly.

The cost of an algorithm is important for employing it in practical problems. The convergence rate is not the only factor, and the method is useful only if it can compete with the most efficient existing solvers of the same type. When we compare (3.5) to (2.9)–(2.10), it is observed that all possess per cycle four matrix products and only one matrix inversion. It is now seen and checked that the proposed methods have wider convergence radii.

To check the global convergence of the presented scheme in contrast to the existing solvers, we may draw attraction basins of the iterative methods when solving the nonlinear equation $g(x) = x^2-1$ on the complex domain $[-2, 2]\times[-2, 2]$. In fact, we divide the domain into a refined mesh and test to what root, each of the mesh points converge. Results of comparisons are brought forward in Figures 1 and 2 by employing the stopping condition

The results are shaded via the number of iterates required to reach the convergence. They too show that (3.4) and (3.5), have wider convergence radii than their competitors of the same order from (2.5).

Here, although the iterative methods (2.1), and (2.9) are globally convergent, the basins for (3.4) and (3.5) contain lighter areas which indicate wider convergence radii.

The larger attraction basins given for the proposed solvers are mainly due to the appropriate selection of the parameters furnished in the second step of the nonlinear equation solver (3.2). In fact, the parameters are selected as limiting values, in which the convergence order of the solver improves to higher order, but it does not select the exact limiting value since it costs further matrix products. In this way, by keeping the number of matrix products fixed, we can derive an iterative method for finding the sign of a nonsingular matrix with larger attraction basins than the competitors from the Padé scheme.

3.4.   Stability

The stability for (3.4) or similarly (3.5) are presented now. Theorem 3 follows from ^[6] on the stability of "pure matrix iterations."

Theorem 3.3. Using (3.5) and having a nonsingular matrix $A$, $\{W_m\}_{m = 0}^{\infty}$ with $W_0 = A$ is stable asymptotically.

Proof. Assume that $\delta_{m}$ is a numerical perturbation which is produced in the $m^{th}$ iteration. So, we can define

In the rest of the proof, it is assumed that $(\delta_{m})^{i} \approx 0, \ i \geq 2$, is used, which is valid based on first order error analysis. This is rational if $\delta_{m}$ is sufficiently small. Now, we obtain

For sufficiently large $m$, viz., at the phase of convergence, it is assumed that

where the following facts are employed (for any the matrix $E$ and any invertible matrix $H$) [15,page 188]:

We used

and

to get that

By further simplifications and using $\delta_{m+1} = \widetilde{W}_{m+1}-W_{m+1}$, we can find

This can lead to the point that the method at the stage $m+1$ is bounded, i.e., we have:

Therefore, the sequence $\{W_{m}\}^{\infty}_{m = 0}$ produced via (3.4) has asymptotical stability. The proof ends now. □

4.   Numerical aspects

Here, the entire implementations were performed by Mathematica 12.0 in machine precision [16,Chapter 1.4]. The iteration methods (2.1), (2.7), (2.9), (3.4) and (3.5) (denoted by NM, HM, Padé, PM1, and PM2, respectively) were compared based on the efficiency and stability on complex test matrices of various sizes. We compared the methods using the following condition:

It is recalled that for (2.1), the process can be accelerated, via calculating an extra parameter per iterate and replacing $W_m$ with $\mu_m W_m$ as follows:

Further acceleration and optimizing of (2.7) were studied in ^[9]. Similar investigations can be performed for the new derived iterative method in order to accelerate the scheme. Anyhow, none of such accelerations will be performed in this section for any of the compared methods in order to have a fair comparison.

Experiment 4.1. Various complex matrices of different sizes are employed here to compare the efficiency of the existing solvers based on both the number iterations as well as the CPU time in seconds. Ten matrices are produced as follows:

${\mathtt{SeedRandom[123];}}$

${\mathtt{number = 10;}}$

${\mathtt{Table[A[l] = RandomComplex[{-5 - 5 I, 5 + 5 I}, {100 l, 100 l}]; , {l, number}];}}$

${\mathtt{tolerance = 10^-4;}}$

Results are provided in Tables 1 and 2 showing that PM1 has the best performance against its competitors. Here, it is important to note that the convergence of the discussed iteration methods relies on an appropriate choice of initial matrices. It is observed PM1 beats all other competitors by converging to the sign as quickly as possible.

This part of the work is ended by stating that after drawing attraction basins of the new methods, such as PM1, in Figures 3 and 4 for higher degree polynomials, it can be stated that like all the other solvers, when one wants to extend NM, HM, Padé or PM1 to find the matrix sector function, then similar conditions can be imposed since they have local convergence (like the Newton's solver) to find the matrix sector function, but PM1 or PM2 enjoy higher rates of convergence, which will reduce the CPU time in practice.

5.   Concluding remarks

This paper has discussed the existing iteration methods to find the matrix sign functions. The most challenging ones arise from (2.5) with arbitrary order of convergence. Such methods are called optimal in some sense in the literature. However, in this article, we have proposed a novel Newton-type solver possessing fourth convergence speed with global convergence behavior to calculate the sign of a matrix. The new method competes with corresponding methods from the Padé family of methods.

A clear application of the sign matrix is in finding the matrix geometric mean of two Hermitian positive definite matrices, which is required mainly in the context of solving a special class of nonlinear matrix equations. For the proposed solvers in this work, (3.4) and (3.5), parallel implementations are simple, employing the available optimized programs for matrix-matrix products. Computational tests in Section 4 have been performed to show the effectiveness of the new iterative methods for a variety of complex matrices of different sizes. The mean CPU times for (3.4) and (3.5) showed better performance than the others.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research work was funded by Institutional Fund Projects under grant no. (IFPIP: 581-130-1443). The authors gratefully acknowledge technical and financial support provided by the Ministry of Education and King Abdulaziz University, DSR, Jeddah, Saudi Arabia.

Conflict of interest

The writers state that there is no conflict of interests in this paper.