Rima Khouja
email: rima.khouja@inria.fr

Bernard Mourrain
email: bernard.mourrain@inria.fr

Jean-Claude Yakoubsohn

Newton-Type Methods For Simultaneous Matrix Diagonalization

Keywords: Simultaneous diagonalization, Newton-type method, eigenproblem, eigenvalues, certification, high precision computation MSC Classification: 65F15, 65H10, 15A18, 65-04 1

This paper proposes a Newton-type method to solve numerically the eigenproblem of several diagonalizable matrices, which pairwise commute. A classical result states that these matrices are simultaneously diagonalizable. From a suitable system of equations associated to this problem, we construct a sequence that converges quadratically towards the solution. This construction is not based on the resolution of a linear system as is the case in the classical Newton method. Moreover, we provide a theoretical analysis of this construction and exhibit a condition to get a quadratic convergence. We also propose numerical experiments, which illustrate the theoretical results.

Introduction 1.Our study

Let us consider p diagonalizable matrices M 1 , • • • , M p in C n×n which pairwise commute. A classical result states that these matrices are simultaneously diagonalizable, i.e., there exists an invertible matrix E and diagonal matrices Σ i , 1 i p, such that EM i E -1 = Σ i , 1 i p, see e.g. [START_REF] Roger | Matrix analysis[END_REF]. The aim of this paper is to compute numerically a solution (E, F, Σ) of the system of equations

f (E, F, Σ) := F E -I n F M E -Σ = 0 (1)
where Σ = (Σ 1 , . . . , Σ p ) and F M E -Σ := (F M 1 E -Σ 1 , . . . , F M p E -Σ p ) = 0. Notice that this system is multi-linear in the unknowns E, F, Σ. We verify that when p = 1 and M 1 is a generic matrix, this system has a solution set of dimension 2 n 2 + n -2n 2 = n (n 2 + n 2 + n unknowns for E, F, Σ and 2 matrix equations corresponding to n 2 + n 2 equations). However, for p > 1 and generic matrices M i , there is no solution. To have a solution, the pencil M must be on the manifold D p of p-tuples of simultaneously diagonalizable matrices.

The system (1) can be generalized to the following system:

f ′ (E, F, Σ ′ ) := F M 0 E -Σ 0 F M E -Σ = 0 (2)
where Σ ′ = (Σ 0 , Σ 1 , . . . , Σ p ), M 0 ∈ C n×n is replacing I n and Σ 0 is a diagonal matrix replacing I n in the first equation of [START_REF] Absil | Joint diagonalization on the oblique manifold for independent component analysis[END_REF]. When the pencil M ′ = (M 0 , M 1 , . . . , M p ) contains an invertible matrix, the solutions of the two systems are closely related. If M 0 is invertible, a solution (E, F, Σ ′ ) of [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF] for M ′ = (M 0 , M 1 , . . . , M p ) gives the solution (F M 0 , EΣ -1 0 , ΣΣ -1 0 ) of (1) for M = (M -1 0 M 1 , . . . , M -1 0 M p ). A similar correspondence between the solution sets can be obtained if a linear combination M ′ 0 = p i=1 λ i M i is invertible. As (2) can be seen as an homogeneisation of (1) and appears in several contexts and applications, we will also study Newton-type methods for this homogenized system.

To solve the system of equations [START_REF] Absil | Joint diagonalization on the oblique manifold for independent component analysis[END_REF], we propose to apply a Newton-like method and to analyze the Newton map associated to an iteration. These ideas have also been developed for instance in [START_REF] Mahony | The constrained newton method on a lie group and the symmetric eigenvalue problem[END_REF] where a Newton method is used for the symmetric eigenvalue problem. A Simultaneous Newton's iteration for ill-conditioned eigenproblem has been introduced in [START_REF] Françoise | Simultaneous newton's iteration for the eigenproblem[END_REF]. For more recent references using the Newton-type approach for eigenproblem see for instance [START_REF] Joho | Joint diagonalization of correlation matrices by using newton methods with application to blind signal separation[END_REF][START_REF] Sato | Riemannian newton-type methods for joint diagonalization on the stiefel manifold with application to independent component analysis[END_REF][START_REF] Joho | Newton method for joint approximate diagonalization of positive definite hermitian matrices[END_REF]. Moreover, similar approach for the fast computation of the singular value decomposition has been presented in a technical report [START_REF] Van Der Hoeven | Certified singular value decomposition[END_REF].

We say that we have a quadratic sequence associated to a system of equations if the sequence converges quadratically towards a solution.

The classical Newton map defines (E + X, F + Y, Σ + S) from (E, F, Σ) in order to cancel the linear part in the Taylor expansion of f (E +X, F +Y, Σ+S). An easy computation shows that the perturbations X, Y and S are solutions of such a Sylvester-type linear system

F E -I n + F X + Y E F M E -Σ -S + XM F + EM Y = 0. (3)
A straight-forward way to solve this linear system is via Kronecker product, see [START_REF] Horn | Topics in Matrix Analysis[END_REF]. This leads to a linear system of size 2n 2 , which can be solved in O(n 6 ) arithmetic operations.

The construction of the methods studied here is based on perturbations of such type (E(I n +X), (I n +Y )F, Σ+S) rather than (E +X, F +Y, Σ+S). More precisely the perturbations X, Y and S that we consider are perturbations which cancel the linear part of the Taylor expansion of f (E(I n + X), (I n + Y )F, Σ+S). In this case, we can produce explicit solutions for the linear system in X, Y and S given by:

Z + X + Y ∆ -S + ΣX + Y Σ = 0. (4)
where Z = F E -I n and ∆ = F M E -Σ. We will see that the linear system (4) admits an explicit solution (X, Y, S) with respect to Z and ∆ for p = 1, 2 in (1). This is because Σ is a diagonal matrix. From these considerations, we define and analyze a sequence that converges quadratically towards a solution of the system (1) without inverting a linear system at each step of this Newton-like method.

Related works

Simultaneous matrix diagonalization is required by many algorithms as it was pointed out in [START_REF] Laub | Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms[END_REF][START_REF] Bunse-Gerstner | A chart of numerical methods for structured eigenvalue problems[END_REF][START_REF] Vollgraf | Quadratic optimization for simultaneous matrix diagonalization[END_REF][START_REF] Jiang | Simultaneous diagonalization of matrices and its applications in quadratically constrained quadratic programming[END_REF][START_REF] Bernard | Simultaneous diagonalization algorithms with applications in multivariate statistics[END_REF]. A numerical analysis for two normal commuting matrices is proposed in [START_REF] Bunse-Gerstner | Numerical methods for simultaneous diagonalization[END_REF] using Jacobi-like methods. Their method adjusts the classical Jacobi method in successively solving n(n-1) 2 two-real-variables optimization problems at each sweep of the algorithm. Their main result states a local quadratic convergence and can be summarized in the following way. Let off 2 (A, B)

2 = i =j |A i,j | 2 + |B i,j | 2 .
Let {α 1 , . . . , α n } (resp. {β 1 , . . . , β n }) be the set of the eigenvalues of A (resp. B). Let A k and B k the matrices obtained at the step k of the Jacobi-like method and

ρ k = off 2 (A k , B k ). If ρ 0 < 1 2 δ := 1 4 min i =j (|α i -α j |, |β i -β j |), then ρ k+1 < 2n(9n -13) ρ 2 k δ .
We will see in Theorems 3 and 5 that the local conditions of the quadratic convergence do not depend on n. Many other papers studied the so-called Jacobi-like methods (see e.g. [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF], [START_REF] Mesloub | Efficient and stable joint eigenvalue decomposition based on generalized givens rotations[END_REF] and references therein).

