# Constructing a High-Order Globally Convergent Iterative Method for Calculating the Matrix Sign Function

Authors: Haifa Bin Jebreen
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This work is concerned with the construction of a new matrix iteration in the form of an iterative method which is globally convergent for finding the sign of a square matrix having no eigenvalues on the axis of imaginary. Toward this goal, a new method is built via an application of a new four-step nonlinear equation solver on a particulate matrix equation. It is discussed that the proposed scheme has global convergence with eighth order of convergence. To illustrate the effectiveness of the theoretical results, several computational experiments are worked out.


---

## Body

## 1. Preliminaries

The sign function for the scalar case is expressed by(1)signω=1,Reω>0,-1,Reω<0,wherein z∈C is not located on the imaginary axis. Roberts in [1] for the first time extended this definition for matrices, which has several important applications in scientific computing, for example see [2–4] and the references cited therein. For example, the off-diagonal decay of the matrix function of sign is also a well-developed area of study in statistics and statistical physics [5].To proceed formally, let us consider thatA∈Cn×n is a square matrix possessing no eigenvalues on the axis of imaginary. We consider(2)A=TJT-1,as a form of Jordan canonical written such that(3)J=diagJ1,J2,and the eigenvalues of J1∈Cp×p are in the open left half-plane, while the eigenvalues of J2∈Cq×q are in the open right half-plane. It is now possible to write the following [6]:(4)S=signA=T-Ip00IqT-1,wherein p+q=n. Noting that sign(A) is definable once A is nonsingular.This procedure takes into account a clear application of the form of Jordan canonical and of the matrixT. Herein none of the matrices T and J are unique. However, it is possible to investigate that sign(A) as provided in (4) does not rely on the special selection of T or J.Here, a simpler interpretation for the sign matrix in the case of Hermitian case (namely, all eigenvalues are real) can be given by(5)S=Udiagsignλ1,…,signλnU∗,wherein(6)U∗AU=diagλ1,…,λn,is a diagonalization of the square matrix A.The significance of calculating and findingS in (4) is because of the point that the function of sign plays a central role in matrix functions theory, specially for principal matrix roots and the polar decomposition; for more one may refer to [7–9].Bini et al. in [10] proved that the principal p-th root of the matrix A could be written as a multiple of the (2,1)-block associated with the sign matrix sign(C), associated with the block companion matrix:(7)C=0I0I⋱⋱⋱IA0∈Cpn×pn.It is requisite to focus of the most general case, i.e., when all the eigenvalues are complex rather than being narrow over a range of specific matrices, such as in (5).An important point goes to the fact that althoughsign(A) is a square root of the unit matrix, it is not equal to I or -I unless the A’s spectrum locates completely in the open right or left half-plane(s), respectively. Thus, the sign function is a nonprimary square root of I.Apart from (4), an efficient way to derive matrix iterations for some matrix functions is to apply the zero-finding iterative methods for solving operator equations which here is a matrix equation. Toward this goal, it is necessary to tackle a nonlinear equation as comes next:(8)FX≔X2-I=0,wherein I is a unit matrix, so as to propose matrix methods for S. The main motivation in this work is to extend the recently published results of the literature [11, 12] in this category by providing a useful novel method for calculating sign matrix. Furthermore, the proposed procedure can be applied for the calculation of polar decomposition, principal matrix square root, and several other scientific computing problems.After this brief introduction about the matrix function of sign in Section1, the remaining sections of this study are given as comes next. Section 2 shortly surveys the existing matrix iterations and their importance for computing S. In Section 3, it is discussed how we construct a new method having global convergence behavior and not belonging to the class of Padé family of iterations for computing S. It too manifests that the proposed scheme is convergent with high order of convergence. Computational reports are furnished to illustrate the higher computational precision of the constructed solvers in Section 4. A final remark of the manuscript is given in Section 5 with some directions for future studies.


## 2. The Literature

