Geometric means of quasi-Toeplitz matrices
GGeometric means of quasi-Toeplitz matrices ∗ Dario A. Bini † Bruno Iannazzo ‡ Jie Meng † Abstract
We study means of geometric type of quasi-Toeplitz matrices, that aresemi-infinite matrices A = ( a i,j ) i,j =1 , ,... of the form A = T ( a ) + E , where E represents a compact operator, and T ( a ) is a semi-infinite Toeplitzmatrix associated with the function a , with Fourier series (cid:80) ∞ (cid:96) = −∞ a (cid:96) e i (cid:96)t ,in the sense that ( T ( a )) i,j = a j − i . If a is real valued and essentiallybounded, then these matrices represent bounded self-adjoint operatorson (cid:96) . We consider the case where a is a continuous function, wherequasi-Toeplitz matrices coincide with a classical Toeplitz algebra, and thecase where a is in the Wiener algebra, that is, has absolutely convergentFourier series.We prove that if a , . . . , a p are continuous and positive functions, orare in the Wiener algebra with some further conditions, then means ofgeometric type, such as the ALM, the NBMP and the Karcher meanof quasi-Toeplitz positive definite matrices associated with a , . . . , a p , arequasi-Toeplitz matrices associated with the geometric mean ( a · · · a p ) /p ,which differ only by the compact correction. We show by numerical teststhat these operator means can be practically approximated. Keywords:
Quasi-Toeplitz matrices, Toeplitz algebra, matrix func-tions, operator mean, Karcher mean, geometric mean, continuous func-tional calculus.
The concept of matrix geometric mean of a set of positive definite matrices,together with its analysis and computation, has received an increasing interestin the past years due to its rich and elegant theoretical properties, the nontrivialalgorithmic issues related to its computation, and also the important role playedin several applications. ∗ The work of the first two authors was partly supported by INdAM (Istituto Nazionale diAlta Matematica) through a GNCS Project. A part of this research has been done during avisit of the third author to the University of Perugia. Version of February 9, 2021 † Dipartimento di Matematica, Universit`a di Pisa, Italy ( [email protected],[email protected] ) ‡ Dipartimento di Matematica e Informatica, Universit`a degli Studi di Perugia, Italy( [email protected] ) a r X i v : . [ m a t h . NA ] F e b fruitful approach to get a matrix geometric mean has been to identify andimpose the right properties required by this function. These properties, knownas Ando-Li-Mathias axioms, are listed in the seminal work [1]. It is interestingto point out that, while in the case of two positive definite matrices A , A the geometric mean G = G ( A , A ) is uniquely defined, in the case of p > A , . . . , A p there are infinitely many functions fulfilling the Ando-Li-Mathias axioms. These axioms are satisfied, in particular, by the ALM mean[1], the NBMP mean, independently introduced in [34] and [10], and by theRiemannian center of mass, the so-called Karcher mean, identified as a matrixgeometric mean in [32]. We refer the reader to the book [4] for an introductionto matrix geometric means with historical remarks.These means have been generalized in a natural way to the infinite dimen-sional case. For instance, the Karcher mean has been extended to self-adjointpositive definite operators and to self-adjoint elements of a C ∗ -algebras in [29]and [27], respectively. Another natural generalization relies on the concept ofweighted mean: for instance, a generalization of the NBMP mean is presentedin [28].Matrix geometric means have played an important role in several appli-cations such as radar detection [26], [42], image processing [37], elasticity [33];while more recent applications include machine learning [40], [19], brain-computerinterface [45], [43], and network analysis [14]. The demand from applicationshas required some effort in the design and analysis of numerical methods formatrix geometric means. See in particular [10], [6], [5], [18], [22], [20], [44], andthe literature cited therein, where algorithmic issues have been investigated.In certain applications, the matrices to be averaged have further structuresoriginated by the peculiar features of the physical model that they describe.In particular, in the radar detection application, the shift invariance propertiesof some quantities involved in the model turn into the Toeplitz structure ofthe matrices. A (possibly infinite) matrix T = ( t i,j ) is said to be Toeplitz if t i,j = a j − i for some { a k } k ∈ Z ⊂ C . The sequence { a k } k ∈ Z may arise as the setof Fourier coefficients of a function a defined on the unit circle T = { z ∈ C : | z | = 1 } , with Fourier series (cid:80) ∞ (cid:96) = −∞ a (cid:96) z (cid:96) . The function a is said to be thesymbol associated with the Toeplitz matrix, and the Toeplitz matrix is denotedby T ( a ).The problem of averaging finite-dimensional Toeplitz matrices has beentreated in some recent papers, see for instance [7], [21], [35], [36], with theaim to provide a definition of matrix geometric mean G = G ( A , . . . , A p ) whichpreserves the matrix structure of the input matrices A k , for k = 1 , . . . , p .In this paper, we consider the problem of analyzing the structure of thematrix G when A k = T ( a k ), for k = 1 , . . . , p , is a self-adjoint positive definiteToeplitz matrix and G is either the ALM, NBMP, and their weighted counter-parts, or the Karcher mean. We start by considering the case where the symbols a k ( z ), k = 1 , . . . , p , are continuous real valued positive functions.We show that if G is the ALM, NBMP, or a weighted mean of the self-adjointpositive definite operators T ( a ) , . . . , T ( a p ), then, even though the Toeplitzstructure is lost in G , a hidden structure is preserved. In fact, we show that G
2s a quasi-Toeplitz matrix, that is, it can be written as G = T ( g ) + K G , where K G is a self-adjoint compact operator on the set (cid:96) formed by infinite vectors v = ( v i ) i ∈ Z + , v i ∈ C such that (cid:80) ∞ i =1 | v i | < ∞ . An interesting feature is thatthe function g ( z ) turns out to be ( a ( z ) a ( z ) · · · a p ( z )) p , as expected. Theseresults are extended to the more general case A k = T ( a k ) + K k , where K k is acompact operator.Among the side results that we have introduced to prove the main propertiesof the geometric mean, it is interesting to point out Theorem 6 which states thatfor any continuous real valued function a ( z ) and for any continuous function f ( t )defined on the spectrum of T ( a ), the difference f ( T ( a )) − T ( f ( a )) is a compactoperator.Algebras of quasi-Toeplitz matrices are also known in the literature as Toeplitzalgebras . See for instance [13] where in the Example 2.8 it is considered the caseof continuous symbols, or in Section 7.3, where the smallest closed subalgebraof B ( (cid:96) p ), 1 < p < ∞ , containing the matrices of the type T ( a ), where a has anabsolute convergent Fourier series, is considered. Here, B ( (cid:96) p ) denotes the setof bounded operators on the set (cid:96) p of infinite vectors v = ( v i ) i ∈ Z + such that (cid:80) ∞ i =1 | v i | p is finite.After analyzing the case of continuous symbols, we identify additional reg-ularity assumptions on the functions a ( z ) , . . . , a p ( z ), such that the summa-tion (cid:80) + ∞ (cid:96) = −∞ | g (cid:96) | is finite, where (cid:80) + ∞ (cid:96) = −∞ g (cid:96) z (cid:96) is the Fourier series of g ( z ) =( a ( z ) · · · a p ( z )) /p . This property is meaningful since it implies that T ( a k ) aswell as G are bounded operators on the set (cid:96) ∞ of infinite vectors with uniformlybounded components.We will apply these ideas also in the case where the matrices A , . . . , A p , areof finite size n and show that, for a sufficiently large value of n , the matrix G is numerically well approximated by a Toeplitz matrix plus a matrix correctionthat is nonzero in the entries close to the top left corner and to the bottom rightcorner of the matrix.Finally we discuss some computational issues. Namely, we examine classicalalgorithms for computing the square root and the p th root of a matrix as CyclicReduction and Newton’s iteration, and analyze their convergence propertieswhen applied to Toeplitz matrices given in terms of the associated symbols.Then we deal with the problem of computing G . We introduce and analyzealgorithms for computing the Toeplitz part T ( g ) and the correction part K G ina very efficient way, both in the case of infinite and of finite Toeplitz matrices. Aset of numerical experiments shows that the geometric means of infinite quasi-Toeplitz matrices can be easily and accurately computed by relying on simpleMATLAB implementations, based on the CQT-Toolbox of [9].The paper is organized as follows: in the next section, after providing somepreliminary results on Banach algebras and C ∗ -algebras, we recall the mostcommon matrix geometric means, introduce the class of quasi-Toeplitz matricesassociated with a continuous symbol and study functions of matrices in thisclass. Moreover, we analyze the case where the symbol of quasi-Toeplitz ma-trices is in the Wiener class, that is, the associated Fourier series is absolutely3onvergent. In Section 3 we prove that the most common definitions of geo-metric mean of operators preserve the structure of quasi-Toeplitz matrix, andthat the symbol associated with the geometric mean is the geometric mean ofthe symbols associated with the input matrices. In the same section, we giveregularity conditions on the symbols in order that the convergence of the ALMsequence of the symbol holds in the Wiener norm, this implies the convergenceof the operator sequences in the infinity norm. In Section 4 we discuss someissues related to the computation of the p th root and of the geometric meanof quasi-Toeplitz matrices, in Section 5 we describe the geometric means of fi-nite quasi-Toeplitz matrices, while in Section 6 we report the results of somenumerical experiments. Finally, Section 7 draws some conclusions and openproblems. Let i be the complex unit, and denote T = { z ∈ C : | z | = 1 } the unit circle inthe complex plane. Given a function f ( z ) : T → C the composition (cid:101) f ( t ) = f ( e i t )is a 2 π -periodic function from R to C which can be restricted to [0 , π ]. For thesake of simplicity, we keep the notation f ( t ) to denote (cid:101) f ( t ), where it should beunderstood that the t variable is real and ranges in the set [0 , π ], while the z variable is complex and ranges in the set T .For a given subset S ⊂ R , and 1 ≤ p < ∞ , we denote by L p ( S ) the setof functions f ( t ) : S → C , such that (cid:107) f (cid:107) p, S = ( (cid:82) S | f ( t ) | p ) p < ∞ . We de-note L ∞ ( S ) the set of measurable functions defined over S with finite essentialsupremum and we set (cid:107) f (cid:107) ∞ , S = ess sup t ∈S | f ( t ) | . If the set S is clear from thecontext we write (cid:107) f (cid:107) L p in place of (cid:107) f (cid:107) p, S for 1 ≤ p < ∞ and (cid:107) f (cid:107) ∞ in place of (cid:107) f (cid:107) ∞ , S .We denote by C ( S ) the set of continuous functions f : S → C . The subsetof f ∈ C ([0 , π ]) such that f (0) = f (2 π ) is denoted by C π , while L p π = L p ([0 , π )).Let Z + be the set of positive integers and denote by (cid:96) p , with 1 ≤ p < ∞ , the set of infinite sequences v = ( v i ) i ∈ Z + , v i ∈ C , such that (cid:107) v (cid:107) p :=( (cid:80) ∞ i =1 | v i | p ) p is finite, while (cid:96) ∞ is the set of sequences such that (cid:107) v (cid:107) ∞ =sup i | v i | < ∞ . If A = ( a i,j ) i,j =1 ,...,n ∈ C n × n is such that a i,j ∈ C , then thefunction (cid:107) A (cid:107) = max (cid:107) x (cid:107) =1 (cid:107) Ax (cid:107) is the matrix norm induced by the vectornorm (cid:107) x (cid:107) = ( (cid:80) ni =1 | x i | ) . Similarly, if A : (cid:96) p → (cid:96) p is a linear operator, thefunction (cid:107) A (cid:107) p = sup (cid:107) v (cid:107) p =1 (cid:107) Av (cid:107) p , is the operator norm induced by the (cid:96) p norm (cid:107) · (cid:107) p .Recall that any linear operator A : (cid:96) → (cid:96) can be represented by a semi-infinite matrix ( a i,j ) i,j ∈ Z + , which we denote with the same symbol A . The setof bounded operators onto (cid:96) is denoted by B ( (cid:96) ), or more simply by B . Werecall that the set K ⊂ B of compact operators is the closure in B of the set ofbounded operators with finite rank. 4iven the matrix (operator) A = ( a i,j ), we define A ∗ = ( a j,i ) the conjugatetranspose of A , where x denotes the complex conjugate of x ∈ C . A matrix(operator) A is Hermitian or self-adjoint if A = A ∗ , moreover, is positive if inaddition v ∗ Av = (cid:80) i v i ( Av ) i > v ∈ C n \{ } ( v ∈ (cid:96) \{ } ). We say that A is positive definite if there exists a constant γ > v ∗ Av ≥ γ (cid:107) v (cid:107) for any v ∈ (cid:96) , v (cid:54) = 0. C ∗ -algebras A Banach algebra B is a Banach space with a product such that (cid:107) ab (cid:107) ≤ (cid:107) a (cid:107)(cid:107) b (cid:107) ,for a, b ∈ B , where (cid:107) · (cid:107) is the norm of B . Trivial examples of Banach algebrasthat we will use are: the set of essentially bounded functions L ∞ ( S ), with S ⊂ R , with the norm (cid:107) f (cid:107) ∞ ,S and pointwise multiplication (as a special casewe have L ∞ π ); the set C ( K ) of continuous functions on a compact K ⊂ C n withpointwise multiplication; and B ( (cid:96) ) with operator composition. If a ( k ) ∈ B is asequence such that lim k (cid:107) a ( k ) − a (cid:107) = 0 for some a ∈ B , we write a ( k ) → B a .An immediate consequence of the definition of Banach algebra is the follow-ing. Lemma 1.
Let B be a Banach algebra and a , . . . , a p ∈ B . If the sequences { a ( k ) i } k ⊂ B are such that a ( k ) i → B a i for i = 1 , . . . , p , then a ( k )1 · · · a ( k ) p → B a · · · a p . Proof.
Let (cid:107) · (cid:107) be the norm in B . We proceed by induction on p . The case p = 1 is trivial and if the property is true for p − (cid:107) a ( k )1 · · · a ( k ) p − a · · · a p (cid:107)≤ (cid:107) a ( k )1 · · · a ( k ) p − a ( k )1 a · · · a p (cid:107) + (cid:107) a ( k )1 a · · · a p − a · · · a p (cid:107)≤ (cid:107) a ( k )1 (cid:107)(cid:107) a ( k )2 · · · a ( k ) p − a · · · a p (cid:107) + (cid:107) a ( k )1 − a (cid:107)(cid:107) a · · · a p (cid:107) , and the latter tends to zero since (cid:107) a ( k )2 · · · a ( k ) p − a · · · a p (cid:107) tends to zero by induc-tive hypothesis, (cid:107) a ( k )1 − a (cid:107) tends to zero by hypothesis, {(cid:107) a ( k )1 (cid:107)} is uniformlybounded since it converges and (cid:107) a · · · a p (cid:107) ≤ (cid:107) a (cid:107) · · · (cid:107) a p (cid:107) that is bounded.To any element A of a Banach algebra, it is assigned the spectrum sp( A )that is the set of λ ∈ C such that A − λI fails to be invertible. The spectralradius ρ ( A ) is defined as ρ ( A ) = sup λ ∈ sp( A ) | λ | .A C ∗ -algebra is a Banach algebra B with an involution A → A ∗ such that( aS + bT ) ∗ = aS ∗ + bT ∗ ; ( ST ) ∗ = T ∗ S ∗ , ( T ∗ ) ∗ = T and (cid:107) T ∗ T (cid:107) = (cid:107) T (cid:107) , for T, S ∈ B and a, b ∈ C . Two important examples of C ∗ -algebras are B ( (cid:96) )with the adjoint operation and C ( K ), for K ⊂ R compact, with the pointwiseconjugation.A C ∗ -subalgebra F ⊂ B is a closed subalgebra that contains the identityand such that A ∈ F implies A ∗ ∈ F . 5n element A of a C ∗ -algebra is self-adjoint if A ∗ = A , normal if A ∗ A = AA ∗ .We need some results about C ∗ -algebras taken from Proposition 4.1.1, andTheorems 4.1.3 and 4.1.6 of [23]. Lemma 2.
Let A be an element of a C ∗ -algebra B . (i) If A is normal then ρ ( A ) = (cid:107) A (cid:107) . Moreover, if A is self-adjoint, then thespectrum sp( A ) of A is a compact subset of R . (ii) If A is self-adjoint, then there exists a unique continuous mapping f → f ( A ) : C (sp( A )) → B such that f ( A ) has its elementary meaning when f is a polynomial. Moreover, f ( A ) is normal, while it is self-adjoint if andonly if f takes real values at the spectrum of A . (iii) If A is self-adjoint and f ∈ C (sp( A )) , then (cid:107) f ( A ) (cid:107) = (cid:107) f (cid:107) ∞ , sp( A ) and sp( f ( A )) = f (sp( A )) . Lemma 2 implies the existence of a continuous functional calculus , that is todefine a function of a self-adjoint element A of a C ∗ -algebra when the functionis continuous at the spectrum of A . A great effort has been done to define properly the geometric mean of two ormore positive definite matrices. The weighted geometric mean of two n × n matrices A and B , with weight t ∈ [0 , A t B = A / ( A − / BA − / ) t A / , while the geometric mean is G ( A, B ) = A B = A / B . Notice that, usingcontinuous functional calculus, the same formula makes sense also for boundedself-adjoint positive definite operators A, B ∈ B ( (cid:96) ).For more than two matrices, several attempts have been done before theright definition were fully understood. For the ease of the reader, here we limitthe treatise to three positive matrices of the same size, while we refer to a latersection for the case of more than three matrices. The ALM mean, proposedby Ando, Li and Mathias [1], has been introduced as the common limit of thesequences A k +1 = B k C k , B k +1 = C k A k , C k +1 = A k B k , with A = A, B = B, C = C . A nice feature of this mean is that the conver-gence of the sequence can be obtained using the Thompson metric [39] d ( A, B ) = log max { ρ ( A − B ) , ρ ( B − A ) } . It is well known that convergence in the Thompson metric implies convergence inthe operator norm induced by the (cid:96) norm, also in the infinite dimensional case.6his fact allows one to define the ALM mean of self-adjoint positive definiteoperators.A variant of this construction, introduced independently by Bini, Meini,Poloni [10] and Nakamura [34], named NBMP mean is obtained by the sequences A k +1 = A k / ( B k C k ) , B k +1 = B k / ( C k A k ) , C k +1 = C k / ( A k B k ) , with A = A , B = B , C = C . The sequences converge to a common limit,different in general from the ALM mean, and with a faster rate than the ALMsequence. The convergence of these sequences in the Thompson metric allowsone to extend them to operators. Weighted versions of the ALM and NBMPmean can be given by introducing suitable parameters, see Section 3 for moredetails.A well-recognized geometric mean of matrices is the unique positive definitesolution of the matrix equationlog( X − / AX − / ) + log( X − / BX − / ) + log( X − / CX − / ) = 0 , the so-called Karcher mean of A, B and C . Its extension to operators hasbeen more complicated, but it has been shown that the Karcher mean can bedefined also for positive definite self-adjoint bounded operators, proving thatthe equation above has a unique positive definite self-adjoint solution [29]. An integrable and 2 π -periodic function a ∈ L π , with Fourier series (cid:80) ∞ k = −∞ a k e i kt can be associated with the semi-infinite Toeplitz matrix T ( a ) = ( t i,j ) i,j ∈ Z + suchthat t i,j = a j − i . The function a is said to be the symbol associated with theToeplitz matrix T ( a ). A classical result states that T ( a ) represents a boundedoperator on (cid:96) if and only if a ∈ L ∞ π [11, Theorem 1.1]. Moreover, if a iscontinuous, then (cid:107) T ( a ) (cid:107) = (cid:107) a (cid:107) ∞ .Toeplitz matrices form a linear subspace of B ( (cid:96) ), closed by involution since T ( a ) ∗ = T ( a ), but not closed under multiplication, that is, by composition ofoperators in B ( (cid:96) ), so they are not a C ∗ -subalgebra of B ( (cid:96) ). Nevertheless, thereis a nice formula for the product of two Toeplitz matrices [11, Propositions 1.10and 1.11]. Theorem 3.
Let a, b ∈ L ∞ π , (1) T ( a ) T ( b ) = T ( ab ) − H ( a − ) H ( b + ) , where H ( a − ) i,j = ( a − i − j +1 ) i,j ∈ Z + and H ( b + ) = ( b i + j − ) i,j ∈ Z + . Moreover, if a, b are continuous functions then H ( a − ) , H ( b + ) are compact operators on (cid:96) . A nice feature of equation (1) is that it relates the symbol of the product tothe product of the symbols and if a, b ∈ C π then T ( a ) T ( b ) − T ( ab ) belongs tothe set K of compact operators on (cid:96) [11].7he smallest C ∗ -subalgebra of the C ∗ -algebra B ( (cid:96) ) containing all Toeplitzoperators with continuous symbols is [13] QT = { T ( a ) + K : a ∈ C π , K ∈ K} , that we call the set of quasi-Toeplitz matrices and it is also known as Toeplitzalgebra. Notice that for A ∈ QT there exists unique a ∈ C π and K ∈ K suchthat A = T ( a ) + K , because the intersection between Toeplitz matrices withcontinuous symbols and compact operators is the zero operator [11, Proposition1.2]. To define a QT matrix, we will use often the notation A = T ( a )+ K ∈ QT ,without saying explicitly that a is a continuous and 2 π -periodic function and K is a compact operator.The following lemma collects and resumes some known results from [11,Sections 1.4 and 1.5], [13, Section 1.1] and from [16], which will be usefulin the sequel. In particular, it states properties of the spectrum sp( A ) of A , of the essential spectrum sp ess ( A ) defined by sp ess ( A ) = { λ ∈ C : A − λI is not Fredholm on B ( (cid:96) ) } , and on the numerical range of A defined asW( A ) = { x ∗ Ax : (cid:107) x (cid:107) = 1 } . Lemma 4.
The following properties hold: (i) T ( a ) represents a bounded operator on (cid:96) if and only if a ∈ L ∞ π , and T ( a ) is self adjoint if and only if a is real valued.Moreover, if a ∈ C π then: (ii) (cid:107) a (cid:107) ∞ = (cid:107) T ( a ) (cid:107) ≤ (cid:107) T ( a ) + K (cid:107) , for K ∈ K ; (iii) sp( T ( a )) = a ([0 , π ]) ∪ { λ ∈ C : wind( a − λ ) (cid:54) = 0 } is compact; (iv) a ([0 , π ]) = sp ess ( T ( a )) = sp ess ( T ( a ) + K ) ⊂ sp( T ( a ) + K ) for K ∈ K ; (v) sp( A ) is contained in the closure of W( A ) for A ∈ B ( (cid:96) ) . The following result will be useful.
Lemma 5.
Let A = T ( a ) + K ∈ QT . If A = A ∗ then T ( a ) = T ( a ) ∗ and K = K ∗ , moreover, a is real valued. If A is positive, then a is nonnegative. Ifin addition A is positive definite then a is strictly positive.Proof. If A = A ∗ , then T ( a ) − T ( a ) ∗ + K − K ∗ = 0, and by the uniquenessof the decomposition of a quasi-Toeplitz matrix, we have that K = K ∗ and T ( a ) = T ( a ) ∗ , that in turn implies that a k = a − k and thus a ( z ) is real. If A ispositive definite, there exists γ > x ∗ Ax ≥ γ (cid:107) x (cid:107) > x (cid:54) = 0.That is, the numerical range W( A ) of A is formed by real numbers greater thanor equal to γ so that its closure contains positive values. Since by Lemma 4, part(v), the spectrum sp( A ) is contained in the closure of W( A ), we have that λ > λ ∈ sp( A ). From Lemma 4 part (iv), it follows that a ( t ) ∈ sp( A ) so that a ( t ) >
0. The case where A is positive is similarly treated since x ∗ Ax ≥ x ∈ (cid:96) so that the closure of W( A ) is formed by nonnegative numbers.8n view of Lemma 2, the fact that QT is a C ∗ -algebra allows one to usecontinuous functional calculus. If A is a self-adjoint quasi-Toeplitz matrix and f is a function continuous on the spectrum of A , then one can define the normalquasi-Toeplitz matrix f ( A ). In particular, the p -th root of a self-adjoint positivedefinite quasi-Toeplitz matrix turns out to be a quasi-Toeplitz matrix. This factimplies that the sequences generated by the ALM and the NBMP constructions,if the initial values are QT matrices, are formed by entries belonging to QT .We will prove that also the limit G of these sequences belongs to QT , and thatit can be written as G = T ( g ) + K G , where g ( t ) is the geometric mean of thesymbols associated with the Toeplitz part of the given matrices. Similarly, wewill show that the Karcher mean of quasi-Toeplitz matrices is a quasi-Toeplitzmatrix (the latter follows also from [27]).In order to prove these properties we use the following results. Theorem 6.
Let a ∈ C π be a real valued function and f a continuous functionon the spectrum of T ( a ) , then (2) f ( T ( a )) = T ( f ◦ a ) + K, where K is a compact operator. If f takes real values on the spectrum of T ( a ) ,then K is self-adjoint and sp( f ( T ( a ))) = sp( T ( f ◦ a )) .Proof. Since a is real valued, by Lemma 4 the operator T ( a ) is self-adjoint andsince a is continuous, its spectrum Σ is a compact set containing the range of a (see Lemma 4). Moreover, one can define f ( T ( a )), for any function f continuouson Σ, using the continuous functional calculus and, since the function f iscontinuous in Σ, by Lemma 2, we have (cid:107) f ( T ( a )) (cid:107) = (cid:107) f (cid:107) ∞ , Σ .From Theorem 3 it follows that if f is a polynomial, then equation (2) holds.By using an approximation argument, we prove that (2) still holds for anycontinuous function f . Indeed, f can be approximated uniformly by a sequenceof polynomials p n such that (cid:107) f − p n (cid:107) ∞ , Σ tends to 0 (by the Weierstrass theorem)and there exist compact operators K k such that p k ( T ( a )) = T ( p k ◦ a ) + K k . Since continuous functional calculus preserves C ∗ -algebras, then f ( T ( a )) ∈ QT ,that is f ( T ( a )) = T ( b ) + K , what we need to prove is that b = f ◦ a .In order to get the result, in view of the uniqueness of the decompositionof a quasi-Toeplitz matrix, it is sufficient to prove that f ( T ( a )) − T ( f ◦ a ) is acompact operator.We have(3) (cid:107) f ( T ( a )) − T ( f ◦ a ) − K k (cid:107) ≤ (cid:107) f ( T ( a )) − p k ( T ( a )) (cid:107) + (cid:107) T ( f ◦ a ) − T ( p k ◦ a ) (cid:107) ≤ (cid:107) f − p k (cid:107) ∞ , Σ , where the last inequality follows by functional calculus, since f − p k is continuousin the spectrum Σ of T ( a ) and from the property (cid:107) T ( a ) (cid:107) = (cid:107) a (cid:107) ∞ (compareLemma 4). Whence we get (cid:107) T ( f ◦ a − p k ◦ a ) (cid:107) = (cid:107) ( f − p k ) ◦ a (cid:107) ∞ , [0 , π ] = (cid:107) f − p k (cid:107) ∞ ,a ([0 , π ]) ≤ (cid:107) f − p k (cid:107) ∞ , Σ . K k converges to f ( T ( a )) − T ( f ◦ a ) in the operatornorm, and since K is closed, we may conclude that f ( T ( a )) − T ( f ◦ a ) is acompact operator, that is what we wanted to prove.Regarding the last statement, by Lemma 2, part (ii) and Lemma 4, we knowthat f ( T ( a )) and T ( f ◦ a ) are self-adjoint and thus K is self-adjoint as well, byLemma 2, part (iii), and Lemma 4, we getsp( f ( T ( a ))) = f (sp( T ( a ))) = f ( a ([0 , π ])) = sp( T ( f ◦ a )) . By slightly modifying the above proof, we can easily arrive at the followinggeneralization.
Theorem 7.
Let a ∈ C π be a real valued function, H a self-adjoint compactoperator and f a continuous function on the spectrum of T ( a )+ H , then f ( T ( a )+ H ) = T ( f ◦ a ) + K , where K is a compact operator. If f takes real values onthe spectrum of T ( a ) , then K is self-adjoint. An immediate consequence of the above result is the following corollaryrelated to the weighted geometric mean.
Corollary 8.
Let
A, B ∈ QT be positive definite operators associated with thecontinuous symbols a, b , respectively. Then, a, b > and for t ∈ [0 , , we have G = A t B ∈ QT and the symbol g associated with G is such that g = a − t b t .Proof. From Lemma 5 we deduce that if A is positive definite then a, b > A t B = A ( A − BA − ) t A so that, by Theorem 7 one has A , ( A − BA − ) t ∈ QT , thus A t B ∈ QT and g = a − t b t .The next lemma shows that if a sequence { A k } k ⊂ QT converges to A ∈B ( (cid:96) ), then A ∈ QT and the Toeplitz part and the compact part of A k convergeto the Toeplitz part and the compact part of A , respectively. Lemma 9.
Let { A k } k be a sequence of quasi-Toeplitz matrices such that A k = T ( a k ) + K k ∈ QT , where a k ∈ C π , and K k ∈ K , for k ∈ Z + . Let A ∈ B ( (cid:96) ) besuch that lim k (cid:107) A − A k (cid:107) = 0 . Then A ∈ QT , that is A = T ( a )+ K with a ∈ C π and K ∈ K and, moreover, lim k (cid:107) a − a k (cid:107) ∞ = 0 and lim k (cid:107) K − K k (cid:107) = 0 .Proof. Since QT is a C ∗ -algebra, then lim k (cid:107) A − A k (cid:107) = 0 implies that A ∈QT . From the inequality (cid:107) T ( a ) (cid:107) ≤ (cid:107) A (cid:107) (see Lemma 4, part (ii)) and from (cid:107) a (cid:107) ∞ = (cid:107) T ( a ) (cid:107) , it follows that (cid:107) a − a k (cid:107) ∞ = (cid:107) T ( a − a k ) (cid:107) ≤ (cid:107) A − A k (cid:107) tendsto 0 and so does (cid:107) K − K k (cid:107) . If a ( z ) = (cid:80) + ∞ (cid:96) = −∞ a (cid:96) z (cid:96) is a continuous function then lim n | a n | = 0, therefore,for an error bound ε > S ε = { n ∈ Z : | a n | < ε } is10nite. This fact allows one to numerically approximate the function a ( z ) to anyprecision in a finite number of operations.Classical results of harmonic analysis relate the regularity of a ( z ) with thedecay of the Fourier coefficients to zero. A better regularity implies a fasterconvergence to zero of { a n } n and, in practice, this allows one to approximate a ( z ) by using fewer coefficients.We will consider Toeplitz matrices associated with functions having an an-alytic or an absolutely convergent Fourier series. The former situation is idealsince it implies an exponential decay to zero of { a n } n .We point out that if the function a is just continuous then (cid:80) + ∞ i = −∞ | a i | isnot necessarily finite so that (cid:107) T ( a ) (cid:107) ∞ = (cid:80) + ∞ i = −∞ | a i | (compare [12, Theorem1.14]) is not bounded in general. This is a reason to determine conditions underwhich a function f ( a ) for a ∈ C π , as well as the geometric mean of a i ∈ C π , i = 1 , . . . , p , have absolutely summable coefficients or are analytic.Another important related issue is to find out under which conditions asequence of matrices A k ∈ QT such that (cid:107) A k (cid:107) ∞ < ∞ converges in the infinitynorm to a limit A ∈ QT such that (cid:107) A (cid:107) ∞ < ∞ . In this section we provide toolsto give an answer to these questions.The set of functions a ∈ L π with absolutely convergent Fourier series,namely W = { f ∈ L π : f = (cid:88) k ∈ Z a k e i kt , (cid:88) k ∈ Z | a k | < ∞} , with the norm (cid:107) f (cid:107) W := (cid:80) k ∈ Z | a k | , is a Banach algebra, also known as theWiener algebra [11]. A function a ∈ W is necessarily continuous, but it mayfail to be differentiable, or even Lipschitz continuous.We consider the convergence of a sequence in three broad classes of functionsincluded in the Wiener algebra.The first class is the set of α -H¨older continuous functions. A continuousfunction ϕ : Ω → C , with Ω subset of R n or C n , is said to be α -H¨older continuouson Ω, with 0 < α ≤
1, if [ ϕ ] α := sup x,y ∈ Ω x (cid:54) = y | ϕ ( x ) − ϕ ( y ) || x − y | α is finite. We denote by C ,α π the set of α -H¨older continuous functions on R withperiod 2 π , that is a Banach algebra with the norm (cid:107) ϕ (cid:107) C ,α := (cid:107) ϕ (cid:107) ∞ + [ ϕ ] α . The second class is the set of functions of bounded variation on [0 , π ], thatis functions f : [0 , π ] → C , such that V ( f ) = sup N =1 , ,... x =0 ≤ x ≤···≤ x N =2 π (cid:110) N − (cid:88) i =0 | f ( x i +1 ) − f ( x i ) | (cid:111)
11s finite.The last class is the set of absolutely continuous functions on C π , that isfunctions f : [0 , π ] → R such that f (2 π ) = f (0) and for any ε >
0, there exists δ > n (cid:88) i =1 | f ( y i ) − f ( x i ) | < ε, whenever { [ x i , y i ] : i = 1 , . . . , n } is a finite collection of mutually disjoint subin-tervals of [0 , π ] with (cid:80) ni =1 ( y i − x i ) < δ .Relations between W and these three classes are stated in the followingsummary of well-known results, that gives estimates of the norm (cid:107) · (cid:107) W . Theorem 10.
Let a ∈ C π , we have (i) If a ∈ C ,α π , with α ∈ (1 / , , then a ∈ W . Moreover, there exists aconstant γ a such that (4) (cid:107) a (cid:107) W ≤ γ a (cid:107) a (cid:107) C ,α . (ii) If a ∈ C ,α π , with α ∈ (0 , and of bounded variation V ( a ) on [0 , π ] , then a ∈ W and there exists a constant γ b such that (5) (cid:107) a (cid:107) W ≤ γ b (cid:0) | a | + V ( a ) / ∞ (cid:88) (cid:96) =1 ( π − (cid:96) ) α/ (cid:1) . (iii) If a ∈ C π is absolutely continuous on [0 , π ] and a (cid:48) ∈ L π , then a ∈ W and there exists a constant γ c such that (6) (cid:107) a (cid:107) W ≤ γ c (cid:0) (cid:107) a (cid:107) L + (cid:107) a (cid:48) (cid:107) L (cid:1) . Proof.
Concerning part (i), the statement about α -H¨older continuous functionswith α ∈ (1 / ,
1] is a classical result by Bernstein [3], while the complete proofwhen 0 < α < α = 1 can be provedanalogously.Concerning part (iii), the statement about the absolute continuous functionwith derivative in L π follows from [24, Theorem I.6.2].It is left to prove part (ii). If a ∈ C ,α π with α ∈ (0 , , π ], then by [24, Theorem I.6.4] (see also [46, p.241]), it followsthat a ∈ W . We just need to prove the inequality (5).Denote by ω ( δ ) = ω ( δ, a ) = sup | a ( t ) − a ( t ) | for t , t ∈ [0 , π ] and | t − t | ≤ δ . Applying the technique as in [46, p.241-242], we get ∞ (cid:88) k =1 | a k | ≤ ∞ (cid:88) (cid:96) =1 (cid:0) ω ( π − (cid:96) ) (cid:1) / V ( a ) / ≤ √ γ ∞ (cid:88) (cid:96) =1 ( π − (cid:96) ) α/ V ( a ) / , ω ( δ ) ≤ r δ α for a constant r independentof δ if a ∈ C ,α π . It can be proved analogously that ∞ (cid:88) k =1 | a − k | ≤ √ γ ∞ (cid:88) (cid:96) =1 ( π − (cid:96) ) α/ V ( a ) / . The proof is completed by choosing γ b = max { , √ γ } .Note that, in particular, Theorem 10 implies that a Lipschitz continuousfunction with period 2 π belongs to W .We are interested in fractional powers and means of functions in the Wieneralgebra. In the case of analytic functions there is an interesting result due toL´evy [31]. Theorem 11.
Let a ∈ W and let f be a complex function, analytic in the rangeof a , then f ◦ a ∈ W . Clearly, if both functions a and f are analytic then f ◦ a is analytic. Asimpler result that we will use in the following is related to H¨older continuousfunctions. Lemma 12.
Let a : Ω → R be α -H¨older continuous and let f ∈ C ( I ) where I is a closed interval containing the range of a . Then f ◦ a is α -H¨older continuousand [ f ◦ a ] α ≤ (cid:107) f (cid:48) (cid:107) ∞ , I [ a ] α . Moreover (cid:107) f ◦ a (cid:107) C ,α ≤ γ ( (cid:107) f (cid:107) ∞ , I + (cid:107) f (cid:48) (cid:107) ∞ , I ) , where γ = max { [ a ] α , } .Proof. If x, y ∈ Ω are such that a ( x ) (cid:54) = a ( y ), then, by the mean value theorem, | f ( a ( x )) − f ( a ( y )) || x − y | α = | f ( a ( x )) − f ( a ( y )) || a ( x ) − a ( y ) | | a ( x ) − a ( y ) || x − y | α ≤ (cid:107) f (cid:48) (cid:107) ∞ , I [ a ] α . The bound holds also when a ( x ) = a ( y ) and we have [ f ◦ a ] α ≤ (cid:107) f (cid:48) (cid:107) ∞ , I [ a ] α .The latter inequality, follows from (cid:107) f ◦ a (cid:107) C ,α = (cid:107) f ◦ a (cid:107) ∞ + [ f ◦ a ] α .As a consequence we have. Corollary 13. If a ( t ) ∈ W is such that a ( t ) > then a ( t ) /p ∈ W . If in addi-tion a ( t ) ∈ C ,α π for some α > then a /p ∈ C ,α π and [ a /p ] α ≤ p [ a ] α min a ( t ) − /p .Proof. Since f ( z ) = z /p is analytic in the range of a , by Theorem 11, a /p ∈ W .Let I = [min a ( t ) , max a ( t )] be a closed interval of positive numbers enclosingthe range of a , then f ∈ C ( I ) and we can apply Lemma 12, using (cid:107) f (cid:48) (cid:107) ∞ , I = p min a ( t ) − /p .Similarly, we have the following result.13 emma 14. Let a : [0 , π ] → R be a real valued function and let f ∈ C ( I ) where I is a closed interval containing the range of a . (i) If a is absolutely continuous and a (cid:48) ∈ L π , then f ◦ a is absolutely contin-uous and with derivative in L π . Moreover, (cid:107) ( f ◦ a ) (cid:48) (cid:107) L ≤ (cid:107) f (cid:48) (cid:107) ∞ , I (cid:107) a (cid:48) (cid:107) L . (ii) If a is of bounded variation on [0 , π ] , then f ◦ a is of bounded variationon [0 , π ] and, moreover, V ( f ◦ a ) ≤ (cid:107) f (cid:48) (cid:107) ∞ , I V ( a ) .Proof. Concerning part (i), it follows from [2, Theorem 5.10, part (d)] that f ◦ a is absolutely continuous. Observe that ( f ◦ a ) (cid:48) = ( f (cid:48) ◦ a ) a (cid:48) and (cid:16)(cid:90) π | ( f ◦ a ) (cid:48) | dt (cid:17) / ≤ (cid:16)(cid:90) π (cid:107) f (cid:48) ◦ a (cid:107) ∞ | a (cid:48) | dt (cid:17) / = (cid:107) f (cid:48) (cid:107) ∞ , I (cid:107) a (cid:48) (cid:107) L < ∞ , which implies that ( f ◦ a ) (cid:48) ∈ L π and (cid:107) ( f ◦ a ) (cid:48) (cid:107) L ≤ (cid:107) f (cid:48) (cid:107) ∞ , I (cid:107) a (cid:48) (cid:107) L .Concerning part (ii), a direct consequence of [2, Theorem 5.10, part (e)]shows that f ◦ a is of bounded variation.Consider a partition x = 0 ≤ x ≤ . . . ≤ x (cid:96) = 2 π , we have (cid:96) − (cid:88) j =0 | f ( a ( x j +1 )) − f ( a ( x j )) | = (cid:96) − (cid:88) j =0 a ( xj ) (cid:54) = a ( xj +1) | f ( a ( x j +1 )) − f ( a ( x j )) || a ( x j +1 ) − a ( x j ) | | a ( x j +1 ) − a ( x j ) |≤ sup x,y ∈I x (cid:54) = y | f ( x ) − f ( y ) || x − y | (cid:96) − (cid:88) j =0 | a ( x j +1 ) − a ( x j ) | ≤ (cid:107) f (cid:48) (cid:107) ∞ , I V ( a ) . Taking the supremum over all partitions, we get the desired inequality. QT matrices We start this section by recalling a general construction of the ALM mean andthe NBMP mean of p ( p ≥
3) positive definite operators.Denote by G t ( A, B ), the weighted geometric mean A t B , with weight t ∈ [0 , p − s , s , . . . , s p − ) with s i ∈ [0 ,
1] and positivedefinite operators A , . . . , A p , the sequences generated by(7) A ( k +1) i = A ( k ) i s G s ,...,s p − ( A ( k )1 , . . . , A ( k ) i − , A ( k ) i +1 , . . . , A ( k ) p ) , i = 1 , . . . , p, with A (0) i = A i , can be recursively defined, and they converge to a commonlimit G s ,...,s p − [30].With the choice ( s , . . . , s p − , s p − ) = (1 , . . . , , /
2) one obtains the ALMmean [1] and with the choice ( s , s , . . . , s p − ) = (( p − /p, ( p − / ( p − , . . . , / w = ( w i ) ∈ R p , i.e., such that w i > (cid:80) pi =1 w i = 1, define thefollowing sequence(8) A ( k +1) i = A ( k ) i − w i G ˆ w ( i ) ( A ( k )1 , . . . , A ( k ) i − , A ( k ) i +1 , . . . , A ( k ) p ) , i = 1 , . . . , p, where ˆ w ( i ) = − w i ( w , . . . , w i − , w i +1 , . . . , w p ) is again a probability vector, and G w ,w ( A , A ) = A w A . If lim k A ( k ) i exists and has the same value for every i , then we denote the common limit by G w ( A , . . . , A p ) and refer to it as theweighted mean. It is proved in [28] that this limit exists in the Thompson metric.Observe that by choosing w = p (1 , , . . . , A , . . . , A p such that A i = T ( a i ) + E i ∈ QT are also QT matrices. In fact, we prove that thesequences { A ( k ) i } k generated by (7) converge to a common limit G ∈ QT withsymbol g = ( a . . . a p ) p and we have lim k (cid:107) a ( k ) i − g (cid:107) ∞ = 0, where a ( k ) i is thesymbol of A ( k ) i . In the case where the symbols a i ∈ W , i = 1 , . . . , p , we providesufficient conditions under which a ( k ) i converges to g in Wiener norm. The ALM sequences { A ( k ) i } ∞ k =0 generated by (7) with ( s , . . . , s p − , s p − ) =(1 , . . . , , / G in the Thompson metric (see [30,Remarks 4.2 and 6.5 and Theorem 4.3]). This implies that lim k (cid:107) A ( k ) i − G (cid:107) = 0since the topology of the Thompson metric agrees with the relative operatornorm topology [39].We will prove that if the positive definite matrices A , . . . , A p , p ≥
3, belongto QT then also the matrices A ( k ) i of the ALM sequence generated by (7) aswell as their limit G belong to QT . Moreover, the symbol associated with theToeplitz part of G is the uniform limit of the symbols a ( k ) i associated with theToeplitz parts of A ( k ) i , which in turn are the functions obtained by applying theALM construction to the symbols a i associated with the Toeplitz parts of thematrices A i . Theorem 15.
Let A i = T ( a i ) + E i ∈ QT be positive definite, for i = 1 , . . . , p ,with p ≥ . The matrices A ( k ) i generated by (7) for the ALM iteration, and theALM mean G = G ( A , . . . , A p ) of A , . . . , A p , satisfy the following properties:1. for any k ≥ and for i = 1 , . . . , p , there exist a ( k ) i ∈ C π and K ( k ) i ∈ K ( (cid:96) ) such that A ( k ) i = T ( a ( k ) i ) + K ( k ) i , that is A ( k ) i ∈ QT ;2. there exist g ∈ C π and K G ∈ K ( (cid:96) ) such that G = T ( g ) + K G , that is, G ∈ QT ; . lim k (cid:107) a ( k ) i − g (cid:107) ∞ = 0 , lim k (cid:107) K ( k ) i − K G (cid:107) = 0 , for i = 1 , . . . , p ;4. the equation a ( k +1) i = a ( k ) i s G s ,...,s p − ( a ( k )1 , . . . , a ( k ) i − , a ( k ) i +1 , . . . , a ( k ) p ) , issatisfied for i = 1 , . . . , p , and for any k . Proof.
It is known from [30] that the sequences { A ( k ) i } ∞ k =0 converge to G inthe Thompson metric and that lim k (cid:107) A ( k ) i − G (cid:107) = 0. We prove parts 1–4 byinduction on p . If p = 3, then(9) A ( k +1)1 = A ( k )1 s G ( A ( k )2 , A ( k )3 ) ,A ( k +1)2 = A ( k )2 s G ( A ( k )3 , A ( k )1 ) ,A ( k +1)3 = A ( k )3 s G ( A ( k )1 , A ( k )2 ) , Recall that if
A, B ∈ QT then the geometric mean G ( A, B ) ∈ QT , and thesymbols a, b, g associated with A, B and G ( A, B ), respectively, are such that g = G ( a, b ) in view of Corollary 8.Using an induction argument on k , we show part 1 of the theorem i.e., A ( k ) i ∈ QT for i ∈ { , , } . We have A (0) i = A i ∈ QT , and assuming A ( k ) i ∈ QT ,for i ∈ { , , } , from (9) and Corollary 8, we deduce that A ( k +1) i ∈ QT , for i ∈ { , , } . Consequently, since lim k (cid:107) A ( k ) i − G (cid:107) = 0, for G = G ( A , A , A ),then from Lemma 9 we deduce part 2, i.e., G ∈ QT and that the symbolassociated with the Toeplitz part of A ( k ) i converges to g uniformly and thatthe compact part of A ( k ) i converges to the compact part of G in norm, i.e.,part 3. Finally, since the symbol associated with A (0) i is a i ( z ), we find thatthe symbol associated with G ( A ( k ) i , A ( k ) j ) is G ( a ( k ) i , a ( k ) j ) in view of Corollary 8.Therefore, from (9) and Corollary 8 we find that a ( k +1) i = a ( k ) i s G ( a ( k ) i − , a ( k ) i +1 )with a ( k )0 := a ( k )3 and a ( k )4 := a ( k )1 . That is, part 4.For the inductive step on p , assume p > p − A i ∈ QT , i = 1 , . . . , p −
1, then they also hold for the sequencegenerated by (7) starting from p matrices A i ∈ QT , i = 1 , . . . , p .To this end, consider equation (7) and use induction on k to prove that A ( k ) i ∈ QT for i = 1 , . . . , p . For k = 0, clearly A (0) i = A i ∈ QT by assumption.Concerning the inductive step on k , assume that A ( k ) i ∈ QT for i = 1 , . . . , p ,and deduce that A ( k +1) i ∈ QT for i = 1 , . . . , p . By the inductive assumptionon k we have A ( k ) i ∈ QT so that by the inductive assumption on p , the matrix G s ,...,s p − ( A ( k )1 , . . . , A ( k ) i − , A ( k ) i +1 , . . . , A ( k ) p ) = G ( A ( k )1 . . . , A ( k ) i − , A ( k ) i +1 , . . . , A ( k ) p ) be-longs to QT . Therefore, in view of Corollary 8 and (7) also A ( k +1) i is in QT .That is part 1 of the theorem. Moreover, since lim k (cid:107) A ( k ) i − G (cid:107) = 0 for Here and in the proof of the theorem, we have s = 1 and thus the notation could besimplified, but we prefer to keep s to let the proof be used also to deal with the NBMP andweighted means. = G ( A , . . . , A p ), then from Lemma 9 we deduce that G ∈ QT and thatthe symbol associated with the Toeplitz part of A ( k ) i uniformly converges to thesymbol g associated with the Toeplitz part of G , and that the compact partof A ( k ) i converges to the compact part of G in norm, i.e., parts 2 and 3 of thetheorem.Concerning part 4, we proceed by induction on p . We have already provedthat for p = 3 the property is satisfied. In order to prove the inductive stepon p we proceed by induction on k . For the initial step, i.e., for k = 0, byCorollary 8 and by the inductive assumption on p , from (9) we find that, a (1) i = a (0) i s G s ,...,s p − ( a (0)1 , . . . , a (0) i − , a (0) i +1 , . . . , a (0) p ). For the induction step on k ,assume that part 4 is satisfied for k and prove it for k + 1. By the inductivehypotheses valid for p − G s ,...,s p − ( A ( k )1 , . . . , A ( k ) i − , A ( k ) i +1 , . . . , A ( k ) p ) ∈ QT and that its symbol is G s ,...,s p − ( a ( k )1 , . . . , a ( k ) i − , a ( k ) i +1 , . . . , a ( k ) p ). Therefore, byCorollary 8 and from (7), the symbol associated with the Toeplitz part of A ( k +1) i is a ( k +1) i = a ( k ) i s G s ,...,s p − ( a ( k )1 , . . . , a ( k ) i − , a ( k ) i +1 , . . . , a ( k ) p ).As a consequence of the above theorem we will show that the symbol g associated with the Toeplitz part of G is such that g ( z ) = ( a ( z ) . . . a p ( z )) p ,where the symbols a i ( z ) take positive values in view of Lemma 5 since A i arepositive definite. In order to prove this representation of g ( z ), consider thesequences a ( k ) i ( z ) defined by the ALM iteration, that is(10) a ( k +1) i = a ( k ) i s G s ,...,s p − ( a ( k )1 , . . . , a ( k ) i − , a ( k ) i +1 , . . . , a ( k ) p ) ,a (0) i ( z ) = a i ( z ) , for i = 1 , . . . , p , k ≥
0. It can be easily verified that(11) a ( k +1) i ( z ) = p (cid:89) j =1 , j (cid:54) = i a ( k ) j ( z ) p − , i = 1 , . . . , p, k = 0 , , . . . We have the following.
Lemma 16.
Let a ( z ) , . . . , a p ( z ) be continuous nonnegative functions and let a ( k ) i , for i = 1 , . . . , p and k = 0 , , . . . , be the sequences defined by the ALMiteration (10) . For k = 0 , , . . . , we have a ( k ) i = a n k − i p (cid:89) j =1 , j (cid:54) = i a n k j , i = 1 , . . . , p, where n k = p (cid:0) ( − k +1 ( p − k (cid:1) . roof. We proceed by induction on k . The case k = 0, follows from n = 0 and n − = 1. For k >
0, from (11), we obtain for the sequences { n k } k the differenceequation ( p − n k +1 = ( p − n k + n k − , with n = 0 and n − = 1, whosesolution is n k = p (1 − ( − (1 / ( p − k )).Observe that lim k n k = p so that lim k a ( k ) i ( z ) = ( (cid:81) pi =1 a i ( z )) p pointwise asexpected. On the other hand, from Theorem 15 it follows that convergence isuniform. We may conclude with the following. Theorem 17. If A i = T ( a i ) + K i ∈ QT , for i = 1 , . . . , p, are positive definiteoperators, then the symbol g ( t ) associated with G = G ( A , . . . , A p ) is such that g ( t ) = ( a ( t ) · · · a p ( t )) p . Now, consider the case where A , A , . . . , A p ∈ QT , are such that theirassociated symbols a i ∈ W and E i ∈ K ( (cid:96) ) for i = 1 , . . . , p . It is clear thatthe sequences { A ( k ) i } generated by the ALM iteration converge to a matrix G ∈ QT with symbol g = ( a · · · a p ) p . Since W is a Banach algebra, we have g = ( a · · · a p ) p ∈ W . Concerning the convergence of the symbols { a ( k ) i } of thesequences { A ( k ) i } , we know that lim k (cid:107) a ( k ) i − g (cid:107) ∞ = 0, but uniform convergencedoes not imply that lim k (cid:107) a ( k ) i − g (cid:107) W = 0 (see [24, Page 34]).Now we show that under some regularity conditions, the sequences { a ( k ) i } k converge to g in Wiener norm, i.e., lim k (cid:107) a ( k ) i − g (cid:107) W = 0. Theorem 18.
Let a , . . . , a p ∈ C π be the symbols of the positive definitematrices A , . . . , A p ∈ QT , and let a ( k ) i , for i = 1 , . . . , p and k = 0 , , , . . . , bethe sequence obtained by the ALM iteration (10) . If one of the conditions(a) a i ∈ C ,α π , with α ∈ ( , ;(b) a i ∈ C ,α π , with α ∈ (0 , and of bounded variation;(c) a i absolutely continuous and with derivative in L π ;for i = 1 , . . . , p , is fulfilled, then a ( k ) i ∈ W , ( a · · · a p ) /p ∈ W and a ( k ) i → W ( a · · · a p ) /p .Proof. Observe that, in all three cases, Theorem 10 implies that a , . . . , a p ∈ W ,so that a ( k ) i , a /pi and ( a · · · a p ) /p belong to W in view of Theorem 11. To show a ( k ) i → W ( a · · · a p ) /p , we can see from Lemmas 1 and 16 that it suffices to show a n k i → W a /pi , where n k is defined in Lemma 16 and lim k n k = p .With the notation of Lemma 16, we have a n k i − a /pi = f k ◦ a i , where f k ( x ) := x n k − x /p = x n k (1 − x ( − / ( p − k )is a sequence of analytic functions on (0 , ∞ ). Observe that a i , for i = 1 , . . . , p ,is a strictly positive function in view of Lemma 5, let I be the set of the range18f a i , then I ⊂ (0 , ∞ ) is a closed interval and the sequences { f k } k and { f (cid:48) k } k converge to 0 uniformly on I by [25, Theorem 1.2]. We show that (cid:107) a n k i − a p i (cid:107) W = (cid:107) f k ◦ a i (cid:107) W → (cid:107) f k ◦ a i (cid:107) W ≤ γ a (cid:107) f k ◦ a i (cid:107) C ,α ≤ γ a γ ( (cid:107) f k (cid:107) ∞ , I + (cid:107) f (cid:48) k (cid:107) ∞ , I ) → . In case (b), by Lemmas 12 and 14, f k ◦ a i ∈ C ,α π and is of bounded variationwith V ( f k ◦ a i ) ≤ (cid:107) f (cid:48) k (cid:107) ∞ , I V ( a i ) → f k ◦ a i )( t ) := h k ( t ) = (cid:80) j ∈ Z h ( k ) j e i jt , then (cid:107) h k (cid:107) ∞ = (cid:107) f k (cid:107) ∞ , I → | h ( k )0 | ≤ π (cid:82) π | h k ( t ) | dt ≤ (cid:107) h k (cid:107) ∞ so that h ( k )0 → k → ∞ .Together with Theorem 10 and (5), it yields (cid:107) f k ◦ a i (cid:107) W ≤ γ b (cid:16) | h ( k )0 | + V ( f k ◦ a i ) ∞ (cid:88) v =1 ( π − v ) α (cid:17) → . In case (c), we use part (i) of Lemma 14 and (6) to get (cid:107) f k ◦ a i (cid:107) W ≤ γ c ( (cid:107) f k ◦ a i (cid:107) L + (cid:107) ( f k ◦ a i ) (cid:48) (cid:107) L ) ≤ γ c ( (cid:107) f k (cid:107) ∞ , I + (cid:107) f (cid:48) k (cid:107) ∞ , I (cid:107) a (cid:48) i (cid:107) L ) → . The analysis performed in the previous section can be repeated here concerningthe NBMP mean. In fact, since the p sequences { A ( k ) i } k =1 , ,... , i = 1 , . . . , p ,generated by (7) for the NBMP iteration converge to a common limit G in theThompson metric [30], then they converge in the operator norm. Following ananalysis similar to the one of Section 3.1, one can see that the NBMP meanof positive definite QT matrices is a QT matrix. Moreover, since the NBMPconstruction applied to scalars converges in just one step, then we have thatfor the symbols a ( k ) i obtained this way it holds that a ( k ) i ( z ) = g ( z ) for k ≥ i = 1 , . . . , p . We may conclude with the following results which can be provedby adapting the proof of Theorem 15. Theorem 19.
Let A i = T ( a i ) + E i ∈ QT , for i = 1 , . . . , p , p ≥ , be positivedefinite. Then the matrices A ( k ) i generated by (7) for the NBMP iteration,and the NBMP mean G = G ( A , . . . , A p ) of A , . . . , A p , satisfy the followingproperties:1. A ( k ) i = T ( g ) + K ( k ) i ∈ QT , i = 1 , . . . , p , for any k ≥ , where g =( a · · · a p ) p ;2. G = T ( g ) + K G ∈ QT ; . lim k (cid:107) K ( k ) i − K G (cid:107) = 0 , for i = 1 , . . . , p . A similar argument can be used for the weighted mean: the p sequences { A ( k ) i } k =1 , ,... , i = 1 , . . . , p , generated by (8) converge to their limit G in theThompson metric [28], and thus then they converge in the operator norm, andas before, the weighted mean of positive definite QT matrices is a QT matrix.The symbols a ( k ) i obtained with this procedure are such that, for k >
1, we have a ( k +1) i ( z ) = a ( k ) i − w i G ˆ w ( i ) ( a ( k )1 , . . . , a ( k ) i − , a ( k ) i +1 , . . . , a ( k ) p ) for i = 1 , . . . , p . Wecan show that a ( k ) i = a w a w · · · a w p p = G w ( a , . . . , a p ) for k ≥ Theorem 20.
Let A i = T ( a i ) + E i ∈ QT , for i = 1 , . . . , p , p ≥ , bepositive definite. Let w = ( w i ) ∈ R p be a probability vector, and A ( k ) i bethe matrix sequences generated by (8) for the weighted iteration. Finally, let G = G w ( A , . . . , A p ) be the weighted mean of A , . . . , A p . Then we have1. A ( k ) i = T ( g ) + K ( k ) i ∈ QT , i = 1 , . . . , p , for any k ≥ , where g = a w a w · · · a w p p ;2. G = T ( g ) + K G ∈ QT ;3. lim k (cid:107) K ( k ) i − K G (cid:107) = 0 , for i = 1 , . . . , p . Let B be a C ∗ -algebra, it has been proved in [27] that the equation(12) p (cid:88) i =1 log( X − / A i X − / ) = 0 , where A , . . . , A p ∈ B has a strictly positive definite solution X ∈ B , unique if B is the C ∗ -algebra of bounded operators over an Hilbert space.This implies that the Karcher mean of quasi-Toeplitz matrices exists and itis unique. Theorem 21.
Let T ( a i ) + K i ∈ QT , with a i > for i = 1 , . . . , p . There existsa unique solution X of the equation (12) and X = T ( g ) + K G ∈ QT where g = ( a · · · a p ) /p .Proof. Since QT is a C ∗ -subalgebra of B ( (cid:96) ) then equation (12) has a uniquesolution X = T ( g ) + K G in QT .Observe that by Theorem 6, and by the elementary arithmetic of quasi-Toeplitz matrices, we have that X − / = T ( g − / ) + K , X − / A i X − / = T ( a i /g ) + K ,i , log( X − / A i X − / ) = T (log( a i /g )) + K ,i and finally0 = p (cid:88) i =1 log( X − / A i X − / ) = T (log( a · · · a p /g p )) + K , K , K ,i , K ,i , K ∈ K , for i = 1 , . . . , p .By the uniqueness of decomposition of quasi-Toeplitz matrices, it followsthat g = ( a · · · a p ) /p . In this section we discuss some issues concerning the effective computation ofthe geometric mean of A , . . . , A p ∈ QT . We rely on the CQT-Toolbox [9] forcomputations in the QT algebra but in order to compute geometric means, weneed to compute some fundamental functions of QT matrices, namely, the p -throot and in particular the square root, that we will discuss in the following.In this section, without loss of generality, we assume that the matrix A = T ( a ) + E A ∈ QT is such that 0 ≤ a ( z ) ≤
1. This condition is satisfied inparticular if A is positive and (cid:107) A (cid:107) ≤ (cid:107) a ( z ) (cid:107) ∞ = (cid:107) T ( a ) (cid:107) ≤ (cid:107) A (cid:107) . The square root has been implemented in [9] in two different ways relying on theDenman and Beavers algorithm, and on the Cyclic Reduction (CR) algorithm,respectively. While the two algorithms are equivalent to the Newton methodfor a matrix A , if the initial value commutes with A , in finite arithmetic theirbehavior differs (see [15, Section 6.3]). In our numerical tests, the CR algorithmprovided better numerical results and thus it appears to be better suited for QT matrices.The CR algorithm for the square root is defined as follows(13) Y k +1 = − Y k W − k Y k , Y = I − A,W k +1 = W k + 2 Y k +1 , W = 2( I + A ) , where we assume that all the matrices W k are invertible. In the finite dimen-sional case it follows that lim k W k = A , lim k Y k = 0, where convergence holdsin any operator norm.The sequences obtained by the CR algorithm are related to the sequencesobtained by the simplified Newton method(14) X k +1 = 12 ( X k + AX − k ) , X = 12 ( I + A ) . Indeed, W k = 4 X k and Y k = 2( X k − X k − ), with X − = A (see [17]).We define the sequence { x k ( x ) } k of rational functions of the variable x , bymeans of x k +1 ( x ) = ( x k ( x ) + xx − k ( x )), for k = 0 , , . . . , with x ( x ) = (1 + x ).Since the sequence x k ( x ) is obtained by applying the Newton method to theequation z = x , customary arguments show that for x ∈ [0 ,
1] the sequence { x k ( x ) } k =0 , ,... is well-defined and monotonically decreasing to √ x . In view ofDini’s theorem we may conclude that convergence is uniform on [0 , X k = x k ( A ) we deduce the following.21 orollary 22. Let A = T ( a ) + K ∈ QT be self-adjoint and such that a is realvalued and v T Av ≥ for any v ∈ (cid:96) , (cid:107) v (cid:107) = 1 . Then the sequences W k and Y k generated by (13) are such that lim k (cid:107) Y k (cid:107) = 0 , lim k (cid:107) W k − A (cid:107) = 0 .Proof. By scaling A , we can assume that v T Av ≤
1. The spectrum sp( A ) of A is a compact set contained in [0 , { v T Av : v T v = 1 } and this closure is contained in [0 , { X k } k be the sequence obtained by (14) and x k ( x ) the corresponding scalarsequence. The function ϕ k ( x ) := x k ( x ) − √ x is continuous in [0 ,
1] and byLemma 2 we have (cid:107) X k − A / (cid:107) = (cid:107) ϕ k ( X k ) (cid:107) = (cid:107) ϕ k ( x ) (cid:107) ∞ , sp( A ) . Since sp( A ) ⊂ [0 , k (cid:107) X k − A / (cid:107) = 0 and using W k = 4 X k and Y k = 2( X k − X k − ) we obtain theproof.From the above results it follows that the symbols associated with the ma-trices in (13) uniformly converge to the symbols of their limit. Moreover thecompact corrections associated with these matrix sequences converge in the B ( (cid:96) ) norm to the compact corrections of their limit.A better convergence of the symbols can be obtained under the assumptionof more regularity of a ( z ) in view of Theorems 10 and 11. p -th root It is well known that if A is an n × n Hermitian matrix, where n is finite,with eigenvalues in [0 ,
1] then the sequence generated by Newton’s iteration Y k +1 = p Y k (( p − I + Y − pk A ), Y = I , is such that lim k Y k = A p [15]. However,this iteration may encounter stability problems when implemented in floatingpoint arithmetic. A stable version is based on the following iteration(15) Y k +1 = Y k ( ( p − I + M k p ) , Y = I,M k +1 = ( ( p − I + M k p ) − p M k , M = A. We prove that the iteration (15) applied to a QT matrix converges in norm. Tothis end we follow the same approach used for the square root. More precisely,we introduce the functional sequences(16) y k +1 ( x ) = y k ( x )( p − m k ( x ) p ) , y ( x ) = 1 ,m k +1 ( x ) = ( p − m k ( x ) p ) − p m k ( x ) , m ( x ) = x. so that we have Y k = y k ( A ), M k = m k ( A ). Now we prove the following result. Theorem 23.
For any x ∈ [0 , , for the sequences (16) we have1. ≤ m k ( x ) ≤ m k +1 ( x ) ≤ , . x p ≤ y k +1 ( x ) ≤ y k ( x ) ,3. lim k (cid:107) m k − (cid:107) ∞ = 0 , lim k (cid:107) y k − x p (cid:107) ∞ = 0 ,where inequalities are strict if x (cid:54) = 0 , .Proof. Concerning Part 1, we prove by induction that 0 ≤ m k ( x ) ≤
1. Clearly,since m ( x ) = x we have 0 ≤ m ( x ) ≤
1. For the inductive step, we observe thatthe function f ( t ) = t ( p − tp ) − p such that m k +1 = f ( m k ), maps monotonicallythe interval [0 ,
1] into itself, so that 0 ≤ m k ≤ ≤ m k +1 ≤
1. Similarly,for the inequality m k ≤ m k +1 we have m ≤ m and the monotonicity of f ( t )inductively implies that m k ≤ m k +1 . Part 2, follows since ( p − m k ( x ) p ) < ≤ m k ( x ) ≤
1, and we know that the limit of thesequence is x p . Finally, Part 3 is proved by applying Dini’s theorem.From the identities Y k = y k ( A ) and M k = m k ( A ) we deduce the following. Corollary 24.
Let A = T ( a ) + E A ∈ QT be self-adjoint and such that a is realvalued and v T Av ≥ for any v ∈ (cid:96) , v (cid:54) = 0 . Then the sequences Y k and M k generated by (15) are such that lim k (cid:107) M k − I (cid:107) = 0 , lim k (cid:107) Y k − A (cid:107) = 0 . As in the square root case, the symbols associated with the matrices in (15)uniformly converge to the symbols of their limit and the compact corrections as-sociated with these matrix sequences converge in the B ( (cid:96) ) norm to the compactcorrections of their limit. Remark 1.
Note that the proofs of Corollaries 22 and 24 rely uniquely on thefact that the spectrum of a self-adjoint operator is contained in the closure ofthe numerical range and thus they can be stated on the milder hypothesis that A is a self-adjoint operator in B ( (cid:96) ). g ( z ) Observe that both the symbol of the ALM mean and the NBMP mean is g ( z ) = ( a ( z ) . . . a p ( z )) /p . The application of the iterations of Section 3 inthe arithmetic of quasi-Toeplitz matrices, provides as a result both the valuesof g ( z ) and of the correction K such that G = T ( g ) + K is the sought geo-metric mean. However, if only the symbol part of G is needed, then it mightbe more convenient to compute the coefficients of g ( z ) by means of the evalua-tion/interpolation technique.In this section, when the symbols are such that a i ∈ W , i = 1 , . . . , p , basedon the evaluation/interpolation at the roots of unity, we provide an algorithm forcomputing the approximation (cid:101) g j , j = − n + 1 , . . . , n , to the Fourier coefficients g j of the symbol g ( z ) = ( a ( z ) · · · a p ( z )) /p .Let n > m = 2 n and let w m = cos πm + i sin πm bethe principal m -th root of 1. There is always a unique Laurent polynomial (cid:101) g ( z ) = (cid:80) n − n +1 (cid:101) g j z j such that g ( w (cid:96)m ) = (cid:101) g ( w (cid:96)m ), that is, (cid:101) g ( z ) interpolates g ( z ) at w (cid:96)m , (cid:96) = − n + 1 , . . . , n . 23he computation of the approximation (cid:101) g j , j = − n + 1 , . . . , n , to the coeffi-cients g j of g ( z ) can proceed by first selecting a positive integer n and evaluating g ( w (cid:96)m ), (cid:96) = − n + 1 , . . . , n ; and then interpolating g ( z ) at w (cid:96)m , (cid:96) = − n + 1 , . . . , n ,by means of the FFT, obtaining the coefficients (cid:101) g j of (cid:101) g ( z ) = (cid:80) nj = − n +1 (cid:101) g j z j .We stop the process if (cid:101) g ( z ) is close enough to g ( z ), otherwise we continue thisprocess by doubling the value of n .Concerning the accuracy of the approximation, we recall the following resultfrom [8, Theorem 3.8](17) (cid:101) g i = g i + ∞ (cid:88) k =1 ( g i + kn + g i − kn ) . If g ∈ W then lim n (cid:80) | j |≥ n | g j | = 0 so that, for a given ε > n such that (cid:80) | j |≥ n | g j | ≤ ε . Thus, from (17) we deducethat | (cid:101) g i − g i | ≤ ε . Under additional assumptions on g ( t ), a guaranteed way fordetermining n such that the above bound is satisfied, can be easily determined(see [8, page 57] for further details). In general, we may adopt the followingheuristic criterion, to halt the evaluation/interpolation procedure (see also [38])where the iteration is terminated if(18) (cid:88) | j | > (cid:100) n (cid:101) | (cid:101) g j | < n (cid:88) j = − n +1 | (cid:101) g j | · ε. We summarize the procedure for computing the coefficients (cid:101) g j , j = − n +1 , . . . , n , as Algorithm 1. The overall computational cost of this algorithm is O ( n log n ) arithmetic operations. Algorithm 1
Approximation of g ( z ) Require:
The coefficients of a i ( z ), i = 1 , . . . , p , and a tolerance ε > Ensure:
Approximations (cid:101) g j , j = − n + 1 , . . . , n , to the coefficients g j of g ( z )such that (cid:80) | j | > (cid:100) n (cid:101) | (cid:101) g j | < (cid:80) nj = − n +1 | (cid:101) g j | · ε . Set n = 4, m = 2 n and w m = cos πm + i sin πm . Evaluate a i ( z ) at z = w (cid:96)m for (cid:96) = − n + 1 , . . . , n and for i = 1 , . . . , p ; For (cid:96) = − n + 1 , . . . , n , compute the p -th root r (cid:96) of ( a ( z ) · · · a p ( z )) at z = w (cid:96)m ; Interpolate the values r (cid:96) , (cid:96) = − n + 1 , . . . , n , by means of FFT and obtainthe coefficients (cid:101) g i of the Laurent polynomial (cid:101) g ( z ) = (cid:80) nj = − n +1 (cid:101) g j z j such that g ( w (cid:96)m ) = (cid:101) g ( w (cid:96)m ), (cid:96) = − n + 1 , . . . , n ; Compute δ m = (cid:80) | j | > (cid:100) n (cid:101) | (cid:101) g j | and κ m = | (cid:101) g | + (cid:80) nj = − n +1 | (cid:101) g j | ; If δ m < κ m ε then exit, else set n = 2 n and compute from Step 2.Observe that if a i > i = 1 , . . . , p , has only real coefficients, then a i ( w (cid:96)m ) = a i ( w − (cid:96)m ) for (cid:96) = − n + 1 , . . . , n and i = 1 , . . . , p , then Step 2 ofAlgorithm 1 can proceed by computing a i ( w (cid:96)m ) for (cid:96) = 0 , . . . , n and setting r (cid:96) = r − (cid:96) for (cid:96) = − , . . . , − n + 1. 24 Geometric means of Finite QT matrices The representation and the arithmetic for QT matrices, up to a certain extent,can be adapted for handing finite dimensional matrices. Accordingly, we showthat iteration (7) for computing the mean G can be applied to finite dimensionalmatrices that can be written as a sum of a Toeplitz matrix and a low-rank matrixcorrection.Given a symbol a ( z ) and m ∈ Z + , we denote by T m ( a ), H m ( a − ) and H m ( a + )the m × m leading principal submatrices of T ( a ), H ( a − ) and H ( a + ), respectively.The following theorem can be seen a finite dimensional version of Theorem 3,which implies that our algorithms for geometric means of QT matrices can beapplied also for finite size matrices. Theorem 25 ([41]) . If a, b ∈ L ∞ , then T m ( a ) T m ( b ) = T m ( ab ) − H m ( a − ) H m ( b + ) − JH m ( a + ) H m ( b − ) J, where J is the m × m flip matrix having 1 on the anti-diagonal and zeros else-where. If we focus on finite size Toeplitz matrices, whose symbols are Laurent poly-nomials of the kind (cid:80) rj = − r a j z j , where the degree r is small compared to thematrix size m , say, r ≤ m/
4, then Theorem 25 shows that the product of ma-trices of this kind can be represented as the sum of a finite Toeplitz matrixand two low-rank matrix corrections with nonzero entries (the support of thecorrection) located in the upper leftmost and in the lower rightmost corners,respectively. Observe also that if the matrices are real symmetric or complexHermitian, then a − = a + and b − = b + so that the compact correction obtainedin the infinite case provides both the corrections in the upper leftmost and thelower rightmost corner.A similar property holds if the matrices can be written as the sum of a bandToeplitz matrix associated with a Laurent polynomial of degree at most r , anda correction with support located in the two opposite corners, provided thatthe values of r and of the size s of the support are suitably smaller than m ,say, r, s ≤ m/
4. Moreover, comparing Theorem 25 with Theorem 3 one cansee that the Toeplitz part T m ( ab ) and the correction H m ( a − ) H m ( b + ) coincidewith the m × m leading principal submatrices of the infinite matrices T ( ab )and H ( a − ) H ( b + ), respectively, so that the infinite QT arithmetic provides thecorresponding result of the finite QT arithmetic. This property holds even inthe case where there is an overlapping of the supports in the two corners ofthe corrections, or if the bandwidth r takes large values. In this situation, theimplementation of the arithmetic is still possible [9], but its not efficient fromthe computational point of view.This allows us to implement the computation of the sequences A ( k ) i of m × m matrices generated by (7) relying on the computation of their infinite counter-parts. This implementation leads to an effective computation if the numeri-cal degree of the symbols as well as the size (the maximum between rows and25olumns) of the supports of the corrections of the matrices A ( k ) i remain boundedfrom above by a constant smaller than or equal to m/ In order to apply the theoretical results of the previous sections, we show bysome tests the effectiveness of computing the ALM and the NBMP means ofthree matrices A , A , A ∈ QT , and the convenience of computing the associ-ated symbol g with the evaluation/interpolation method.The algorithms of Section 4 have been implemented in MATLAB by follow-ing the lines of the corresponding implementations for finite positive definite ma-trices of the Matrix Means Toolbox ( http://bezout.dm.unipi.it/software/mmtoolbox/ ). They rely on the package CQT-Toolbox of [9] for the storage andarithmetic of QT matrices. The software can be provided by the authors uponrequest.The tests have been run on a cluster with 128GB of RAM and 24 cores. Theinternal precision of the CQT toolbox has been set to threshold = 1.e-15 ,while the value of ε = ∗ threshold = . e − has been used to truncatethe values of the symbol and of the correction. That is, only the values g i , i = 0 , . . . , k are computed, where k is such that | g j | ≤ ε (cid:107) g (cid:107) ∞ for j > k .Similarly, for the correction E G , written in the form E G = U V T , where U and V are thin matrices, we computed only the values u i,j , i = 1 , . . . , m and v i,j , i = 1 , . . . , n , where m, n are such that | u i,j | < ε max i | u i,j | for any i > m and j and | v i,j | < ε max i | v i,j | for any i > n and j , respectively.The test examples have been constructed by relying on trigonometric sym-bols in the class f ( t ) = f + 2 f cos( t ) + 2 f cos(2 t ), for different values of thecoefficients f , f , f . A Toeplitz matrix associated with a symbol of this typeturns out to be pentadiagonal.With the choices of ( f , f , f ) in the set { (2 , , , (3 , , , (9 , , } , thesymbol f ( t ) takes values in the intervals [0 , , , f into f + ϑ where ϑ >
0, we are able to tune the ratio max t f ( t ) / min t f ( t )which for finite matrices is related to the condition number of the associatedToeplitz matrix. More specifically, with the choices ϑ = 1 , . , .
01 we havethree sets of matrices with increasing condition numbers.In Table 1, for each value of ϑ we report the numerical length of the sym-bol g ( z ) together with the CPU time needed by the evaluation/interpolationtechnique to compute the coefficients of g ( z ). The CPU time needed for thiscomputation is quite negligible and the numerical length of the symbol, as wellas the number of interpolation points, grow as the ratio min t g ( t ) / max t g ( t ) getscloser to 0, as expected.In Table 2, for each value of ϑ we report the number of iterations and theCPU time needed to compute the ALM mean G ALM and the NBMP mean G NBMP , together with the numerical size and the numerical rank of the compactcorrection. 26 length n CPU1 110 512
Table 1: Length of the symbol, number of interpolation points and CPU timein seconds for computing the coefficients of g ( z ).ALM NBMP ϑ iter. CPU size rank iter. CPU size rank1 43
95 17 3
108 160.1 43
290 30 3 3.2e1 345 280.01 43 g ( t ). A closer analysis, performed by using the MATLAB profiler, shows thatthe most part of time is spent to perform compression operations in the CQT-Toolbox. By compressing a matrix E with respect to a threshold value ε , wemean to find a matrix (cid:101) E of the lowest rank such that (cid:107) E − (cid:101) E (cid:107) ≤ ε .In order to give an idea of the structure of the geometric mean, for ϑ = 1,in Figure 1 we show in logarithmic scale the graph of the symbol g ( t ), i.e.,the plot of the pairs ( j, log | g j | ) for j ≥
0, and of the compact corrections E = ( e i,j ), i.e., the plot of the triples ( i, j, log | e i,j | ) for G ALM . While in Figure2 we show the graph of the correction of G NBMP together with the modulus ofcomponentwise difference G ALM − G NBMP .We may see that the geometric mean, even though is an infinite dimensionalmatrix, can be effectively approximated by a finite number of parameters onceit is decomposed into its Toeplitz part and its compact correction. From thegraphs in Figures 1 and 2 we may appreciate some rounding errors in the com-pact corrections. For the ALM mean, errors are located in the components ofmodulus less than 10 − , while for the NBMP mean the affected componentshave modulus less than 10 − along an edge of the domain.Since the ALM and the NBMP means differ only in the compact correction,an interesting remark is that this difference has values of modulus less than 10 − ,27igure 1: Absolute value of the symbol and of the correction of the ALM meanin log scale.Figure 2: Absolute value of the symbol of the NBMP mean (left) and of thedifference between the ALM and the NBMP mean (right) in log scale.see Figure 2. This reflects also in the infinite dimensional case, the relativelysmall difference that has been usually observed between the two mean in thefinite dimensional case [5].Finally, we have compared the CPU time needed by the ALM and NBMPiterations applied to QT matrices and to their finite truncation to size n . Thevalue of n is chosen equal to 3 times the size of the compact correction. Thischoice is motivated by the goal to keep well separated the compact correction inthe top leftmost corner from that in the bottom rightmost corner of the finitesize matrix. The CPU time is reported in Table 3. Observe that, while in thecase ϑ = 1 of well conditioned matrices, the application of the ALM and NBMPiterations is faster for finite size matrices than for QT matrices, in the slightlymore ill conditioned cases the time needed for finite size matrices gets largerthan the time taken for QT matrices. In particular, for ϑ = 0 .
01 the speed-upreached by the computation based on the quasi-Toeplitz technology is roughly3 for the ALM mean and 9.7 for the NBMP mean.28
ALM NBMP1
Table 3: CPU time in seconds, needed to compute the ALM and the NBMPmeans in the case of finite matrices whose size is 3 times the size of the cor-responding compact correction. For comparison, the CPU time needed in theinfinite case is written between bracket.
The common definitions of geometric means of positive definite matrices, namelythe ALM, NBMP, weighted and Karcher mean, have been extended to the caseof infinite matrices which are the sum of a Toeplitz matrix, associated witha continuous symbol, plus a compact correction. We have shown that thesemeans still keep the same structure of the input matrices, i.e., they can bewritten as the sum of a Toeplitz matrix and a compact correction. Moreover, thesymbol associated with the Toeplitz part of the mean is the mean of the symbolsassociated with the averaged matrices. We have given, moreover, conditionsunder which the geometric mean belongs to B ( (cid:96) ∞ ) and under which the ALMsequence of the symbols converges in Wiener norm.Numerical computations show that the ALM and the NBMP means can becomputed also for infinite matrices in the QT class and that these means canbe represented with good precision in terms of a finite number of parameters.Perhaps surprising, these ideas can be useful also when computing means offinite matrices.An open issue, which we would like to investigate further, is exploiting theavailability of the symbol g ( t ) of the geometric mean G in order to acceleratethe computation of the compact correction. In fact, even though g ( t ) can becomputed separately at a much lower cost, our current implementations of theALM and NBMP iterations do not take advantage of the availability of g ( t )and compute simultaneously approximation to the symbol and to the compactcorrection of the mean. Moreover, the large CPU time needed for computingthe correction is mainly spent for the operation of low-rank approximations oflarge matrices. We believe that this part can be much improved by means ofrandomized techniques both for the task of computing general matrix functionsand for computing the geometric mean. References [1] T. Ando, C.-K. Li, and R. Mathias. Geometric means.
Linear AlgebraAppl. , 385:305–334, 2004. 292] J. Appell, J. Bana´ s , and N. Merentes. Bounded Variation and Around . DeGruyter, Berlin, 2013.[3] M. S. Bernstein. Sur la convergence absolue des s´eries trigonom´etriques.
Comptes rendu , t. 158:1161–1163, 1914.[4] R. Bhatia.
Positive definite matrices . Princeton Series in Applied Mathe-matics. Princeton University Press, Princeton, NJ, 2007. [2015] paperbackedition of the 2007 original [ MR2284176].[5] D. A. Bini and B. Iannazzo. A note on computing matrix geometric means.
Adv. Comput. Math. , 35(2-4):175–192, 2011.[6] D. A. Bini and B. Iannazzo. Computing the Karcher mean of symmetricpositive definite matrices.
Linear Algebra Appl. , 438(4):1700–1710, 2013.[7] D. A. Bini, B. Iannazzo, B. Jeuris, and R. Vandebril. Geometric means ofstructured matrices.
BIT , 54(1):55–83, 2014.[8] D. A. Bini, G. Latouche, and B. Meini.
Numerical methods for structuredMarkov chains . Numerical Mathematics and Scientific Computation. Ox-ford University Press, New York, 2005. Oxford Science Publications.[9] D. A. Bini, S. Massei, and L. Robol. Quasi-Toeplitz matrix arithmetic: amatlab toolbox.
Numerical Algorithms , 81:741–769, 2019.[10] D. A. Bini, B. Meini, and F. Poloni. An effective matrix geometric meansatisfying the Ando-Li-Mathias properties.
Math. Comp. , 79(269):437–452,2010.[11] A. B¨ottcher and S. M. Grudsky.
Toeplitz matrices, asymptotic linear alge-bra, and functional analysis . Birkh¨auser Verlag, Basel, 2000.[12] A. B¨ottcher and S. M. Grudsky.
Spectral properties of banded Toeplitz ma-trices . Society for Industrial and Applied Mathematics (SIAM), Philadel-phia, PA, 2005.[13] A. B¨ottcher and B. Silbermann.
Introduction to large truncated Toeplitzmatrices . Universitext. Springer-Verlag, New York, 1999.[14] M. Fasi and B. Iannazzo. Computing the weighted geometric mean of twolarge-scale matrices and its inverse times a vector.
SIAM Journal on MatrixAnalysis and Applications , 39(1):178–203, 2018.[15] N. J. Higham.
Functions of Matrices: Theory and Computation . Societyfor Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008.[16] S. Hildebrandt. The closure of the numerical range of an operator as spec-tral set.
Comm. Pure Appl. Math. , 17:415–421, 1964.3017] B. Iannazzo. A note on computing the matrix square root.
Calcolo ,40(4):273–283, 2003.[18] B. Iannazzo. The geometric mean of two matrices from a computationalviewpoint.
Numer. Linear Algebra Appl. , 23(2):208–229, 2016.[19] B. Iannazzo, B. Jeuris, and F. Pompili. The Derivative of the MatrixGeometric Mean with an Application to the Nonnegative Decompositionof Tensor Grids. In
Structured Matrices in Numerical Linear Algebra , pages107–128. Springer, 2019.[20] B. Iannazzo and M. Porcelli. The Riemannian Barzilai–Borwein methodwith nonmonotone line search and the matrix geometric mean computation.
IMA Journal of Numerical Analysis , 38(1):495–517, 04 2017.[21] B. Jeuris and R. Vandebril. The K¨ahler mean of block-Toeplitz matriceswith Toeplitz structured blocks.
SIAM J. Matrix Anal. Appl. , 37(3):1151–1175, 2016.[22] B. Jeuris, R. Vandebril, and B. Vandereycken. A survey and compari-son of contemporary algorithms for computing the matrix geometric mean.
Electron. Trans. Numer. Anal. , 39:379–402, 2012.[23] R. V. Kadison and J. R. Ringrose.
Fundamentals of the theory of operatoralgebras. Vol. I , volume 100 of
Pure and Applied Mathematics . AcademicPress, Inc. [Harcourt Brace Jovanovich, Publishers], New York, 1983. Ele-mentary theory.[24] Y. Katznelson.
An Introduction to Harmonic Analysis . Cambridge Univer-sity Press, New York, 3rd edition, 2004.[25] S. Lang. Complex Analysis.
Springer-Verlag , 1999.[26] J. Lapuyade-Lahorgue and F. Barbaresco. Radar detection using siegeldistance between autoregressive processes, application to hf and x-bandradar. In , pages 1–6, 2008.[27] J. Lawson. Existence and uniqueness of the Karcher mean on unital C ∗ -algebras. J. Math. Anal. Appl. , 483(2):123625, 16, 2020.[28] J. Lawson, H. Lee, and Y. Lim. Weighted geometric means.
Forum Math. ,24(5):1067–1090, 2012.[29] J. Lawson and Y. Lim. Karcher means and Karcher equations of positivedefinite operators.
Trans. Amer. Math. Soc. Ser. B , 1:1–22, 2014.[30] H. Lee, Y. Lim, and T. Yamazaki. Multi-variable weighted geometric meansof positive definite matrices.
Linear Algebra Appl. , 435(2):307–322, 2011.[31] P. L´evy. Sur la convergence absolue des s´eries de Fourier.
CompositioMath. , 1:1–14, 1935. 3132] M. Moakher. A differential geometric approach to the geometric meanof symmetric positive-definite matrices.
SIAM J. Matrix Anal. Appl. ,26(3):735–747, 2005.[33] M. Moakher. On the averaging of symmetric positive-definite tensors.
J.Elasticity , 82(3):273–296, 2006.[34] N. Nakamura. Geometric means of positive operators.
Kyungpook Math.J. , 49(1):167–181, 2009.[35] E. Nobari. A monotone geometric mean for a class of Toeplitz matrices.
Linear Algebra Appl. , 511:1–18, 2016.[36] E. Nobari and B. Ahmadi Kakavandi. A geometric mean for Toeplitz andToeplitz-block block-Toeplitz matrices.
Linear Algebra Appl. , 548:189–202,2018.[37] Y. Rathi, A. Tannenbaum, and O. Michailovich. Segmenting images onthe tensor manifold. In , pages 1–8, 2007.[38] L. Robol. Rational Krylov and ADI iteration for infinite size quasi-Toeplitzmatrix equations.
Linear Algebra Appl. , 604:210–235, 2020.[39] A. C. Thompson. On certain contraction mappings in a partially orderedvector space.
Proc. Amer. Math. soc. , 14:438–443, 1963.[40] Y. Wang, S. Qiu, X. Ma, and H. He. A prototype-based spd matrix net-work for domain adaptation eeg emotion recognition.
Pattern Recognition ,110:107626, 2021.[41] H. Widom. Asymptotic behavior of block Toeplitz matrices and determi-nants. II.
Advances in Math. , 21:1–29, 1976.[42] L. Yang, M. Arnaudon, and F. Barbaresco. Geometry of covariance ma-trices and computation of median. In
Bayesian inference and maximumentropy methods in science and engineering , volume 1305 of
AIP Conf.Proc. , pages 479–486. Amer. Inst. Phys., Melville, NY, 2010.[43] F. Yger, M. Berar, and F. Lotte. Riemannian approaches in brain-computerinterfaces: A review.
IEEE Transactions on Neural Systems and Rehabili-tation Engineering , 25(10):1753–1762, 2017.[44] X. Yuan, W. Huang, P.-A. Absil, and K. A. Gallivan. Computing thematrix geometric mean: Riemannian versus Euclidean conditioning, im-plementation techniques, and a Riemannian BFGS method.
NumericalLinear Algebra with Applications , 27(5):e2321, 2020.3245] P. Zanini, M. Congedo, C. Jutten, S. Said, and Y. Berthoumieu. Trans-fer Learning: A Riemannian Geometry Framework With Applications toBrain–Computer Interfaces.
IEEE Transactions on Biomedical Engineer-ing , 65(5):1107–1116, 2018.[46] A. Zygmund.