In [START_REF] Hoeven | Efficient certification of numeric solutions to eigenproblems[END_REF] a sequence with proof of its convergence towards a numerical solution of the system (1) when p = 1, i.e., for M 1 , with the assumption of M 1 being a diagonalizable matrix, is presented. It requires matrix inversion. Furthermore, under some extra assumptions, its quadratic convergence is established.

For a pencil of real symmetric matrices C = (C 1 , . . . , C s ), several algorithms based on Riemannian optimization methods (see [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF]) have been developed in order to find an approximate joint diagonalizer (see e.g. [START_REF] Bouchard | Approximate joint diagonalization with Riemannian optimization on the general linear group[END_REF][START_REF] Absil | Joint diagonalization on the oblique manifold for independent component analysis[END_REF][START_REF] Rahbar | Geometric optimization methods for blind source separation of signals[END_REF][START_REF] Joho | Joint diagonalization of correlation matrices by using Newton methods with application to blind signal separation[END_REF]). The idea is to find a local minimizer B ∈ R n×n of an objective function f which measures the degree of non-diagonality of the pencil (BC 1 B T , . . . , BC s B T ) over a Riemannian manifold (see [START_REF] Wang | Penalty functionbased joint diagonalization approach for convolutive blind separation of nonstationary sources[END_REF][START_REF] Bouchard | Approximate joint diagonalization with Riemannian optimization on the general linear group[END_REF][START_REF] Afsari | Sensitivity analysis for the problem of matrix joint diagonalization[END_REF] for some examples of objective functions). This Riemannian manifold is defined according to the geometric constraints considered on B. For instance, the diagonalizer is supposed to be orthogonal in some of these algorithms after a pre-whitening step (see e.g. [START_REF] Cardoso | Blind beamforming for non gaussian signals[END_REF][START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF][START_REF] Flury | An algorithm for simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly diagonal form[END_REF][START_REF] Rahbar | Geometric optimization methods for blind source separation of signals[END_REF][START_REF] Douglas | Self-stabilized gradient algorithms for blind source separation with orthogonality constraints[END_REF][START_REF] Joho | Joint diagonalization of correlation matrices by using Newton methods with application to blind signal separation[END_REF][START_REF] Nikpour | Algorithms on the Stiefel manifold for joint diagonalisation[END_REF][START_REF] Nishimori | Learning algorithms utilizing quasi-geodesic flows on the Stiefel manifold[END_REF]). Due to inaccuracies in the computation of the diagonalizer with orthogonality constraints (see. [START_REF] Yeredor | Non-orthogonal joint diagonalization in the least-squares sense with application in blind source separation[END_REF]), oblique constraints, i.e., all the rows of the diagonalizer have unit Euclidean norm, have also been considered instead of the former constraints in more recent works (see e.g. [START_REF] Absil | Joint diagonalization on the oblique manifold for independent component analysis[END_REF][START_REF] Bouchard | Approximate joint diagonalization with Riemannian optimization on the general linear group[END_REF]). These algorithms can be used when the pencil of symmetric matrices is simultaneously diagonalizable. In this case, we aim to find a zero of the objective function f . However, these algorithms have a computation complexity higher than the Newton-type algorithm that we propose (see Proposition 4). For instance, most of them combine line search [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF]Ch4] or trust region [2, Ch7] methods, and matrix inversions at each iteration (see the exact Riemannian Newton iteration in [START_REF] Absil | Joint diagonalization on the oblique manifold for independent component analysis[END_REF]). Moreover, the points on the Riemannian manifold are updated using a retraction operator (see [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF]Ch4] or [START_REF] Bouchard | Approximate joint diagonalization with Riemannian optimization on the general linear group[END_REF] for an example of a retraction operator on the oblique manifold). In the Newtontype method described in Sections 3 and 4 the points are updated by using direct and explicit formulas. They have lower complexity than the Riemannian optimization-based algorithms and they are well-adapted to computation with high precision.

Simultaneous matrix diagonalization appears in many applications. For instance, in the solution of multivariate polynomial equations by algebraic methods, the isolated roots of the system are obtained from the computation of common eigenvectors of commuting operators of multiplication in the quotient ring and from their eigenvalues [START_REF] Cox | Using Algebraic Geometry[END_REF], [START_REF] Elkadi | Introduction à la résolution des systèmes polynomiaux[END_REF]. In the case of simple roots, this reduces to simultaneous diagonalization of a pencil of matrices.

The approach of approximate joint diagonalizer for a pencil of real symmetric matrices is used to solve Blind Source Separation (BSS) problem, with potential applications in wide domains of engineering (see e.g. [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF]).

Simultaneous matrix diagonalization of pencils of general matrices also appears in the rank (or canonical) decomposition of tensors [START_REF] De | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF]. Under certain conditions this rank decomposition is unique [START_REF] Nicholas | On the uniqueness of multilinear decomposition of n-way arrays[END_REF]. In this case simultaneous matrix diagonalization allows to compute this rank decomposition which plays a crucial role in numerous applications such that Psychometric [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of "Eckart-Young" decomposition[END_REF], Signal Processing and Machine Learning [START_REF] Cichocki | Tensor decompositions for signal processing applications: From two-way to multiway component analysis[END_REF], [START_REF] Nicholas D Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF], Sensor array processing [START_REF] Sørensen | Coupled canonical polyadic decompositions and multiple shift invariance in array processing[END_REF], Arithmetic Complexity [START_REF] Bürgisser | Algebraic complexity theory[END_REF], wireless communications [START_REF] Sørensen | Blind multichannel deconvolution and convolutive extensions of canonical REFERENCES polyadic and block term decompositions[END_REF], multidimensional harmonic retrieval [START_REF] Sørensen | Multidimensional harmonic retrieval via coupled canonical polyadic decomposition-part i: Model and identifiability[END_REF], [START_REF] Sørensen | Multidimensional harmonic retrieval via coupled canonical polyadic decomposition-part ii: Algorithm and multirate sampling[END_REF], Chemometrics [START_REF] Bro | Parafac tutorial and applications[END_REF], and Principal components analysis [START_REF] Jolliffe | Principal Component Analysis[END_REF].

Outline

Our contributions are a new iteration for the simultaneous diagonalization of matrices, with a local quadratic convergence and its analysis. The iteration is different from a Newton iteration. It does not require to invert a large linear system, but performs simple matrix operations. We analyse the numerical behavior of the method and provide a certification test for the convergence. Sections 2, 3, 4, and 5 are devoted to respectively constructing a sequence to solve numerically:

• F E -I n = 0,

• the system (1) when p = 1,

• the system (2) when p = 1,

• the system (1) for any p.

Moreover, we provide for these cases, a certification that the sequence converges to a nearby solution, and a test to detect when this convergence is quadratic from an initial point. More precisely, in Section 3 we show that a triplet (E 0 , F 0 , Σ 0 ) must satisfy a property depending on the quantity ε 0 := max(κ 2 0 K 2 0 Z 0 , κ 2 0 K 0 ∆ 0 ) to get a quadratic convergence where

1-Z 0 = F 0 E 0 -I n , 2-∆ 0 = F 0 M E 0 -Σ 0 , 3-κ 0 = max 1, max 1 j<k n 1 |σ 0,k -σ 0,j | , 4-K 0 = max k 1, |σ 0,k | ,
where σ 0,1 , . . . , σ 0,n denote the diagonal entries of Σ 0 . The quantity κ is the condition number of the studied methods. Based on the same methodology as in Section 3, Sections 4 and 5 exhibit a certification of the convergence of the sequence constructed to the studied case towards the solution with a sufficient condition on the initial point. In Section 6 we perform numerical experimentation. The final section is for our conclusions and future works.

Notation and preliminaries

Throughout this work, we will use the infinity vector norm and the corresponding matrix norm. For a given vector v ∈ C n and matrix M ∈ C n×n , they are respectively given by:

v = max{|v 1 |, . . . , |v n |} M = max v =1 M v . Explicitly, M = max{|m i,1 | + . . . + |m i,n | : 1 ≤ i ≤ n}.
For a second matrix N ∈ C n×n , we have