In the current research work, iterative methods are the main focus for calculatingS. As a matter of fact, such matrix iteration methods are Newton-type schemes that are in essence fixed-point schemes by providing a convergent matrix sequence by imposing a sharp initial value.The connection of matrix iterative expressions with the function of sign is not that straightforward, but in practice, such methods could be constructed by considering a suitable root-finding method to the nonlinear matrix equation (8). Noting that sign(A) is a solution of this equation (refer to [12] and the references cited therein).Applying the classic Newton’s method (NM) to (8) yields(9)Xk+1=12Xk+Xk-1.An inversion-free version of (9), called Newton-Schultz iteration [6], is defined by(10)Xk+1=12Xk3I-Xk2,by applying the well-known Schulz inverse-finder in order to remove the computation of the inverse matrix per computing step.The Newton-Schulz scheme is a second-order convergent, inversion-free method in calculating the sign matrix, but it suffers from the drawback that its convergence unlike the Newton’s method (9) is local.Analogously, the third-order convergent Halley’s method (HM) [13] for calculating the sign matrix is defined by(11)Xk+1=I+3Xk2Xk3I+Xk2-1.It is noted that all the above-mentioned schemes are particular cases of the Padé family presented and discussed in [13, 14]. The Padé approximation belongs to a broader category of rational approximations. Coincidentally, the best uniform approximation of the sign function on a pair of symmetric but disjoint intervals can be expressed as a rational function.Recently a fourth-order iterative method was furnished in [15] as follows:(12)Xk+1=I+18Xk2+13Xk4Xk7I+22Xk2+3Xk4-1.


## 3. Construction of a New Matrix Method

Assume the following nonlinear equation:(13)gx=0,wherein g:D⊆C→C is a scalar function. In what follows, let us first present a new scheme in nonlinear equation solving. The idea of increasing the order (see, e.g., [16, 17]) is to consider several substeps, while the newly appearing first derivatives are approximated via a secant-like approximation. Thus, we may write(14)yk=xk-sk,zk=xk-1+gykgxk-5/3gykgxkg′xk,wk=zk-gzkgyk,zk,xk+1=wk-gwkgwk,zk,whereas the first-order divided difference operator is defined by(15)gl,ħ≔għ-glħ-l.Here the point is that the new method should not aim at having theoptimal convergence order (in the sense of Kung-Traub, see, e.g., [16]) since such schemes lose the global convergence behavior when applied to (13). Since the final objective is to propose a method for matrix sign, we should take two things into account, which are having global convergence behavior and the novelty. In fact, the optimality conjecture of Kung-Traub is not useful once we extend iterative methods for calculating the sign of a matrix. Because, the optimality conjecture here ruins the final structure of the matrix method.One may also ask that why method (14) has been selected and what are the other ways for improving it. To answer these, we recall that the best way to improve the performance of (14) is to add one Newton’s substep at the end of its fourth step, which is costly since it includes the computation of the first derivative per cycle. In a similar way as in (14), we can add a secant-like substep and increase the convergence order. Generally speaking, a family of iterations can now be constructed in this way. On the other hand, since very higher order methods may not be useful in double precision arithmetic, namely, the practical environment that most researchers work in, we here only provide (14) and discuss its application and extension for matrix sign.Theorem 1.
Assume thatα∈D is a simple root of a function g:D⊆C→C, which is sufficiently differentiable and contains x0 as an initial value. Accordingly, the iterative expression (14) reads(16)ek+1=19c27ek8+Oek9,where cj=g(j)(α)/j!g′(α), and ek=xk-α.Proof.
The steps of proving the convergence order for this iterative method are via Taylor expansion, which is straightforward.Applying (14) on the matrix equation (8) will yield a novel matrix scheme to calculate (4) in its reciprocal form as follows:(17)Xk+1=Xk12I+200Xk2+560Xk4+344Xk6+36Xk8Ξk-1,wherein(18)Ξk=I+64Xk2+406Xk4+532Xk6+145Xk8+4Xk10,and the initial approximation is(19)X0=A.Applying (14) on the nonlinear equation g(x)=x2-1 contributes a global convergence in the complex plane (excluding the values locating on the axis of imaginary). The basins of attraction for (10) (locally convergent) and (9) (globally convergent) are portrayed in Figure 1. This global behavior of the proposed scheme, that is kept for matrix case, has been shown in Figures 2-3Figure 1
Basins of attractions for (9) in left and (10) in right (shading is done based on the number of cycles to achieve the convergence).Figure 2
Basins of attractions for (17) in left and its density plot on the same domain in right (shading is done based on the number of cycles to achieve the convergence).Figure 3
The density plot of the basins of attractions for (17) on [-50,50]×[-50,50] in left and on [-100,100]×[-100,100] in right (shading is done based on the number of cycles to achieve the convergence).To draw the basins of attractions, we consider a squareΓ=[-2,2]×[-2,2]∈C and allocate a color to any point x0∈Γ according to the simple zero, at which the new methods (or the existing methods for comparisons) converge. Subsequently, we highlight the point as black once the scheme diverge. Herein, we consider the stopping termination for convergence as gxk≤10-2. Using a different stopping criterion |xk+1-xk|≤10-2, the density plot of the basins of attraction only for the new high order of convergence method (17) are brought forward in Figures 2-3 on different domains.


