Convergence theory of efficient parametric iterative methods for solving the Yang-Baxter-like matrix equation

Abstract

Purpose

In this study, we present a novel parametric iterative method for computing the polar decomposition and determining the matrix sign function.

Design/methodology/approach

This method demonstrates exceptional efficiency, requiring only two matrix-by-matrix multiplications and one matrix inversion per iteration. Additionally, we establish that the convergence order of the proposed method is three and four, and confirm that it is asymptotically stable.

Findings

In conclusion, we extend the iterative method to solve the Yang-Baxter-like matrix equation. The efficiency indices of the proposed methods are shown to be superior compared to previous approaches.

Originality/value

The efficiency and accuracy of our proposed methods are demonstrated through various high-dimensional numerical examples, highlighting their superiority over established methods.

Full text

Translate

Turn on search term navigation

1. Introduction

In recent years, much attention has been given in the literature to solve matrix equations (Erfanifar and Hajarian, 2023, 2024a, b, d; Erfanifar et al., 2023). Consider the quadratic matrix equation(1.1) $X A X = A X A,$ where $A, X \in C^{n \times n}$ , is called the nonlinear Yang-Baxter-like matrix equation (YBLME). The term “Yang-Baxter equation” actually stems from the work of statistical physicists C.N. Yang and R.J. Baxter in the 1970s. The equation that they introduced is a fundamental concept in theoretical physics, especially in the context of integrable systems, statistical mechanics, and quantum groups. The Yang-Baxter equation in mathematics and physics has many applications, including in statistical mechanics, quantum physics, knot theory, and quantum computation (Tian, 2016; Dong and Ding, 2016; Mansour et al., 2017).

The YBLME possesses significant relevance in the realms of completely integrable quantum and classical systems, as well as in exactly solvable models within statistical physics. In recent years, this equation has garnered substantial attention, undergoing intensive study due to its implications and applications. Additionally, its profound connections with areas of mathematics such as group theory and algebraic geometry have become increasingly evident (Kumar et al., 2022; Baxter, 1972).

Integrable quantum theory is a broad field with applications spanning condensed matter and statistical physics, establishing numerous connections to various branches of mathematics. It presents a distinctive avenue to explore the intricacies of quantum field theory, enabling a more profound understanding of relativistic particle theories by treating them as scaling limits of quantum chains and classical lattice systems. This approach offers deeper insights beyond conventional methods. Acquiring knowledge in integrable quantum theory equips theoretical physicists with valuable tools and contributes to a broader comprehension of subjects that are of interest to mathematicians.

Over the last 3 decades, the study of YBLME has extended beyond statistical physics and garnered interest from mathematicians due to its profound connections to various branches of pure mathematics. Notably, recent advancements in knot theory have revealed the formulation of knot invariants using statistical mechanics methodologies. While this novel approach to knot invariants can be explored in diverse contexts, a particularly insightful perspective is provided by statistical mechanics, with YBLME at its core. YBLME’s significance manifests in knot invariants, where the cabling approach contributes to constructing new explicit solutions of YBLME derived from existing ones. Solutions to these equations have been uncovered as commuting objects arising from second tensor products of deformations of Lie algebras, specifically quantum groups, in their standard representations. Additionally, YBLME plays a pivotal role in analyzing integrable systems, theoretical physics, quantum and statistical mechanics, knot theory, and the theory of quantum groups. Moreover, the theory of integrable Hamiltonian systems greatly benefit from solutions to the one-parameter form of the Yang-Baxter equation, as coefficients from the power series expansion of such solutions facilitate the computation of integrals of motion. YBLME finds applications in various domains, including differential equations, braid groups, and numerous other disciplines (Bayoumi, 2023; Roberts, 1980; Alhayani and Abdallah, 2021; Ballester-Bolinches et al., 2021; Castelli et al., 2021; Bachiller and Cedó, 2014).

The YBLME possesses at least two trivial solutions: X = 0 and X = A. The main focus lies in discovering nontrivial solutions to (1.1). This pursuit has proven to be challenging due to the intricate nature of characterizing the entire set of solutions for a general matrix A (Ding and Rhee, 2012, 2013).

In this paper, we introduce innovative parametric iterative methods aimed at determining the polar decomposition of a matrix. Additionally, we present a new approach for computing the sign function of a matrix. Finally, we extend the iterative methods to compute solutions for the YBLME.

The sign function of a nonsingular square matrix holds significant mathematical value and finds applications across various mathematical domains (Soheili et al., 2015; Nakatsukasa and Freund, 2016; Gomilko et al., 2012; Erfanifar and Hajarian, 2024c).

Suppose that A has the following Jordan canonical as follows:(1.2) $A = W Λ W^{- 1},$ where Λ = diag(Λ₁, Λ₂), the eigenvalues of $Λ_{1} \in C^{p \times p}$ lie in open left half-plane and those of $Λ_{2} \in C^{q \times q}$ lie in open right half-plane, then from Roberts (1980), we have $S = sign (A) = W (\begin{matrix} - I_{p} & 0 \\ 0 & I_{q} \end{matrix}) W^{- 1} .$ The matrix S is the matrix sign of A having no pure imaginary eigenvalue.

In the following, some properties of sign(A) are collected.