M + N M + N (sub-additivity) M N M N (sub-multiplicativity).
Moreover, for a given matrix M ∈ C n×n , we denote by M L,Tri and M F rob the following:

M L,Tri := max 1≤i≤n 1≤j≤i-1 |m i,j |,
i.e the max matrix norm of the lower triangular part of M,

M F rob := n i=1 n j=1 |m i,j | 2 ,
i.e., the Frobenius norm of M . Furthermore, we consider in this paper the regular case of diagonalizable matrices, that is, the matrices are diagonalizable with simple eigenvalues. Thus we will use the following notation

W n := {M ∈ C n×n | M with pairwise distinct eigenvalues}. It is well-known that W n is dense in C n×n .
The Lie group of n×n invertible matrices, denoted by GL n , is the so-called general linear group [START_REF] Andruchow | The left invariant metric in the general linear group[END_REF]. We denote by D n the vector space of diagonal matrices of size n and D ′ n denotes the subset of D n in which the diagonal matrices are

of n distinct diagonal entries. Let E, F ∈ GL n and Σ ∈ D ′ n . The tangent space of GL n at E (resp. F ) is denoted by T E GL n (resp. T F GL n ) and the tangent space of D ′ n at Σ is denoted by T Σ D ′ n .
The perturbation of respectively E, F and Σ that we consider in this paper are of the following form: E + Ė, F + Ḟ and Σ + Σ, where Ė and Ḟ are respectively in T E GL n and T F GL n and Σ is in T Σ D ′ n . As GL n is a Lie group, Ė and Ḟ can be written as EX and Y F such that X, Y are in the Lie algebra of GL n which is equal to C n×n (since this Lie algebra is

T In GL n and GL n is an open subset in C n×n ). As D ′ n is open in D n then T Σ D ′ n = D n , herein Σ = S ∈ D n .
Finally, the perturbations of E, F and Σ that we consider are as follows: E + EX, F + Y F and Σ + S, such that X and Y are in C n×n and S is a diagonal matrix in C n×n .

For a matrix M ∈ C n×n , let diag(M ) be the diagonal matrix with the same diagonal as M and let off(M ) be the matrix where the diagonal term of M are replaced by 0. We have M = diag(M ) + off(M ). We say that M is an off-matrix if M = off(M ). In addition, let (λ 1 , . . . , λ n ) ∈ C n , diag(λ 1 , . . . , λ n ) is the diagonal matrix in C n×n of diagonal entries λ 1 , . . . , λ n .

The superscripts . t , . * and .

-1 are used respectively for the transpose, Hermitian conjugate, and the inverse matrix.

We state the following lemma which will be used in some of the proofs.

Lemma 1 Let ϕ(ε, u) = j 0 (1+uε 2 j )-1 εu . Given ε 1 2 , u 1
, and i 0, we have

j 0 (1 + uε 2 j+i ) 1 + 2uε 2 i ( 5
)
Proof Modulo taking ε 2 i instead of ε, it suffices to consider the case when i = 0. Now ϕ(ε, u) is an increasing function in ε and u, since its power series expansion in ε and u admits only positive coefficients. Consequently, ϕ(ε, u) ϕ( 12 , 1) = 2.

2 Newton-type method for the system

F E -I n = 0. Let f : GL n × GL n → C n×n , (E, F ) → F E -I n .
We consider the following perturbations E + EX, F + Y F of respectively E and F where X, Y ∈ C n×n .

To define the Newton sequence we have to solve the linear system obtained by canceling the linear part in the Taylor expansion of f (E + EX, F + Y F ).

The same methodology will be adopted in the next sections for the other considered systems. Hereafter, we detail the computation of the Newton sequence associated to the system F E -I n = 0. Moreover, a sufficient condition on the initial point for the quadratic convergence of this Newton sequence will be established.

Let Z = F E -I n . We observe that

f (E + EX, F + Y F ) = (F + Y F )(E + EX) -I n (6) = Z + (Z + I n )X + Y (Z + I n ) + Y (Z + I n )X. ( 7
)
We assume here that Z is of small norm, i.e., we start from an initial point (E 0 , F 0 ) close from the solution of the system F E -I n = 0. Consequently, the linear system of first order terms to solve is

Z + X + Y = 0. ( 8
)
Hence X = Y = -Z 2 is a solution of Equation ( 8). Moreover we get, by substituting in Equation ( 7) X and Y by -Z 2 ,

(F + Y F )(E + EX) -I n = Z 2 - 3 4 I n + Z 4 . (9)
Proposition

1 Let Z 0 = F 0 E 0 -In. Define X 0 = -Z0 2 , E 1 = E 0 (In + X 0 ), F 1 = (In + X 0 )F 0 and Z 1 = F 1 E 1 -In. Assume that Z 0 1. Then Z 1 Z 0 2 ( 10
)
Proof It follows easily from [START_REF] Bürgisser | Algebraic complexity theory[END_REF].

Theorem 2 Let E 0 and F 0 two complex square matrices of size n. Let Z 0 = F 0 E 0 -In and assume that ε = Z 0 < 1 2 . The sequences defined for i 0

Z i = F i E i -In X i = - Z i 2 E i+1 = E i (In + X i ) F i+1 = (In + X i )F i
converge quadratically towards the solution of F E -In = 0. Each E i , respectively F i are invertible and, if E∞ and F∞ are respectively the limits of sequences (E i ) i 0 and (F i ) i 0 we have for i 0,

E i -E∞ (1 + 2ε)2 -2 i+1 +1 ε E 0 , F i -F∞ (1 + 2ε)2 -2 i+1 +1 ε F 0 .
Proof First, by the assumption

F 0 E 0 -In = Z 0 < 1 2 , we have E 0 and F 0 are invertible. In fact, E 0 F 0 = In + E 0 F 0 -In = In + Z 0 is invertible when Z 0 < 1 which is the case since we suppose Z 0 < 1 2 . Let us prove by induction that Z k 2 -2 k +1 ε. Since ε < 1 2 , we have Z k+1 Z k 2 from (10) ε2 -2 k+1 +2 ε 2 -2 k+1 +1 ε. Consequently Z∞ = 0. Since X k = -Z k 2 we deduce X k 2 -2 k ε.
It follows X∞ = 0. We have

E k = E k-1 (In + X k-1 ) = E 0 (In + X 0 ) • • • (In + X k-1 )
.

Denoting W i = 0 k i (In + X k ), W∞ = k 0 (In + X k ) we compute W∞ -In k 0 (1 + 2 -2 k ε) -1 2ε by using Lemma 1.
Then W∞ is invertible and

W -1 ∞ 1 1 -2ε . Let E∞ = E 0 W∞. Hence E 0 = E∞W -1 ∞ .
In the same way

F 0 = W -1 ∞ F∞.
Finally, the identity F∞E∞ -In = 0 permits to conclude that E 0 and F 0 are invertible. In the same way we prove easily that

W i -In 2ε. It follows that W i is invertible. Since E i = E 0 W i we deduce that E i is invertible. Moreover W i -W∞ W i 1 - k i+1 (1 + X k ) (1 + W i -In ) k 0 (1 + 2 -2 k+i+1 ε) -1 (1 + 2ε)2 -2 i+1 +1 ε by using Lemma 1.
We deduce that

E i -E∞ (1 + 2ε)2 -2 i+1 +1 ε E 0 .
These properties also hold for the F i 's. The theorem is proved.

3 Newton-like method for diagonalizable matrices.

Let M ∈ W n , Σ ∈ D ′ n , E, F ∈ GL n . We aim to construct Newton sequences which converge towards the numerical solution of f (E, F, Σ) = 0 where f :