### 3.1. Convergence Study

Now, it is shown that the proposed schemes are convergent, under standard conditions, namely, when there are no pure imaginary eigenvalues in the absence of rounding errors.Theorem 2.
Assume thatA∈Cn×n possess no eigenvalues on the axis of imaginary. Accordingly, the iterations {Xk}k=0k=∞ expressed by (17) is convergent to S, using (19).Proof.
Assume thatR is a rational operator in accordance with (17). Since X∈Cn×n has a form of Jordan canonical, there exists a matrix Z, so that(20)X=ZJZ-1.It is recalled that diagonalizable and nondiagonalizable matrices have a Jordan normal form A=TJT-1, whereas J includes the Jordan blocks. So,(21)RX=ZRJZ-1.An eigenvalue λ of Xk is transferred into an eigenvalue of R(λ) of Xk+1 via the iterative expression (17). This relevance and relation among the eigenvalues show that it is required to search how R(λ) transforms the complex plane into itself. It is recalled that R has the feature of sign preservation, namely,(22)signRx=signxx∈C.Moreover, it should have the global convergence, that is, the sequence defined by(23)xk+1=Rxk,with x0=x converges to sign(x) while x is not located on the axis of imaginary. At this moment, assume that the square matrix A have a form of canonical Jordan considered just like [6, p. 107]:(24)Z-1AZ=Λ=C00N,wherein Z is a not singular and C,N are square Jordan blocks in association with eigenvalues locating in C- and C+, respectively. We show by λ1,…,λp and λp+1,…,λn locating on the main diagonals of blocks C and N, respectively. Applying (24), it is possible to write(25)signA=Z-Ip00In-pZ-1.Taking (25) into account, it is easy to deduce(26)signΛ=signZ-1AZ=Z-1signAZ=signλ10⋱signλpsignλp+1⋱0signλn.From D0=Z-1AZ, we define(27)Dk=Z-1XkZ,in order to obtain a convergent sequence to sign(Λ). Thence, from the scheme (17), we simply could obtain(28)Dk+1=Dk12I+200Dk2+560Dk4+344Dk6+36Dk8×I+64Dk2+406Dk4+532Dk6+145Dk8+4Dk10-1.When D0 is a diagonal matrix, it is possible to show that all successive Dk are diagonal matrices, via mathematical induction. The other case when D0 is not diagonal and will be handled in the remaining part of the proof.
By rearranging (28) as n uncoupled scalar iterations as comes next:(29)dk+1i=12dki+200dki3+560dki5+344dki7+36dki91+64dki2+406dki4+532dki6+145dki8+4dki10, where(30)dki=Dki,i,1≤i≤n.Using (28) and (29), we should investigate the convergence of {dki} to sign(λi), for all 1≤i≤n. From (29) and because the eigenvalues of A are not pure imaginary, it is possible to write(31)signλi=si=±1.Thus, we attain(32)dk+1i-sidk+1i+si=--si+dkisi+dki8si-2dkisi+2dki2.Noting that the factor si-2dki/si+2dki, is bounded due to choosing an appropriate initial matrix (19). Since(33)d0i-sid0i+si<1,we attain(34)limk→∞dk+1i-sidk+1i+si=0,and, therefore,(35)limk→∞dki=si=signλi.At this point, it possible to derive that(36)limk→∞Dk=signΛ.Noting that if D0 is not diagonal, the relation among the iterates’ eigenvalues must be dealt with for (17). The eigenvalues of Xk are transformed from the iteration k to the iteration k+1 via(37)λk+1i=12λki+200λki3+560λki5+344λki7+36λki9×1+64λki2+406λki4+532λki6+145λki8+4λki10-1,1≤i≤n.The relation (37) manifests that the eigenvalues generally are convergent to si=±1, namely,(38)limk→∞λk+1i-siλk+1i+si=0.Ultimately, it would be easy to write that(39)limk→∞Xk=Zlimk→∞DkZ-1=ZsignΛZ-1=signA.This ends the convergence proof for (17) to calculate S.Now by considering (18) and the facts that Xk are rational functions of A, so, similar to A, commute with S, and S2=I, S-1=S, S2j=I, and S2j+1=S, j≥1. It can be shown that the new scheme reads the following error inequality:(40)Xk+1-S≤Ξk-1I-2Xk2Xk-S8.Inequality (40) shows the eighth order of convergence.A scaling technique to accelerate the initial phase of convergence is normally requisite since the convergence rate cannot be seen in the initial iterates. Such an idea was discussed fully in [18] for Newton’s method. A robust procedure to improve the initial convergence speed is to scale the iterations before each iteration; i.e., Xk should be moved to μkXk.If the scaling parameter (for the Newton’s method) is defined by [18],(41)ζk=Xk-1Xk,normscaling,ρXk-1ρXk,spectralscaling,detXk-1/n,determinantalscaling,then the accelerated forms of the proposed matrix iteration for S is defined as follows:(42)X0=A,ζk=isthescalingparametercomputedby41,Xk+1=ζkXk12I+200ζk2Xk2+560ζk4Xk4+344ζk6Xk6+36ζk8Xk8×I+64ζk2Xk2+406ζk4Xk4+532ζk6Xk6+145ζk8Xk8+4ζk10Xk10-1,where(43)limk→∞ζk=1,(44)limk→∞Xk=S.