Theorem 1.1. Higham (2008) Let $A \in C^{n \times n}$ have no eigenvalues on the axis of imaginary and suppose that S = sign(A). Then

The equation g(x) = 1 has unique solution in the interval [0, 1) as follows:

T = \frac{3 - 9 b - \sqrt{9 - 22 b - 15 b^{2}}}{2 (3 b - 1)} .

The function g is increasing on (0, T). So, if x_k ∈ (0, T), then there is an index k₀ ≥ k such that $x_{k_{0}} > T$ , in fact $x_{k_{0} + 1} \in (T, 1)$ .
If x_k ∈ (T, 1), then we have x_k+1 ∈ (1, λ).
If x_k+1 ∈ (1, λ), then the sequence ${x_{k + m}}_{m = 1}^{\infty}$ converges to x = 1.

According to the mentioned notes, it is concluded that the sequence x_k+1 = g(x_k) is convergent to x = 1, and since g′(1) = g″(1) = 0, the convergence order is three.

Since $\frac{1}{5} \leq b$ , the critical point of the function g is obtained as follows

g^{'} (x) = 0 \Rightarrow x = 1 .

We know that x_critical = 1, g(x_critical) = 1, and the function g is strictly increasing in [0, 1]. So, if x₀ ∈ (0, 1), then we obtain x₁ = g(x₀) ≤ 1. Therefore, after some iterations, the sequence

{x_{k}}_{k \geq 0}

falls in (0, 1) and to x = 1 is monotonically convergent. The proof for b ≤ −1 is similar to this process.

The proof (iii) is easily obtained similar to previous parts. □

In the following, we propose a method for finding the factor U of the polar decomposition of a matrix. By using the sequence (2.5), we can obtain(3.2) $U_{k + 1} = U_{k} \{\frac{3}{8} {(1 + b)}^{3} {(b I + U_{k}^{*} U_{k})}^{- 1} - \frac{3}{8} (b^{2} + 3 b - 2) I + \frac{1}{8} (3 b - 1) (U_{k}^{*} U_{k})\} .$

Regarding the iterative method (3.2), the following results are obtained.

Theorem 3.2. Assume that $A \in C^{m \times n}$ . Then, the matrix sequence ${U_{k}}_{k \geq 0}$ generated by (3.2) with U₀ = αA converges to the factor U.

Proof. By using the singular value decomposition of A in the form(3.3) $A = P (\begin{matrix} Σ \\ 0 \end{matrix}) Q^{*},$ where $P \in C^{m \times m}$ and $Q \in C^{n \times n}$ are unitary matrices and $Σ = diag (σ_{1}, σ_{2}, \dots, σ_{r}), σ_{1} \geq σ_{2} \geq \dots \geq σ_{r} \geq 0 .$ Note that if Range(U*) = Range(H), then U and H are unique. Hence, one can use a result from Ben-Israel and Greville (2003) about the uniqueness U and H as follows: $U = U_{r} V_{r}^{*}, H = U^{*} A,$ where U_r which is made of the first r columns of P and V_r is made of the first r columns of Q. Define(3.4) $D_{k} = P^{*} U_{k} Q .$ Subsequently, by substituting (3.4) in (3.2), we have(3.5) $\begin{matrix} D_{0} = α Σ, \\ D_{k + 1} = D_{k} \{\frac{3}{8} {(1 + b)}^{3} {(b I + D_{k}^{2})}^{- 1} - \frac{3}{8} (b^{2} + 3 b - 2) I + \frac{1}{8} (3 b - 1) D_{k}^{2}\} . \end{matrix}$ Since $D_{0} \in R^{m \times n}$ is diagonal with nonnegative entries, therefore, it is defined by induction on k $D_{k} = (\begin{matrix} diag (d_{i}^{(k)}) & 0 \\ 0 & 0 \end{matrix}), d_{i}^{(k)} > 0, 1 \leq i \leq r,$ where zeros may not exist. Accordingly, (3.5) represents r iterations as follows:(3.6) $\begin{matrix} d_{i}^{(0)} = α σ_{i}, \\ d_{i}^{(k + 1)} = d_{i}^{(k)} \{\frac{3}{8} {(1 + b)}^{3} {(b + d_{i}^{{(k)}^{2}})}^{- 1} - \frac{3}{8} (b^{2} + 3 b - 2) I + \frac{1}{8} (3 b - 1) d_{i}^{{(k)}^{2}}\}, \end{matrix} 1 \leq i \leq r .$

With simple manipulations, and using (3.6), the following relationship is obtained:(3.7) $\frac{d_{i}^{(k + 1)} - 1}{d_{i}^{(k + 1)} + 1} = \frac{- 8 b + d_{i}^{(k)} (3 + 15 b) - 8 d_{i}^{{(k)}^{2}} + d_{i}^{{(k)}^{3}} (6 - 10 b) + d_{i}^{{(k)}^{5}} (3 b - 1)}{8 b + d_{i}^{(k)} (3 + 15 b) + 8 d_{i}^{{(k)}^{2}} + d_{i}^{{(k)}^{3}} (6 - 10 b) + d_{i}^{{(k)}^{5}} (3 b - 1)} .$ Since σ_i is positive, and (3.7) holds for every i, therefore, we have $\lim_{k \to \infty} |\frac{d_{i}^{(k + 1)} - 1}{d_{i}^{(k + 1)} + 1}| = 1,$ and we can obtain for every i: $D_{k} \to (\begin{matrix} I_{r} & 0 \\ 0 & 0 \end{matrix}) .$ Thus, as $\lim_{k \to \infty} U_{k} = U_{r} V_{r}^{*} = U$ . The proof is complete. □

In the following, we prove that the method (3.2) is third-order convergent.

Theorem 3.3. Let $A \in C^{m \times n}$ , then the method (3.2) has third order for finding the unitary matrix U for $b \neq \frac{1}{5}$ .

Proof. By applying Theorem 3.2, the method (3.2) transforms the singular values of U_k according to(3.8) $σ_{i}^{(k + 1)} = σ_{i}^{(k)} \{\frac{3}{8} {(1 + b)}^{3} {(b + σ_{i}^{{(k)}^{2}})}^{- 1} - \frac{3}{8} (b^{2} + 3 b - 2) + \frac{1}{8} (3 b - 1) σ_{i}^{{(k)}^{2}}\}, 1 \leq i \leq r .$

From (3.8), we indicate that convergence of $σ_{i}^{(k)}$ is third order for k ≥ 1:(3.9) $\frac{σ_{i}^{(k + 1)} - 1}{σ_{i}^{(k + 1)} + 1} = \frac{{(σ_{i}^{(k)} - 1)}^{3} (- σ_{i}^{(k)} (3 + σ_{i}^{(k)}) + b (8 + 3 σ_{i}^{(k)} (3 + σ_{i}^{(k)})))}{{(σ_{i}^{(k)} + 1)}^{3} (- σ_{i}^{(k)} (- 3 + σ_{i}^{(k)}) + b (8 + 3 σ_{i}^{(k)} (- 3 + σ_{i}^{(k)})))} .$ Then, we have(3.10) $|\frac{σ_{i}^{(k + 1)} - 1}{σ_{i}^{(k + 1)} + 1}| \leq {|\frac{σ_{i}^{(k)} - 1}{σ_{i}^{(k)} + 1}|}^{3} |\frac{(- σ_{i}^{(k)} (3 + σ_{i}^{(k)}) + b (8 + 3 σ_{i}^{(k)} (3 + σ_{i}^{(k)})))}{(- σ_{i}^{(k)} (- 3 + σ_{i}^{(k)}) + b (8 + 3 σ_{i}^{(k)} (- 3 + σ_{i}^{(k)})))}| .$ Thus the method (3.2) has third-order convergence.□

In the following, we propose two third-order convergence iterative methods. For example, if $b = \frac{1}{20}$ , then we have(3.11) $\begin{matrix} V_{k} = U_{k}^{*} U_{k}, \\ U_{k + 1} = U_{k} (\frac{27783}{64000} {(\frac{1}{20} I + V_{k})}^{- 1} + \frac{2217}{3200} I - \frac{17}{160} V_{k}), & k \geq 0, \end{matrix}$ and if $b = \frac{1}{30}$ , then we obtain(3.12) $\begin{matrix} V_{k} = U_{k}^{*} U_{k}, \\ U_{k + 1} = U_{k} (\frac{29791}{72000} {(\frac{1}{30} I + V_{k})}^{- 1} + \frac{1709}{2400} I - \frac{9}{80} V_{k}), & k \geq 0 . \end{matrix}$

Corollary 3.1. The novel methods (3.11) and (3.12) have third order to obtain the matrix U.

Theorem 3.4. Let $b = \frac{1}{5}$ , then the method (3.2) has fourth order to obtain the matrix U.

Proof. According to (3.9), if $b = \frac{1}{5}$ , then we obtain(3.13) $\frac{σ_{i}^{(k + 1)} - 1}{σ_{i}^{(k + 1)} + 1} = \frac{{(σ_{i}^{(k)} - 1)}^{4} (4 + σ_{i}^{(k)})}{{(σ_{i}^{(k)} + 1)}^{4} (- 4 + σ_{i}^{(k)})} .$ Therefore,(3.14) $|\frac{σ_{i}^{(k + 1)} - 1}{σ_{i}^{(k + 1)} + 1}| \leq {|\frac{σ_{i}^{(k)} - 1}{σ_{i}^{(k)} + 1}|}^{4} |\frac{(4 + σ_{i}^{(k)})}{(- 4 + σ_{i}^{(k)})}| .$

Finally, the order of convergence method (3.2) for $b = \frac{1}{5}$ is four.□

Here is the provided method for finding the unitary matrix A:(3.15) $U_{k + 1} = U_{k} (\frac{81}{125} {(\frac{1}{5} I + U_{k}^{*} U_{k})}^{- 1} + \frac{51}{100} I - \frac{1}{20} U_{k}^{*} U_{k}), k \geq 0,$ or equivalently:(3.16) $\begin{matrix} V_{k} = U_{k}^{*} U_{k}, \\ U_{k + 1} = U_{k} (\frac{81}{125} {(\frac{1}{5} I + V_{k})}^{- 1} + \frac{51}{100} I - \frac{1}{20} V_{k}), & k \geq 0 . \end{matrix}$

Note that the methods (3.16), (3.11) and (3.12) are denoted by R8, R9 and R10, respectively.

3.1 Extension of the methods for matrix sign

In the following, we explore an application of the iterative method (3.2) to compute the matrix sign function. A connection exists between the polar decomposition and the matrix sign function. Table 2 showcases established methods utilized for determining the matrix sign function.

Now, we develop the proposed methods in this study to find the matrix sign function. In fact, the iterative methods (3.16), (3.11) and (3.12) are(3.17) $Y_{k + 1} = Y_{k} (\frac{81}{125} {(\frac{1}{5} I + Y_{k}^{2})}^{- 1} + \frac{51}{100} I - \frac{1}{20} Y_{k}^{2}), k \geq 0,$

(3.18) $Y_{k + 1} = Y_{k} (\frac{27783}{64000} {(\frac{1}{20} I + Y_{k}^{2})}^{- 1} + \frac{2217}{3200} I - \frac{17}{160} Y_{k}^{2}), k \geq 0,$ and(3.19) $Y_{k + 1} = Y_{k} (\frac{29791}{72000} {(\frac{1}{30} I + Y_{k}^{2})}^{- 1} + \frac{1709}{2400} I - \frac{9}{80} Y_{k}^{2}), k \geq 0,$ respectively.

Note that the methods (3.17), (3.18) and (3.19) are denoted by E6, E7 and E8, respectively.

Remark 3.1. The simple definitions for each consisted of the following (Birkhoff, 1927)

An equilibrium point is said to be stable if for some initial value close to the equilibrium point, the solution will eventually stay close to the equilibrium point.
An equilibrium point is said to be asymptotically stable if for some initial value close to the equilibrium point, the solution will converge to the equilibrium point.

Theorem 3.5. The sequence ${Y_{k}}_{k = 0}^{k = \infty}$ generated by (3.17) is asymptotically stable.

Proof. Let Δ_k be a numerical perturbation introduced at the k-th iterate of (3.17). Then, we have ${\hat{Y}}_{k} = Y_{k} + Δ_{k} .$ Now, according to Soleymani et al. (2015), for any nonsingular matrix E and the matrix F, we can get ${(E + F)}^{- 1} \approx E^{- 1} - E^{- 1} F E^{- 1},$ and assuming Y_k ≈ sign(A) = S for enough large k, S² = I, S⁻¹ = S, and also ${(Δ_{k})}^{i} \approx 0$ , i ≥ 2, we have $\begin{matrix} {\hat{Y}}_{k + 1} & = \frac{81}{125} {\hat{Y}}_{k} {(\frac{1}{5} I + {\hat{Y}}_{k}^{2})}^{- 1} + \frac{51}{100} {\hat{Y}}_{k} - \frac{1}{20} {\hat{Y}}_{k}^{3} \\ = \frac{81}{125} (Y_{k} + Δ_{k}) {(\frac{1}{5} I + {(Y_{k} + Δ_{k})}^{2})}^{- 1} + \frac{51}{100} (Y_{k} + Δ_{k}) - \frac{1}{20} {(Y_{k} + Δ_{k})}^{3} \\ \approx \frac{81}{125} (S + Δ_{k}) {(\frac{6}{5} I + 2 S Δ_{k})}^{- 1} + \frac{46}{100} (S + Δ_{k}) \\ \approx (S - \frac{4}{5} Δ_{k} - \frac{9}{5} Δ_{k} S Δ_{k}) . \end{matrix}$ since $Δ_{k + 1} = {\hat{Y}}_{k + 1} - Y_{k + 1} = {\hat{Y}}_{k + 1} - S,$ therefore, we can write $Δ_{k + 1} \approx - \frac{4}{5} Δ_{k} - \frac{9}{5} Δ_{k} S Δ_{k} .$ now, we can conclude that the perturbation at the iterate k + 1 is bounded, and we have $‖ Δ_{k + 1} ‖ \approx \frac{1}{5^{k + 1}} ‖ 4 Δ_{0} + 9 Δ_{0} S Δ_{0} ‖ .$ Then, the iterative method (3.17) is asymptotically stable. This ends the proof. □

Corollary 3.2. The sequences ${Y_{k}}_{k = 0}^{k = \infty}$ generated by (3.18) and (3.19) are asymptotically stable.

Theorem 3.6. Suppose that $A \in C^{n \times n}$ has no pure imaginary eigenvalues. Then, the iterative method (3.17) converges to the matrix sign function S, by choosing Y₀ = αA.

Proof. Since all matrices have a Jordan canonical form, like, A = WΛW⁻¹, where Λ includes Jordan blocks and W is a nonsingular matrix. Therefore, let A have a Jordan canonical form arranged as $W^{- 1} A W = Λ = (\begin{matrix} C & 0 \\ 0 & N \end{matrix}),$ here C and N are square Jordan blocks corresponding to eigenvalues lying in $C^{-}$ and $C^{+}$ , where $C^{-}$ and $C^{+}$ are open left-half complex plane and open right-half complex plane, respectively.

As we know A is invertible, according to Higham (2008), we have(3.20) $sign (A) = W (\begin{matrix} - I_{q} & 0 \\ 0 & I_{n - q} \end{matrix}) W^{- 1} .$ Then we can get(3.21) $\begin{matrix} sign (Λ) & = sign (W^{- 1} A W) = W^{- 1} sign (A) W \\ = (\begin{matrix} sign (μ_{1}) \\ ⋱ \\ sign (μ_{q}) \\ sign (μ_{q + 1}) \\ ⋱ \\ sign (μ_{n}) \end{matrix}), \end{matrix}$ where μ₁, μ₂, …, μ_q and μ_q+1, μ_q+2, …, μ_n eigenvalues lying on the main diagonals of C and N, respectively.

In the following, we define H_k = W⁻¹Y_kW, and from the method (3.17), we obtain(3.22) $H_{k + 1} = H_{k} (\frac{81}{125} {(\frac{1}{5} I + H_{k}^{2})}^{- 1} + \frac{51}{100} I - \frac{1}{20} H_{k}^{2}) .$ note that if H₀ is a diagonal matrix, therefore all the other H_k are diagonal. Now it is enough to prove that from (3.22), the iterative method (3.17) converges to sign(A).

We can write (3.22) to solve h(y) = y² − 1 = 0, given by(3.23) $δ_{k + 1}^{i} = δ_{k}^{i} (\frac{81}{125 (\frac{1}{5} I + δ_{k}^{i^{2}})} + \frac{51}{100} I - \frac{1}{20} δ_{k}^{i^{2}}),$ where $δ_{k}^{i} = {(H_{k})}_{i}$ , 1 ≤ i ≤ n. By using (3.22) and (3.23), it suffices to prove that $δ_{k}^{i}$ converges to sign(μ_i) = s_i = ±1. $\frac{δ_{k + 1}^{i} - 1}{δ_{k + 1}^{i} + 1} = \frac{{(δ_{k + 1}^{i} - 1)}^{4} (4 + δ_{k}^{i})}{{(δ_{k}^{i} + 1)}^{4} (- 4 + δ_{k}^{i})} .$ Now since $δ_{0}^{i} = | μ_{i} | > 0$ , we obtain $\lim_{k \to \infty} |\frac{δ_{k + 1}^{i} - 1}{δ_{k + 1}^{i} + 1}| \leq {|\frac{δ_{k}^{i} - 1}{δ_{k}^{i} + 1}|}^{4} |\frac{(4 + δ_{k}^{i})}{(- 4 + δ_{k}^{i})}| = 0,$ so $\lim_{k \to \infty} | δ_{k}^{i} | = 1 = | sign (μ_{i}) | .$ As a result lim_k→∞H_k = sign(Λ), and from H_k = W⁻¹Y_kW, we can get $\lim_{k \to \infty} Y_{k} = W (\lim_{k \to \infty} H_{k}) W^{- 1} = W sign (Λ) W^{- 1} = sign (A) .$

Finally, the proof is completed.□

Corollary 3.3. Suppose that $A \in C^{n \times n}$ has no pure imaginary eigenvalues. Then, the iterative methods (3.18) and (3.19) converge to the matrix sign function S.

In the following, we give some theorems for the relationship between the matrix sign function and the solutions of YBLME.

Theorem 3.7. Soleymani and Kumar (2019) Let A be the square matrix and S be its matrix sign function, then $\frac{A (I + S)}{2}$ and $\frac{A (I - S)}{2}$ are solutions of (1.1).

Theorem 3.8. Soleymani and Kumar (2019) Suppose that A is the square matrix, and $β \in C$ , if

(1)S₁ is the sign function of G = A + βI, then $\frac{A (I + S_{1})}{2}$ and $\frac{A (I - S_{1})}{2}$ are solutions of (1.1).
(2)S₂ is the sign function of G = βA, then $\frac{A (I + S_{2})}{2}$ and $\frac{A (I - S_{2})}{2}$ are solutions of (1.1).

Algorithm 1 proposes the method (3.19) to solve YBLME.

Algorithm 1. The method (3.19) to compute the solution of YBLME.

4. Efficiency index

The concept of efficiency stands as a cornerstone in numerical methods. Enhancing efficiency, often synonymous with effectiveness in engineering contexts, involves minimizing resource consumption in an operation. As high-speed methods may incur substantial costs, establishing an efficiency index for numerical techniques becomes crucial. Hence, prioritizing a theoretical discourse on practical implications becomes pivotal to achieve faster implementation times in problem-solving endeavors.

Note that to obtain a fair comparison, if the cost of m is β, then the cost of c would be more than 2β.

Now, the approximate value EI for the proposed methods would be in Table 1 wherein s_i, is all the number iterations required for the convergence of the methods, respectively in the same environment. For example, if s_i for method R1 is s, then for method R4 is $\frac{2}{3} s$ . Table 3 shows s, m, c and EI of the proposed methods. Figure 2 shows EI = 1 + ɛ of each method for s.

EI supports the superiority of the new method for solving some numerical examples in the next section.

5. Numerical examples

We compare the iterative methods for several examples. All programming runs with MATLAB software and the criterion for stopping is: $‖ U_{k + 1} - U_{k} ‖_{1} < 1 0^{- 10} ‖ U_{k + 1} ‖_{1} .$ Furthermore, we set $α = \frac{1}{\sqrt{‖ A ‖_{1} ‖ A ‖_{\infty} + 1}}$ .

Example 5.1. The performance evaluation of various iterative methods employed to compute the polar decomposition of matrices is outlined in Table 4. The obtained outcomes are visually depicted in Figures 3–5 and detailed in Tables 5–7.

Example 5.2. The performance analysis of various iterative methods utilized to compute the sign function of matrices is presented in Table 8. The obtained results are visually represented in Figures 6–8 and detailed in Tables 9–11.

We investigate square random matrices of various dimensions, and the outcomes are visually illustrated in Figure 9.

Figures 5–8 along with Tables 6–11 illustrate the iteration counts and computational cost of each method. Notably, the iterative methods (3.16), (3.11), and (3.12) for determining the polar decomposition of matrices, as well as the iterative methods (3.17), (3.18), and (3.19) for computing the sign function of matrices, outperform the others significantly. Tables 6 through 11 exhibit ‖AXA − XAX‖, thereby showcasing that the solution of YBLME is efficiently achieved at an optimal computational expense.

To conclude, we investigated several random real matrices and random complex matrices with varying dimensions. The obtained results are visually depicted in Figure 9.

6. Conclusions

The solution to the quadratic matrix equation $X A X = A X A .$ has been instrumental in solving numerous problems in knot theory, differential equations, braid groups, and quantum groups. One approach to solving this equation involves utilizing the sign function. In this study, we extensively explored new iterative methods for computing the matrix polar decomposition and sign function. The investigation demonstrated that the proposed methods exhibits global convergence of orders three and four. Through comprehensive testing with various examples, it was evident that the proposed methods showcased significant superiority over other well-established methods.

The authors wish to express their gratitude to the editor and anonymous reviewers for helpful remarks and suggestions.

Conflicts of interest: The authors have no conflict of interest to declare.

Data availability: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Figure 1

Plots of the line y = x and the function g for different values of b

[Figure omitted. See PDF]

Figure 2

EI of mentioned methods for s

[Figure omitted. See PDF]

Figure 3

The residual of the results of the several iterative methods for matrix A with n = 20, 30

[Figure omitted. See PDF]

Figure 4

The residual of the results of the several iterative methods for matrix B with n = 20, 30

[Figure omitted. See PDF]

Figure 5

The residual of the results of the several iterative methods for matrix C with n = 20, 30

[Figure omitted. See PDF]

Figure 6

The residual of the results of the several iterative methods for matrix D with n = 20, 30

[Figure omitted. See PDF]

Figure 7

The residual of the results of the several iterative methods for matrix E with n = 20, 30

[Figure omitted. See PDF]

Figure 8

The residual of the results of the several iterative methods for matrix F with n = 20, 30

[Figure omitted. See PDF]

Figure 9

Numerical results for random real and random complex matrices with different dimensions

[Figure omitted. See PDF]

Table 1

Methods and their computational efficiency

Method	Form	s_i	m	c	p	EI
R1 (Kovarik, 1970)	$V_{k} = U_{k}^{*} U_{k}, K_{k} = (I - V_{k}) {(I + V_{k})}^{- 1}, U_{k + 1} = U_{k} (I + K_{k})$	s	3	1	2	$2^{\frac{1}{4 s}}$
R2 (Higham, 1986)	$U_{k + 1} = \frac{1}{2} (U_{k} + U_{k}^{† *})$	s	0	1	2	$2^{\frac{1}{s}}$
R3 (Gander, 1990)	$U_{k + 1} = \frac{1}{2} U_{k} (I + V_{k}^{- 1})$	s	2	1	2	$2^{\frac{1}{4 s}}$
R4 (Gander, 1985)	$U_{k + 1} = U_{k} (V_{k} + 3 I) {(3 V_{k} + I)}^{- 1}$	$\frac{2}{3} s$	3	1	3	$3^{\frac{3}{8 s}}$
R5 (Haghani and Soleymani, 2015)	$Z_{k} = V_{k}^{2}, U_{k + 1} = U_{k} (7 I + V_{k}) (I + 3 V_{k}) {(I + 18 V_{k} + 13 Z_{k})}^{- 1}$	$\frac{1}{2} s$	5	1	4	$4^{\frac{1}{3 s}}$
R6 (Khaksar Haghani, 2014)	$U_{k + 1} = U_{k} (38 I + 42 V_{k}) {(9 I + 60 V_{k} + 11 Z_{k})}^{- 1}$	$\frac{2}{3} s$	4	1	3	$3^{\frac{3}{10 s}}$
R7 (Soleymani et al., 2016)	$\begin{matrix} W_{k} = Z_{k} V_{k}, T_{k} = W_{k} V_{k}, Y_{k} = 81 I + 2524 V_{k} + 6990 Z_{k} + 3084 W_{k} + 121 T_{k} \\ U_{k + 1} = U_{k} (684 I + 5316 V_{k} + 5876 Z_{k} + 924 W_{k}) Y_{k}^{- 1} \end{matrix}$	$\frac{1}{3} s$	6	1	6	$6^{\frac{3}{7 s}}$

Source(s): Created by authors

Table 2

The methods to find the matrix sign function

Method	Form
E1	$Y_{k + 1} = \frac{1}{2} (Y_{k} + Y_{k}^{- 1})$
E2	$Y_{k + 1} = \frac{1}{2} Y_{k} (3 I + Y_{k}^{2})$
E3	$Y_{k + 1} = Y_{k} (3 I + Y_{k}^{2}) {(I + 3 Y_{k})}^{- 1}$
E4	$Y_{k + 1} = Y_{k} (18 I - 20 Y_{k}^{2} - 30 Y_{k}^{4}) {(5 I + 15 Y_{k}^{2} - 45 Y_{k}^{4} - 7 Y_{k}^{6})}^{- 1}$
E5	$Y_{k + 1} = Y_{k} (105 I - 252 Y_{k}^{2} - 210 Y_{k}^{4} + 564 Y_{k}^{6} + 49 Y_{k}^{8}) {(25 I + 21 (4 Y_{k}^{2} - 26 Y_{k}^{4} + 20 Y_{k}^{6} + 13 Y_{k}^{8}))}^{- 1}$

Source(s): Created by authors

Table 3

Methods and their computational efficiency

Method	Form	s_i	m	c	p	EI
R8	$U_{k + 1} = U_{k} (\frac{81}{125} {(\frac{1}{5} I + V_{k})}^{- 1} + \frac{51}{100} I - \frac{1}{20} V_{k})$	$\frac{1}{2} s$	2	1	4	$4^{\frac{2}{3 s}}$
R9	$U_{k + 1} = U_{k} (\frac{27783}{64000} {(\frac{1}{20} I + V_{k})}^{- 1} + \frac{2217}{3200} I - \frac{17}{160} V_{k})$	$\frac{2}{3} s$	2	1	3	$3^{\frac{1}{2 s}}$
R10	$U_{k + 1} = U_{k} (\frac{29791}{72000} {(\frac{1}{30} I + V_{k})}^{- 1} + \frac{1709}{2400} I - \frac{9}{80} V_{k})$	$\frac{1}{3} s$	2	1	3	$3^{\frac{1}{2 s}}$

Source(s): Created by authors

Table 4

Several matrices for Example 5.1

Matrix	Entry
A (Pascal)	a_i1 = a_1j = 1, a_ij = a_i−1,j + a_i,j−1
B (Lotkin)	a_1j = 1, $a_{i j} = \frac{1}{i + j - 1}$
C (Hilbert,Hilbert)	$a_{i j} = (\frac{1}{i + j - 1}, \frac{1}{i + j - 1})$

Source(s): Created by authors

Table 5

Results of iterations for matrix A

n	Method	R1	R2	R3	R4	R5	R6	R7	R8	R9	R10
10	m	111	0	–	69	65	72	72	38	26	24
	c	37	37	–	23	13	18	12	19	13	12
	m + 2c	185	74	–	115	91	108	96	76	52	48
	CPU	0.0200	0.0162	–	0.0127	0.0139	0.0117	0.0111	0.0091	0.0083	0.0055
20	m	225	0	–	147	130	148	144	76	48	44
	c	75	123	–	49	26	37	24	38	24	22
	m + 2c	375	246	–	445	182	222	192	152	96	88
	CPU	0.0210	0.0202	–	0.0142	0.0143	0.0149	0.0135	0.0131	0.0105	0.0084
30	m	285	–	–	171	170	176	180	98	60	54
	c	95	–	–	57	34	44	30	49	30	27
	m + 2c	475	–	–	285	238	264	240	196	120	108
	CPU	0.0538	–	–	0.0363	0.0246	0.0265	0.0258	0.0236	0.0185	0.0175

Source(s): Created by authors

Table 6

Results of iterations for matrix B

n	Method	R1	R2	R3	R4	R5	R6	R7	R8	R9	R10
10	m	147	0	–	93	90	96	96	52	34	30
	c	49	49	–	31	18	24	16	26	17	15
	m + 2c	245	98	–	155	126	144	128	104	68	60
	CPU	0.0200	0.0176	–	0.0132	0.0115	0.0120	0.0125	0.0099	0.0077	0.0052
20	m	216	126	–	135	115	128	132	70	42	40
	c	72	72	–	45	23	32	22	35	21	20
	m + 2c	360	270	–	225	161	192	176	140	84	80
	CPU	0.0228	0.0112	–	0.0170	0.0189	0.0149	0.0171	0.0134	0.0089	0.0080
30	m	198	–	–	126	130	132	126	70	44	40
	c	66	–	–	42	26	33	21	35	22	20
	m + 2c	332	–	–	110	182	198	168	140	88	80
	CPU	0.0324	–	–	0.0204	0.0169	0.0171	0.0177	0.0158	0.0108	0.0093

Source(s): Created by authors

Table 7

Results of iterations for matrix C

n	Method	R1	R2	R3	R4	R5	R6	R7	R8	R9	R10
10	m	147	0	–	93	90	96	96	52	32	30
	c	49	49	–	31	18	24	16	26	16	15
	m + 2c	245	98	–	155	126	144	128	104	64	60
	CPU	0.0196	0.0159	–	0.0159	0.0129	0.0120	0.0123	0.0114	0.0087	0.0069
20	m	195	214	–	129	115	124	126	74	42	40
	c	65	107	–	43	23	31	21	37	21	20
	m + 2c	325	328	–	215	161	186	168	148	84	80
	CPU	0.0241	0.0254	–	0.0160	0.0151	0.0181	0.0172	0.0144	0.0109	0.0098
30	m	210	–	–	123	120	128	126	70	46	40
	c	67	–	–	41	24	32	21	35	23	20
	m + 2c	335	–	–	205	168	192	168	140	92	80
	CPU	0.0582	–	–	0.0335	0.0267	0.0288	0.0265	0.0253	0.0187	0.0156

Source(s): Created by authors

Table 8

Several matrices for Example 5.2

Matrix	Entry
D (Hilbert)	$a_{i j} = \frac{1}{i + j - 1}$
E (Cauchy)	$a_{i j} = \frac{1}{i - j + 0.5}$
F(Lehmer,Lotkin)	$a_{i j} = (a_{i j} = a_{j i} = i / j, a_{i j} = 1, a_{i j} = \frac{1}{i + j - 1})$

Source(s): Created by authors

Table 9

Results of iterations for matrix D

n	Method	E1	E2	E3	E4	E5	E6	E7	E8
10	m	0	162	93	135	144	52	32	30
	c	49	0	31	27	24	26	16	15
	m + 2c	98	98	155	189	192	104	648	60
	CPU	0.0075	0.0123	0.0097	0.0104	0.0115	0.0104	0.0067	0.0052
	‖AXA − XAX‖	1.02E-15	1.25E-15	2.68E-14	8.54E-15	3.60E-15	1.69E-15	8.54E-15	1.05E-15
20	m	0	210	126	170	180	70	44	40
	c	69	0	42	34	30	35	22	20
	m + 2c	138	210	210	238	240	140	88	80
	CPU	0.0148	0.0183	0.0122	0.0134	0.0152	0.0130	0.0104	0.0069
	‖AXA − XAX‖	2.6105	2.30E-13	1.09E-14	5.73E-14	1.83E-12	1.81E-14	1.30E-14	6.61E-14
30	m	0	–	126	180	186	72	44	44
	c	66	–	42	36	31	36	22	22
	m + 2c	132	–	210	252	248	144	88	88
	CPU	0.0162	–	0.0137	0.0180	0.0191	0.0162	0.0099	0.0086
	‖AXA − XAX‖	3.2969	–	1.88E-14	8.22E-14	1.54E-11	3.26E-14	2.49E-14	4.45E-14

Source(s): Created by authors

Table 10

Results of iterations for matrix E

n	Method	E1	E2	E3	E4	E5	E6	E7	E8
10	m	0	0	96	140	150	54	34	32
	c	51	84	32	28	25	27	17	16
	m + 2c	102	168	160	196	200	108	68	64
	CPU	0.0054	0.0053	0.0034	0.0043	0.0049	0.0031	0.0021	0.0019
	‖AXA − XAX‖	6.98E-15	1.25E-15	3.05E-16	5.64E-14	8.23E-13	1.45E-15	8.69E-15	3.27E-16
20	m	0	–	-132	170	180	72	42	38
	c	64	–	44	34	30	36	21	19
	m + 2c	128	–	220	238	240	144	84	76
	CPU	0.0105	–	0.0080	0.0093	0.0105	0.0095	0.0072	0.0048
	‖AXA − XAX‖	0.5635	–	4.51E-14	5.26E-14	2.77E-13	4.70E-15	1.33E-14	1.26E-14
30	m	0	–	123	175	186	68	44	40
	c	69	–	41	35	31	34	22	20
	m + 2c	138	–	205	245	248	136	88	80
	CPU	0.0162	–	0.0093	0.0088	0.0099	0.0061	0.0041	0.0037
	‖AXA − XAX‖	8.54E-15	–	6.58E-15	5.32E-15	8.54E-15	8.02E-15	4.56E-15	9.34E-15

Source(s): Created by authors

Table 11

Results of iterations for matrix F

n	Method	E1	E2	E3	E4	E5	E6	E7	E8
10	m	0	34	21	30	30	12	10	10
	c	11	0	7	6	5	6	5	5
	m + 2c	22	34	35	42	40	24	20	20
	CPU	0.0053	0.055	0.0060	0.0082	0.0097	0.0062	0.0055	0.0038
	‖AXA − XAX‖	8.09E-14	9.49E-14	8.04E-14	6.25E-14	6.44E-14	8.69E-14	3.54E-14	4.52E-14
20	m	0	40	27	35	54	14	12	12
	c	13	0	9	7	9	7	6	6
	m + 2c	26	40	45	49	72	28	24	24
	CPU	0.0093	0.0104	0.0106	0.0108	0.0129	0.0109	0.0098	0.0060
	‖AXA − XAX‖	8.54E-15	5.98E-13	9.58E-13	9.54E-13	4.56E-13	1.42E-13	8.04E-13	3.50E-13
30	m	0	44	27	40	48	16	12	12
	c	14	0	9	8	8	8	6	6
	m + 2c	28	44	45	56	64	32	24	24
	CPU	0.0096	0.0094	0.0096	0.0121	0.0129	0.0106	0.0094	0.0071
	‖AXA − XAX‖	8.54E-15	8.59E-13	6.04E-13	8.24E-13	1.24E-12	5.98E-12	1.54E-13	9.66E-13

Source(s): Created by authors

Word count: 2751

Show less

Convergence theory of efficient parametric iterative methods for solving the Yang-Baxter-like matrix equation

Content area

Abstract

Full text

1. Introduction

3.1 Extension of the methods for matrix sign

4. Efficiency index

5. Numerical examples

6. Conclusions