GL n × GL n × D ′ n → C n×n × C n×n , (E, F, Σ) → (F E -I n , F M E -Σ).
We consider in the same way as before the perturbations E + EX and F + Y F of respectively E and F and in addition the perturbation Σ + S of Σ such that S ∈ D n . We get with Z = F E -I n and ∆ = F M E -Σ :

(F + Y F )(E + EX) -I n =Z + (Z + I n )X + Y (Z + I n ) + Y (Z + I n )X ( 11
) (F + Y F )M (E + EX) -Σ -S =F M E -Σ -S + F M EX + Y F M E + Y F M EX =∆ -S + ΣX + Y Σ + ∆X + Y ∆ + Y (∆ + Σ)X (12)
As in the previous section we assume that (E, F, Σ) is sufficiently close to the solution of f (E, F, Σ) = 0, thus the linear system that we obtain from [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF] and ( 12) is

Z + X + Y = 0 ∆ -S + ΣX + Y Σ = 0
The following lemma gives a solution of this linear system.

Lemma 2 Let Σ = diag(σ 1 , • • • , σn), Z = (z i,j ) 1≤i,j≤n and ∆ = (δ i,j ) 1≤i,j≤n be given matrices in C n×n . Assume that σ i = σ j for i = j. Let S, X and Y be matrices defined by

S = diag(∆ -ZΣ) (13)
x i,i = 0 (14)

x i,j = -δ i,j + z i,j σ j σ i -σ j , i = j ( 15
)
y i,i = -z i,i (16)
y i,j = δ i,j -z i,j σ i σ i -σ j , i = j. ( 17
)
Then we have

Z + X + Y = 0 (18) ∆ -S + ΣX + Y Σ = 0 (19) Moreover X , Y κε(K + 1) (20)
where

ε max( Z , ∆ ), κ = max 1, max i =j 1 |σ i -σ j | and K = max(1, max i |σ i |).
Proof It is easy to verify that X + Y + Z = 0. In this way the equation ( 19) is equivalent to

∆ -S -ZΣ + ΣX -XΣ = 0.
Since diag(∆ -S -ZΣ) = diag(ΣX -XΣ) = 0 the formulas which define X follow easily. The bounds (20) also are obvious to establish.

In the next theorem we introduce the Newton sequences associated to the system f (E, F, Σ) = 0 with a sufficient condition on the initial point for its quadratic convergence.

Theorem 3 Let E 0 , F 0 ∈ GLn and Σ 0 ∈ D ′ n be given such that they define the sequences for i 0,

Z i = F i E i -In ∆ i = F i M E i -Σ i S i = diag(∆ i -Z i Σ i ) E i+1 = E i (In + X i ) F i+1 = (In + Y i )F i Σ i+1 = Σ i + S i ,
where S i , X i and Y i are defined by the formulas [START_REF] Cichocki | Tensor decompositions for signal processing applications: From two-way to multiway component analysis[END_REF][START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF][START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] De | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF][START_REF] Douglas | Self-stabilized gradient algorithms for blind source separation with orthogonality constraints[END_REF]. Let us define

κ 0 = max 1, max i =j 1 |σ 0,i -σ 0,j | , K 0 = max(1, max i |σ 0,i |) and ε 0 = max(κ 2 0 K 2 0 Z 0 , κ 2 0 K 0 ∆ 0 ). Assume that ε 0 0.033. ( 21
)
Then the sequences (Σ i, E i , F i ) i 0 converge quadratically to the solution of (F E -In, F M E -Σ) = 0. More precisely E 0 and F 0 are invertible and

E i -E∞ 8.1 × 2 1-2 i+1 E 0 ε 0 κK F i -F∞ 8.1 × 2 1-2 i+1 F 0 ε 0 κK . Σ i -Σ∞ 1.85 × 2 1-2 i ε 0 κ 2 K .
Proof Let us denote for each i 0,

ε = ε 0 ε i = max(κ 2 i K 2 i Z i , κ 2 i K i ∆ i ) κ = κ 0 κ i = max 1, max 1 j<k n 1 |σ i,k -σ i,j | K = K 0 K i = max 1≤k≤n 1, |σ i,k | ,
where σ i,1 , . . . , σ i,n denote the diagonal entries of Σ i . Let us show by induction on i that

ε i 2 1-2 i ε ( 22
) Σ i -Σ 0 (2 -2 2-2 i ) 2a κ ε (23)
with a = 1 1 -8ε

. These inequalities clearly hold for i = 0. Assuming that the induction hypothesis holds for a given i and let us prove it for i+1. We first prove that

Σ i+1 -Σ 0 (2 -2 2-2 i+1 ) 2a κ ε under the assumption Σ i -Σ 0 (2 -2 2-2 i ) 2a κ ε.
To do this, at the first step we show that this implies K -

4a κ ε K i ≤ K + 4a κ ε and 1 1 + 8aε κ κ i ≤ κ 1 -8aε . Let us prove K - 4a κ ε K i K + 4a κ ε.
We have

K i := Σ i ≤ Σ 0 + Σ i -Σ 0 K + (2 -2 2-2 i ) 2a κ ε K + 4a κ ε K(1 + 4aε).
This implies simultaneously

K i K -|K -K i | K - 4a κ ε and K i K(1 -4aε). Let us show that κ i ≤ κ 1 -8aε
. In fact, if the σ i,j 's are the diagonal values of Σ i , the Weyl's bound [START_REF] Weyl | Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung)[END_REF] implies that

|σ i,j -σ 0,j | Σ i -Σ 0 4a κ ε for 1 j n.
So that for 1 j < k n, we obtain using 1 -8aε 0 :

|σ i,k -σ i,j | |σ 0,k -σ 0,j | -|σ i,k -σ 0,k | -|σ i,j -σ 0,j | |σ 0,k -σ 0,j |(1 -κ|σ i,k -σ 0,k | -κ|σ i,j -σ 0,j |) |σ 0,j -σ 0,k |(1 -8aε) 0.
Finally, we get :

κ i κ 1 -8aε .
On the other hand the inequality

|σ i,k -σ i,j | |σ 0,k -σ 0,j | + |σ i,k -σ 0,k | + |σ i,j -σ 0,j |
implies in the same way that above

κ i 1 1 + 8aε κ.
Next we prove [START_REF] Van Der Hoeven | Certified singular value decomposition[END_REF] for i + 1. We know

S i = diag(∆ i -Z i Σ i ). Since ε i = max(κ 2 i K 2 i Z i , κ 2 i K i ∆ i ) and κ i , K i 1 then S i ≤ 2 κ i ε i 2(1 + 8aε) κ 2 1-2 i ε.
It follows :

Σ i+1 -Σ 0 S i + Σ i -Σ 0 2(1 + 8aε) κ 2 1-2 i ε + (2 -2 2-2 i ) 2a κ ε 2 -2 1-2 i (2 -1) 2a κ ε since 1 + 8aε = a 2 -2 1-2 i 2a κ ε But it is eay to see that 2 1-2 i 2 2-2 i+1 . Finally we get Σ i+1 -Σ 0 2 -2 2-2 i+1 2a κ ε.
Hence we can also write

K i - 2a κ i ε Σ i -Σ i+1 -Σ i K i+1 Σ i + Σ i+1 -Σ i K i + 2a κ i ε
Using more the Weyl's bound we can easily get that

κ i 1 + 4aε κ i+1 κ i 1 -4aε . Now we bound κ 2 i+1 K 2 i+1 Z i+1 . We have Z i+1 = Z i X i + Y i Z i + Y i (Z i + In)X i . Since X i , Y i ≤ κ i ( ∆ i + K i Z i ) 2 κ i K i ε i , we can write κ 2 i+1 K 2 i+1 Z i+1 κ 2 i+1 K 2 i+1 κ 3 i K 3 i 4ε 2 i + κ 2 i+1 K 2 i+1 κ 4 i K 4 i 4ε 3 i + κ 2 i+1 K 2 i+1 κ 2 i K 2 i 4ε 2 i 4 (2 + ε i ) κ i+1 K i+1 κ i K i 2 ε 2 i 4 (2 + ε i ) 1 + 2aε 1 -4aε 2 ε 2 i
On the other hand