## 3.1. Convergence Study

Now, it is shown that the proposed schemes are convergent, under standard conditions, namely, when there are no pure imaginary eigenvalues in the absence of rounding errors.Theorem 2.
Assume thatA∈Cn×n possess no eigenvalues on the axis of imaginary. Accordingly, the iterations {Xk}k=0k=∞ expressed by (17) is convergent to S, using (19).Proof.
Assume thatR is a rational operator in accordance with (17). Since X∈Cn×n has a form of Jordan canonical, there exists a matrix Z, so that(20)X=ZJZ-1.It is recalled that diagonalizable and nondiagonalizable matrices have a Jordan normal form A=TJT-1, whereas J includes the Jordan blocks. So,(21)RX=ZRJZ-1.An eigenvalue λ of Xk is transferred into an eigenvalue of R(λ) of Xk+1 via the iterative expression (17). This relevance and relation among the eigenvalues show that it is required to search how R(λ) transforms the complex plane into itself. It is recalled that R has the feature of sign preservation, namely,(22)signRx=signxx∈C.Moreover, it should have the global convergence, that is, the sequence defined by(23)xk+1=Rxk,with x0=x converges to sign(x) while x is not located on the axis of imaginary. At this moment, assume that the square matrix A have a form of canonical Jordan considered just like [6, p. 107]:(24)Z-1AZ=Λ=C00N,wherein Z is a not singular and C,N are square Jordan blocks in association with eigenvalues locating in C- and C+, respectively. We show by λ1,…,λp and λp+1,…,λn locating on the main diagonals of blocks C and N, respectively. Applying (24), it is possible to write(25)signA=Z-Ip00In-pZ-1.Taking (25) into account, it is easy to deduce(26)signΛ=signZ-1AZ=Z-1signAZ=signλ10⋱signλpsignλp+1⋱0signλn.From D0=Z-1AZ, we define(27)Dk=Z-1XkZ,in order to obtain a convergent sequence to sign(Λ). Thence, from the scheme (17), we simply could obtain(28)Dk+1=Dk12I+200Dk2+560Dk4+344Dk6+36Dk8×I+64Dk2+406Dk4+532Dk6+145Dk8+4Dk10-1.When D0 is a diagonal matrix, it is possible to show that all successive Dk are diagonal matrices, via mathematical induction. The other case when D0 is not diagonal and will be handled in the remaining part of the proof.
By rearranging (28) as n uncoupled scalar iterations as comes next:(29)dk+1i=12dki+200dki3+560dki5+344dki7+36dki91+64dki2+406dki4+532dki6+145dki8+4dki10, where(30)dki=Dki,i,1≤i≤n.Using (28) and (29), we should investigate the convergence of {dki} to sign(λi), for all 1≤i≤n. From (29) and because the eigenvalues of A are not pure imaginary, it is possible to write(31)signλi=si=±1.Thus, we attain(32)dk+1i-sidk+1i+si=--si+dkisi+dki8si-2dkisi+2dki2.Noting that the factor si-2dki/si+2dki, is bounded due to choosing an appropriate initial matrix (19). Since(33)d0i-sid0i+si<1,we attain(34)limk→∞dk+1i-sidk+1i+si=0,and, therefore,(35)limk→∞dki=si=signλi.At this point, it possible to derive that(36)limk→∞Dk=signΛ.Noting that if D0 is not diagonal, the relation among the iterates’ eigenvalues must be dealt with for (17). The eigenvalues of Xk are transformed from the iteration k to the iteration k+1 via(37)λk+1i=12λki+200λki3+560λki5+344λki7+36λki9×1+64λki2+406λki4+532λki6+145λki8+4λki10-1,1≤i≤n.The relation (37) manifests that the eigenvalues generally are convergent to si=±1, namely,(38)limk→∞λk+1i-siλk+1i+si=0.Ultimately, it would be easy to write that(39)limk→∞Xk=Zlimk→∞DkZ-1=ZsignΛZ-1=signA.This ends the convergence proof for (17) to calculate S.Now by considering (18) and the facts that Xk are rational functions of A, so, similar to A, commute with S, and S2=I, S-1=S, S2j=I, and S2j+1=S, j≥1. It can be shown that the new scheme reads the following error inequality:(40)Xk+1-S≤Ξk-1I-2Xk2Xk-S8.Inequality (40) shows the eighth order of convergence.A scaling technique to accelerate the initial phase of convergence is normally requisite since the convergence rate cannot be seen in the initial iterates. Such an idea was discussed fully in [18] for Newton’s method. A robust procedure to improve the initial convergence speed is to scale the iterations before each iteration; i.e., Xk should be moved to μkXk.If the scaling parameter (for the Newton’s method) is defined by [18],(41)ζk=Xk-1Xk,normscaling,ρXk-1ρXk,spectralscaling,detXk-1/n,determinantalscaling,then the accelerated forms of the proposed matrix iteration for S is defined as follows:(42)X0=A,ζk=isthescalingparametercomputedby41,Xk+1=ζkXk12I+200ζk2Xk2+560ζk4Xk4+344ζk6Xk6+36ζk8Xk8×I+64ζk2Xk2+406ζk4Xk4+532ζk6Xk6+145ζk8Xk8+4ζk10Xk10-1,where(43)limk→∞ζk=1,(44)limk→∞Xk=S.