∆ i+1 = ∆ i X i + Y i ∆ i + Y i (∆ i + Σ i )X i .
Hence

κ 2 i+1 K i+1 ∆ i+1 κ 2 i+1 K i+1 κ 3 i K 2 i 4ε 2 i + κ 2 i+1 K i+1 κ 4 i K 3 i 4ε 3 i + κ 2 i+1 K i+1 κ 2 i K i 4ε 2 i 4 (2 + ε i ) κ 2 i+1 K i+1 κ 2 i K i ε 2 i 4 (2 + ε i ) 1 + 2aε (1 -4aε) 2 ε 2 i It follows ε i+1 4 (2 + ε) 1 + 2aε 1 -4aε 2 ε 2 i 8 (2 + ε) 1 -6ε 1 -12ε 2 ε 2 1-2 i+1 ε 2 1-2 i+1 ε since 8 (2 + ε) 1 -6ε 1 -12ε 2 ε 1 for ε 0.033.
This completes the proof of the two induction hypothesis [START_REF] Hoeven | Efficient certification of numeric solutions to eigenproblems[END_REF][START_REF] Van Der Hoeven | Certified singular value decomposition[END_REF] at order i + 1. Let

W i = i k=0 (In + X k ). Since X k 2 κ k K k ε k 2(1 + 8aε) κK(1 -4aε) ε2 1-2 k 2 κK(1 -12ε) ε2 1-2 k Consequently, W∞ -In i 0 (1 + 2 κK(1 -12ε) ε2 1-2 i ) -1 4 κK(1 -12ε) ε from Lemma 1 0.22 κK since ε 0.033..
Hence W∞ is invertible and E 0 = E∞W -1 ∞ . This implies that E 0 is invertible. Moreover,

W i -W∞ W i 1 - k i+1 (1 + X k ) (1 + W i -In ) k 0 (1 + 2 κK(1 -12ε) ε × 2 1-2 k+i+1 ) -1 (1 + 0.22) × 4 κK(1 -12ε) × 2 1-2 i+1 ε from Lemma 1 8.1 κK × 2 1-2 i+1 ε.
We deduce that

E i -E∞ 8.1 κK × 2 1-2 i+1 E 0 ε.
In the same way we show that F 0 is invertible and

F i -F∞ 8.1 κK × 2 1-2 i+1 F 0 ε. Finally Σ i -Σ∞ k i Σ k+1 -Σ k k i 2 κ 2 k K k ε k   k 0 2 -2 k   2 1-2 i 2 κ 2 K(1 -12ε)(1 -8ε) ε 0.82 × 2.25 × 2 1-2 i ε κK since k 0 2 -2 k 0.82 and ε 0.033. 1.85 × 2 1-2 i ε 0 .
The theorem is proved.

Proposition 4

The complexity of one Newton iteration in Theorem 3 is in O(n 3 ).

Proof The computation of all the entries x i,j , y i,j of X i and Y i by the formulas [START_REF] Cichocki | Tensor decompositions for signal processing applications: From two-way to multiway component analysis[END_REF][START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF][START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] De | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF][START_REF] Douglas | Self-stabilized gradient algorithms for blind source separation with orthogonality constraints[END_REF] requires in total O(n 2 ) arithmetic operations. The computation of Z i , ∆ i , S i , E i+1 , F i+1 , which requires 6 backward stable matrix multiplications and diagonal matrix operations, has a complexity in O(n 3 ). Consequently, the complexity of each iteration is in O(n 3 ).

Remark 1 It is possible to generalize this approach to the case where the diagonal matrices are replaced by Jordan matrices.

4 Newton-like method for two simultaneously diagonalizable matrices.

Let M 1 , M 2 be two commuting matrices in W n , thus M 1 and M 2 are simultaneously diagonalizable. We aim to find E, F ∈ GL n which diagonalize simultaneously M 1 , M 2 so that:

F M k E = Σ k | k ∈ {1, 2}, and Σ 1 , Σ 2 ∈ D ′ n . This equivalent to find the numerical solution of f (E, F, Σ 1 , Σ 2 ) = 0 such that f : (E, F, Σ 1 , Σ 2 ) → (F M 1 E -Σ 1 , F M 1 E -Σ 1 )
We consider as before the perturbations E + EX, F + Y F and Σ k + S k of respectively E, F and Σ k for k∈ {1, 2} . Letting Z k = FM k E -Σ k for k = 1, 2, we have:

(F + YF)M k (E + EX) -(Σ k + S k ) = Z k -S k + Σ k X + Y Σ k + Z k X + YZ k +Y (Z k + Σ k )X (24)
By assuming Z 1 , Z 2 are of small norm, the linear system to solve from Equation ( 24) is the following

Z k -S k + Σ k X + Y Σ k = 0, k = 1, 2 (25)
A solution of ( 25) is given by the following lemma.

Lemma 3 Let Σ k = diag(σ k 1 , • • • , σ k n ), Z k = (z k i,j ) 1≤i,j≤n be given matrices in C n×n for k ∈ {1, 2}. Assume that σ 1 j σ 2 j σ 1 i σ 2 i = 0 for i = j.
Let X, Y , and S k be the matrices defined by

x i,i = 0 (26)

x i,j = σ 1 j z 1 i,j σ 2 j z 2 i,j σ 1 i σ 1 j σ 2 i σ 2 j , i = j ( 27
)
y i,i = 0 ( 28
)
y i,j = - σ 1 i z 1 i,j σ 2 i z 2 i,j σ 1 i σ 1 j σ 2 i σ 2 j , i = j ( 29
)
S k = diag(Z k ), k = 1, 2. ( 30
)
Then we have

Z k -S k + Σ k X + Y Σ k = 0, k = 1, 2 (31)
Moreover

X , Y 2κεK (32)
where

ε = max( Z 1 , Z 2 ), κ = max       1, max i =j 1 σ 1 i σ 1 j σ 2 i σ 2 j       , K = max(1, max i,k |σ k i |).
Proof It is easy to verify that the equation [START_REF] Laub | Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms[END_REF] implies that for i = j, σ k i x i,j + σ k j y i,j + z k i,j = 0 and that the solution of these equations is given by the formula ( 27), [START_REF] Joho | Joint diagonalization of correlation matrices by using newton methods with application to blind signal separation[END_REF]. Choosing

x i,i = y i,i =0, we take S k = diag(Z k + Σ k X + Y Σ k ) = diag(Z k ) since Σ k X + Y Σ k
is an off-matrix, to satisfy the equation [START_REF] Laub | Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms[END_REF]. The bounds [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF] follows easily from ( 27), [START_REF] Joho | Joint diagonalization of correlation matrices by using newton methods with application to blind signal separation[END_REF].

Theorem 5 Let E 0 , F 0 ∈ GLn and Σ 0,k = diag(σ k 0,1 , . . . , σ k 0,n ) ∈ D ′ n , k = 1, 2
, be given and let define the sequences for i 0 and k = 1, 2 by:

Z i,k = F i M k E i -Σ i,k S i,k = diag(Z i,k ) E i+1 = E i (In + X i ) F i+1 = (In + Y i )F i Σ i+1,k = Σ i,k + S i,k ,
where X i , Y i are defined by the formulas [START_REF] Jiang | Simultaneous diagonalization of matrices and its applications in quadratically constrained quadratic programming[END_REF][START_REF] Joho | Joint diagonalization of correlation matrices by using Newton methods with application to blind signal separation[END_REF][START_REF] Joho | Newton method for joint approximate diagonalization of positive definite hermitian matrices[END_REF][START_REF] Joho | Joint diagonalization of correlation matrices by using newton methods with application to blind signal separation[END_REF]. Let ε 0 = max( Z 0,1 , Z 0,2 ),

κ 0 = max     1, max i =j 1 σ 1 0,i σ 1 0,j σ 2 0,i σ 2 0,j     and K 0 = max(1, max j,k |σ k 0,j |). Assume that u := 4ε 0 κ 2 0 K 3 0 0.094. ( 33
)
Then the sequences

(Σ i,k, E i , F i ) i 0 converge quadratically to the solution of F M k E -Σ k for k = 1, 2.
More precisely E 0 and F 0 are invertible and

E i -E ∞ 1.46 × 2 1-2 i+1 E 0 u F i -F ∞ 1.46 × 2 1-2 i+1 F 0 u.
Proof Let us denote for each i 0,

ε = ε 0 ε i = max( Z i,1 , Z i,2 ) κ = κ 0 κ i = max       1, max 1 j<k n 1 σ 1 i,j σ 1 i,k σ 2 i,j σ 2 i,k       K = K 0 K i = max(1, max j,k (|σ k i,j |)),
where σ k i,1 , . . . , σ k i,n are the diagonal entries of Σ i,k . Let us show by induction on i that

ε i 2 1-2 i ε ( 34
) Σ i,k -Σ 0,k (2 -2 2-2 i )ε (35)
These inequalities clearly hold for i = 0. Assuming that the induction hypothesis holds for a given i and let us prove it for i + 1. We can notice that ε i ≤ 1. In fact by induction hypothesis, we have ε i ≤ 2 1-2 i ε 0 and from ( 33)

ε 0 = u 4κ 2 0 K 3 0 ≤ 1, since u ≤ 1 and κ 0 , K 0 ≥ 1. As 2 1-2 i ≤ 1, ∀i ≥ 0, we have ε i ≤ 1. We first prove that Σ i+1,k -Σ 0,k (2 -2 2-2 i+1 )ε under the assump- tion Σ i,k -Σ 0,k (2 -2 2-2 i )ε.
To do this, at the first step we show that this implies

K i ≤ K + 2ε and κ i ≤ κ 1 -8κε(K + ε) . Let us prove K i K + 2ε.
We have

K i := Σ i ≤ Σ 0 + Σ i -Σ 0 K + (2 -2 2-2 i )ε K + 2ε. Let us show that κ i ≤ κ 1 -8κε(K + ε)
. In fact, if the σ i,j k 's are the diagonal values of Σ k i , we have |σ k i,jσ k 0,j | Σ i,k -Σ 0,k 2ε for 1 j n and k = 1, 2. It follows :

|σ 1 i,j σ 2 i,k -σ 1 0,j σ 2 0,k | = |σ 1 i,j σ 2 i,k -σ 1 0,j σ 2 i,k + σ 1 0,j σ 2 i,k -σ 1 0,j σ 2 0,k | = |σ 2 i,k (σ 1 i,j -σ 1 0,j ) + σ 1 0,j (σ 2 i,k -σ 2 0,k )| 2ε|σ 2 i,k | + 2ε|σ 1 0,j | 2ε(K + 2ε) + 2εK = 4ε(K + ε). Now, |σ 1 i,j σ 2 i,k -σ 1 i,k σ 2 i+1,j | |σ 1 0,j σ 2 0,k -σ 1 0,k σ 2 0,j | -|σ 1 0,j σ 2 0,k -σ 1 i+1,j σ 2 i,k | -|σ 1 i,k σ 2 i,j -σ 1 0,k σ 2 0,j | |σ 1 0,j σ 2 0,k -σ 1 0,k σ 2 0,j |(1 -8kε(K + ε)).
Finally, we get :

κ i κ 1 -8κε(K + ε) .
To prove [START_REF] Nikpour | Algorithms on the Stiefel manifold for joint diagonalisation[END_REF] it is sufficient to write

Σ i+1,k -Σ 0,k S i,k + Σ i+1,k -Σ 0,k ε i + (2 -2 2-2 i )ε (2 1-2 i + 2 -2 2-2 i )ε (2 -2 2-2 i+1 )ε.
Let us prove [START_REF] Mesloub | Efficient and stable joint eigenvalue decomposition based on generalized givens rotations[END_REF]. Since we have

Z i+1,k = Z i,k X i + Y i Z i,k + Y i (Z i,k + Σ i,k )X i . we deduce Z i+1,k 2ε 2 i κ i K i + 2ε 2 i κ i K i + 4ε 2 i κ 2 i K 2 i (ε i + K i ) 4ε 2 i κ 2 i K i + 4ε 2 i κ 2 i K 2 i (1 + K i ) since ε i 1and κ i 1 3 × 4ε 2 i κ 2 i K 3 i = 12ε 2 i κ 2 i K 3 i since K i 1.
It follows

ε i+1 12κ 2 (K + 2ε) 3 (1 -8κε(K + ε)) 2 ε 2 i 12εκ 2 (K + 2ε) 3 (1 -8κε(K + ε)) 2 2 2-2 i+1 ε 3 1 + u 2 3 1 -2u 1 + u 4 2 u2 2-2 i+1 ε since ε K u 4 , κε u 4 2 1-2 i+1 ε since 3 1 + u 2 3 1 -2u 1 + u 4 2
2 -1 for u 0.094.

Let W i = i k=0 (I n + X k ). Since X l 2κ l K l ε l 2 κ 1 -8κε(K + ε) (K + 2ε)ε2 1-2 l 1 + u 2 u 2 1 -2u 1 + u 4 2 1-2 l 0.65 × 2 1-2 l u since u 0.094. Consequently, W ∞ -I n i 0 (1 + 0.65 × 2 1-2 i u) -1 1.3u from Lemma 1 1.3 × 0.094 = 0.1222 Hence W ∞ is invertible and E 0 = E ∞ W -1 ∞ . This implies that E 0 is invertible. Moreover, W i -W ∞ W i 1 - k i+1 (1 + X k ) (1 + W i -I n ) k 0 (1 + 0.059 × 2 1-2 k+i+1 ) -1 (1 + 0.1222) × 1.3 × 2 1-2 i+1 u 1.46 × 2 1-2 i+1 u.
We deduce that

E i -E ∞ 1.46 × 2 1-2 i+1 E 0 u.
In the same way we show that F 0 is invertible and

F i -F ∞ 1.46 × 2 1-2 i+1 F 0 u.
The theorem is proved.

5 Convergence of a pencil of simultaneously diagonalizable matrices.

In this section we present two strategies to solve the system (1) of a pencil of commuting matrices (M i ) 1≤i≤p in W n . The first strategy is trivial and consists of finding the common diagonalizers E and F of the pencil by numerically solving one of the systems (F E -

I n , F M 1 E -Σ 1 ) = 0 or (F M 1 E -Σ 1 , F M 2 E -Σ 1 )
= 0 using Theorem 3 or Theorem 5. Next we deduce the remaining diagonal matrices Σ i using the formulas

Σ i,k = E(:, k) * M i E(:, k) E(:, k) * E(:, k) 1 k n, 2 or 3 i p,
where E(:, k) is the k-th column in E.

In this strategy we use that a diagonalizer of one or two matrices of the pencil can diagonalize the other matrices of the pencil. We note that, in general, we don't have this property for simultaneously diagonalizable matrices, where, for instance, it is posssible to find a diagonalizer of M 1 which is not a common diagonalizer for the other matrices of the pencil. Nevertheless, this property holds here since we suppose that the matrices M i have simple eigenvalues. Another strategy is to find a "good" linear combination of the M i 's. This is based on Lemma 4 and Theorem 6.

Lemma 4 Let us suppose that the M i commute pairwise and they are linearly independent, i.e., p i=1 a i M i = 0 ⇒ a i = 0, i = 1 : p. Let E ∈ GLn and Σ i ∈ D ′ n be such that