## 4. Experiments

Herein, several experiments are discussed for the calculation of the sign matrix. The direct application of the new formulas for findingS is given below, though the application for computing the polar decomposition, finding the Yang-Baxter matrix equation can be given similarly. The simulations are run on an office laptop with Windows 7 Ultimate equipped Intel(R) Core(TM) i5-2430M CPU 315 2.40GHz processor and 16.00 GB of RAM on a 64-bit operating system. In addition, the simulations are done in Mathematica 11.0 [19].Various schemes are compared in respect to the iteration numbers and the elapsed CPU time. Globally convergent schemes are only included for comparison. The compared matrix methods are NM, HM, ANM, and PM1 (i.e., (17)). We do not include comparisons with methods having local convergence behavior such as the Newton-Schulz method (10) or (computationally expensive) methods from different categories such as the ones based on the computation of the Cauchy integral(45)signA=S=2π∫0∞t2I+A2-1dt.The stopping criterion for our simulations is defined by(46)Xk2-I2≤10-4.The reason of choosing (46) lies in the fact that, at each iterate, the obtained approximation should satisfy the main matrix equation. Thus, this criterion is much more trustful than other Cauchy-like terminations when calculating the sign of a matrix.Experiment 1.
Here, we calculate the sign matrix of 20 generated randomlycomplex matrices (with uniform distributions via the following piece of codes in the Mathematica environment) as comes nextSeedRandom[123]; number = 20;Table[A[l] = RandomComplex[-5 - 5 I, 5 + 5 I,50 l, 50 l];,l,number];
Noting that hereI=-1.The numerical reports are provided in Tables1-2 for various sizes of the input matrices based on the required number of steps and the elapsed CPU times. The results uphold the analytical parts and discussions of Sections 2-3. They manifest that there is a clear improvement in the iterations’ numbers and the total elapsed CPU time by applying (17). As a matter of fact, the mean of number of iterations and the CPU times listed in the last rows of each table indicate that the scheme (17) has the best performance in general.Table 1
Comparison of iterations’ numbers for Experiment1.
Matrix No. NM ANM HM PM1 A50×50 12 11 8 4 A100×100 15 14 10 5 A150×150 18 17 12 6 A200×200 15 14 10 5 A250×250 18 18 12 6 A300×300 23 21 14 6 A350×350 19 17 12 8 A400×400 17 16 11 6 A450×450 16 16 10 6 A500×500 20 20 13 7 A550×550 17 18 11 6 A600×600 20 19 13 6 A650×650 21 21 13 7 A700×700 17 18 11 6 A750×750 20 21 13 6 A800×800 20 22 13 6 A850×850 20 21 13 9 A900×900 20 21 13 8 A950×950 20 20 13 6 A1000×1000 21 24 14 8  Mean 18.45 18.45 11.95 6.35Table 2
Comparison of the elapsed (CPU) times for Experiment1.
Matrix No. NM ANM HM PM1 A50×50 0.148498 0.0230943 0.00861964 0.0102797 A100×100 0.0943384 0.109786 0.0727258 0.0392688 A150×150 0.232303 0.506882 0.189903 0.137378 A200×200 0.330597 0.709762 0.280812 0.21293 A250×250 0.657373 1.44232 0.555998 0.458002 A300×300 1.29408 2.56028 1.13594 0.801164 A350×350 1.60244 2.90168 1.36108 1.53705 A400×400 2.06766 3.83108 1.84136 1.8719 A450×450 2.6057 5.23433 2.41827 2.14952 A500×500 4.38739 9.09064 4.03174 3.20028 A550×550 4.89925 10.7225 4.83132 3.62766 A600×600 7.07534 14.5142 7.22801 4.66437 A650×650 9.72633 19.8393 8.51814 6.87333 A700×700 9.34249 21.6787 8.03064 7.05715 A750×750 13.9398 31.1319 11.5503 8.65203 A800×800 16.7687 39.0702 14.9441 10.4224 A850×850 19.9251 45.4667 18.0246 18.9427 A900×900 24.1888 52.8144 20.6969 19.2772 A950×950 28.187 58.5243 24.4635 17.3038 A1000×1000 34.3487 82.3412 29.5465 26.1742  Mean 9.09109 20.1257 7.98652 6.67063It is pointed out that the calculation ofXk2 per iteration for computing the stopping termination introduces one matrix product for NM, while the HM and the presented scheme calculate this matrix in the middle of each cycle.Experiment 2.
The aim of this test is to check the convergence history of different methods for a randomly1000×1000 generated complex matrix as follows:SeedRandom[1];A= RandomComplex[-100 - 100 I, 100 + 100 I,1000, 1000];
using the stopping criterion(47)Error=Xk2-IF.The results for Experiment2 are given in Figure 4 showing a stable and consistent behavior of the proposed scheme for finding matrix sign function.Figure 4
The convergence history of different methods in Experiment2.The numerical reports and evidences in Section4 improve themean of the CPU time, clearly. This was the main target of this paper in order to propose an efficient method.


## 5. Discussion

In various fields of numerical linear algebra and scientific computing, the theory and computation of matrix functions are very much useful. In the modern numerical linear algebra, it underlies an effective way introducing one to resolve the topical problems of the control theory. The function of a matrix can be defined in several ways, of which the following three are generally the most useful: Jordan canonical form, polynomial interpolation, and finally Cauchy integral.In this research work, we have focused on iterative methods for this purpose. Hence, a high-order nonlinear equation solver has been employed for constructing a novel scheme in calculating the sign matrix, which does not have pure imaginary eigenvalues.It was shown that the convergence is global via attraction basins in the complex plane and the rate of convergence is eight. Finally, some numerical experiments in double precision arithmetic were performed to manifest the superiority and applicability of (17). Outlines for future works can be forced to extend the discussed matrix iterations for calculating polar decompositions in future studies based on the application of the new schemes.



---

*Source: 8973867-2018-06-21.xml*