E -1 M i E -Σ i = 0, i = 1 : p.
Let S ∈ C n×p and the column i of S is the diagonal of Σ i . Let σ = (σ 1, . . . , σn) and Σ = diag(σ). Then the matrix S has a full rank and α = (S * S) -1 S * σ satisfies p i=1

α i E -1 M i E -Σ = 0.
Proof Since the matrices M i are simultaneously diagonalizable there exists E be such that

E -1 M i E -Σ i = 0. The condition p i=1 α i Σ i -Σ = 0
is written as Sα = σ where S ∈ C n×p . The assumption p i=1 a i M i = 0 ⇒ a i = 0, i = 1 : p implies that the matrix has a full rank. Consequently, α = (S * S) -1 S * σ.

The lemma follows.

Theorem 6 Let M 1 , . . . , Mp ∈ C n×n be p simultaneously diagonalizable matrices and verify the assumption of linearly independent. Let us consider matrices E 0 , F 0 and Σ 0,i = diag(F 0 M E 0 ), i = 1 : p. Let us define the matrix S ∈ C n×p in which the column i is the diagonal of Σ 0,i . Let σ = 1, e 2iπ n , . . . , e 2i(n-1)π n

, Σ = diag(σ)
and α = (S * S) -1 S * σ. We consider the system

EF -I n F M E -Σ = 0 ( 36
)
where

M = p i=1 α i M i . If n 2 max( Z 0 , ∆ 0 ) 16 × 0.033
then (F 0 , E 0 , Σ) satisfies the condition ( 21) of Theorem 3.

Proof In this case the quantity κ defined in the Theorem 3 is equal to

κ = 1 2 |sin π n | n 4 since |sin π n | 2 n for n 2.
Since K 0 = 1 we get

ε 0 = max(κ 2 0 K 2 0 Z 0 , κ 2 0 K 0 ∆ 0 ) ≤ n 2 16 max( Z 0 , ∆ 0 ).
The condition max( Z 0 , ∆ 0 ) ≤ 0.033 16 n 2 , gives the result.

Numerical illustration

We use a Julia implementation of the Newton sequences in the numerical experiments. The experimentation has been done on a Dell Windows desktop with 8 GB memory and Intel 2.3 GHz CPU. We use the Julia package ArbNumerics for the computation in high precision.

Simulation

In this section we apply the Newton iterations presented in Theorem 3 (resp. Theorem 5) on examples of diagonalizable matrices (resp. of two simultaneously diagonalizable matrices). We validate experimentally the sufficiency of the condition established in Theorem 3 (resp. Theorem 5) to have a quadratic sequence (Tables 1, 2, 6 and7). On the other hand, as this condition is sufficient but not necessary, we show through some other examples how this Newton sequence starting from an initial point which is not verifying this condition could converge quadratically (Tables 3, 4, 8 and9). We note that the the computation in the aforementioned tables is done in high precision. Nevertheless, we test also the two Newton-type sequences using machine precision (Tables 5 and10) and this to show that these sequences have the same numerical behavior of a classical Newton method, i.e., if the solution is in the neighborhood of the initial point the Newton-type iterations will converge towards this solution with a few number of iterations and the residual error obtained at the end is in double precision.

This allows us to have an heuristic estimation on the numerical dependency of the Newton sequences from this condition to converge. Furthermore, these examples reveal the possibility of achieving computation in such problem with high precision. For example, in the case of a diagonalizable matrix of simple eigenvalues, we can compute its eigenvalues using one of the solvers which works with a double precision. Then we take this point as an initial point for the Newton sequence of Theorem 3 in order to increase the precision. Hereafter, we give some details about the tests: Test-1 for Theorem 3 and Test-2 for Theorem 5, considered in this section.

Test-1. Let K = R or C, M = EΣE -1 + 10 -e A, where e ∈ {3, 6}. The matrices E, Σ, and A ∈ K n×n are chosen randomly following standard normal distributions such that E is invertible, Σ is diagonal with n different diagonal entries and A is a random square matrix obeying normal distribution of size n and Frobenius norm equal to 1. Since M is a small perturbation of EΣE -1 , more precisely M -EΣE -1 F rob = 10 -e , M is a diagonalizable matrix of simple eigenvalues. Herein, we apply the Newton iteration of Theorem 3 on M with initial point E 0 = E, F 0 = E -1 and Σ 0 = Σ. The residual error reported in this test at iteration k is given by:

err res = max( F k E k -I n , F k M E k -Σ k ). Test-2. Let K = R or C, M 1 = F -1 Σ 1 E -1 , M 2 = F -1 Σ 2 E -1
, where E, F , Σ 1 and Σ 2 ∈ K n×n are randomly sampled according to standard normal distributions, such that E and F are invertible, Σ 1 and Σ 2 are diagonal with n different diagonal entries. The Newton iteration in Theorem 5 is applied on M 1 and M 2 with initial point E 0 , F 0 , Σ 0,1 and Σ 0,2 , such that these matrices are obtained by applying a small perturbation on respectively E, F , Σ 1 and Σ 2 as follows: E 0 = E + 10 -e A, F 0 = F + 10 -e B, Σ 0,1 = Σ 1 + 10 -e C, Σ 0,2 = Σ 2 + 10 -e D, where e ∈ {3, 6}, A and B (resp. C and D) are random square matrices (resp. random diagonal matrices with different diagonal entries) of size n and Frobenius norm equal to 1, with entries in K following standard normal distributions. The residual error reported in this test at iteration k is given by:

err res = max( F k M 1 E k -Σ k,1 , F k M 2 E k -Σ k,2
). We notice that the condition established in Theorem 3 (resp. Theorem 5) is reached in Test-1 (resp. Test-2 ) for matrices of size 10 with order of perturbation equal to 10 -6 , and we can see in Tables 1, 2, 6 and 7 that the Newton sequences with initial point verifying the condition in the associated theorem converge quadratically. We can notice also that by increasing the perturbation up to 10 -3 (the initial point does not verify the condition in the associated theorem), the Newton sequences converge quadratically for different sizes of matrices n = 10, 50, 100 (see Tables 3,4, 8 and 9). Moreover, we can notice in Table 5 the Newton-type iteration of Theorem 3 applied in double precision converges with a few number of iterations ∼ 5 and the final residual error measured with the Frobenius norm is of order machine precision ∼ 10 -14 and it is of the same order obtained by the standard Julia method eigen to compute the eigen decomposition. The same remarks are valid for Table 10 where the Newton-type sequence of Theorem 5 needs, in double precision, a few iterations to converges towards the solution given by using the Frobenius norm a residual error of order machine precision.

Table 1 The computational results throughout 7 iterations of an example of implementation of Test-1 with K = R, n = 10 and e = 6 in precision 1024.

Iteration ε := max(κ 2 0 K 2 0 Z0 , κ 2 0 K0 ∆0 ) ≤ 0.

Cauchy matrix

In this section we present an example for a Cauchy matrix of size n = 13 of entries a i,j = 1 i+j , ∀1 ≤ i, j ≤ 13, that illustrates how the Newton-type iteration can be used to increase the accuracy of the eigenvalues. We take the eigen decomposition given by the standard Julia method eigen from the package Table 11 The relative error between σeigen from the method eigen and σnewton from the Newton-type method for the Cauchy matrix

Sub-matrix iterations

It is possible to adapt the proposed method, taking into account the condition of the eigenvalue σ i given by the quantity

κ(σ i ) = max i =j 1, 1 |σ i -σ j |
Theoretical results imply that the computation of clusters of eigenvalues is illconditioned. However, one can apply Theorem 3 on sub-matrices to improve the well-conditioned eigenvalues. We denote δ = K ∆ 0 0.033 and p the index such that Σ = Σ p Σ n-p , Σ p = diag(σ 1 , . . . , σ p ), Σ n-p = diag(σ p+1 , . . . , σ n ) and |σ iσ j | > δ for all 1 i ≤ p and i < j n. We adapt Newton iteration to the block associated with the well-conditioned eigenvalues by defining the matrices X, Y and S as follows:

x i,i = 0

x i,j =    -δ i,j + z i,j σ j σ i -σj if |σ i -σ j | > δ 0 otherwise Y = -Z -X S = diag(-∆ + ZΣ).
Table 12 (resp. Table 13) shows the residual error err res as in Test-1 for the Cauchy matrix of size 200 (resp. the Rosser matrix of size 256 [START_REF] Rosser | Separation of close eigen-values of a real symmetric matrix[END_REF]) by applying the aforementioned sequences, the initial point is given by the Julia method eigen. The computation is done in precision 1024.

Conclusion

Taking a Newton approach towards systems of equations describing the simultaneous diagonalization problem of diagonalizable matrices, leads us to new algorithmic insights. We exhibit a Newton-type method without solving a linear system at each step as is the case of a classical Newton method. The numerical experiments corroborate the quadratic convergence predicted by the theoretical analysis. We focused on the regular case. Some improvements and extensions can be considered, such as the treatment of clusters of eigenvalues. Another direction that can be explored, is the construction of higher-order methods.

Table 2

 2 The computational results throughout 7 iterations of an example of implementation of Test-1 with K = C, n = 10 and e = 6 in precision 1024.

	033	errres

Table 3

 3 The residual error throughout 7 iterations given by the implementation of Test-1 with K = R, e = 3 and n = 10, 50, 100 in precision 1024.

	Iteration	n = 10	n = 50	n = 100
	1	8.57e -3	7.93e -2	3.22e -2
	2	1.91e -4	5.76e -2	1.38e -2
	3	1.58e -8	6.19e -3	6.12e -4
	4	4.79e -16	8.74e -5	5.42e -7
	5	3.56e -31	1.31e -8	3.83e -13
	6	1.39e -61	2.39e -16	1.80e -25
	7	1.91e -122	7.03e -32	3.81e -50

Table 4

 4 The residual error throughout 7 iterations given by the implementation of Test-1 with K = C, e = 3 and n = 10, 50, 100 in precision 1024.

	Iteration	n = 10	n = 50	n = 100
	1	8.84e -3	9.75e -2	1.61e -2
	2	8.59e -6	6.39e -5	1.03e -4
	3	3.91e -11	3.99e -9	4.68e -9
	4	9.87e -22	1.87e -17	3.13e -17
	5	7.60e -43	4.42e -34	8.84e -34
	6	5.14e -85	2.50e -67	9.45e -67
	7	2.64e -169 8.28e -134	1.05e -132

Table 5

 5 The residual error throughout 5 iterations given by the implementation of Test-1 with K = R, e = 3 and n = 10, 20, 30, in double precision.

	Iteration	n = 10	n = 20	n = 30
	1	4.78e -3	1.01e -2	1.01e -2
	2	4.71e -3	2.55e -3	1.14e -3
	3	2.29e -5	1.97e -5	4.08e -7
	4	1.43e -9	2.36e -10	2.26e -13
	5	4.06e -15	1.23e -14	5.04e -14
	M -EeigenΣeigenE -1 eigen F rob M -EnewtonΣnewtonE -1 newton F rob	9.49e -15 2.96e -15	2.83e -14 1.01e -14	7.45e -14 3.42e -14

Table 6

 6 The computational results throughout 7 iterations of an example of implementation of Test-2 with K = R, n = 10 and e = 6 in precision 1024.

	Iteration 4κ 2 K 3 ε ≤ 0.094	errres
	1	7.65e -2	6.72e -6
	2	1.73e -7	1.52e -11
	3	5.58e -18	4.90e -22
	4	5.49e -39	4.82e -43
	5	3.10e -81	2.73e -85
	6	2.28e -165	2.01e -169
	7	2.20e -279	1.94e -283

Table 7

 7 The computational results throughout 7 iterations of an example of implementation of Test-2 with K = C, n = 10 and e = 6 in precision 1024.

	Iteration 4κ 2 K 3 ε ≤ 0.094	errres
	1	6.86e -3	9.16e -6
	2	7.14e -9	9.53e -12
	3	9.51e -21	1.26e -23
	4	6.69e -44	8.92e -47
	5	3.77e -90	5.04e -93
	6	2.59e -182	3.45e -185
	7	1.65e -281	2.20e -284

  1 i+j 1≤i,j≤13 .

	Eigenvalue	σeigen	σnewton	|σnewton-σeigen| σnewton
	1	2.4030587641505818e-17	5.958203769841865e-19	39.33
	2	1.8824087522342697e-16	1.7156976132548192e-16	0.09716
	3	2.3152722725223998e-14	2.3178576801522747e-14	0.00111
	4	1.9513972147589434e-12	1.951356013568409e-12	2.11e -5
	5	1.1466969172503778e-10	1.1466967568738049e-10	1.39e -7
	6	4.991788233415145e-9	4.991788235245136e-9	3.66e -10
	7	1.6668681228080362e-7	1.666868122813953e-7	3.54e -12
	8	4.360227301207107e-6	4.360227301206033e-6	2.46e -13
	9	9.040674871074817e-5	9.040674871075823e-5	1.11e -13
	10	0.0014925044272821445	0.0014925044272821172	1.83e -14
	11	0.01955788569925287	0.01955788569925287	4.81e -17
	12	0.19958813407010345	0.19958813407010337	4.64e -16
	13	1.3693334145989837	1.3693334145989824	9.98e -16

Table 12

 12 The residual error throughout 6 iterations with the Cauchy matrix of size 200.

	Iteration p = 12, δ = 4.51e -7 p = 5, δ = 4.51e -7
	1	2.45e -15	2.35e -15
	2	9.63e -26	3.75e -29
	3	1.56e -36	1.21e -53
	4	1.54e -45	1.81e -83
	5	1.15e -54	3.49e -110
	6	5.08e -64	8.67e -137

Table 13

 13 The residual error throughout 6 iterations with the Rosser matrix of size 256.Iteration p = 11, δ = 1.11e -3 p = 5, δ = 1.11e -

	1	7.15e -12	1.65e -12
	2	7.18e -20	7.18e -20
	3	1.42e -40	1.81e -41
	4	1.73e -53	1.56e -85
	5	7.17e -66	1.75e -119
	6	8.79e -79	8.11e -153

Newton-Type Methods For Simultaneous Matrix Diagonalization

Table 8 The residual error throughout 7 iterations given by the implementation of Test-2 with K = R, e = 3 and n = 10, 50, 100 in precision 1024. LinearAlgebra as an initial point of Newton sequences in Theorem 3 with 5 iterations. The computation is done with the precision 1024 using ArbNumerics package. The initial point given by eigen is in double precision. It is converted to the precision 1024 using ArbNumerics package, in order to apply Newtons iterations with this precision of 1024 bits. In Table 11 we report the eigenvalues given by eigen (σ eigen ) and the eigenvalues rounded to the double precision given by Newton-type sequence (σ newton ) initialized with eigen.

We also report the relative error |σnewton-σeigen| σnewton in order to show the refinement amount realized by the Newton method. As we can see the matrix of this example is ill-conditioned (Cauchy matrices are in general ill-conditioned). There is a cluster of eigenvalues nearby zero. The accuracy enhancement obtained by applying Newton-type iterations can be clearly seen in Table 11, in particular for the first four smallest eigenvalues. For instance, the smallest eigenvalue returned by eigen is of order 10 -17 close to the second smallest eigenvalues of order 10 -16 . Newton-type method shows that the smallest eigenvalue of the order 10 -19 yields a large relative error ∼ 39.33. This also shows that all the eigenvalues are well-separated.

Conflict of Interest

The authors declare no conflicting interest that is directly or indirectly related to the work submitted for publication.